[PWGJE] [PWGJE]Add cell time smearing in MC (#8697)

Author: nstrangm

PWGJE/TableProducer/emcalCorrectionTask.cxx

@@ -22,6 +22,7 @@
 #include <string>
 #include <tuple>
 #include <vector>
+#include <random>
 
 #include "CCDB/BasicCCDBManager.h"
 #include "Framework/runDataProcessing.h"
@@ -90,7 +91,7 @@ struct EmcalCorrectionTask {
   Configurable<float> exoticCellInCrossMinAmplitude{"exoticCellInCrossMinAmplitude", 0.1, "Minimum energy of cells in cross, if lower not considered in cross"};
   Configurable<bool> useWeightExotic{"useWeightExotic", false, "States if weights should be used for exotic cell cut"};
   Configurable<bool> isMC{"isMC", false, "States if run over MC"};
-  Configurable<int> applyCellTimeShift{"applyCellTimeShift", 0, "apply shift to the cell time for data and MC; For data: 0 = off; non-zero = log function extracted from data - For MC: 0 = off; 1 = const shift; 2 = eta-dependent shift"};
+  Configurable<bool> applyCellTimeCorrection{"applyCellTimeCorrection", true, "apply a correction to the cell time for data and MC: Shift both average cell times to 0 and smear MC time distribution to fit data better"};
 
   // Require EMCAL cells (CALO type 1)
   Filter emccellfilter = aod::calo::caloType == selectedCellType;
@@ -112,6 +113,11 @@ struct EmcalCorrectionTask {
   // QA
   o2::framework::HistogramRegistry mHistManager{"EMCALCorrectionTaskQAHistograms"};
 
+  // Random number generator to draw cell time smearing for MC
+  std::random_device rd{};
+  std::mt19937_64 rdgen{rd()};
+  std::normal_distribution<> normalgaus{0, 1}; // mean = 0, stddev = 1 (apply amplitude of smearing after drawing random for performance reasons)
+
   // EMCal geometry
   o2::emcal::Geometry* geometry;
 
@@ -189,8 +195,9 @@ struct EmcalCorrectionTask {
     mHistManager.add("hCellE", "hCellE", o2HistType::kTH1F, {energyAxis});
     mHistManager.add("hCellTowerID", "hCellTowerID", o2HistType::kTH1D, {{20000, 0, 20000}});
     mHistManager.add("hCellEtaPhi", "hCellEtaPhi", o2HistType::kTH2F, {etaAxis, phiAxis});
+    mHistManager.add("hCellTimeEnergy", "hCellTime", o2HistType::kTH2F, {{300, -30, 30}, {200, 0., 20.}});
     // NOTE: Reversed column and row because it's more natural for presentation.
-    mHistManager.add("hCellRowCol", "hCellRowCol;Column;Row", o2HistType::kTH2D, {{97, 0, 97}, {600, 0, 600}});
+    mHistManager.add("hCellRowCol", "hCellRowCol;Column;Row", o2HistType::kTH2D, {{96, -0.5, 95.5}, {208, -0.5, 207.5}});
     mHistManager.add("hClusterE", "hClusterE", o2HistType::kTH1F, {energyAxis});
     mHistManager.add("hClusterNLM", "hClusterNLM", o2HistType::kTH1F, {nlmAxis});
     mHistManager.add("hClusterEtaPhi", "hClusterEtaPhi", o2HistType::kTH2F, {etaAxis, phiAxis});
@@ -770,6 +777,7 @@ struct EmcalCorrectionTask {
     // For convenience, use the clusterizer stored geometry to get the eta-phi
     for (auto& cell : cellsBC) {
       mHistManager.fill(HIST("hCellE"), cell.getEnergy());
+      mHistManager.fill(HIST("hCellTimeEnergy"), cell.getTimeStamp(), cell.getEnergy());
       mHistManager.fill(HIST("hCellTowerID"), cell.getTower());
       auto res = mClusterizers.at(0)->getGeometry()->EtaPhiFromIndex(cell.getTower());
       mHistManager.fill(HIST("hCellEtaPhi"), std::get<0>(res), TVector2::Phi_0_2pi(std::get<1>(res)));
@@ -801,33 +809,29 @@ struct EmcalCorrectionTask {
   // In data this is done to correct for the time walk effect
   float getCellTimeShift(const int16_t cellID, const float cellEnergy)
   {
+    if (!applyCellTimeCorrection) {
+      return 0.f;
+    }
+    float timeshift = 0.f;
+    float timesmear = 0.f;
     if (isMC) {
-      if (applyCellTimeShift == 1) { // constant shift
-        LOG(debug) << "shift the cell time by 15ns";
-        return -15.f;                       // roughly calculated by assuming particles travel with v=c (photons) and EMCal is 4.4m away from vertex
-      } else if (applyCellTimeShift == 2) { // eta dependent shift ( as larger eta values are further away from collision point)
-        // Use distance between vertex and EMCal (at eta = 0) and distance on EMCal surface (cell size times column) to calculate distance to cell
-        // 0.2 is cell size in m (0.06) divided by the speed of light in m/ns (0.3)
-        // 47.5 is the "middle" of the EMCal (2*48 cells in one column)
-        float timeCol = 0.2f * (geometry->GlobalCol(cellID) - 47.5f); // calculate time to get to specific column
-        float time = -sqrt(215.f + timeCol * timeCol);                // 215 is 14.67ns^2 (time it takes to get the cell at eta = 0)
-        LOG(debug) << "shift the cell time by " << time << "  applyCellTimeShift " << applyCellTimeShift;
-        return time;
-      } else {
-        return 0.f;
-      }
-    } else { // data
-      if (applyCellTimeShift != 0) {
-        if (cellEnergy < 0.3) // Cells with tless than 300 MeV cannot be the leading cell in the cluster, so their time does not require precise calibration
-          return 0.f;
-        else if (cellEnergy < 4.)                                         // Low energy regime
-          return (0.57284 + 0.82194 * TMath::Log(1.30651 * cellEnergy));  // Parameters extracted from LHC22o (pp), but also usable for other periods
-        else                                                              // High energy regime
-          return (-0.05858 + 1.50593 * TMath::Log(0.97591 * cellEnergy)); // Parameters extracted from LHC22o (pp), but also usable for other periods
-      } else {                                                            // Dont apply cell time shift if applyCellTimeShift == 0
-        return 0.f;
-      }
+      // Shift the time to 0, as the TOF was simulated -> eta dependent shift (as larger eta values are further away from collision point)
+      // Use distance between vertex and EMCal (at eta = 0) and distance on EMCal surface (cell size times column) to calculate distance to cell
+      // 0.2 is cell size in m (0.06) divided by the speed of light in m/ns (0.3) - 47.5 is the "middle" of the EMCal (2*48 cells in one column)
+      float timeCol = 0.2f * (geometry->GlobalCol(cellID) - 47.5f); // calculate time to get to specific column
+      timeshift = -sqrt(215.f + timeCol * timeCol);                 // 215 is 14.67ns^2 (time it takes to get the cell at eta = 0)
+      // Also smear the time to account for the broader time resolution in data than in MC
+      timesmear = normalgaus(rdgen) * (1.6 + 9.5 * TMath::Exp(-3. * cellEnergy)); // Parameters extracted from LHC22o (pp), but also usable for other periods
+    } else {                                                                      // data
+      if (cellEnergy < 0.3)                                                       // Cells with tless than 300 MeV cannot be the leading cell in the cluster, so their time does not require precise calibration
+        timeshift = 0.;
+      else if (cellEnergy < 4.)                                            // Low energy regime
+        timeshift = 0.57284 + 0.82194 * TMath::Log(1.30651 * cellEnergy);  // Parameters extracted from LHC22o (pp), but also usable for other periods
+      else                                                                 // High energy regime
+        timeshift = -0.05858 + 1.50593 * TMath::Log(0.97591 * cellEnergy); // Parameters extracted from LHC22o (pp), but also usable for other periods
     }
+    LOG(debug) << "Shift the cell time by " << timeshift << " + " << timesmear << " ns";
+    return timeshift + timesmear;
   }
 };