[ALICE3] Add digitization for forward disks

Author: marcovanleeuwen

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

@@ -25,7 +25,8 @@ namespace trk
 {
 
 struct AlmiraParam : public o2::conf::ConfigurableParamHelper<AlmiraParam> {
-  static constexpr size_t kNLayers = constants::VD::petal::nLayers + constants::ML::nLayers + constants::OT::nLayers; // TODO: include disks
+  // This should be part of geometryTGeo
+  static constexpr size_t kNLayers = constants::VD::petal::nLayers + constants::ML::nLayers + constants::OT::nLayers + constants::MLOTDisks::nLayers;
   static constexpr size_t getNLayers() { return kNLayers; }
 
   int roFrameLengthInBCPerLayer[kNLayers] = {0};  ///< ROF length in BC per layer

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

@@ -89,6 +89,8 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache
   int getNumberOfActivePartsVD() const { return mNumberOfActivePartsVD; }
   int getNumberOfHalfStaves(int lay) const { return mNumberOfHalfStaves[lay]; }
 
+  int getNumberOfDisksMLOT() const { return mNumberOfDisksMLOT; }
+  int getNumberOfStavesInDisk(int lay) const { return mFirstStaveIndexDisc[lay+1] - mFirstStaveIndexDisc[lay]; }
   bool isOwner() const { return mOwner; }
   void setOwner(bool v) { mOwner = v; }
 
@@ -259,8 +261,8 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache
   std::vector<unsigned short> mLastChipIndexVD;   ///< max ID of the detector in the layer for the VD
   // std::vector<unsigned short> mLastChipIndexMLOT; ///< max ID of the detector in the layer for the MLOT
   std::vector<unsigned short> mFirstChipIndexMLOTDisc; ///< ID of the first sensor chip in the layer for the MLOT; array size is one larger than the number of disks; last element equals nChips+1
-  std::vector<int> mFirstStaveIndexDisc;               ///< Index of first stave (abs ID) in each MLOT Disc
-  std::vector<int> mFirstChipIndexStave;               ///< Index of first chip on stave (Discs)
+  std::vector<unsigned short> mFirstStaveIndexDisc;               ///< Index of first stave (abs ID) in each MLOT Disc
+  std::vector<unsigned short> mFirstChipIndexStave;               ///< Index of first chip on stave (Discs)
   std::array<char, MAXLAYERS> mLayerToWrapper; ///< Layer to wrapper correspondence, not implemented yet
 
   bool mOwner = true; //! is it owned by the singleton?

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

@@ -116,7 +116,7 @@ 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) {
+    } else if (subDetID == 1 || subDetID == 2) {
       pitchRow = PitchRowMLOT;
       pitchCol = PitchColMLOT;
       maxWidth = constants::moduleMLOT::chip::width - constants::moduleMLOT::chip::passiveEdgeReadOut;
@@ -135,7 +135,7 @@ 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) { // ML/OT
+    } else if (subDetID == 1 || subDetID == 2) { // ML/OT
       maxWidth = constants::moduleMLOT::chip::width - constants::moduleMLOT::chip::passiveEdgeReadOut;
       maxLength = constants::moduleMLOT::chip::length;
     }
@@ -151,7 +151,7 @@ class SegmentationChip
       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) {
+    } else if (subDetID == 1 || subDetID == 2) {
       nRows = constants::moduleMLOT::chip::nRows;
       nCols = constants::moduleMLOT::chip::nCols;
     }
@@ -196,7 +196,7 @@ 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) { // ML/OT
+    } else if (subDetID == 1 || subDetID == 2) { // ML/OT
       xRow = 0.5 * (constants::moduleMLOT::chip::width - constants::moduleMLOT::chip::passiveEdgeReadOut - PitchRowMLOT) - (row * PitchRowMLOT);
       zCol = col * PitchColMLOT + 0.5 * (PitchColMLOT - constants::moduleMLOT::chip::length);
     }

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

@@ -142,6 +142,12 @@ constexpr double pitchZ{10.0 * mu};
 constexpr double responseYShift{5 * mu}; /// center of the epitaxial layer
 constexpr double thickness{20 * mu};
 } // namespace alice3resp
+
+namespace MLOTDisks
+{
+constexpr int nLayers{12};  // number of disks in the ML and OT
+}
+
 } // namespace o2::trk::constants
 
 #endif

Detectors/Upgrades/ALICE3/TRK/base/src/GeometryTGeo.cxx

@@ -302,7 +302,19 @@ int GeometryTGeo::getLayerTRK(int index) const
     return -1; /// disks do not have a global layer index
   }
   int subDetID = getSubDetID(index);
-  return subDetID * o2::trk::constants::VD::petal::nLayers + getLayer(index); // MLOT: offset by number of VD layers
+  int firstDetLayer = 0;
+  
+  // NOTE: taking these from o2::trk::constants, instead 
+  // of the geometry that is constructed in the 'Build' function
+  // risks inconsistencies...
+
+  if (subDetID == 1) {
+    firstDetLayer = o2::trk::constants::VD::petal::nLayers;
+  }
+  else if (subDetID == 2) {
+    firstDetLayer = o2::trk::constants::VD::petal::nLayers + o2::trk::constants::ML::nLayers + o2::trk::constants::OT::nLayers;
+  }
+  return firstDetLayer + getLayer(index);
 }
 //__________________________________________________________________________
 int GeometryTGeo::getStave(int index) const

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

