Geometry stability fixes: remove nanometre gaps and overlaps caused by float rounding

Detectors/Base/include/DetectorsBase/Detector.h

@@ -87,8 +87,8 @@ class Detector : public FairDetector
   /// \param phi1 azimuthal angle for axis I
   /// \param phi2 azimuthal angle for axis II
   /// \param phi3 azimuthal angle for axis III
-  void Matrix(Int_t& nmat, Float_t theta1, Float_t phi1, Float_t theta2, Float_t phi2, Float_t theta3,
-              Float_t phi3) const;
+  void Matrix(Int_t& nmat, Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2, Double_t theta3,
+              Double_t phi3) const;
 
   static void setDensityFactor(Float_t density)
   {

Detectors/Base/src/Detector.cxx

@@ -101,8 +101,8 @@ void Detector::SpecialProcess(Int_t numed, EProc parID, int val)
   mgr.SpecialProcess(GetName(), numed, parID, val);
 }
 
-void Detector::Matrix(Int_t& nmat, Float_t theta1, Float_t phi1, Float_t theta2, Float_t phi2, Float_t theta3,
-                      Float_t phi3) const
+void Detector::Matrix(Int_t& nmat, Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2, Double_t theta3,
+                      Double_t phi3) const
 {
   TVirtualMC::GetMC()->Matrix(nmat, theta1, phi1, theta2, phi2, theta3, phi3);
 }

Detectors/CPV/simulation/include/CPVSimulation/GeometryParams.h

@@ -40,7 +40,7 @@ class GeometryParams final : public TNamed
     return sGeomParam;
   }
 
