IOTOF: add in-pixel efficiency (#15786)

* add pixel efficiency map to iotof digitization

* add path to the available file in a comment

* change binning of the histograms

* clang-format

* revert size of the histograms

Author: GiorgioAlbertoLucia

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

@@ -30,6 +30,8 @@ struct DPLDigitizerParam : public o2::conf::ConfigurableParamHelper<DPLDigitizer
   float timeResolution = 0.020f;          ///< time resolution sigma in ns (20 ps default)
   float tdcBin = 0.010f;                  ///< TDC time bin (10 ps default)
   float efficiency = 0.98f;               ///< detection efficiency
+  std::string efficiencyFilePath{};       ///< optional efficiency map file path.
+                                          ///< The efficiency map is currently available at /alice/cern.ch/user/g/glucia/ALICE3/IOTOF/pixelEfficiency/PixelEfficiencyMap_TH2.root. FIXME to be removed once switch to CCDBFetcher
   int chargeThreshold = 100;              ///< charge threshold in Nelectrons
   int minChargeToAccount = 7;             ///< minimum charge contribution to account
   int nSimSteps = 1;                      ///< number of steps in response simulation

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

@@ -23,8 +23,9 @@
 #include <deque>
 #include <memory>
 
-#include "Rtypes.h"  // for Digitizer::Class
-#include "TObject.h" // for TObject
+#include <TH2D.h>
+#include <Rtypes.h>  // for Digitizer::Class
+#include <TObject.h> // for TObject
 
 #include "ITSMFTSimulation/Hit.h"
 #include "DataFormatsIOTOF/Digit.h"
@@ -90,8 +91,13 @@ class Digitizer : public TObject
   /// Convert energy loss to charge
   int energyToCharge(float energyLoss) const;
 
+  /// Load the efficiency map from a file
+  void loadEfficiencyMap(const std::string& filePath);
+
   /// Check if the hit passes efficiency cut
-  bool isEfficient() const;
+  /// \param x Detector local coordinate x in cm with respect to the center of the sensitive volume.
+  /// \param z Detector local coordinate z in cm with respect to the center of the sensitive volume.
+  bool isEfficient(const float x, const float z) const;
 
   std::vector<o2::iotof::McLabelRef>* getExtraLabelBuffer(uint32_t roFrame)
   {
@@ -108,8 +114,10 @@ class Digitizer : public TObject
   }
 
   static constexpr float sec2ns = 1e9f; ///< seconds to nanoseconds conversion
+  static constexpr float cm2um = 1e4f;  ///< centimeters to micrometers conversion
 
   const o2::iotof::GeometryTGeo* mGeometry = nullptr; ///< IOTOF geometry
+  TH2D* mEfficiencyMap = nullptr;                     ///< Efficiency map for the detector
 
   std::vector<o2::iotof::Chip> mChips;                                               //! Chips in the detector, indexed by chip ID
   std::deque<std::unique_ptr<std::vector<o2::iotof::McLabelRef>>> mExtraLabelBuffer; //! buffer for multiple mc labels to the same pixel

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

@@ -20,7 +20,11 @@
 #include "IOTOFSimulation/DPLDigitizerParam.h"
 #include "DetectorsRaw/HBFUtils.h"
 
+#include <TCollection.h>
+#include <TFile.h>
+#include <TKey.h>
 #include <TRandom.h>
+
 #include <vector>
 #include <iostream>
 #include <numeric>
@@ -51,6 +55,9 @@ void Digitizer::init()
   }
 
   const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance();
+  if (!digitizerParams.efficiencyFilePath.empty()) {
+    loadEfficiencyMap(digitizerParams.efficiencyFilePath);
+  }
 
   LOG(info) << "Initializing IOTOF digitizer";
   LOG(info) << "  Time resolution: " << digitizerParams.timeResolution * 1e3 << " ps";
@@ -95,12 +102,6 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
 {
   // Process a single hit and create a digit if it passes all cuts
 
-  // Apply efficiency cut
-  if (!isEfficient()) {
-    LOG(debug) << "Hit rejected by efficiency cut";
-    return;
-  }
-
   // Get detector element ID
   const int chipID = hit.GetDetectorID();
   auto& chip = mChips[chipID];
@@ -109,6 +110,26 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
     return;
   }
 
