PWGJE/EMCal: Add option to shift cluster time in MC (#8153)

Co-authored-by: jokonig <jokonig@cern.ch>

Author: jokonig

PWGJE/TableProducer/emcalCorrectionTask.cxx

@@ -87,6 +87,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; 0 = off; 1 = const shift; 2 = eta-dependent shift"};
 
   // Require EMCAL cells (CALO type 1)
   Filter emccellfilter = aod::calo::caloType == selectedCellType;
@@ -108,6 +109,9 @@ struct EmcalCorrectionTask {
   // QA
   o2::framework::HistogramRegistry mHistManager{"EMCALCorrectionTaskQAHistograms"};
 
+  // EMCal geometry
+  o2::emcal::Geometry* geometry;
+
   void init(InitContext const&)
   {
     LOG(debug) << "Start init!";
@@ -124,7 +128,7 @@ struct EmcalCorrectionTask {
       mCcdbManager->get<TGeoManager>("GLO/Config/Geometry");
     }
     LOG(debug) << "After load geometry!";
-    o2::emcal::Geometry* geometry = o2::emcal::Geometry::GetInstanceFromRunNumber(223409);
+    geometry = o2::emcal::Geometry::GetInstanceFromRunNumber(223409);
     if (!geometry) {
       LOG(error) << "Failure accessing geometry";
     }
@@ -175,6 +179,7 @@ struct EmcalCorrectionTask {
     using o2HistType = o2::framework::HistType;
     using o2Axis = o2::framework::AxisSpec;
     o2Axis energyAxis{200, 0., 100., "E (GeV)"},
+      timeAxis{300, -100, 200., "t (ns)"},
       etaAxis{160, -0.8, 0.8, "#eta"},
       phiAxis{72, 0, 2 * 3.14159, "phi"};
     mHistManager.add("hCellE", "hCellE", o2HistType::kTH1F, {energyAxis});
@@ -184,6 +189,7 @@ struct EmcalCorrectionTask {
     mHistManager.add("hCellRowCol", "hCellRowCol;Column;Row", o2HistType::kTH2D, {{97, 0, 97}, {600, 0, 600}});
     mHistManager.add("hClusterE", "hClusterE", o2HistType::kTH1F, {energyAxis});
     mHistManager.add("hClusterEtaPhi", "hClusterEtaPhi", o2HistType::kTH2F, {etaAxis, phiAxis});
+    mHistManager.add("hClusterTime", "hClusterTime", o2HistType::kTH1F, {timeAxis});
     mHistManager.add("hGlobalTrackEtaPhi", "hGlobalTrackEtaPhi", o2HistType::kTH2F, {etaAxis, phiAxis});
     mHistManager.add("hGlobalTrackMult", "hGlobalTrackMult", o2HistType::kTH1D, {{200, -0.5, 199.5, "N_{trk}"}});
     mHistManager.add("hCollisionType", "hCollisionType;;#it{count}", o2HistType::kTH1D, {{3, -0.5, 2.5}});
@@ -258,7 +264,7 @@ struct EmcalCorrectionTask {
         }
         cellsBC.emplace_back(cell.cellNumber(),
                              amplitude,
-                             cell.time(),
+                             cell.time() + getCellTimeShift(cell.cellNumber()),
                              o2::emcal::intToChannelType(cell.cellType()));
         cellIndicesBC.emplace_back(cell.globalIndex());
       }
@@ -377,7 +383,7 @@ struct EmcalCorrectionTask {
         }
         cellsBC.emplace_back(cell.cellNumber(),
                              amplitude,
-                             cell.time(),
+                             cell.time() + getCellTimeShift(cell.cellNumber()),
                              o2::emcal::intToChannelType(cell.cellType()));
         cellIndicesBC.emplace_back(cell.globalIndex());
         cellLabels.emplace_back(cell.mcParticleIds(), cell.amplitudeA());
@@ -479,7 +485,7 @@ struct EmcalCorrectionTask {
       for (auto& cell : cellsInBC) {
         cellsBC.emplace_back(cell.cellNumber(),
                              cell.amplitude(),
-                             cell.time(),
+                             cell.time() + getCellTimeShift(cell.cellNumber()),
                              o2::emcal::intToChannelType(cell.cellType()));
         cellIndicesBC.emplace_back(cell.globalIndex());
       }
@@ -618,6 +624,7 @@ struct EmcalCorrectionTask {
       } // end of cells of cluser loop
       // fill histograms
       mHistManager.fill(HIST("hClusterE"), cluster.E());
+      mHistManager.fill(HIST("hClusterTime"), cluster.getClusterTime());
       mHistManager.fill(HIST("hClusterEtaPhi"), pos.Eta(), TVector2::Phi_0_2pi(pos.Phi()));
       if (IndexMapPair && trackGlobalIndex) {
         for (unsigned int iTrack = 0; iTrack < std::get<0>(*IndexMapPair)[iCluster].size(); iTrack++) {
@@ -670,7 +677,7 @@ struct EmcalCorrectionTask {
         clustercellsambiguous(clustersAmbiguous.lastIndex(),
                               cellIndicesBC[cellindex]);
       } // end of cells of cluster loop
-    }   // end of cluster loop
+    } // end of cluster loop
   }
 
   template <typename Collision>
@@ -782,6 +789,30 @@ struct EmcalCorrectionTask {
       return 1.f;
     }
   }
+
+  // Apply shift of the cell time
+  // This has to be done to shift the cell time in MC (which is not calibrated to 0 due to the flight time of the particles to the EMCal surface (~15ns))
+  float getCellTimeShift(const int16_t cellID)
+  {
+    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
+      return 0.f;
+    }
+  }
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)