ALICE3-TRK: adapting to new segmentation in chips in ML/OT and using o2::trk::Hit instead of o2::itsmft::Hit

Author: atriolo

Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/SegmentationChip.h

@@ -88,7 +88,7 @@ class SegmentationChip
   /// \param int disk Disk number (0 to 5 for VD)
   static bool localToDetector(float xRow, float zCol, int& iRow, int& iCol, int subDetID, int layer, int disk) noexcept
   {
-    if (!isValidGlob(xRow, zCol, subDetID, layer)) {
+    if (!isValidLoc(xRow, zCol, subDetID, layer)) {
       LOGP(debug, "Local coordinates not valid: row = {} cm, col = {} cm", xRow, zCol);
       return false;
     }
@@ -116,68 +116,46 @@ class SegmentationChip
       maxWidth = constants::VD::petal::layer::width[layer];
       maxLength = constants::VD::petal::layer::length;
       // TODO: change this to use the layer and disk
-    } else if (subDetID == 1 && layer <= 3) { // ML
+    } else if (subDetID == 1) {
       pitchRow = PitchRowMLOT;
       pitchCol = PitchColMLOT;
-      maxWidth = constants::ML::width;
-      maxLength = constants::ML::length;
-    } else if (subDetID == 1 && layer == 4) { // ML/OT (mixed layer, length = ML but staggered as OT)
-      pitchRow = PitchRowMLOT;
-      pitchCol = PitchColMLOT;
-      maxWidth = constants::OT::halfstave::width;
-      maxLength = constants::ML::length;
-    } else if (subDetID == 1 && layer > 4) { // OT
-      pitchRow = PitchRowMLOT;
-      pitchCol = PitchColMLOT;
-      maxWidth = constants::OT::halfstave::width;
-      maxLength = constants::OT::halfstave::length;
+      maxWidth = constants::moduleMLOT::chip::width - constants::moduleMLOT::chip::passiveEdgeReadOut;
+      maxLength = constants::moduleMLOT::chip::length;
     }
     // convert to row/col
-    iRow = static_cast<int>(((maxWidth / 2 - xRow) / pitchRow));
-    iCol = static_cast<int>(((zCol + maxLength / 2) / pitchCol));
+    iRow = static_cast<int>(std::floor((maxWidth / 2 - xRow) / pitchRow));
+    iCol = static_cast<int>(std::floor((zCol + maxLength / 2) / pitchCol));
   };
 
   // Check local coordinates (cm) validity.
-  static constexpr bool isValidGlob(float x, float z, int subDetID, int layer) noexcept
+  static constexpr bool isValidLoc(float x, float z, int subDetID, int layer) noexcept
   {
     float maxWidth(0), maxLength(0);
     if (subDetID == 0) {
       maxWidth = constants::VD::petal::layer::width[layer];
       maxLength = constants::VD::petal::layer::length;
       // TODO: change this to use the layer and disk
-    } else if (subDetID == 1 && layer <= 3) { // ML
-      maxWidth = constants::ML::width;
-      maxLength = constants::ML::length;
-    } else if (subDetID == 1 && layer == 4) { // ML/OT (mixed layer, length = ML but staggered as OT)
-      maxWidth = constants::OT::halfstave::width;
-      maxLength = constants::ML::length;
-    } else if (subDetID == 1 && layer > 4) { // OT
-      maxWidth = constants::OT::halfstave::width;
-      maxLength = constants::OT::halfstave::length;
+    } else if (subDetID == 1) { // ML/OT
+      maxWidth = constants::moduleMLOT::chip::width - constants::moduleMLOT::chip::passiveEdgeReadOut;
+      maxLength = constants::moduleMLOT::chip::length;
     }
     return (-maxWidth / 2 < x && x < maxWidth / 2 && -maxLength / 2 < z && z < maxLength / 2);
   }
 
   // Check detector coordinates validity.
