Merge 01309a0 into sapling-pr-archive-ktf

Author: ktf

Detectors/FOCAL/base/include/FOCALBase/Geometry.h

@@ -135,6 +135,9 @@ class Geometry
   bool getInsertFrontPadLayers() const { return mInsertFrontPadLayers; }
   bool getInsertHCalReadoutMaterial() const { return mInsertFrontHCalReadoutMaterial; }
 
+  float getDetectorOpeningRight() const { return mGlobal_DetectorOpening_Right; }
+  float getDetectorOpeningLeft() const { return mGlobal_DetectorOpening_Left; }
+
   std::vector<const Composition*> getFOCALMicroModule(int layer) const;
   const Composition* getComposition(int layer, int stack) const;
   std::string_view getTowerGapMaterial() const { return mGlobal_Gap_Material; }
@@ -175,6 +178,9 @@ class Geometry
   float mWaferSizeX = 0.0;      // Wafer X size
   float mWaferSizeY = 0.0;      // Wafer Y size
 
+  float mGlobal_DetectorOpening_Right = 0.0; // detector opening in X
+  float mGlobal_DetectorOpening_Left = 0.0;  // detector opening in Y
+
   // PIX setup
   float mGlobal_Pixel_Size = 0.0;  // pixel size
   float mGlobal_PIX_SizeX = 0.0;   // sensor size X

Detectors/FOCAL/base/src/Geometry.cxx

@@ -351,6 +351,16 @@ void Geometry::setParameters(std::string geometryfile)
         LOG(debug) << "Z-Location of the FoCAL is set to : " << mGlobal_FOCAL_Z0;
       }
 
+      if (command.find("DetectorOpen_Right") != std::string::npos) {
+        mGlobal_DetectorOpening_Right = std::stof(tokens[1]);
+        LOG(debug) << "Detector opening on the right : " << mGlobal_DetectorOpening_Right;
+      }
+
+      if (command.find("DetectorOpen_Left") != std::string::npos) {
+        mGlobal_DetectorOpening_Left = std::stof(tokens[1]);
+        LOG(debug) << "Detector opening on the left : " << mGlobal_DetectorOpening_Left;
+      }
+
       if (command.find("HCAL_TOWER_SIZE") != std::string::npos) {
         mGlobal_HCAL_Tower_Size = std::stof(tokens[1]);
         LOG(debug) << "The size of the HCAL readout tower will be : " << mGlobal_HCAL_Tower_Size;
@@ -578,8 +588,8 @@ void Geometry::setParameters(std::string geometryfile)
             }
           }
         } // end for itowerY
-      }   // end for itowerX
-    }     // end else
+      } // end for itowerX
+    } // end else
     center_z += tmpComp.getThickness();
   } // end loop over pad layer compositions
   LOG(debug) << "============ Created all pad layer compositions (" << mPadCompositionBase.size() << " volumes)";

Detectors/FOCAL/simulation/geometryFiles/geometry_Spaghetti.txt

@@ -59,6 +59,10 @@ COMMAND_INSERT_PIX_AT_L9
   GLOBAL_TOWER_TOLX       0.02  Air
   GLOBAL_TOWER_TOLY       0.8   Al
   GLOBAL_FOCAL_Z          763.5
+# Open the detector on the right and left in cm,
+# can only work if the GLOBAL_HCAL_TOWER_NY is odd number and GLOBAL_HCAL_TOWER_NX is even number
+#  GLOBAL_DetectorOpen_Right 5
+#  GLOBAL_DetectorOpen_Left  5
   GLOBAL_Tower_NX      2
   GLOBAL_Tower_NY      11
   GLOBAL_MIDDLE_TOWER_OFFSET 5

Detectors/FOCAL/simulation/src/Detector.cxx

@@ -451,7 +451,7 @@ void Detector::ConstructGeometry()
   }
 
   float pars[4];
-  pars[0] = (mGeometry->getFOCALSizeX() + 2 * mGeometry->getMiddleTowerOffset()) / 2;
+  pars[0] = (mGeometry->getFOCALSizeX() + 2 * mGeometry->getMiddleTowerOffset() + mGeometry->getDetectorOpeningRight() + mGeometry->getDetectorOpeningLeft()) / 2;
   pars[1] = mGeometry->getFOCALSizeY() / 2;
   pars[2] = mGeometry->getFOCALSizeZ() / 2;
   // Add space to place 2 SiPad layers in front of ECAL
