diff --git a/SimFastTiming/FastTimingCommon/interface/BTLDeviceSim.h b/SimFastTiming/FastTimingCommon/interface/BTLDeviceSim.h index 475b6117d3f80..08aa9f125324a 100644 --- a/SimFastTiming/FastTimingCommon/interface/BTLDeviceSim.h +++ b/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_; }; diff --git a/SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h b/SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h index 2d07d14d22b04..ec8f0b1a032c0 100644 --- a/SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h +++ b/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_; diff --git a/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py b/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py index f437c447f6d3c..6b3a96ec68ba5 100644 --- a/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py +++ b/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] diff --git a/SimFastTiming/FastTimingCommon/src/BTLDeviceSim.cc b/SimFastTiming/FastTimingCommon/src/BTLDeviceSim.cc index ca6a7ca585505..28191c0d30dcc 100644 --- a/SimFastTiming/FastTimingCommon/src/BTLDeviceSim.cc +++ b/SimFastTiming/FastTimingCommon/src/BTLDeviceSim.cc @@ -20,8 +20,7 @@ BTLDeviceSim::BTLDeviceSim(const edm::ParameterSet& pset, edm::ConsumesCollector bxTime_(pset.getParameter("bxTime")), LightYield_(pset.getParameter("LightYield")), LightCollEff_(pset.getParameter("LightCollectionEff")), - LightCollSlopeR_(pset.getParameter("LightCollectionSlopeR")), - LightCollSlopeL_(pset.getParameter("LightCollectionSlopeL")), + LightCollSlope_(pset.getParameter("LightCollectionSlope")), PDE_(pset.getParameter("PhotonDetectionEff")), LCEpositionSlope_(pset.getParameter("LCEpositionSlope")) {} @@ -88,8 +87,8 @@ void BTLDeviceSim::getHitsResponse(const std::vector(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; diff --git a/SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc b/SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc index 54dae0e5fbc1c..ac40c6829dea2 100644 --- a/SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc +++ b/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("LightYield") * pset.getParameter("LightCollectionEff") * + pset.getParameter("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 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; +}