-  static constexpr bool isValidDet(float row, float col, int subDetID, int layer) noexcept
+  static constexpr bool isValidDet(int row, int col, int subDetID, int layer) noexcept
   {
     // Check if the row and column are within the valid range
     int nRows(0), nCols(0);
     if (subDetID == 0) {
       nRows = constants::VD::petal::layer::nRows[layer];
       nCols = constants::VD::petal::layer::nCols;
       // TODO: change this to use the layer and disk
-    } else if (subDetID == 1 && layer <= 3) { // ML
-      nRows = constants::ML::nRows;
-      nCols = constants::ML::nCols;
-    } else if (subDetID == 1 && layer == 4) { // ML/OT (mixed layer, length = ML but staggered as OT)
-      nRows = constants::OT::halfstave::nRows;
-      nCols = constants::ML::nCols;
-    } else if (subDetID == 1 && layer > 4) { // OT
-      nRows = constants::OT::halfstave::nRows;
-      nCols = constants::OT::halfstave::nCols;
+    } else if (subDetID == 1) {
+      nRows = constants::moduleMLOT::chip::nRows;
+      nCols = constants::moduleMLOT::chip::nCols;
     }
-    return (row >= 0 && row < static_cast<float>(nRows) && col >= 0 && col < static_cast<float>(nCols));
+    return (row >= 0 && row < nRows && col >= 0 && col < nCols);
   }
 
   /// Transformation from Detector cell coordinates to Geant detector centered
@@ -202,7 +180,7 @@ class SegmentationChip
     detectorToLocalUnchecked(iRow, iCol, xRow, zCol, subDetID, layer, disk);
     LOG(debug) << "Result from detectorToLocalUnchecked: iRow " << iRow << " -> xRow " << xRow << ", iCol " << iCol << " -> zCol " << zCol << " on subDetID, layer, disk: " << subDetID << " " << layer << " " << disk;
 
-    if (!isValidGlob(xRow, zCol, subDetID, layer)) {
+    if (!isValidLoc(xRow, zCol, subDetID, layer)) {
       LOGP(debug, "Local coordinates not valid: row = {} cm, col = {} cm", xRow, zCol);
       return false;
     }
@@ -218,15 +196,9 @@ class SegmentationChip
     if (subDetID == 0) {
       xRow = 0.5 * (constants::VD::petal::layer::width[layer] - PitchRowVD) - (row * PitchRowVD);
       zCol = col * PitchColVD + 0.5 * (PitchColVD - constants::VD::petal::layer::length);
-    } else if (subDetID == 1 && layer <= 3) { // ML
-      xRow = 0.5 * (constants::ML::width - PitchRowMLOT) - (row * PitchRowMLOT);
-      zCol = col * PitchRowMLOT + 0.5 * (PitchRowMLOT - constants::ML::length);
-    } else if (subDetID == 1 && layer == 4) { // ML/OT (mixed layer, length = ML but staggered as OT)
-      xRow = 0.5 * (constants::OT::halfstave::width - PitchRowMLOT) - (row * PitchRowMLOT);
-      zCol = col * PitchRowMLOT + 0.5 * (PitchRowMLOT - constants::ML::length);
-    } else if (subDetID == 1 && layer > 4) { // OT
-      xRow = 0.5 * (constants::OT::halfstave::width - PitchRowMLOT) - (row * PitchRowMLOT);
-      zCol = col * PitchColMLOT + 0.5 * (PitchColMLOT - constants::OT::halfstave::length);
+    } else if (subDetID == 1) { // ML/OT
+      xRow = 0.5 * (constants::moduleMLOT::chip::width - constants::moduleMLOT::chip::passiveEdgeReadOut - PitchRowMLOT) - (row * PitchRowMLOT);
+      zCol = col * PitchRowMLOT + 0.5 * (PitchRowMLOT - constants::moduleMLOT::chip::length);
     }
   }
 
@@ -282,20 +254,13 @@ class SegmentationChip
     LOG(info) << "Number of rows:\nVD L0: " << constants::VD::petal::layer::nRows[0]
               << "\nVD L1: " << constants::VD::petal::layer::nRows[1]
               << "\nVD L2: " << constants::VD::petal::layer::nRows[2]
-              << "\nML stave: " << constants::ML::nRows
-              << "\nOT half stave: " << constants::OT::halfstave::nRows;
+              << "\nML/OT chip: " << constants::moduleMLOT::chip::nRows;
 
     LOG(info) << "Number of cols:\nVD: " << constants::VD::petal::layer::nCols