@@ -539,6 +539,8 @@ void Detector::CreateHCALSpaghetti()
     }
   }
 
+  bool splitDet = mGeometry->getDetectorOpeningRight() > 0.0 || mGeometry->getDetectorOpeningLeft() > 0.0;
+
   double TowerSize = mGeometry->getHCALTowerSize();
   double CuBoxThickness = 0.3; // Thickness of the Cu box carrying capillary tubes
 
@@ -598,25 +600,57 @@ void Detector::CreateHCALSpaghetti()
   Columns = 0;
   RowPos = 0.;
   Int_t NumTowers = 1;
-  for (Rows = 0; Rows < nTowersY; Rows++) {
 
-    float ColumnPos = 0.;
-    RowPos = Rows * TowerSize;
-    for (Columns = 0; Columns < nTowersX; Columns++) {
-      ColumnPos = Columns * TowerSize;
-      TGeoTranslation* trans = new TGeoTranslation(ColumnPos - SizeXHCAL / 2 + TowerSize / 2, RowPos - SizeYHCAL / 2 + TowerSize / 2, 0.);
+  if (splitDet) {
+    SizeXHCAL = SizeXHCAL / 2;
 
-      // Remove the Towers that overlaps with the beam pipe
-      Double_t RadialDistance = TMath::Power(trans->GetTranslation()[0], 2) + TMath::Power(trans->GetTranslation()[1], 2);
+    TGeoVolumeAssembly* volHalfHCAL = new TGeoVolumeAssembly("HalfHCAL");
 
-      if (RadialDistance < MinRadius * MinRadius || TMath::Abs(trans->GetTranslation()[0]) > SizeXHCAL / 2) {
-        continue;
+    for (Rows = 0; Rows < nTowersY; Rows++) {
+
+      float ColumnPos = 0.;
+      RowPos = Rows * TowerSize;
+      for (Columns = 0; Columns < nTowersX / 2; Columns++) {
+        ColumnPos = Columns * TowerSize;
+        TGeoTranslation* trans = new TGeoTranslation(ColumnPos - SizeXHCAL / 2 + TowerSize / 2, RowPos - SizeYHCAL / 2 + TowerSize / 2, 0.);
+
+        // Shit the beampipe towers by TowerSize/2
+        if (Rows == nTowersY / 2) {
+          trans->SetDx(trans->GetTranslation()[0] + TowerSize / 2);
+        }
+
+        // Adding the Tower to the HCAL
+        volHalfHCAL->AddNode(volTowerHCAL, NumTowers, trans);
+
+        NumTowers++;
       }
+      volHCAL->AddNode(volHalfHCAL, 1, new TGeoTranslation(SizeXHCAL / 2 + mGeometry->getDetectorOpeningRight(), 0, 0));
+      TGeoRotation* rotFlipZ = new TGeoRotation();
+      rotFlipZ->RotateY(180); // Flip around Y to reverse Z
+      TGeoCombiTrans* combHalf = new TGeoCombiTrans(-SizeXHCAL / 2 - mGeometry->getDetectorOpeningLeft(), 0., 0., rotFlipZ);
+      volHCAL->AddNode(volHalfHCAL, 2, combHalf);
+    }
+  } else {
+    for (Rows = 0; Rows < nTowersY; Rows++) {
 
-      // Adding the Tower to the HCAL
-      volHCAL->AddNode(volTowerHCAL, NumTowers, trans);
+      float ColumnPos = 0.;
+      RowPos = Rows * TowerSize;
+      for (Columns = 0; Columns < nTowersX; Columns++) {
+        ColumnPos = Columns * TowerSize;
+        TGeoTranslation* trans = new TGeoTranslation(ColumnPos - SizeXHCAL / 2 + TowerSize / 2, RowPos - SizeYHCAL / 2 + TowerSize / 2, 0.);
 
-      NumTowers++;
+        // Remove the Towers that overlaps with the beam pipe
+        Double_t RadialDistance = TMath::Power(trans->GetTranslation()[0], 2) + TMath::Power(trans->GetTranslation()[1], 2);
+
+        if (RadialDistance < MinRadius * MinRadius || TMath::Abs(trans->GetTranslation()[0]) > SizeXHCAL / 2) {
+          continue;
+        }
+
+        // Adding the Tower to the HCAL
+        volHCAL->AddNode(volTowerHCAL, NumTowers, trans);
+
+        NumTowers++;
+      }
     }
   }
 
@@ -791,6 +825,8 @@ void Detector::CreateECALGeometry()
   // this shifts all the pixel layers to the center near the beampipe
   double pixshift = geom->getTowerSizeX() - (geom->getGlobalPixelWaferSizeX() * geom->getNumberOfPIXsInX());
 
+  bool splitDet = mGeometry->getDetectorOpeningRight() > 0.0 || mGeometry->getDetectorOpeningLeft() > 0.0;
+
   float offset = pars[2];
   // gMC->Gsvolu("EMSC1", "BOX", idtmed[3698], pars, 4);//Left towers (pixels shifted right)
   // gMC->Gsvolu("EMSC2", "BOX", idtmed[3698], pars, 4);//Right towers (pixels shifted left)
@@ -948,7 +984,7 @@ void Detector::CreateECALGeometry()
   // Place the towers in the ECAL
   // --- Place the ECAL in FOCAL
   float fcal_pars[4];
-  fcal_pars[0] = (geom->getFOCALSizeX() + 2. * geom->getMiddleTowerOffset()) / 2.;
+  fcal_pars[0] = (geom->getFOCALSizeX() + 2. * geom->getMiddleTowerOffset() + mGeometry->getDetectorOpeningRight() + mGeometry->getDetectorOpeningLeft()) / 2.;
   fcal_pars[1] = geom->getFOCALSizeY() / 2.;
   fcal_pars[2] = geom->getECALSizeZ() / 2.;
   fcal_pars[3] = 0.;
@@ -977,9 +1013,13 @@ void Detector::CreateECALGeometry()
     // const auto towerCenter = geom->getGeoTowerCenter(number); //only ECAL part, second parameter = -1 by default
     // xp = std::get<0>towerCenter;
     // std::tie(xp, yp, zp) = geom->getGeoTowerCenter(number);
-    const auto [xp, yp, zp] = geom->getGeoTowerCenter(number); // only ECAL part, second parameter = -1 by default
+    auto [xp, yp, zp] = geom->getGeoTowerCenter(number); // only ECAL part, second parameter = -1 by default
 
     if (itowerx == 0) {
+      if (splitDet) {
+        xp -= geom->getDetectorOpeningLeft();
+      }
+
       TVirtualMC::GetMC()->Gspos("EMSC1", number + 1, "ECAL", xp, yp, 0, 0, "ONLY");
       // Add the SiPad front volumes directly under the FOCAL volume
       if (geom->getInsertFrontPadLayers()) {
@@ -992,6 +1032,10 @@ void Detector::CreateECALGeometry()
       }
     }
     if (itowerx == 1) {
+      if (splitDet) {
+        xp += geom->getDetectorOpeningRight();
+      }
+
       TVirtualMC::GetMC()->Gspos("EMSC2", number + 1, "ECAL", xp, yp, 0, 0, "ONLY");
       // Add the SiPad front volumes directly under the FOCAL volume
       if (geom->getInsertFrontPadLayers()) {

Detectors/ITSMFT/ITS/workflow/src/.ThresholdCalibratorSpec.cxx.swo

@@ -38,7 +38,7 @@ class AlpideRespSimMat
   static int constexpr getNPix() { return NPix; }
 
   AlpideRespSimMat() = default;
-  ~AlpideRespSimMat() = default;
+  virtual ~AlpideRespSimMat() = default;
 
   void adopt(const AlpideRespSimMat& src, bool flipRow = false, bool flipCol = false)
   {
@@ -69,7 +69,7 @@ class AlpideRespSimMat
  private:
   std::array<float, MatSize> data;
 
-  ClassDefNV(AlpideRespSimMat, 1);
+  ClassDef(AlpideRespSimMat, 1);
 };
 
 /*
@@ -91,6 +91,7 @@ class AlpideSimResponse
   int getDepthBin(float pos) const;
   std::string composeDataName(int colBin, int rowBin);
 
+ protected:
   int mNBinCol = 0;                                   /// number of bins in X(col direction)
   int mNBinRow = 0;                                   /// number of bins in Y(row direction)
   int mNBinDpt = 0;                                   /// number of bins in Z(sensor dept)

Detectors/ITSMFT/common/simulation/include/ITSMFTSimulation/AlpideSimResponse.h

@@ -12,12 +12,11 @@
 /// \file SegmentationMosaix.h
 /// \brief Definition of the SegmentationMosaix class
 /// \author [email protected]
+/// \author [email protected]
 
 #ifndef ALICEO2_ITS3_SEGMENTATIONMOSAIX_H_
 #define ALICEO2_ITS3_SEGMENTATIONMOSAIX_H_
 
-#include <type_traits>
-
 #include "MathUtils/Cartesian.h"
 #include "ITS3Base/SpecsV2.h"
 
@@ -43,24 +42,22 @@ class SegmentationMosaix
   // 3. The detector coordinate system. Defined by the row and column segmentation
   //    defined at the upper edge in the flat coord.
 
-  // row,col=0
-  // |
-  // v
-  // x----------------------x
-  // |           |          |
-  // |           |          |
-  // |           |          |                        ^ x
-  // |           |          |                        |
-  // |           |          |                        |
-  // |           |          |                        |
-  // |-----------X----------|  X marks (x,z)=(0,0)   X----> z
-  // |           |          |
+  // O----------------------|
   // |           |          |
+  // |           |          |  ^ x
+  // |           |          |  |
+  // |           |          |  |
+  // |           |          |  |
+  // |           |          |  X----> z   X marks (x,z)=(0,0)
+  // |-----------X----------|
+  // |           |          |  O----> col O marks (row,col)=(0,0)
+  // |           |          |  |
+  // |           |          |  |
+  // |           |          |  v
+  // |           |          |  row
   // |           |          |
-  // |           |          |
-  // |           |          |
-  // |           |          |
-  // x----------------------x
+  // |----------------------|
+
  public:
   constexpr SegmentationMosaix(int layer) : mRadius(static_cast<float>(constants::radiiMiddle[layer])) {}
   constexpr ~SegmentationMosaix() = default;
@@ -79,7 +76,6 @@ class SegmentationMosaix
   static constexpr float PitchCol{constants::pixelarray::pixels::mosaix::pitchZ};
   static constexpr float PitchRow{constants::pixelarray::pixels::mosaix::pitchX};
   static constexpr float SensorLayerThickness{constants::totalThickness};
-  static constexpr float NominalYShift{constants::nominalYShift};
 
   /// Transformation from the curved surface to a flat surface.
   /// Additionally a shift in the flat coordinates must be applied because
@@ -102,10 +98,10 @@ class SegmentationMosaix
     // stack
     float dist = std::hypot(xCurved, yCurved);
     float phi = std::atan2(yCurved, xCurved);
-    xFlat = (mRadius * phi) - WidthH;
     // the y position is in the silicon volume however we need the chip volume (silicon+metalstack)
     // this is accounted by a y shift
-    yFlat = dist - mRadius + NominalYShift;
+    xFlat = WidthH - mRadius * phi;
+    yFlat = dist - mRadius;
   }
 
   /// Transformation from the flat surface to a curved surface
@@ -122,11 +118,12 @@ class SegmentationMosaix
   {
     // MUST align the flat surface with the curved surface with the original pixel array is on and account for metal
     // stack
+    float dist = yFlat + mRadius;
+    float phi = (WidthH - xFlat) / mRadius;
     // the y position is in the chip volume however we need the silicon volume
     // this is accounted by a -y shift
-    float dist = yFlat - NominalYShift + mRadius;
-    xCurved = dist * std::cos((xFlat + WidthH) / mRadius);
-    yCurved = dist * std::sin((xFlat + WidthH) / mRadius);
+    xCurved = dist * std::cos(phi);
+    yCurved = dist * std::sin(phi);
   }
 
   /// Transformation from Geant detector centered local coordinates (cm) to
@@ -142,8 +139,11 @@ class SegmentationMosaix
   /// \param int iCol Detector z cell coordinate.
   constexpr bool localToDetector(float const xRow, float const zCol, int& iRow, int& iCol) const noexcept
   {
+    if (!isValidLoc(xRow, zCol)) {
+      return false;
+    }
     localToDetectorUnchecked(xRow, zCol, iRow, iCol);
-    if (!isValid(iRow, iCol)) {
+    if (!isValidDet(iRow, iCol)) {
       iRow = iCol = -1;
       return false;
     }
@@ -167,49 +167,54 @@ class SegmentationMosaix
   /// center of the sensitive volume.
   /// If iRow and or iCol is outside of the segmentation range a value of -0.5*Dx()
   /// or -0.5*Dz() is returned.
-  constexpr bool detectorToLocal(int const iRow, int const iCol, float& xRow, float& zCol) const noexcept
+  bool detectorToLocal(float const row, float const col, float& xRow, float& zCol) const noexcept
   {
-    if (!isValid(iRow, iCol)) {
+    if (!isValidDet(row, col)) {
       return false;
     }
-    detectorToLocalUnchecked(iRow, iCol, xRow, zCol);
-    return isValid(xRow, zCol);
+    detectorToLocalUnchecked(row, col, xRow, zCol);
+    return isValidLoc(xRow, zCol);
   }
 
   // Same as detectorToLocal w.o. checks.
   // We position ourself in the middle of the pixel.
-  constexpr void detectorToLocalUnchecked(int const iRow, int const iCol, float& xRow, float& zCol) const noexcept
+  void detectorToLocalUnchecked(float const row, float const col, float& xRow, float& zCol) const noexcept
   {
-    xRow = -(static_cast<float>(iRow) + 0.5f) * PitchRow + WidthH;
-    zCol = (static_cast<float>(iCol) + 0.5f) * PitchCol - LengthH;
+    xRow = -(row + 0.5f) * PitchRow + WidthH;
+    zCol = (col + 0.5f) * PitchCol - LengthH;
   }
 
-  bool detectorToLocal(int const row, int const col, math_utils::Point3D<float>& loc) const noexcept
+  bool detectorToLocal(float const row, float const col, math_utils::Point3D<float>& loc) const noexcept
   {
     float xRow{0.}, zCol{0.};
     if (!detectorToLocal(row, col, xRow, zCol)) {
       return false;
     }
-    loc.SetCoordinates(xRow, NominalYShift, zCol);
+    loc.SetCoordinates(xRow, 0.0f, zCol);
     return true;
   }
 
-  void detectorToLocalUnchecked(int const row, int const col, math_utils::Point3D<float>& loc) const noexcept
+  void detectorToLocalUnchecked(float const row, float const col, math_utils::Point3D<float>& loc) const noexcept
   {
     float xRow{0.}, zCol{0.};
     detectorToLocalUnchecked(row, col, xRow, zCol);
-    loc.SetCoordinates(xRow, NominalYShift, zCol);
+    loc.SetCoordinates(xRow, 0.0f, zCol);
   }
 
  private:
+  // Check local coordinates (cm) validity.
   template <typename T>
-  [[nodiscard]] constexpr bool isValid(T const row, T const col) const noexcept
+  constexpr bool isValidLoc(T const x, T const z) const noexcept
   {
-    if constexpr (std::is_floating_point_v<T>) { // compares in local coord.
-      return (-WidthH < row && row < WidthH && -LengthH < col && col < LengthH);
-    } else { // compares in rows/cols
-      return !static_cast<bool>(row < 0 || row >= static_cast<int>(NRows) || col < 0 || col >= static_cast<int>(NCols));
-    }
+    return (-WidthH < x && x < WidthH && -LengthH < z && z < LengthH);
+  }
+
+  // Check detector coordinates validity.
+  template <typename T>
+  constexpr bool isValidDet(T const row, T const col) const noexcept
+  {
+    return (row >= 0 && row < static_cast<T>(NRows) &&
+            col >= 0 && col < static_cast<T>(NCols));
   }
 
   float mRadius;