ALICE3-TRK: Detector.cxx modified to retrieve the correct Detector ID following the chosen indexing scheme

Author: atriolo

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

@@ -80,14 +80,20 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache
   int getDisk(int index) const;
 
   /// This routine computes the chip index number from the subDetID, petal, disk, layer, stave /// TODO: retrieve also from chip when chips will be available
-  /// in substave
   /// \param int subDetID The subdetector ID, 0 for VD, 1 for MLOT
   /// \param int petalcase The petal case number for VD, from 0 to 3
   /// \param int disk The disk number for VD, from 0 to 5
   /// \param int lay The layer number. Starting from 0 both for VD and MLOT
   /// \param int stave The stave number for MLOT. Starting from 0
   int getChipIndex(int subDetID, int petalcase, int disk, int lay, int stave) const;
 
+  /// This routine computes the chip index number from the subDetID, volume, layer, stave /// TODO: retrieve also from chip when chips will be available
+  /// \param int subDetID The subdetector ID, 0 for VD, 1 for MLOT
+  /// \param int volume is needed only with the current configuration for VD where each single element is a volume. // TODO: when the geometry naming scheme will be changed, change this method
+  /// \param int lay The layer number for the MLOT. In the current configuration for VD this is not needed. // TODO: when the geometry naming scheme will be changed, change this method
+  /// \param int stave The stave number in each layer for MLOT. Starting from 0.
+  int getChipIndex(int subDetID, int volume, int lay, int stave) const;
+
   /// This routine computes subDetID, petal, disk, layer, stave given the chip index number /// TODO: copute also from chip when chips will be available
   /// \param int index The chip index number, starting from 0
   /// \param int subDetID The subdetector ID, 0 for VD, 1 for MLOT

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

@@ -155,7 +155,7 @@ int GeometryTGeo::getLayer(int index) const
     while (index > mLastChipIndex[lay]) {
       lay++;
     }
-    return lay - mNumberOfPetalsVD; /// numeration of MLOT layesrs  starting from 1
+    return lay - mNumberOfPetalsVD; /// numeration of MLOT layesrs  starting from 0
   }
   return -1; /// -1 if not found
 }
@@ -207,6 +207,18 @@ int GeometryTGeo::getChipIndex(int subDetID, int petalcase, int disk, int lay, i
   return -1; // not found
 }
 