-              << "\nML stave: " << constants::ML::nCols
-              << "\nOT half stave: " << constants::OT::halfstave::nCols;
-
-    LOG(info) << "Pitch rows [cm]:\nVD: " << PitchRowVD
-              << "\nML stave: " << PitchRowMLOT
-              << "\nOT stave: " << PitchRowMLOT;
+              << "\nML/OT chip: " << constants::moduleMLOT::chip::nCols;
 
-    LOG(info) << "Pitch cols [cm]:\nVD: " << PitchColVD
-              << "\nML stave: " << PitchColMLOT
-              << "\nOT stave: " << PitchColMLOT;
+    LOG(info) << "Pitch rows x cols [um]:\nVD: " << PitchRowVD * 1e4 << "x" << PitchColVD * 1e4
+              << "\nML/OT chip: " << PitchRowMLOT * 1e4 << "x" << PitchColMLOT * 1e4;
   }
 };
 

Detectors/Upgrades/ALICE3/TRK/base/include/TRKBase/Specs.h

@@ -76,15 +76,14 @@ constexpr double thickness{0 * mu}; // thickness of the copper metal stack - for
 } // namespace metalstack
 namespace chip
 {
-constexpr double width{25 * mm};                                             // width of the chip
-constexpr double length{32 * mm};                                            // length of the chip
-constexpr double pitchX{50 * mu};                                            // pitch of the row
-constexpr double pitchZ{50 * mu};                                            // pitch of the column
-constexpr int nRows{static_cast<int>(width / pitchX)};                       // number of columns in the chip
-constexpr int nCols{static_cast<int>(length / pitchZ)};                      // number of rows in the chip
-constexpr double totalThickness{silicon::thickness + metalstack::thickness}; // total thickness of the chip
-/// Set to 0 for the moment, to be adjusted with the actual design of the chip if needed
-static constexpr float passiveEdgeReadOut{1.5 * mm}; // width of the readout edge
+constexpr double width{25 * mm};                                              // width of the chip
+constexpr double length{32 * mm};                                             // length of the chip
+constexpr double pitchX{50 * mu};                                             // pitch of the row
+constexpr double pitchZ{50 * mu};                                             // pitch of the column
+constexpr double totalThickness{silicon::thickness + metalstack::thickness};  // total thickness of the chip
+static constexpr double passiveEdgeReadOut{1.5 * mm};                         // width of the readout edge -> dead zone
+constexpr int nRows{static_cast<int>((width - passiveEdgeReadOut) / pitchX)}; // number of rows in the chip
+constexpr int nCols{static_cast<int>(length / pitchZ)};                       // number of columns in the chip
 } // namespace chip
 namespace gaps
 {

Detectors/Upgrades/ALICE3/TRK/simulation/include/TRKSimulation/Digitizer.h

@@ -25,7 +25,7 @@
 #include "TRKSimulation/ChipDigitsContainer.h"
 
 #include "TRKSimulation/DigiParams.h"
-#include "ITSMFTSimulation/Hit.h"
+#include "TRKSimulation/Hit.h"
 #include "TRKBase/GeometryTGeo.h"
 #include "DataFormatsITSMFT/Digit.h"
 #include "DataFormatsITSMFT/ROFRecord.h"
@@ -54,7 +54,7 @@ class Digitizer
   o2::trk::ChipSimResponse* getChipResponse(int chipID);
 
   /// Steer conversion of hits to digits
-  void process(const std::vector<itsmft::Hit>* hits, int evID, int srcID);
+  void process(const std::vector<o2::trk::Hit>* hits, int evID, int srcID);
   void setEventTime(const o2::InteractionTimeRecord& irt);
   double getEndTimeOfROFMax() const
   {
@@ -83,7 +83,7 @@ class Digitizer
   void setDeadChannelsMap(const o2::itsmft::NoiseMap* mp) { mDeadChanMap = mp; }
 
  private:
-  void processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID);
+  void processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, int srcID);
   void registerDigits(o2::trk::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF,
                       uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl);
 
@@ -102,31 +102,27 @@ class Digitizer
   /// Get the number of columns according to the subdetector
   /// \param subDetID 0 for VD, 1 for ML/OT
   /// \param layer 0 to 2 for VD, 0 to 7 for ML/OT
