[EMCAL-518] Handling of TOF in time sampling

- Add templates for different TOF to the EMCAL surface
- Limit accepted TOF to the EMCAL readout window
- Resolution 1 ns for 1.5 mus (1500 templates per phase)
- Do not do SDigitization when sampling for different
  Time-of-flight, as delays in the same tower can be
  sufficiently different

Author: mfasDa

Detectors/EMCAL/simulation/include/EMCALSimulation/Digitizer.h

@@ -94,21 +94,28 @@ class Digitizer : public TObject
   const std::vector<o2::emcal::TriggerRecord>& getTriggerRecords() const { return mDigits.getTriggerRecords(); }
   const o2::dataformats::MCTruthContainer<o2::emcal::MCLabel>& getMCLabels() const { return mDigits.getMCLabels(); }
 
+  static constexpr int getTOFSamplingBins() { return EMC_TOF_BINS; }
+
  private:
-  static constexpr int EMC_PHASES = 4; ///< Number of phases
-  short mEventTimeOffset = 0;          ///< event time difference from trigger time (in number of bins)
-  short mPhase = 0;                    ///< event phase
-  UInt_t mROFrameMin = 0;              ///< lowest RO frame of current digits
-  UInt_t mROFrameMax = 0;              ///< highest RO frame of current digits
-  bool mSmearEnergy = true;            ///< do time and energy smearing
-  bool mSimulateTimeResponse = true;   ///< simulate time response
-  const SimParam* mSimParam = nullptr; ///< SimParam object
+  using TimeSampleContainer = std::array<double, constants::EMCAL_MAXTIMEBINS>;
+  static constexpr int EMC_PHASES = 4;                                                  ///< Number of phases
+  static constexpr int EMC_TOF_BINS = 1500;                                             ///< Number of bins in TOF sampling of the time response
+  static constexpr double EMC_TOF_MIN = 0;                                              ///< Min TOF
+  static constexpr double EMC_TOF_MAX = 1500.;                                          ///< Max TOF
+  static constexpr double EMC_TOF_BINWITH = (EMC_TOF_MAX - EMC_TOF_MIN) / EMC_TOF_BINS; ///< Number time samples simulated
+  short mEventTimeOffset = 0;                                                           ///< event time difference from trigger time (in number of bins)
+  short mPhase = 0;                                                                     ///< event phase
+  UInt_t mROFrameMin = 0;                                                               ///< lowest RO frame of current digits
+  UInt_t mROFrameMax = 0;                                                               ///< highest RO frame of current digits
+  bool mSmearEnergy = true;                                                             ///< do time and energy smearing
+  bool mSimulateTimeResponse = true;                                                    ///< simulate time response
+  const SimParam* mSimParam = nullptr;                                                  ///< SimParam object
 
   std::vector<Digit> mTempDigitVector;     ///< temporary digit storage
   o2::emcal::DigitsWriteoutBuffer mDigits; ///< used to sort digits and labels by tower
 
   TRandom3* mRandomGenerator = nullptr; ///< random number generator
-  std::array<std::array<double, constants::EMCAL_MAXTIMEBINS>, EMC_PHASES>
+  std::array<std::array<TimeSampleContainer, EMC_TOF_BINS>, EMC_PHASES>
     mAmplitudeInTimeBins; ///< template of the sampled time response function: amplitude of signal for each time bin (per phase)
 
   std::unique_ptr<o2::utils::TreeStreamRedirector> mDebugStream = nullptr;

Detectors/EMCAL/simulation/src/Digitizer.cxx

@@ -64,9 +64,12 @@ void Digitizer::init()
       // parameter 1: Handling phase + delay
       // phase: 25 ns * phase index (-4)
       // delay: Average signal delay