-  void GetModuleAngle(int module, float angle[3][2]) const
+  void GetModuleAngle(int module, double angle[3][2]) const
   {
     for (int i = 0; i < 3; i++) {
       for (int ian = 0; ian < 2; ian++) {
@@ -49,9 +49,9 @@ class GeometryParams final : public TNamed
     }
   }
 
-  float GetCPVAngle(Int_t index) const { return mCPVAngle[index - 1]; }
+  double GetCPVAngle(Int_t index) const { return mCPVAngle[index - 1]; }
 
-  void GetModuleCenter(int module, float* pos) const
+  void GetModuleCenter(int module, double* pos) const
   {
     for (int i = 0; i < 3; i++) {
       pos[i] = mModuleCenter[module][i];
@@ -61,18 +61,18 @@ class GeometryParams final : public TNamed
   int GetNModules() const { return mNModules; }
   int GetNumberOfCPVPadsPhi() const { return mNumberOfCPVPadsPhi; }
   int GetNumberOfCPVPadsZ() const { return mNumberOfCPVPadsZ; }
-  float GetCPVPadSizePhi() const { return mCPVPadSizePhi; }
-  float GetCPVPadSizeZ() const { return mCPVPadSizeZ; }
-  float GetCPVBoxSize(int index) const { return mCPVBoxSize[index]; }
-  float GetCPVActiveSize(int index) const { return mCPVActiveSize[index]; }
+  double GetCPVPadSizePhi() const { return mCPVPadSizePhi; }
+  double GetCPVPadSizeZ() const { return mCPVPadSizeZ; }
+  double GetCPVBoxSize(int index) const { return mCPVBoxSize[index]; }
+  double GetCPVActiveSize(int index) const { return mCPVActiveSize[index]; }
   int GetNumberOfCPVChipsPhi() const { return mNumberOfCPVChipsPhi; }
   int GetNumberOfCPVChipsZ() const { return mNumberOfCPVChipsZ; }
-  float GetGassiplexChipSize(int index) const { return mGassiplexChipSize[index]; }
-  float GetCPVGasThickness() const { return mCPVGasThickness; }
-  float GetCPVTextoliteThickness() const { return mCPVTextoliteThickness; }
-  float GetCPVCuNiFoilThickness() const { return mCPVCuNiFoilThickness; }
-  float GetFTPosition(int index) const { return mFTPosition[index]; }
-  float GetCPVFrameSize(int index) const { return mCPVFrameSize[index]; }
+  double GetGassiplexChipSize(int index) const { return mGassiplexChipSize[index]; }
+  double GetCPVGasThickness() const { return mCPVGasThickness; }
+  double GetCPVTextoliteThickness() const { return mCPVTextoliteThickness; }
+  double GetCPVCuNiFoilThickness() const { return mCPVCuNiFoilThickness; }
+  double GetFTPosition(int index) const { return mFTPosition[index]; }
+  double GetCPVFrameSize(int index) const { return mCPVFrameSize[index]; }
 
  private:
   ///
@@ -86,23 +86,23 @@ class GeometryParams final : public TNamed
   int mNModules;                // Number of CPV modules
   int mNumberOfCPVPadsPhi;      // Number of CPV pads in phi
   int mNumberOfCPVPadsZ;        // Number of CPV pads in z
-  float mCPVPadSizePhi;         // CPV pad size in phi
-  float mCPVPadSizeZ;           // CPV pad size in z
-  float mCPVBoxSize[3];         // Outer size of CPV box
-  float mCPVActiveSize[2];      // Active size of CPV box (x,z)
+  double mCPVPadSizePhi;        // CPV pad size in phi
+  double mCPVPadSizeZ;          // CPV pad size in z
+  double mCPVBoxSize[3];        // Outer size of CPV box
+  double mCPVActiveSize[2];     // Active size of CPV box (x,z)
   int mNumberOfCPVChipsPhi;     // Number of CPV Gassiplex chips in phi
   int mNumberOfCPVChipsZ;       // Number of CPV Gassiplex chips in z
-  float mGassiplexChipSize[3];  // Size of a Gassiplex chip (0 - in z, 1 - in phi, 2 - thickness (in ALICE radius))
-  float mCPVGasThickness;       // Thickness of CPV gas volume
-  float mCPVTextoliteThickness; // Thickness of CPV textolite PCB (without moil)
-  float mCPVCuNiFoilThickness;  // Thickness of CPV Copper-Nickel moil of PCB
-  float mFTPosition[4];         // Positions of the 4 PCB vs the CPV box center
-  float mCPVFrameSize[3];       // CPV frame size (0 - in phi, 1 - in z, 2 - thickness (along ALICE radius))
-  float mIPtoCPVSurface;        // Distance from IP to CPV front cover
-  float mModuleAngle[5][3][2];  // Orientation angles of CPV modules
-  float mCPVAngle[5];           // Direction to the center of CPV modules in phi
-  float mModuleCenter[5][3];    // Coordunates of modules centra in ALICE system
-  ClassDefOverride(GeometryParams, 1);
+  double mGassiplexChipSize[3]; // Size of a Gassiplex chip (0 - in z, 1 - in phi, 2 - thickness (in ALICE radius))
+  double mCPVGasThickness;      // Thickness of CPV gas volume
+  double mCPVTextoliteThickness; // Thickness of CPV textolite PCB (without moil)
+  double mCPVCuNiFoilThickness;  // Thickness of CPV Copper-Nickel moil of PCB
+  double mFTPosition[4];         // Positions of the 4 PCB vs the CPV box center
+  double mCPVFrameSize[3];       // CPV frame size (0 - in phi, 1 - in z, 2 - thickness (along ALICE radius))
+  double mIPtoCPVSurface;        // Distance from IP to CPV front cover
+  double mModuleAngle[5][3][2];  // Orientation angles of CPV modules
+  double mCPVAngle[5];           // Direction to the center of CPV modules in phi
+  double mModuleCenter[5][3];    // Coordunates of modules centra in ALICE system
+  ClassDefOverride(GeometryParams, 2);
 };
 } // namespace cpv
 } // namespace o2

Detectors/CPV/simulation/src/Detector.cxx