-  /// \return Number of columns (for the moment, in the entire layer(VD) or stave (ML/OT)
+  /// \return Number of columns (In the entire layer(VD) or chip (ML/OT)
   int getNCols(int subDetID, int layer)
   {
     if (subDetID == 0) { // VD
       return constants::VD::petal::layer::nCols;
-    } else if (subDetID == 1 && layer <= 3) { // ML
-      return constants::ML::nCols;
-    } else if (subDetID == 1 && layer >= 4) { // OT
-      return constants::OT::nCols;
+    } else if (subDetID == 1) { // ML/OT: the smallest element is a chip of 500 rows and 640 cols
+      return constants::moduleMLOT::chip::nCols;
     }
     return 0;
   }
 
   /// Get the number of rows according to the subdetector
   /// \param subDetID 0 for VD, 1 for ML/OT
   /// \param layer 0 to 2 for VD, 0 to 7 for ML/OT
-  /// \return Number of rows (for the moment, in the entire layer(VD) or stave (ML/OT)
+  /// \return Number of rows (In the entire layer(VD) or chip (ML/OT)
   int getNRows(int subDetID, int layer)
   {
     if (subDetID == 0) { // VD
       return constants::VD::petal::layer::nRows[layer];
-    } else if (subDetID == 1 && layer <= 3) { // ML
-      return constants::ML::nRows;
-    } else if (subDetID == 1 && layer >= 4) { // OT
-      return constants::OT::nRows;
+    } else if (subDetID == 1) { // ML/OT
+      return constants::moduleMLOT::chip::nRows;
     }
     return 0;
   }

Detectors/Upgrades/ALICE3/TRK/simulation/src/Digitizer.cxx

@@ -26,7 +26,7 @@
 #include <fairlogger/Logger.h> // for LOG
 
 using o2::itsmft::Digit;
-using o2::itsmft::Hit;
+using o2::trk::Hit;
 using Segmentation = o2::trk::SegmentationChip;
 
 using namespace o2::trk;
@@ -68,8 +68,8 @@ void Digitizer::init()
   LOG(info) << " Depth max MLOT: " << mChipSimRespMLOT->getDepthMax();
   LOG(info) << " Depth min MLOT: " << mChipSimRespMLOT->getDepthMin();
 
-  float thicknessVD = 0.0095; // cm --- hardcoded based on geometry currently present
-  float thicknessMLOT = 0.1;  // cm --- hardcoded based on geometry currently present
+  float thicknessVD = 0.0095;                                            // cm --- hardcoded based on geometry currently present
+  float thicknessMLOT = o2::trk::SegmentationChip::SiliconThicknessMLOT; // 0.01 cm = 100 um --- based on geometry currently present
 
   mSimRespVDScaleX = o2::trk::constants::apts::pitchX / o2::trk::SegmentationChip::PitchRowVD;
   mSimRespVDScaleZ = o2::trk::constants::apts::pitchZ / o2::trk::SegmentationChip::PitchColVD;
@@ -115,8 +115,8 @@ void Digitizer::process(const std::vector<Hit>* hits, int evID, int srcID)
             << " cont.mode: " << isContinuous()
             << " Min/Max ROFrames " << mROFrameMin << "/" << mROFrameMax;
 
-  // std::cout << "Printing segmentation info: " << std::endl;
-  // SegmentationChip::Print();
+  std::cout << "Printing segmentation info: " << std::endl;
+  SegmentationChip::Print();
 
   // // is there something to flush ?
   if (mNewROFrame > mROFrameMin) {
@@ -252,7 +252,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast)
 }
 
 //_______________________________________________________________________