-      RawResponse.SetParameter(1, 0.25 * phase + mSimParam->getSignalDelay() / constants::EMCAL_TIMESAMPLE);
-      for (int sample = 0; sample < constants::EMCAL_MAXTIMEBINS; sample++) {
-        mAmplitudeInTimeBins[phase][sample] = RawResponse.Eval(sample);
+      for (int itofbin = 0; itofbin < EMC_TOF_BINS; itofbin++) {
+        double tofbincenter = itofbin * EMC_TOF_BINWITH + 0.5 * EMC_TOF_BINWITH;
+        RawResponse.SetParameter(1, 0.25 * phase + (tofbincenter + mSimParam->getSignalDelay()) / constants::EMCAL_TIMESAMPLE);
+        for (int sample = 0; sample < constants::EMCAL_MAXTIMEBINS; sample++) {
+          mAmplitudeInTimeBins[phase][itofbin][sample] = RawResponse.Eval(sample);
+        }
       }
     }
   } else {
@@ -158,15 +161,29 @@ void Digitizer::sampleSDigit(const Digit& sDigit)
 
   Double_t energies[15];
   if (mSimulateTimeResponse) {
-    for (int sample = 0; sample < mAmplitudeInTimeBins[mPhase].size(); sample++) {
+    if (sDigit.getTimeStamp() + mSimParam->getSignalDelay() + mPhase * 25 > EMC_TOF_MAX) {
+      // Digit time larger than sampling window, will not be sampled
+      // For time response simulation take also signal delay and phase into account
+      return;
+    }
+    int tofbin = static_cast<int>(sDigit.getTimeStamp() / EMC_TOF_BINWITH);
+    if (tofbin >= EMC_TOF_BINS) {
+      tofbin = EMC_TOF_BINS - 1;
+    }
+    for (int sample = 0; sample < mAmplitudeInTimeBins[mPhase][tofbin].size(); sample++) {
 
-      double val = energy * (mAmplitudeInTimeBins[mPhase][sample]);
+      double val = energy * (mAmplitudeInTimeBins[mPhase][tofbin][sample]);
       energies[sample] = val;
       double digitTime = mEventTimeOffset * constants::EMCAL_TIMESAMPLE;
       Digit digit(tower, val, digitTime);
       mTempDigitVector.push_back(digit);
     }
   } else {
+    if (sDigit.getTimeStamp() > EMC_TOF_MAX) {
+      // Digit time larger than sampling window, will not be sampled
+      // In non-sampled mode only apply the max. time window
+      return;
+    }
     Digit digit(tower, energy, smearTime(sDigit.getTimeStamp(), energy));
     mTempDigitVector.push_back(digit);
   }

Detectors/EMCAL/workflow/include/EMCALWorkflow/EMCALDigitizerSpec.h

@@ -70,6 +70,7 @@ class DigitizerSpec final : public o2::base::BaseDPLDigitizer, public o2::framew
  private:
   Bool_t mFinished = false;                   ///< Flag for digitization finished
   bool mIsConfigured = false;                 ///< Initialization status of the digitizer
+  bool mRunSDitizer = false;                  ///< Run SDigitization
   Digitizer mDigitizer;                       ///< Digitizer object
   o2::emcal::SDigitizer mSumDigitizer;        ///< Summed digitizer
   std::shared_ptr<CalibLoader> mCalibHandler; ///< Handler of calibration objects

Detectors/EMCAL/workflow/src/EMCALDigitizerSpec.cxx

@@ -119,7 +119,18 @@ void DigitizerSpec::run(framework::ProcessingContext& ctx)
 
       LOG(info) << "For collision " << collID << " eventID " << part.entryID << " found " << mHits.size() << " hits ";
 
-      std::vector<o2::emcal::LabeledDigit> summedDigits = mSumDigitizer.process(mHits);
+      std::vector<o2::emcal::LabeledDigit> summedDigits;
+      if (mRunSDitizer) {
+        summedDigits = mSumDigitizer.process(mHits);
+      } else {
+        for (auto& hit : mHits) {
+          o2::emcal::MCLabel digitlabel(hit.GetTrackID(), part.entryID, part.sourceID, false, 1.);
+          if (hit.GetEnergyLoss() < __DBL_EPSILON__) {
+            digitlabel.setAmplitudeFraction(0);
+          }
+          summedDigits.emplace_back(hit.GetDetectorID(), hit.GetEnergyLoss(), hit.GetTime(), digitlabel);
+        }
+      }
 
       // call actual digitization procedure
       mDigitizer.process(summedDigits);