@@ -113,7 +113,7 @@ class Digitizer
   {
     if (subDetID == 0) { // VD
       return constants::VD::petal::layer::nCols;
-    } else if (subDetID == 1) { // ML/OT: the smallest element is a chip of 470 rows and 640 cols
+    } else if (subDetID == 1 || subDetID == 2) { // ML/OT: the smallest element is a chip of 470 rows and 640 cols
       return constants::moduleMLOT::chip::nCols;
     }
     return 0;
@@ -127,7 +127,7 @@ class Digitizer
   {
     if (subDetID == 0) { // VD
       return constants::VD::petal::layer::nRows[layer];
-    } else if (subDetID == 1) { // ML/OT
+    } else if (subDetID == 1  || subDetID == 2) { // ML/OT
       return constants::moduleMLOT::chip::nRows;
     }
     return 0;

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

@@ -107,7 +107,7 @@ const o2::trk::ChipSimResponse* Digitizer::getChipResponse(int chipID)
     return mChipSimRespVD;
   }
 
-  else if (mGeometry->getSubDetID(chipID) == 1) { /// ML/OT
+  else if (mGeometry->getSubDetID(chipID) == 1 || mGeometry->getSubDetID(chipID) == 2) { /// ML/OT
     return mChipSimRespMLOT;
   }
   return nullptr;
@@ -121,12 +121,12 @@ void Digitizer::process(const std::vector<Hit>* hits, int evID, int srcID, int l
   LOG(info) << " Digitizing " << mGeometry->getName() << " (ID: " << mGeometry->getDetID()
             << ") hits of event " << evID << " from source " << srcID
             << " at time " << mEventTime.getTimeNS() << " ROFrame = " << mNewROFrame
-            << " Min/Max ROFrames " << mROFrameMin << "/" << mROFrameMax;
+            << " Min/Max ROFrames " << mROFrameMin << "/" << mROFrameMax << " layer " << layer;
 
-  std::cout << "Printing segmentation info: " << std::endl;
-  SegmentationChip::Print();
+  //std::cout << "Printing segmentation info: " << std::endl;
+  //SegmentationChip::Print();
 
-  // // is there something to flush ?
+  // is there something to flush ?
   if (mNewROFrame > mROFrameMin) {
     fillOutputContainer(mNewROFrame - 1, layer); // flush out all frames preceding the new one
   }
@@ -257,11 +257,11 @@ void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, i
   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);
 
-  int layer = mGeometry->getLayer(chipID);
+  int layer = mGeometry->getLayer(chipID); // local layer nr for response
   int disk = mGeometry->getDisk(chipID);
 
   if (disk != -1) {
-    LOG(debug) << "Skipping disk " << disk;
+    LOG(debug) << "Skipping VD disk " << disk;
     return; // skipping hits on disks for the moment
   }
 
@@ -310,7 +310,7 @@ void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, i
 
   const auto& matrix = mGeometry->getMatrixL2G(hit.GetDetectorID());
   // matrix.print();
-
+ 
   /// transorm from the global detector coordinates to the local detector coordinates
   math_utils::Vector3D<float> xyzLocS(matrix ^ (hit.GetPosStart())); // start position in sensor frame
   math_utils::Vector3D<float> xyzLocE(matrix ^ (hit.GetPos()));      // end position in sensor frame
@@ -403,8 +403,8 @@ void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, i
   float cRowPix = 0.f, cColPix = 0.f; // local coordinate of the current pixel center
 
   const o2::trk::ChipSimResponse* resp = getChipResponse(chipID);
-  // std::cout << "Printing chip response:" << std::endl;
-  // resp->print();
+  //std::cout << "Printing chip response:" << std::endl;
+  //resp->print();
 
   // 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
@@ -463,7 +463,7 @@ void Digitizer::processHit(const o2::trk::Hit& hit, uint32_t& maxFr, int evID, i
       }
     }
   }
-
+  LOG(info) << "Response done; adding labels; making digits";
   // fire the pixels assuming Poisson(n_response_electrons)
   o2::MCCompLabel lbl(hit.GetTrackID(), evID, srcID, false);
   auto roFrameAbs = mNewROFrame + roFrameRel;

Steer/DigitizerWorkflow/src/TRKDigitizerSpec.cxx

@@ -279,7 +279,7 @@ class TRKDPLDigitizerTask : BaseDPLDigitizer
       const auto& dopt = o2::trk::DPLDigitizerParam<o2::detectors::DetID::TRK>::Instance();
       // pc.inputs().get<o2::trk::AlmiraParam*>("TRK_almiraparam");
       const auto& aopt = o2::trk::AlmiraParam::Instance();
-      mLayers = constants::VD::petal::nLayers + geom->getNumberOfLayersMLOT();
+      mLayers = constants::VD::petal::nLayers + geom->getNumberOfLayersMLOT() + geom->getNumberOfDisksMLOT();
       mDigits.resize(mLayers);
       mROFRecords.resize(mLayers);
       mROFRecordsAccum.resize(mLayers);