-void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID)
+void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, int srcID)
 {
   int chipID = hit.GetDetectorID(); //// the chip ID at the moment is not referred to the chip but to a wider detector element (e.g. quarter of layer or disk in VD, stave in ML, half stave in OT)
   int subDetID = mGeometry->getSubDetID(chipID);
@@ -333,13 +333,6 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
   // std::cout<< "Example hit in local frame: " << exampleLoc << std::endl;
   // std::cout<<"Going back to glob coordinates: " << (matrix * exampleLoc) << std::endl;
 
-  //// adapting the depth (Y) of the chip to the APTS response maximum depth
-  LOG(debug) << "local original: startPos = " << xyzLocS << ", endPos = " << xyzLocE << std::endl;
-  xyzLocS.SetY(xyzLocS.Y());
-  xyzLocE.SetY(xyzLocE.Y());
-
-  LOG(debug) << "rescaled Y: startPos = " << xyzLocS << ", endPos = " << xyzLocE << std::endl;
-
   math_utils::Vector3D<float> step(xyzLocE);
   step -= xyzLocS;
   step *= nStepsInv; // position increment at each step
@@ -418,7 +411,10 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
   // take into account that the ChipSimResponse depth defintion has different min/max boundaries
   // although the max should coincide with the surface of the epitaxial layer, which in the chip
   // local coordinates has Y = +SensorLayerThickness/2
+  // LOG(info)<<"SubdetID = " << subDetID<< " shift: "<<mSimRespVDShift<<" or "<<mSimRespMLOTShift;
+  // LOG(info)<< " Before shift: S = " << xyzLocS.Y()*1e4 << "  E = " << xyzLocE.Y()*1e4;
   xyzLocS.SetY(xyzLocS.Y() + ((subDetID == 0) ? mSimRespVDShift : mSimRespMLOTShift));
+  // LOG(info)<< " After shift: S = " << xyzLocS.Y()*1e4 << "  E = " << xyzLocE.Y()*1e4;
 
   // collect charge in every pixel which might be affected by the hit
   for (int iStep = nSteps; iStep--;) {
@@ -451,7 +447,8 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
       LOG(debug) << "Error in rspmat for step " << iStep << " / " << nSteps;
       continue;
     }
-    LOG(debug) << "rspmat valid! for step " << iStep << " / " << nSteps << ", (row,col) = (" << row << "," << col << ")";
+    // LOG(info) << "rspmat valid! for step " << iStep << " / " << nSteps << ", (row,col) = (" << row << "," << col << ")";
+    // LOG(info) << "rspmat valid! for step " << iStep << " / " << nSteps << " Y= " << xyzLocS.Y()*1e4 << " , (row,col) = (" << row << "," << col << ")";
     // rspmat->print(); // print the response matrix for debugging
 
     for (int irow = AlpideRespSimMat::NPix; irow--;) {
@@ -472,17 +469,17 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID
   // fire the pixels assuming Poisson(n_response_electrons)
   o2::MCCompLabel lbl(hit.GetTrackID(), evID, srcID, false);
   auto roFrameAbs = mNewROFrame + roFrameRel;
-  LOG(debug) << "Spanning through rows and columns; rowspan = " << rowSpan << " colspan = " << colSpan << " = " << colE << " - " << colS << " +1 " << std::endl;
+  LOG(debug) << "\nSpanning through rows and columns; rowspan = " << rowSpan << " colspan = " << colSpan << " = " << colE << " - " << colS << " +1 ";
   for (int irow = rowSpan; irow--;) {          // irow ranging from 4 to 0
     uint16_t rowIS = irow + rowS;              // row distant irow from the row of the hit start
     for (int icol = colSpan; icol--;) {        // icol ranging from 4 to 0
       float nEleResp = respMatrix[irow][icol]; // value of the probability of the response in this pixel
       if (nEleResp <= 1.e-36) {
         continue;
       }
-      LOG(debug) << "nEleResp: value " << nEleResp << " for pixel " << irow << " " << icol << std::endl;
+      LOG(debug) << "nEleResp: value " << nEleResp << " for pixel " << irow << " " << icol;
       int nEle = gRandom->Poisson(nElectrons * nEleResp); // total charge in given pixel = number of electrons generated in the hit multiplied by the probability of being detected in their position
-      LOG(debug) << "Charge detected in the pixel: " << nEle << " for pixel " << irow << " " << icol << std::endl;
+      LOG(debug) << "Charge detected in the pixel: " << nEle << " for pixel " << irow << " " << icol;
       // ignore charge which have no chance to fire the pixel
       if (nEle < mParams.getMinChargeToAccount()) { /// TODO: substitute with the threshold?
         LOG(debug) << "Ignoring pixel with nEle = " << nEle << " < min charge to account "