[ALICE3] added segmentation to the particle propagation in TOF3 (#15591)

* adding stepping in digitization

* added segmentation to the digitization

* added non-segmented version, with configurable choice

* clang format

* removed redundant hit response in the middle, nsteps=1 is now default

* Charge threshold 0 by defult, using efficiency set by the digitizer

Changed chargeThreshold from 75 to 0.

* Please consider the following formatting changes

---------

Co-authored-by: ALICE Action Bot <alibuild@cern.ch>

Author: GiorgioAlbertoLucia

Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/GeometryTGeo.h

@@ -39,29 +39,29 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache
   static const char* getIOTOFVolPattern() { return sIOTOFVolumeName.c_str(); }
 
   // Inner TOF
-  const int getITOFNumberOfChips() { return mNumberOfChipsIOTOF[0]; }
+  const int getITOFNumberOfChips() const { return mNumberOfChipsIOTOF[0]; }
   static const char* getITOFLayerPattern() { return sITOFLayerName.c_str(); }
   static const char* getITOFStavePattern() { return sITOFStaveName.c_str(); }
   static const char* getITOFModulePattern() { return sITOFModuleName.c_str(); }
   static const char* getITOFChipPattern() { return sITOFChipName.c_str(); }
   static const char* getITOFSensorPattern() { return sITOFSensorName.c_str(); }
 
   // Outer TOF
-  const int getOTOFNumberOfChips() { return mNumberOfChipsIOTOF[1]; }
+  const int getOTOFNumberOfChips() const { return mNumberOfChipsIOTOF[1]; }
   static const char* getOTOFLayerPattern() { return sOTOFLayerName.c_str(); }
   static const char* getOTOFStavePattern() { return sOTOFStaveName.c_str(); }
   static const char* getOTOFModulePattern() { return sOTOFModuleName.c_str(); }
   static const char* getOTOFChipPattern() { return sOTOFChipName.c_str(); }
   static const char* getOTOFSensorPattern() { return sOTOFSensorName.c_str(); }
 
   // Forward TOF
-  const int getFTOFNumberOfChips() { return mNumberOfChipsFTOF; }
+  const int getFTOFNumberOfChips() const { return mNumberOfChipsFTOF; }
   static const char* getFTOFLayerPattern() { return sFTOFLayerName.c_str(); }
   static const char* getFTOFChipPattern() { return sFTOFChipName.c_str(); }
   static const char* getFTOFSensorPattern() { return sFTOFSensorName.c_str(); }
 
   // Backward TOF
-  const int getBTOFNumberOfChips() { return mNumberOfChipsBTOF; }
+  const int getBTOFNumberOfChips() const { return mNumberOfChipsBTOF; }
   static const char* getBTOFLayerPattern() { return sBTOFLayerName.c_str(); }
   static const char* getBTOFChipPattern() { return sBTOFChipName.c_str(); }
   static const char* getBTOFSensorPattern() { return sBTOFSensorName.c_str(); }

Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h

@@ -27,10 +27,11 @@ struct DPLDigitizerParam : public o2::conf::ConfigurableParamHelper<DPLDigitizer
   float noisePerPixel = DEFNoisePerPixel(); ///< ALPIDE Noise per channel
 
   double timeOffset = 0.;                 ///< time offset (in seconds!) to calculate ROFrame from hit time
-  int chargeThreshold = 75;               ///< charge threshold in Nelectrons
+  int chargeThreshold = 0;                ///< charge threshold in Nelectrons
   int minChargeToAccount = 7;             ///< minimum charge contribution to account
-  int nSimSteps = 475;                    ///< number of steps in response simulation
+  int nSimSteps = 1;                      ///< number of steps in response simulation
   float energyToNElectrons = 1. / 3.6e-9; // conversion of eloss to Nelectrons
+  int responseMatrixSize = 1;             ///< size of the response matrix (odd number)
 
   std::string noiseFilePath{}; ///< optional noise masks file path. FIXME to be removed once switch to CCDBFetcher
 

Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h

@@ -93,6 +93,8 @@ class Digitizer : public TObject
   void registerDigits(Chip& chip, uint32_t roFrame, double time, int nROF,
                       uint16_t row, uint16_t col, int nElectrons, o2::MCCompLabel& label);
 
+  void stepping(const o2::itsmft::Hit& hit, float**& respMatrix, int& rowStart, int& colStart, int& rowSpan, int& colSpan);
+
   /// Apply time smearing to simulate detector resolution
   double smearTime(double time) const;
 

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

@@ -17,6 +17,7 @@
 ///
 
 #include "IOTOFSimulation/Digitizer.h"
+#include "IOTOFSimulation/DPLDigitizerParam.h"
 #include "DetectorsRaw/HBFUtils.h"
 
 #include <TRandom.h>
@@ -30,7 +31,6 @@ namespace o2::iotof
 {
 
 o2::iotof::Segmentation* Digitizer::sSegmentation = nullptr;
-
 //_______________________________________________________________________
 void Digitizer::init()
 {
@@ -100,7 +100,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
   }
 
   // Get detector element ID
-  int chipID = hit.GetDetectorID();
+  const int chipID = hit.GetDetectorID();
   auto& chip = mChips[chipID];
   if (chip.isDisabled()) {
     LOG(debug) << "Hit rejected because chip " << chipID << " is disabled";
@@ -110,6 +110,8 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
   // Convert energy loss to charge (number of electrons)
   float energyLoss = hit.GetEnergyLoss(); // in GeV
   int charge = energyToCharge(energyLoss);
+  const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance();
+  int electronsPerStep = static_cast<int>(charge / digitizerParams.nSimSteps);
 
   // Apply charge threshold
   if (charge < mChargeThreshold) {
@@ -128,25 +130,126 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
     LOG(debug) << "Invalid detector ID: " << chipID << ", geometry size: " << mGeometry->getSize();
     return; // invalid detector ID
   }
+
+  // Create the digit with time information
+  o2::MCCompLabel label(hit.GetTrackID(), evID, srcID, false);
+  const int roFrameAbs = 0; // For now, we can set this to 0 or calculate based on time if needed
+  const int nROF = 1;       // For now, we can assume the signal is contained in one ROF, this can be extended to multiple ROFs based on the time
+
+  float** respMatrix = nullptr;
+  int rowStart = 0, colStart = 0, rowSpan = 0, colSpan = 0;
+  stepping(hit, respMatrix, rowStart, colStart, rowSpan, colSpan);
+
+  for (int irow = rowSpan; irow--;) {
+    uint16_t rowIS = irow + rowStart;
+    for (int icol = colSpan; icol--;) {
+      uint16_t colIS = icol + colStart;
+      float nEleResp = respMatrix[irow][icol];
+      if (!nEleResp) {
+        continue;
+      }
+      const int nElectronsSampled = gRandom->Poisson(electronsPerStep * nEleResp);
+      // Noise can be added here if needed
+
+      registerDigits(chip, roFrameAbs, smearedTime, nROF,
+                     static_cast<uint16_t>(rowIS), static_cast<uint16_t>(colIS), nElectronsSampled, label);
+    }
+  }
+
+  for (int irow = 0; irow < rowSpan; ++irow) {
+    delete[] respMatrix[irow];
+  }
+  delete[] respMatrix;
+}
+
+void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, int& rowStart, int& colStart, int& rowSpan, int& colSpan)
+{
   const auto& matrix = mGeometry->getMatrixL2G(hit.GetDetectorID());
+  const int chipID = hit.GetDetectorID();
+  const int subdetectorID = mGeometry->getIOTOFLayer(chipID);
+
+  auto xyzPositionStart(matrix ^ (hit.GetPosStart())); // start position in sensor frame
+  auto xyzPositionEnd(matrix ^ (hit.GetPos()));        // end position in sensor frame
+
+  const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance();
+  const auto stepVector = (xyzPositionEnd - xyzPositionStart) / digitizerParams.nSimSteps;
+  xyzPositionStart = xyzPositionStart + stepVector * 0.5f; // center the start position in the middle of the step
+  xyzPositionEnd = xyzPositionEnd - stepVector * 0.5f;     // center the end position in the middle of the step
+
+  rowStart = -1;
+  colStart = -1;
+  int rowEnd = -1, colEnd = -1, nSkip = 0, nSteps = digitizerParams.nSimSteps;
+  while (!sSegmentation->localToDetector(xyzPositionStart.X(), xyzPositionStart.Z(), rowStart, colStart, mGeometry->getIOTOFLayer(chipID))) {
+    if (++nSkip > digitizerParams.nSimSteps) { // additional check to add: should we exclude something?
+      LOG(debug) << "Hit position out of bounds for detector ID " << chipID;
+      return; // hit is outside the active area
+    }
+    xyzPositionStart += stepVector;
+  }
+
+  while (!sSegmentation->localToDetector(xyzPositionEnd.X(), xyzPositionEnd.Z(), rowEnd, colEnd, mGeometry->getIOTOFLayer(chipID))) {
+    if (++nSkip > digitizerParams.nSimSteps) { // additional check to add: should we exclude something?
+      LOG(debug) << "Hit position out of bounds for detector ID " << chipID;
+      return; // hit is outside the active area
+    }
+    xyzPositionEnd += stepVector;
+  }
 
-  math_utils::Vector3D<float> xyzPositionStart(matrix ^ (hit.GetPosStart())); // start position in sensor frame
-  // math_utils::Vector3D<float> xyzPositionEnd(matrix ^ (hit.GetPos()));      // end position in sensor frame
+  if (rowStart > rowEnd) {
+    std::swap(rowStart, rowEnd);
+  }
+  if (colStart > colEnd) {
+    std::swap(colStart, colEnd);
+  }
 
-  int row = 0; // Will be determined from start hit position
-  int col = 0; // Will be determined from start hit position
+  // Expand the range to take into account the effects of charge sharing
+  rowStart -= digitizerParams.responseMatrixSize / 2;
+  rowEnd += digitizerParams.responseMatrixSize / 2;
+  rowStart = std::max(rowStart, 0);
+  colStart = std::max(colStart, 0);
 
-  if (!sSegmentation->localToDetector(xyzPositionStart.X(), xyzPositionStart.Z(), row, col, mGeometry->getIOTOFLayer(chipID))) {
-    LOG(debug) << "Hit position out of bounds for detector ID " << chipID;
-    return; // hit is outside the active area
+  rowEnd = std::min(rowEnd, (subdetectorID == 0 ? sSegmentation->mITofSpecsConfig.NRows : sSegmentation->mOTofSpecsConfig.NRows) - 1);
+  colEnd = std::min(colEnd, (subdetectorID == 0 ? sSegmentation->mITofSpecsConfig.NCols : sSegmentation->mOTofSpecsConfig.NCols) - 1);
+  rowSpan = rowEnd - rowStart + 1;
+  colSpan = colEnd - colStart + 1;
+
+  respMatrix = new float*[rowSpan];
+  for (int i = 0; i < rowSpan; ++i) {
+    respMatrix[i] = new float[colSpan]();
   }
 
-  // Create the digit with time information
-  o2::MCCompLabel label(hit.GetTrackID(), evID, srcID, false);
-  const int roFrameAbs = 0;  // For now, we can set this to 0 or calculate based on time if needed
-  const int nROF = 1;        // For now, we can assume the signal is contained in one ROF, this can be extended to multiple ROFs based on the time
+  int rowPrev = -1, colPrev = -1, row = 0, col = 0;
+  if (!respMatrix || rowSpan <= 0 || colSpan <= 0) {
+    return;
+  }
+  if (nSkip) {
+    nSteps -= nSkip;
+  }
 
-  registerDigits(chip, roFrameAbs, smearedTime, nROF, static_cast<uint16_t>(row), static_cast<uint16_t>(col), charge, label);
+  auto& currentPosLocal = xyzPositionStart;
+  for (int iStep = nSteps; iStep--;) {
+    sSegmentation->localToDetector(currentPosLocal.X(), currentPosLocal.Z(), row, col, subdetectorID);
+    if (row != rowPrev || col != colPrev) {
+      rowPrev = row;
+      colPrev = col;
+    }
+
+    currentPosLocal += stepVector; // Move to the next step position
+
+    for (int irow = digitizerParams.responseMatrixSize; irow--;) {
+      int rowDest = row + irow - (digitizerParams.responseMatrixSize / 2) - rowStart; // destination row in the respMatrix
+      if (rowDest < 0 || rowDest >= rowSpan) {
+        continue;
+      }
+      for (int icol = digitizerParams.responseMatrixSize; icol--;) {
+        int colDest = col + icol - (digitizerParams.responseMatrixSize / 2) - colStart; // destination column in the respMatrix
+        if (colDest < 0 || colDest >= colSpan) {
+          continue;
+        }
+        respMatrix[rowDest][colDest] += 1.;
+      }
+    }
+  }
 }
 
 //_______________________________________________________________________
@@ -200,10 +303,9 @@ void Digitizer::fillOutputContainer()
     auto& chipDigits = chip.getDigits();
     for (const auto& [key, digit] : chipDigits) {
 
-      /// Charge threshold not implemented yet
-      /// if (digit.getCharge() < mChargeThreshold) {
-      ///   continue; // skip digits below threshold
-      /// }
+      if (digit.getCharge() < mChargeThreshold) {
+        continue; // skip digits below threshold
+      }
 
       int digitID = mDigits->size();
       mDigits->emplace_back(digit.getChipIndex(), digit.getRow(), digit.getColumn(), digit.getCharge(), digit.getTime());