Merge pull request #41131 from fabiocos/fc-btlresmodel131X

MTD digitization: Add calculation of per SiPM time resolution to BtlElectronicsSim

SimFastTiming/FastTimingCommon/interface/BTLDeviceSim.h

@@ -43,8 +43,7 @@ class BTLDeviceSim {
   const float LightYield_;
   const float LightCollEff_;
 
-  const float LightCollSlopeR_;
-  const float LightCollSlopeL_;
+  const float LightCollSlope_;
   const float PDE_;
   const float LCEpositionSlope_;
 };

SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h

@@ -42,6 +42,12 @@ class BTLElectronicsSim {
 private:
   float sigma2_pe(const float& Q, const float& R) const;
 
+  float sigma_stochastic(const float& npe) const;
+
+  float sigma2_DCR(const float& npe) const;
+
+  float sigma2_electronics(const float npe) const;
+
   const bool debug_;
 
   const float bxTime_;

SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py

@@ -1,5 +1,12 @@
 import FWCore.ParameterSet.Config as cms
 
+_common_BTLparameters = cms.PSet(
+    bxTime                   = cms.double(25),      # [ns]
+    LightYield               = cms.double(40000.),  # [photons/MeV]
+    LightCollectionEff       = cms.double(0.25),
+    PhotonDetectionEff       = cms.double(0.20),
+)
+
 _barrel_MTDDigitizer = cms.PSet(
     digitizerName     = cms.string("BTLDigitizer"),
     inputSimHits      = cms.InputTag("g4SimHits:FastTimerHitsBarrel"),
@@ -9,16 +16,12 @@
     premixStage1MinCharge = cms.double(1e-4),
     premixStage1MaxCharge = cms.double(1e6),
     DeviceSimulation = cms.PSet(
-        bxTime                   = cms.double(25),      # [ns]
-        LightYield               = cms.double(40000.),  # [photons/MeV]
-        LightCollectionEff       = cms.double(0.25),
-        LightCollectionSlopeR    = cms.double(0.075),   # [ns/cm]
-        LightCollectionSlopeL    = cms.double(0.075),   # [ns/cm]
-        PhotonDetectionEff       = cms.double(0.20),
+        _common_BTLparameters,
+        LightCollectionSlope     = cms.double(0.075),   # [ns/cm]
         LCEpositionSlope         = cms.double(0.071),   # [1/cm] LCE variation vs longitudinal position shift
         ),
     ElectronicsSimulation = cms.PSet(
-        bxTime                     = cms.double(25),    # [ns]
+        _common_BTLparameters,
         TestBeamMIPTimeRes         = cms.double(4.293), # This is given by 0.048[ns]*sqrt(8000.), in order to
                                                         # rescale the time resolution of 1 MIP = 8000 p.e.
         ScintillatorRiseTime       = cms.double(1.1),   # [ns]

SimFastTiming/FastTimingCommon/src/BTLDeviceSim.cc

@@ -20,8 +20,7 @@ BTLDeviceSim::BTLDeviceSim(const edm::ParameterSet& pset, edm::ConsumesCollector
       bxTime_(pset.getParameter<double>("bxTime")),
       LightYield_(pset.getParameter<double>("LightYield")),
       LightCollEff_(pset.getParameter<double>("LightCollectionEff")),
-      LightCollSlopeR_(pset.getParameter<double>("LightCollectionSlopeR")),
-      LightCollSlopeL_(pset.getParameter<double>("LightCollectionSlopeL")),
+      LightCollSlope_(pset.getParameter<double>("LightCollectionSlope")),
       PDE_(pset.getParameter<double>("PhotonDetectionEff")),
       LCEpositionSlope_(pset.getParameter<double>("LCEpositionSlope")) {}
 
@@ -88,8 +87,8 @@ void BTLDeviceSim::getHitsResponse(const std::vector<std::tuple<int, uint32_t, f
     double distR = 0.5 * topo.pitch().first - convertMmToCm(hit.localPosition().x());
     double distL = 0.5 * topo.pitch().first + convertMmToCm(hit.localPosition().x());
 
-    double tR = std::get<2>(hitRef) + LightCollSlopeR_ * distR;
-    double tL = std::get<2>(hitRef) + LightCollSlopeL_ * distL;
+    double tR = std::get<2>(hitRef) + LightCollSlope_ * distR;
+    double tL = std::get<2>(hitRef) + LightCollSlope_ * distL;
 
     // --- Accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
     const int iBXR = std::floor(tR / bxTime_) + mtd_digitizer::kInTimeBX;

SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc

@@ -46,7 +46,26 @@ BTLElectronicsSim::BTLElectronicsSim(const edm::ParameterSet& pset, edm::Consume
       SPTR2_(SinglePhotonTimeResolution_ * SinglePhotonTimeResolution_),
       DCRxRiseTime_(DarkCountRate_ * ScintillatorRiseTime_),
       SigmaElectronicNoise2_(SigmaElectronicNoise_ * SigmaElectronicNoise_),
-      SigmaClock2_(SigmaClock_ * SigmaClock_) {}
+      SigmaClock2_(SigmaClock_ * SigmaClock_) {
+#ifdef EDM_ML_DEBUG
+  float lightOutput = 4.2f * pset.getParameter<double>("LightYield") * pset.getParameter<double>("LightCollectionEff") *
+                      pset.getParameter<double>("PhotonDetectionEff");  // average npe for 4.2 MeV
+  float s1 = sigma_stochastic(lightOutput);
+  float s2 = std::sqrt(sigma2_DCR(lightOutput));
+  float s3 = std::sqrt(sigma2_electronics(lightOutput));
+  float s4 = SinglePhotonTimeResolution_ / std::sqrt(TimeThreshold1_);
+  float s5 = SigmaClock_;
+  LogDebug("BTLElectronicsSim") << " BTL time resolution model (ns, 1 SiPM), for an average light output of "
+                                << std::fixed << std::setw(14) << lightOutput << " N_pe:"
+                                << "\n sigma stochastic   = " << std::setw(14) << s1
+                                << "\n sigma DCR          = " << std::setw(14) << s2
+                                << "\n sigma electronics  = " << std::setw(14) << s3
+                                << "\n sigma sptr         = " << std::setw(14) << s4
+                                << "\n sigma clock        = " << std::setw(14) << s5 << "\n ---------------------"
+                                << "\n sigma total        = " << std::setw(14)
+                                << std::sqrt(s1 * s1 + s2 * s2 + s3 * s3 + s4 * s4 + s5 * s5);
+#endif
+}
 
 void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
                             BTLDigiCollection& output,
@@ -62,8 +81,8 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
     toa2.fill(0.f);
     for (size_t iside = 0; iside < 2; iside++) {
       // --- Fluctuate the total number of photo-electrons
-      float Npe = CLHEP::RandPoissonQ::shoot(hre, (it->second).hit_info[2 * iside][iBX]);
-      if (Npe < EnergyThreshold_)
+      float npe = CLHEP::RandPoissonQ::shoot(hre, (it->second).hit_info[2 * iside][iBX]);
+      if (npe < EnergyThreshold_)
         continue;
 
       // --- Get the time of arrival and add a channel time offset
@@ -77,7 +96,7 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
       // --- Calculate and add the time walk: the time of arrival is read in correspondence
       //                                      with two thresholds on the signal pulse
       std::array<float, 3> times =
-          btlPulseShape_.timeAtThr(Npe / ReferencePulseNpe_, TimeThreshold1_ * Npe_to_V_, TimeThreshold2_ * Npe_to_V_);
+          btlPulseShape_.timeAtThr(npe / ReferencePulseNpe_, TimeThreshold1_ * Npe_to_V_, TimeThreshold2_ * Npe_to_V_);
 
       // --- If the pulse amplitude is smaller than TimeThreshold2, the trigger does not fire
       if (times[1] == 0.)
@@ -99,7 +118,7 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
         }  // ibx loop
 
         if (rate_oot > 0.) {
-          float sigma_oot = sqrt(rate_oot * ScintillatorRiseTime_) * ScintillatorDecayTime_ / Npe;
+          float sigma_oot = sqrt(rate_oot * ScintillatorRiseTime_) * ScintillatorDecayTime_ / npe;
           float smearing_oot = CLHEP::RandGaussQ::shoot(hre, 0., sigma_oot);
           finalToA1 += smearing_oot;
           finalToA2 += smearing_oot;
@@ -110,7 +129,7 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
       if (testBeamMIPTimeRes_ > 0.) {
         // In this case the time resolution is parametrized from the testbeam.
         // The same parameterization is used for both thresholds.
-        float sigma = testBeamMIPTimeRes_ / sqrt(Npe);
+        float sigma = sigma_stochastic(npe);
         float smearing_stat_thr1 = CLHEP::RandGaussQ::shoot(hre, 0., sigma);
         float smearing_stat_thr2 = CLHEP::RandGaussQ::shoot(hre, 0., sigma);
 
@@ -121,21 +140,17 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
         // In this case the time resolution is taken from the literature.
         // The fluctuations due to the first TimeThreshold1_ p.e. are common to both times
         float smearing_stat_thr1 =
-            CLHEP::RandGaussQ::shoot(hre, 0., ScintillatorDecayTime_ * sqrt(sigma2_pe(TimeThreshold1_, Npe)));
+            CLHEP::RandGaussQ::shoot(hre, 0., ScintillatorDecayTime_ * sqrt(sigma2_pe(TimeThreshold1_, npe)));
         float smearing_stat_thr2 = CLHEP::RandGaussQ::shoot(
-            hre, 0., ScintillatorDecayTime_ * sqrt(sigma2_pe(TimeThreshold2_ - TimeThreshold1_, Npe)));
+            hre, 0., ScintillatorDecayTime_ * sqrt(sigma2_pe(TimeThreshold2_ - TimeThreshold1_, npe)));
         finalToA1 += smearing_stat_thr1;
         finalToA2 += smearing_stat_thr1 + smearing_stat_thr2;
       }
 
       // --- Add in quadrature the uncertainties due to the SiPM timing resolution, the SiPM DCR,
       //     the electronic noise and the clock distribution:
-      float slew2 = ScintillatorDecayTime2_ / Npe / Npe;
-
-      float sigma2_tot_thr1 =
-          SPTR2_ / TimeThreshold1_ + (DCRxRiseTime_ + SigmaElectronicNoise2_) * slew2 + SigmaClock2_;
-      float sigma2_tot_thr2 =
-          SPTR2_ / TimeThreshold2_ + (DCRxRiseTime_ + SigmaElectronicNoise2_) * slew2 + SigmaClock2_;
+      float sigma2_tot_thr1 = SPTR2_ / TimeThreshold1_ + sigma2_DCR(npe) + sigma2_electronics(npe) + SigmaClock2_;
+      float sigma2_tot_thr2 = SPTR2_ / TimeThreshold2_ + sigma2_DCR(npe) + sigma2_electronics(npe) + SigmaClock2_;
 
       // --- Smear the arrival times using the correlated uncertainties:
       float smearing_thr1_uncorr = CLHEP::RandGaussQ::shoot(hre, 0., sqrt(sigma2_tot_thr1));
@@ -145,7 +160,7 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
       finalToA2 += sinPhi_ * smearing_thr1_uncorr + cosPhi_ * smearing_thr2_uncorr;
 
       //Smear the energy according to TOFHIR energy branch measured resolution
-      float tofhir_ampnoise_relsigma = ROOT::Math::Chebyshev4(Npe,
+      float tofhir_ampnoise_relsigma = ROOT::Math::Chebyshev4(npe,
                                                               sigmaRelTOFHIRenergy_[0],
                                                               sigmaRelTOFHIRenergy_[1],
                                                               sigmaRelTOFHIRenergy_[2],
@@ -235,3 +250,23 @@ float BTLElectronicsSim::sigma2_pe(const float& Q, const float& R) const {
 
   return sigma2;
 }
+
+float BTLElectronicsSim::sigma_stochastic(const float& npe) const { return testBeamMIPTimeRes_ / std::sqrt(npe); }
+
+float BTLElectronicsSim::sigma2_DCR(const float& npe) const {
+  // trick to safely switch off the electronics contribution for resolution studies
+
+  if (DCRxRiseTime_ == 0.) {
+    return 0.;
+  }
+  return DCRxRiseTime_ * ScintillatorDecayTime2_ / npe / npe;
+}
+
+float BTLElectronicsSim::sigma2_electronics(const float npe) const {
+  // trick to safely switch off the electronics contribution for resolution studies
+
+  if (SigmaElectronicNoise2_ == 0.) {
+    return 0.;
+  }
+  return SigmaElectronicNoise2_ * ScintillatorDecayTime2_ / npe / npe;
+}