+  // middle position of the hit in the sensor frame
+  const auto& matrix = mGeometry->getMatrixL2G(chipID);
+  auto xyzPositionStart = matrix ^ hit.GetPosStart();
+  auto xyzPositionEnd = matrix ^ hit.GetPos();
+  const auto xMid = 0.5f * (xyzPositionStart.X() + xyzPositionEnd.X());
+  const auto zMid = 0.5f * (xyzPositionStart.Z() + xyzPositionEnd.Z());
+  // move this to the local pixel coordinates for the efficiency map
+  int row, col;
+  float xPixelCenter, zPixelCenter;
+  if (!sSegmentation->localToDetector(xMid, zMid, row, col, mGeometry->getIOTOFLayer(chipID))) {
+    LOG(debug) << "Hit rejected because position (" << xMid << ", " << zMid << ") is outside the active area of chip " << chipID;
+    return; // hit is outside the active area
+  }
+  sSegmentation->detectorToLocalUnchecked(row, col, xPixelCenter, zPixelCenter, mGeometry->getIOTOFLayer(chipID));
+
+  if (!isEfficient(xMid - xPixelCenter, zMid - zPixelCenter)) {
+    LOG(debug) << "Hit rejected by efficiency cut";
+    return;
+  }
+
   // Convert energy loss to charge (number of electrons)
   float energyLoss = hit.GetEnergyLoss(); // in GeV
   int charge = energyToCharge(energyLoss);
@@ -126,7 +147,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
   double hitTime = hit.GetTime() * sec2ns;                // convert to ns
   double eventTimeInBC = mEventTime.getTimeOffsetWrtBC(); // event time wrt bc
   double hitTimeWrtBC = hitTime + eventTimeInBC;          // hit time wrt bc
-  double smearedTime = smearTime(hitTimeWrtBC);           // apply detector resolution
+  double smearedTime = smearTime(hitTimeWrtBC);
 
   if (chipID < 0 || chipID >= mGeometry->getSize() || mGeometry->getSize() < 1) {
     LOG(debug) << "Invalid detector ID: " << chipID << ", geometry size: " << mGeometry->getSize();
@@ -166,8 +187,8 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
 
 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 auto& matrix = mGeometry->getMatrixL2G(chipID);
   const int subdetectorID = mGeometry->getIOTOFLayer(chipID);
 
   auto xyzPositionStart(matrix ^ (hit.GetPosStart())); // start position in sensor frame
@@ -277,10 +298,45 @@ int Digitizer::energyToCharge(float energyLoss) const
 }
 
 //_______________________________________________________________________
-bool Digitizer::isEfficient() const
+void Digitizer::loadEfficiencyMap(const std::string& filePath)
+{
+  // Load the efficiency map from a file
+  TFile* file = TFile::Open(filePath.c_str());
+  if (!file || !file->IsOpen()) {
+    LOG(error) << "Failed to open efficiency map file: " << filePath;
+    return;
+  }
+
+  auto* rawMap = dynamic_cast<TH2D*>(file->Get("hEfficiencyMap"));
+  if (!rawMap) {
+    LOG(error) << "Failed to retrieve efficiency map from file: " << filePath;
+    LOG(error) << "Available keys in the file:";
+    TIter next(file->GetListOfKeys());
+    TKey* key;
+    while ((key = dynamic_cast<TKey*>(next()))) {
+      LOG(error) << "  " << key->GetName() << " (" << key->GetClassName() << ")";
+    }
+    file->Close();
+    return;
+  }
+  mEfficiencyMap = dynamic_cast<TH2D*>(rawMap->Clone("mEfficiencyMap"));
+  mEfficiencyMap->SetDirectory(nullptr); // Detach from file to avoid deletion when file is closed
+
+  file->Close();
+}
+
+//_______________________________________________________________________
+bool Digitizer::isEfficient(const float x, const float z) const
 {
   // Apply efficiency cut using random number
   const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance();
+  if (mEfficiencyMap) {
+    // int bin = mEfficiencyMap->FindBin(x * o2::iotof::Digitizer::cm2um, z * o2::iotof::Digitizer::cm2um);
+    int bin = mEfficiencyMap->FindBin(x * o2::iotof::Digitizer::cm2um, z * o2::iotof::Digitizer::cm2um);
+    float efficiency = mEfficiencyMap->GetBinContent(bin);
+    LOG(debug) << "Efficiency map check: x=" << x * o2::iotof::Digitizer::cm2um << ", z=" << z * o2::iotof::Digitizer::cm2um << ", bin=" << bin << ", efficiency=" << efficiency;
+    return gRandom->Uniform() < efficiency;
+  }
   return gRandom->Uniform() < digitizerParams.efficiency;
 }