+//__________________________________________________________________________
+int GeometryTGeo::getChipIndex(int subDetID, int volume, int lay, int stave) const
+{
+  if (subDetID == 0) { // VD
+    return volume;     /// In the current configuration for VD, each volume is the sensor element = chip. // TODO: when the geometry naming scheme will be changed, change this method
+
+  } else if (subDetID == 1) { // MLOT
+    return getFirstChipIndex(lay, -1, subDetID) + stave;
+  }
+  return -1; // not found
+}
+
 //__________________________________________________________________________
 bool GeometryTGeo::getChipID(int index, int& subDetID, int& petalcase, int& disk, int& lay, int& stave) const
 {
@@ -291,7 +303,7 @@ TGeoHMatrix* GeometryTGeo::extractMatrixSensor(int index) const
   auto path = getMatrixPath(index);
 
   static TGeoHMatrix matTmp;
-  gGeoManager->PushPath();
+  // gGeoManager->PushPath(); // Preserve the modeler state.
 
   // if (!gGeoManager->cd(path.Data())) {
   //   gGeoManager->PopPath();
@@ -440,7 +452,6 @@ int GeometryTGeo::extractNumberOfActivePartsVD() const
   if (vdV == nullptr) {
     LOG(fatal) << getName() << " volume " << getTRKVolPattern() << " is not in the geometry";
   }
-  LOG(info) << "Volume name: " << getTRKVolPattern();
 
   // Loop on all TRKV nodes, count Layer volumes by checking names
   TObjArray* nodes = vdV->GetNodes();
@@ -470,7 +481,6 @@ int GeometryTGeo::extractNumberOfDisksVD() const
   if (vdV == nullptr) {
     LOG(fatal) << getName() << " volume " << getTRKVolPattern() << " is not in the geometry";
   }
-  LOG(info) << "Volume name: " << getTRKVolPattern();
 
   // Loop on all TRKV nodes, count Layer volumes by checking names
   TObjArray* nodes = vdV->GetNodes();
@@ -500,7 +510,6 @@ int GeometryTGeo::extractNumberOfPetalsVD() const
   if (vdV == nullptr) {
     LOG(fatal) << getName() << " volume " << getTRKVolPattern() << " is not in the geometry";
   }
-  LOG(info) << "Volume name: " << getTRKVolPattern();
 
   // Loop on all TRKV nodes, count Layer volumes by checking names
   TObjArray* nodes = vdV->GetNodes();
@@ -530,7 +539,6 @@ int GeometryTGeo::extractNumberOfLayersVD() const
   if (vdV == nullptr) {
     LOG(fatal) << getName() << " volume " << getTRKVolPattern() << " is not in the geometry";
   }
-  LOG(info) << "Volume name: " << getTRKVolPattern();
 
   // Loop on all TRKV nodes, count Layer volumes by checking names
   TObjArray* nodes = vdV->GetNodes();
@@ -560,7 +568,6 @@ int GeometryTGeo::extractNumberOfChipsPerPetalVD() const
   if (vdV == nullptr) {
     LOG(fatal) << getName() << " volume " << getTRKVolPattern() << " is not in the geometry";
   }
-  LOG(info) << "Volume name: " << getTRKVolPattern();
 
   // Loop on all TRKV nodes, count Layer volumes by checking names
   TObjArray* nodes = vdV->GetNodes();

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

@@ -31,6 +31,9 @@ namespace trk
 class Detector : public o2::base::DetImpl<Detector>
 {
  public:
+  static constexpr Int_t mNumberOfVolumes = 44;   /// hardcoded for the current geometry = 8 MLOT layers + 36 volumes in the VD. TODO: automatize or change according to the current geometry
+  static constexpr Int_t mNumberOfVolumesVD = 36; /// hardcoded for the current geometry = 36 volumes in the VD. TODO: automatize or change according to the current geometry
+
   Detector(bool active);
   Detector();
   ~Detector();
@@ -96,6 +99,18 @@ class Detector : public o2::base::DetImpl<Detector>
 
   void defineSensitiveVolumes();
 
+ protected:
+  std::vector<int> mSensorID;       //! layer identifiers
+  std::vector<TString> mSensorName; //! layer names
+
+ public:
+  static constexpr Int_t sNumberVDPetalCases = 4;          //! Number of VD petals
+  int getNumberOfLayers() const { return mLayers.size(); } //! Number of TRK layers
+  int getNumberOfLayersVD() const { return mPetalCases[0].mPetalLayers.size(); }
+  int getNumberOfDisksVD() const { return mPetalCases[0].mPetalDisks.size(); }
+
+  void Print(FairVolume* vol, int volume, int subDetID, int layer, int stave, int halfstave, int chipID) const;
+
   template <typename Det>
   friend class o2::base::DetImpl;
   ClassDefOverride(Detector, 1);

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

@@ -53,6 +53,28 @@ Detector::Detector(bool active)
     configServices();
   }
 
+  mSensorName.resize(mNumberOfVolumes); // hardcoded. TODO: change size when a different naming scheme for VD is in place. Ideally could be 4 petals + 8 layers = 12
+  int VDvolume = 0;
+  for (int i = 0; i < 4; i++) { /// VD
+    for (int j = 0; j < 3; j++) {
+      mSensorName[VDvolume].Form("%s%d_%s%d_%s%d", GeometryTGeo::getTRKPetalPattern(), i, GeometryTGeo::getTRKPetalLayerPattern(), j, GeometryTGeo::getTRKSensorPattern(), j);
+      VDvolume++;
+    }
+    for (int j = 0; j < 6; j++) {
+      mSensorName[VDvolume].Form("%s%d_%s%d_%s%d", GeometryTGeo::getTRKPetalPattern(), i, GeometryTGeo::getTRKPetalDiskPattern(), j, GeometryTGeo::getTRKSensorPattern(), j);
+      VDvolume++;
+    }
+  }
+
+  for (int i = 0; i < 8; i++) { /// MLOT
+    mSensorName[VDvolume].Form("%s%d", GeometryTGeo::getTRKSensorPattern(), i);
+    VDvolume++;
+  }
+
+  for (auto vd : mSensorName) {
+    std::cout << "Volume name: " << vd << std::endl;
+  }
+
   LOGP(info, "Summary of TRK configuration:");
   for (auto& layer : mLayers) {
     LOGP(info, "Layer: {} name: {} r: {} cm | z: {} cm | thickness: {} cm", layer.getNumber(), layer.getName(), layer.getInnerRadius(), layer.getZ(), layer.getChipThickness());
@@ -241,7 +263,7 @@ void Detector::createGeometry()
   mPetalCases.clear();
   // Add petal cases (the sensitive layers inside the petal cases get constructed here too)
   auto& trkPars = TRKBaseParam::Instance();
-  for (Int_t petalCaseNumber = 0; petalCaseNumber < 4; ++petalCaseNumber) {
+  for (Int_t petalCaseNumber = 0; petalCaseNumber < sNumberVDPetalCases; ++petalCaseNumber) {
     mPetalCases.emplace_back(petalCaseNumber, vTRK, trkPars.irisOpen);
     mServices.excavateFromVacuum(mPetalCases[petalCaseNumber].getFullName());
   }
@@ -253,6 +275,12 @@ void Detector::InitializeO2Detector()
   LOG(info) << "Initialize TRK O2Detector";
   mGeometryTGeo = GeometryTGeo::Instance();
   defineSensitiveVolumes();
+
+  mSensorID.resize(mNumberOfVolumes); // hardcoded. TODO: change size when a different namingh scheme for VD is in place. Ideally could be 4 petals + 8 layers = 12
+  for (int i = 0; i < mNumberOfVolumes; i++) {
+    mSensorID[i] = gMC ? TVirtualMC::GetMC()->VolId(mSensorName[i]) : 0; // Volume ID from the Geant geometry
+    LOGP(info, "{}: mSensorID={}", i, mSensorID[i]);
+  }
 }
 
 void Detector::defineSensitiveVolumes()
@@ -264,7 +292,7 @@ void Detector::defineSensitiveVolumes()
   LOGP(info, "Adding TRK Sensitive Volumes");
 
   // Add petal case sensitive volumes
-  for (int petalCase = 0; petalCase < 4; ++petalCase) {
+  for (int petalCase = 0; petalCase < sNumberVDPetalCases; ++petalCase) {
     // Petal layers
     for (int petalLayer = 0; petalLayer < mPetalCases[petalCase].mPetalLayers.size(); ++petalLayer) {
       volumeName = mPetalCases[petalCase].mPetalLayers[petalLayer].getSensorName();
@@ -338,17 +366,36 @@ bool Detector::ProcessHits(FairVolume* vol)
     return false;
   }
 
-  int lay = vol->getVolumeId();
+  int subDetID = -1;
+  int layer = -1;
+  int volume = 0;
+  int stave = -1;
   int volID = vol->getMCid();
 
+  bool notSens = false;
+  while ((volume < mNumberOfVolumes) && (notSens = (volID != mSensorID[volume]))) {
+    ++volume; /// there are 44 volumes, 36 for the VD (1 for each sensing element) and 8 for the MLOT (1 for each layer)
+  }
+
+  if (notSens) {
+    return kFALSE; // RS: can this happen? This method must be called for sensors only?
+  }
+
+  if (volume < mNumberOfVolumesVD) {
+    subDetID = 0; // VD. For the moment each "chip" is a volume./// TODO: change this logic once the naming scheme is changed
+  } else {
+    subDetID = 1; // MLOT
+    layer = volume - mNumberOfVolumesVD;
+  }
+
   // Is it needed to keep a track reference when the outer ITS volume is encountered?
   auto stack = (o2::data::Stack*)fMC->GetStack();
   // if (fMC->IsTrackExiting() && (lay == 0 || lay == mLayers.size() - 1)) {
   if (fMC->IsTrackExiting() && InsideFirstOrLastLayer(vol->GetName())) {
     // Keep the track refs for the innermost and outermost layers only
     o2::TrackReference tr(*fMC, GetDetId());
     tr.setTrackID(stack->GetCurrentTrackNumber());
-    tr.setUserId(lay);
+    tr.setUserId(volume);
     stack->addTrackReference(tr);
   }
   bool startHit = false, stopHit = false;
@@ -398,13 +445,18 @@ bool Detector::ProcessHits(FairVolume* vol)
     TLorentzVector positionStop;
     fMC->TrackPosition(positionStop);
     // Retrieve the indices with the volume path
-    int stave(0), halfstave(0), chipinmodule(0), module;
-    fMC->CurrentVolOffID(1, chipinmodule);
-    fMC->CurrentVolOffID(2, module);
-    fMC->CurrentVolOffID(3, halfstave);
-    fMC->CurrentVolOffID(4, stave);
+    // int subDetID, int petalcase, int disk, int lay, int stave
+    int petal(0), ll(0), stave(0), halfstave(0);
+    if (subDetID == 1) {
+      fMC->CurrentVolOffID(1, halfstave);
+      fMC->CurrentVolOffID(2, stave);
+    } /// if VD, for the moment the volume is the "chipID" so no need to retrieve other elments
 
-    Hit* p = addHit(stack->GetCurrentTrackNumber(), lay, mTrackData.mPositionStart.Vect(), positionStop.Vect(),
+    int chipID = mGeometryTGeo->getChipIndex(subDetID, volume, layer, stave);
+
+    Print(vol, volume, subDetID, layer, stave, halfstave, chipID);
+
+    Hit* p = addHit(stack->GetCurrentTrackNumber(), chipID, mTrackData.mPositionStart.Vect(), positionStop.Vect(),
                     mTrackData.mMomentumStart.Vect(), mTrackData.mMomentumStart.E(), positionStop.T(),
                     mTrackData.mEnergyLoss, mTrackData.mTrkStatusStart, status);
     // p->SetTotalEnergy(vmc->Etot());
@@ -424,6 +476,23 @@ o2::itsmft::Hit* Detector::addHit(int trackID, int detID, const TVector3& startP
   mHits->emplace_back(trackID, detID, startPos, endPos, startMom, startE, endTime, eLoss, startStatus, endStatus);
   return &(mHits->back());
 }
+
+void Detector::Print(FairVolume* vol, int volume, int subDetID, int layer, int stave, int halfstave, int chipID) const
+{
+  int currentVol(0);
+  LOG(INFO) << "Current volume name: " << fMC->CurrentVolName() << " and ID " << fMC->CurrentVolID(currentVol);
+  LOG(INFO) << "volume: " << volume << "/" << mNumberOfVolumes - 1;
+  if (subDetID == 1) {
+    LOG(INFO) << "off volume name 1 " << fMC->CurrentVolOffName(1) << "  halfstave: " << halfstave;
+    LOG(INFO) << "off volume name 2  " << fMC->CurrentVolOffName(2) << "  stave: " << stave;
+    LOG(INFO) << "SubDetector ID: " << subDetID << "  Layer: " << layer << "  staveinLayer: " << stave << "  Chip ID: " << chipID;
+
+  } else {
+    LOG(INFO) << "SubDetector ID: " << subDetID << "  Chip ID: " << chipID;
+  }
+  LOG(INFO);
+}
+
 } // namespace trk
 } // namespace o2