In NA62 we use a large number of lead glass blocks recovered from the OPAL experiment. These blocks have all trapezoidal shape: the two parallel end faces are isosceles trapezoids with different dimensions and the side face connecting the two major bases of the trapezoids is ortogonal to the two end faces.
The blocks belong to 9 different categories according to the dimensions of the two end faces: these dimensions were defined by OPAL when the blocks were built and were validated by NA62.
To model these blocks with GEANT4 we used the G4Trap class with the block dimensions obtained from OPAL. This did not create a problem up to GEANT4 version 10.3. Since version 10.4 the G4Trap creator crashes reporting a planarity problem on one of the side faces of the trapezoid:
-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : GeomSolids0002
issued by : G4Trap::MakePlanes()
Side face ~-X is not planar for solid: pbglblock
Discrepancy: 0.000332441 mm
-----------------------------------------------------------
*** Dump for solid: pbglblock ***
===================================================
Solid type: G4Trap
Parameters:
half length Z: 185 mm
half length Y, face -Dz: 54.976 mm
half length X, face -Dz, side -Dy1: 54.763 mm
half length X, face -Dz, side +Dy1: 53.5015 mm
half length Y, face +Dz: 49.125 mm
half length X, face +Dz, side -Dy2: 48.7875 mm
half length X, face +Dz, side +Dy2: 47.659 mm
theta: 0.9059720706634017°
phi: -90°
alpha, face -Dz: 0°
alpha, face +Dz: 0°
-----------------------------------------------------------
Checking the code, we found that the planarity tolerance check changed completely from version 10.3 to 10.4: in the older versions it was
if (std::fabs(Vcross.dot(v14)/(Vcross.mag()*v14.mag())) > kCoplanar_Tolerance)
{
good = false;
}
with kCoplanar_Tolerance set to 1E-4, while in the new version the check is
// compute distances and check planarity
G4double d1 = std::abs(normal.dot(p1) + plane.d);
G4double d2 = std::abs(normal.dot(p2) + plane.d);
G4double d3 = std::abs(normal.dot(p3) + plane.d);
G4double d4 = std::abs(normal.dot(p4) + plane.d);
G4double dmax = std::max(std::max(std::max(d1,d2),d3),d4);
return (dmax > 1000 * kCarTolerance) ? false : true;
where kCarTolerance is set by G4GeometryTolerance to 1E-9 mm, thus requiring a planarity of 1E-6 mm (1 nm).
The dimensions of the block are correct but they are known with a precision of O(1um): we could try changing them a bit to remove the planarity problem but this would require a lengthy trial and error procedure we would rather avoid and which would be completely arbitrary.
Is a there a more reasonable way to fix this problem, which also does not require creating a private version of GEANT4 with a looser planarity tolerance in G4Trap?
Thank you
Emanuele Leonardi
P.S. this is the code we use to create the trapezoids. It comes from the DetectorConstruction class.
// Block OPAL 15
G4double blockZLength = 370.000*mm; // Z length (PbGl only)
G4double blockL1Length = 109.952*mm; // L1 length (long base of back face)
G4double blockL2Length = 109.918*mm; // L2 length (short base of back face)
G4double blockL3Length = 95.420*mm; // L3 length (long base of front face)
G4double blockL4Length = 95.391*mm; // L4 length (short base of front face)
G4double blockW1Length = 109.952*mm; // W1 length (height of back face)
G4double blockW2Length = 95.461*mm; // W2 length (height of front face)
// Create solid for PbGl block
G4Trap* solidPbGlBlock = PbGlShape("pbglblock",blockL1Length,blockL2Length,blockL3Length,
blockL4Length,blockW1Length,blockW2Length,blockZLength);
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4Trap* B1DetectorConstruction::PbGlShape
(const G4String& name,
G4double l1, G4double l2, // bases of back face
G4double l3, G4double l4, // bases of front face
G4double w1, G4double w2, // heights of back and front faces
G4double z // Z-length of trapezoid
)
{
G4Trap* thisTrap;
//--------------------------------------
// Create trapezoid solid for PbGl block
//--------------------------------------
// Get half length of sides
G4double hl1 = 0.5*l1;
G4double hl2 = 0.5*l2;
G4double hl3 = 0.5*l3;
G4double hl4 = 0.5*l4;
G4double hw1 = 0.5*w1;
G4double hw2 = 0.5*w2;
G4double hz = 0.5*z;
G4double theta = atan2(0.5*(w1-w2),z);
G4double phi = -90.*deg;
G4cout << "Theta " << theta << " Phi " << phi << G4endl;
thisTrap = new G4Trap(name,
hz, // Half z length
theta, // Polar angle of the line joining the centres of the faces at -/+pDz
phi, // Azimuthal angle of the line joining the centre of the face at -pDz to the centre of the face at +pDz
hw1, // Half y length at -pDz
hl1, // Half x length of the side at y=-pDy1 of the face at -pDz
hl2, // Half x length of the side at y=+pDy1 of the face at -pDz
0., // Angle with respect to the y axis from the centre of the side (lower endcap)
hw2, // Half y length at +pDz
hl3, // Half x length of the side at y=-pDy2 of the face at +pDz
hl4, // Half x length of the side at y=+pDy2 of the face at +pDz
0. // Angle with respect to the y axis from the centre of the side (upper endcap)
);
// Trapezoid
return thisTrap;
}