@@ -358,7 +358,7 @@ void Detector::ConstructGeometry()
   CreateMaterials();
 
   // Create a CPV modules-containers which will be filled with the stuff later.
-  float par[3], x, y, z;
+  double par[3], x, y, z;
 
   // The box containing all CPV filled with air
   par[0] = geomParams->GetCPVBoxSize(0) / 2.0;
@@ -372,11 +372,11 @@ void Detector::ConstructGeometry()
     if (!mActiveModule[iModule]) {
       continue;
     }
-    float angle[3][2] = {0};
+    double angle[3][2] = {0};
     geomParams->GetModuleAngle(iModule, angle);
     Matrix(idrotm[iModule], angle[0][0], angle[0][1], angle[1][0], angle[1][1], angle[2][0], angle[2][1]);
 
-    float pos[3] = {0};
+    double pos[3] = {0};
     geomParams->GetModuleCenter(iModule, pos);
 
     fMC->Gspos("CPV", iModule, "barrel", pos[0], pos[1] + 30., pos[2], idrotm[iModule], "ONLY");
@@ -396,8 +396,8 @@ void Detector::ConstructGeometry()
   fMC->Gspos("CPVC", 1, "CPVG", 0, y, 0, 0, "ONLY");
 
   // Position of the chip inside CPV
-  float xStep = geomParams->GetCPVActiveSize(0) / (geomParams->GetNumberOfCPVChipsPhi() + 1);
-  float zStep = geomParams->GetCPVActiveSize(1) / (geomParams->GetNumberOfCPVChipsZ() + 1);
+  double xStep = geomParams->GetCPVActiveSize(0) / (geomParams->GetNumberOfCPVChipsPhi() + 1);
+  double zStep = geomParams->GetCPVActiveSize(1) / (geomParams->GetNumberOfCPVChipsZ() + 1);
   int copy = 0;
   y = geomParams->GetCPVFrameSize(1) / 2 - geomParams->GetFTPosition(0) +
       geomParams->GetCPVTextoliteThickness() / 2 + geomParams->GetGassiplexChipSize(1) / 2 + 0.1;
@@ -451,10 +451,20 @@ void Detector::ConstructGeometry()
   par[2] = geomParams->GetCPVFrameSize(2) / 2;
   fMC->Gsvolu("CPVF2", "BOX ", getMediumID(ID_AL), par, 3);
 
+  // The frame's inner faces are the edges of the active plane. Deriving them from the box and
+  // frame sizes instead of from the active size left CPVF2 38 nm away from CPVF and CPVAr:
+  // algebraically the same plane, different float sums. TGeo relocates from a point pushed far
+  // past a boundary and never notices; VecGeom answers DistanceToOut = -1 and gives the volume
+  // back to its mother. Place the frame against the halves the active volumes were actually
+  // given, so the faces land on the same double.
+  const double activeHalfX = geomParams->GetCPVActiveSize(0) / 2;
+  const double activeHalfZ = geomParams->GetCPVActiveSize(1) / 2;
+  const double frameHalfX = geomParams->GetCPVFrameSize(0) / 2;
+  const double frameHalfZ = geomParams->GetCPVFrameSize(2) / 2;
   for (int j = 0; j <= 1; j++) {
-    x = TMath::Sign(1, 2 * j - 1) * (geomParams->GetCPVBoxSize(0) - geomParams->GetCPVFrameSize(0)) / 2;
+    x = TMath::Sign(1, 2 * j - 1) * (activeHalfX + frameHalfX);
     fMC->Gspos("CPVF1", j + 1, "CPV", x, 0, 0, 0, "ONLY");
-    z = TMath::Sign(1, 2 * j - 1) * (geomParams->GetCPVBoxSize(2) - geomParams->GetCPVFrameSize(2)) / 2;
+    z = TMath::Sign(1, 2 * j - 1) * (activeHalfZ + frameHalfZ);
     fMC->Gspos("CPVF2", j + 1, "CPV", 0, 0, z, 0, "ONLY");
   }
 
@@ -590,7 +600,7 @@ void Detector::addAlignableVolumes() const
     LOG(debug) << "Got TGeoPNEntry " << alignableEntry;
 
     if (alignableEntry) {
-      Float_t angle = geom->GetCPVAngle(iModule);
+      double angle = geom->GetCPVAngle(iModule);
       TGeoHMatrix* globMatrix = alignableEntry->GetGlobalOrig();
 
       TGeoHMatrix* matTtoL = new TGeoHMatrix;

Detectors/FIT/FT0/simulation/src/Detector.cxx

@@ -203,10 +203,17 @@ void Detector::SetOneMCP(TGeoVolume* ins)
 
   Double_t x, y, z;
 
-  Float_t ptop[3] = {1.324, 1.324, 1.};      // Cherenkov radiator
-  Float_t ptopref[3] = {1.3241, 1.3241, 1.}; // Cherenkov radiator wrapped with reflector
-  Double_t prfv[3] = {0.0002, 1.323, 1.};    // Vertical refracting layer bettwen radiators and between radiator and not optical Air
-  Double_t prfh[3] = {1.323, 0.0002, 1.};    // Horizontal refracting layer bettwen radiators and ...
+  Float_t ptop[3] = {1.324, 1.324, 1.};   // Cherenkov radiator
+  Double_t prfv[3] = {0.0002, 1.323, 1.}; // Vertical refracting layer bettwen radiators and between radiator and not optical Air
+  Double_t prfh[3] = {1.323, 0.0002, 1.}; // Horizontal refracting layer bettwen radiators and ...
+  // Cherenkov radiator wrapped with reflector. The container has to hold the radiator plus one
+  // reflector strip on each side, so its half-width is the radiator's plus the full strip
+  // thickness. It was written as the radiator plus 1 um while a strip is 4 um thick, which left
+  // each strip sticking 3 um out of its own mother. The clearance keeps the strips off the
+  // container wall, so no daughter face sits exactly on the mother's.
+  const double kWrapClearance = 1.e-4; // cm
+  Float_t ptopref[3] = {static_cast<Float_t>(ptop[0] + 2 * prfv[0] + kWrapClearance),
+                        static_cast<Float_t>(ptop[1] + 2 * prfh[1] + kWrapClearance), 1.};
   Float_t pmcp[3] = {2.949, 2.949, 0.66};    // MCP
   Float_t pmcpinner[3] = {2.749, 2.749, 0.1};
   Float_t pmcpbase[3] = {2.949, 2.949, 0.675};
@@ -231,39 +238,52 @@ void Detector::SetOneMCP(TGeoVolume* ins)
   Int_t ntops = 0, nrfvs = 0, nrfhs = 0;
   x = y = z = 0;
   topref->AddNode(top, 1, new TGeoTranslation(0, 0, 0));
-  float xinv = -ptop[0] - prfv[0];
+  // Each reflector strip has to sit against the radiator with no gap: SetBorderSurface below
+  // makes 0TOP/0RFV and 0TOP/0RFH mirrors, and a border surface only acts where the two
+  // volumes actually touch. Rounding the offset to float left 0.33 nm of air between them, so
+  // the mirrors never applied and the light left through the side.
+  auto seat = [](Float_t half, Double_t strip) { return double(half) + strip; };
+  double xinv = -seat(ptop[0], prfv[0]);
   topref->AddNode(rfv, 1, new TGeoTranslation(xinv, 0, 0));
-  xinv = ptop[0] + prfv[0];
+  xinv = seat(ptop[0], prfv[0]);
   topref->AddNode(rfv, 2, new TGeoTranslation(xinv, 0, 0));
-  float yinv = -ptop[1] - prfh[1];
+  double yinv = -seat(ptop[1], prfh[1]);
   topref->AddNode(rfh, 1, new TGeoTranslation(0, yinv, 0));
-  yinv = ptop[1] + prfh[1];
+  yinv = seat(ptop[1], prfh[1]);
   topref->AddNode(rfh, 2, new TGeoTranslation(0, yinv, 0));
 
+  // The wrapped radiator, the MCP top glass, the photocathode and the MCP are stacked along z
+  // with no gap between them. Each layer's position is built here from the back face of the one
+  // before it, using the same double addition the box shape uses for its own half-length, so
+  // adjacent faces land on the same double. Summing the float thicknesses in a different order
+  // for each layer, as this code did, left the faces up to 0.24 nm apart.
+  const double zRadiatorCentre = -mInStart[2] + double(ptopref[2]);
+  const double zRadiatorBack = zRadiatorCentre + double(ptopref[2]);
+  const double zTopGlassCentre = zRadiatorBack + double(pmcptopglass[2]);
+  const double zTopGlassBack = zTopGlassCentre + double(pmcptopglass[2]);
+  const double zCathodeCentre = zTopGlassBack + double(preg[2]);
+  const double zCathodeBack = zCathodeCentre + double(preg[2]);
+
   // container for radiator, cathode
   for (Int_t ix = 0; ix < 2; ix++) {
     float xin = -mInStart[0] + 0.3 + (ix + 0.5) * 2 * ptopref[0];
     for (Int_t iy = 0; iy < 2; iy++) {
       float yin = -mInStart[1] + 0.3 + (iy + 0.5) * 2 * ptopref[1];
       ntops++;
-      z = -mInStart[2] + ptopref[2];
-      ins->AddNode(topref, ntops, new TGeoTranslation(xin, yin, z));
-      LOG(debug) << " n " << ntops << " x " << xin << " y " << yin << " z radiator " << z;
-      z += ptopref[2] + 2. * pmcptopglass[2] + preg[2];
-      ins->AddNode(cat, ntops, new TGeoTranslation(xin, yin, z));
-      LOG(debug) << " n " << ntops << " x " << xin << " y " << yin << " z cathod " << z;
+      ins->AddNode(topref, ntops, new TGeoTranslation(xin, yin, zRadiatorCentre));
+      LOG(debug) << " n " << ntops << " x " << xin << " y " << yin << " z radiator " << zRadiatorCentre;
+      ins->AddNode(cat, ntops, new TGeoTranslation(xin, yin, zCathodeCentre));
+      LOG(debug) << " n " << ntops << " x " << xin << " y " << yin << " z cathod " << zCathodeCentre;
     }
   }
   // MCP
   TVirtualMC::GetMC()->Gsvolu("0MTO", "BOX", getMediumID(kOpGlass), pmcptopglass, 3); // Op  Glass
   TGeoVolume* mcptop = gGeoManager->GetVolume("0MTO");
-  z = -mInStart[2] + 2 * ptopref[2] + pmcptopglass[2];
-  ins->AddNode(mcptop, 1, new TGeoTranslation(0, 0, z));
+  ins->AddNode(mcptop, 1, new TGeoTranslation(0, 0, zTopGlassCentre));
 
   TVirtualMC::GetMC()->Gsvolu("0MCP", "BOX", getMediumID(kAir), pmcp, 3); // glass
   TGeoVolume* mcp = gGeoManager->GetVolume("0MCP");
-  z = -mInStart[2] + 2 * ptopref[2] + 2 * pmcptopglass[2] + 2 * preg[2] + pmcp[2];
-  ins->AddNode(mcp, 1, new TGeoTranslation(0, 0, z));
+  ins->AddNode(mcp, 1, new TGeoTranslation(0, 0, zCathodeBack + double(pmcp[2])));
 
   TVirtualMC::GetMC()->Gsvolu("0MSI", "BOX", getMediumID(kMCPwalls), pmcpside, 3); // glass
   TGeoVolume* mcpside = gGeoManager->GetVolume("0MSI");