FastTimingCommon/src/BTLElectronicsSim.cc

0001 #include "SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h"
0002
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "CLHEP/Random/RandPoissonQ.h"
0007 #include "CLHEP/Random/RandGaussQ.h"
0008
0009 #include "Math/ChebyshevPol.h"
0010
0011 using namespace mtd;
0012
0013 BTLElectronicsSim::BTLElectronicsSim(const edm::ParameterSet& pset, edm::ConsumesCollector iC)
0014     : debug_(pset.getUntrackedParameter<bool>("debug", false)),
0015       bxTime_(pset.getParameter<double>("bxTime")),
0016       testBeamMIPTimeRes_(pset.getParameter<double>("TestBeamMIPTimeRes")),
0017       ScintillatorRiseTime_(pset.getParameter<double>("ScintillatorRiseTime")),
0018       ScintillatorDecayTime_(pset.getParameter<double>("ScintillatorDecayTime")),
0019       ChannelTimeOffset_(pset.getParameter<double>("ChannelTimeOffset")),
0020       smearChannelTimeOffset_(pset.getParameter<double>("smearChannelTimeOffset")),
0021       EnergyThreshold_(pset.getParameter<double>("EnergyThreshold")),
0022       TimeThreshold1_(pset.getParameter<double>("TimeThreshold1")),
0023       TimeThreshold2_(pset.getParameter<double>("TimeThreshold2")),
0024       ReferencePulseNpe_(pset.getParameter<double>("ReferencePulseNpe")),
0025       SinglePhotonTimeResolution_(pset.getParameter<double>("SinglePhotonTimeResolution")),
0026       DarkCountRate_(pset.getParameter<double>("DarkCountRate")),
0027       SigmaElectronicNoise_(pset.getParameter<double>("SigmaElectronicNoise")),
0028       SigmaClock_(pset.getParameter<double>("SigmaClock")),
0029       smearTimeForOOTtails_(pset.getParameter<bool>("SmearTimeForOOTtails")),
0030       Npe_to_pC_(pset.getParameter<double>("Npe_to_pC")),
0031       Npe_to_V_(pset.getParameter<double>("Npe_to_V")),
0032       sigmaRelTOFHIRenergy_(pset.getParameter<std::vector<double>>("SigmaRelTOFHIRenergy")),
0033       adcNbits_(pset.getParameter<uint32_t>("adcNbits")),
0034       tdcNbits_(pset.getParameter<uint32_t>("tdcNbits")),
0035       adcSaturation_MIP_(pset.getParameter<double>("adcSaturation_MIP")),
0036       adcBitSaturation_(std::pow(2, adcNbits_) - 1),
0037       adcLSB_MIP_(adcSaturation_MIP_ / adcBitSaturation_),
0038       adcThreshold_MIP_(pset.getParameter<double>("adcThreshold_MIP")),
0039       toaLSB_ns_(pset.getParameter<double>("toaLSB_ns")),
0040       tdcBitSaturation_(std::pow(2, tdcNbits_) - 1),
0041       CorrCoeff_(pset.getParameter<double>("CorrelationCoefficient")),
0042       cosPhi_(0.5 * (sqrt(1. + CorrCoeff_) + sqrt(1. - CorrCoeff_))),
0043       sinPhi_(0.5 * CorrCoeff_ / cosPhi_),
0044       ScintillatorDecayTime2_(ScintillatorDecayTime_ * ScintillatorDecayTime_),
0045       ScintillatorDecayTimeInv_(1. / ScintillatorDecayTime_),
0046       SPTR2_(SinglePhotonTimeResolution_ * SinglePhotonTimeResolution_),
0047       DCRxRiseTime_(DarkCountRate_ * ScintillatorRiseTime_),
0048       SigmaElectronicNoise2_(SigmaElectronicNoise_ * SigmaElectronicNoise_),
0049       SigmaClock2_(SigmaClock_ * SigmaClock_) {
0050 #ifdef EDM_ML_DEBUG
0051   float lightOutput = 4.2f * pset.getParameter<double>("LightYield") * pset.getParameter<double>("LightCollectionEff") *
0052                       pset.getParameter<double>("PhotonDetectionEff");  // average npe for 4.2 MeV
0053   float s1 = sigma_stochastic(lightOutput);
0054   float s2 = std::sqrt(sigma2_DCR(lightOutput));
0055   float s3 = std::sqrt(sigma2_electronics(lightOutput));
0056   float s4 = SinglePhotonTimeResolution_ / std::sqrt(TimeThreshold1_);
0057   float s5 = SigmaClock_;
0058   LogDebug("BTLElectronicsSim") << " BTL time resolution model (ns, 1 SiPM), for an average light output of "
0059                                 << std::fixed << std::setw(14) << lightOutput << " N_pe:"
0060                                 << "\n sigma stochastic   = " << std::setw(14) << s1
0061                                 << "\n sigma DCR          = " << std::setw(14) << s2
0062                                 << "\n sigma electronics  = " << std::setw(14) << s3
0063                                 << "\n sigma sptr         = " << std::setw(14) << s4
0064                                 << "\n sigma clock        = " << std::setw(14) << s5 << "\n ---------------------"
0065                                 << "\n sigma total        = " << std::setw(14)
0066                                 << std::sqrt(s1 * s1 + s2 * s2 + s3 * s3 + s4 * s4 + s5 * s5);
0067 #endif
0068 }
0069
0070 void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
0071                             BTLDigiCollection& output,
0072                             CLHEP::HepRandomEngine* hre) const {
0073   MTDSimHitData chargeColl, toa1, toa2;
0074
0075   for (MTDSimHitDataAccumulator::const_iterator it = input.begin(); it != input.end(); it++) {
0076     // --- Digitize only the in-time bucket:
0077     const unsigned int iBX = mtd_digitizer::kInTimeBX;
0078
0079     chargeColl.fill(0.f);
0080     toa1.fill(0.f);
0081     toa2.fill(0.f);
0082     for (size_t iside = 0; iside < 2; iside++) {
0083       // --- Fluctuate the total number of photo-electrons
0084       float npe = CLHEP::RandPoissonQ::shoot(hre, (it->second).hit_info[2 * iside][iBX]);
0085       if (npe < EnergyThreshold_)
0086         continue;
0087
0088       // --- Get the time of arrival and add a channel time offset
0089       float finalToA1 = (it->second).hit_info[1 + 2 * iside][iBX] + ChannelTimeOffset_;
0090
0091       if (smearChannelTimeOffset_ > 0.) {
0092         float timeSmearing = CLHEP::RandGaussQ::shoot(hre, 0., smearChannelTimeOffset_);
0093         finalToA1 += timeSmearing;
0094       }
0095
0096       // --- Calculate and add the time walk: the time of arrival is read in correspondence
0097       //                                      with two thresholds on the signal pulse
0098       std::array<float, 3> times =
0099           btlPulseShape_.timeAtThr(npe / ReferencePulseNpe_, TimeThreshold1_ * Npe_to_V_, TimeThreshold2_ * Npe_to_V_);
0100
0101       // --- If the pulse amplitude is smaller than TimeThreshold2, the trigger does not fire
0102       if (times[1] == 0.)
0103         continue;
0104
0105       float finalToA2 = finalToA1 + times[1];
0106       finalToA1 += times[0];
0107
0108       // --- Estimate the time uncertainty due to photons from earlier OOT hits in the current BTL cell
0109       if (smearTimeForOOTtails_) {
0110         float rate_oot = 0.;
0111         // Loop on earlier OOT hits
0112         for (int ibx = 0; ibx < mtd_digitizer::kInTimeBX; ++ibx) {
0113           if ((it->second).hit_info[2 * iside][ibx] > 0.) {
0114             float hit_time = (it->second).hit_info[1 + 2 * iside][ibx] + bxTime_ * (ibx - mtd_digitizer::kInTimeBX);
0115             float npe_oot = CLHEP::RandPoissonQ::shoot(hre, (it->second).hit_info[2 * iside][ibx]);
0116             rate_oot += npe_oot * exp(hit_time * ScintillatorDecayTimeInv_) * ScintillatorDecayTimeInv_;
0117           }
0118         }  // ibx loop
0119
0120         if (rate_oot > 0.) {
0121           float sigma_oot = sqrt(rate_oot * ScintillatorRiseTime_) * ScintillatorDecayTime_ / npe;
0122           float smearing_oot = CLHEP::RandGaussQ::shoot(hre, 0., sigma_oot);
0123           finalToA1 += smearing_oot;
0124           finalToA2 += smearing_oot;
0125         }
0126       }  // if smearTimeForOOTtails_
0127
0128       // --- Uncertainty due to the fluctuations of the n-th photon arrival time:
0129       if (testBeamMIPTimeRes_ > 0.) {
0130         // In this case the time resolution is parametrized from the testbeam.
0131         // The same parameterization is used for both thresholds.
0132         float sigma = sigma_stochastic(npe);
0133         float smearing_stat_thr1 = CLHEP::RandGaussQ::shoot(hre, 0., sigma);
0134         float smearing_stat_thr2 = CLHEP::RandGaussQ::shoot(hre, 0., sigma);
0135
0136         finalToA1 += smearing_stat_thr1;
0137         finalToA2 += smearing_stat_thr2;
0138
0139       } else {
0140         // In this case the time resolution is taken from the literature.
0141         // The fluctuations due to the first TimeThreshold1_ p.e. are common to both times
0142         float smearing_stat_thr1 =
0143             CLHEP::RandGaussQ::shoot(hre, 0., ScintillatorDecayTime_ * sqrt(sigma2_pe(TimeThreshold1_, npe)));
0144         float smearing_stat_thr2 = CLHEP::RandGaussQ::shoot(
0145             hre, 0., ScintillatorDecayTime_ * sqrt(sigma2_pe(TimeThreshold2_ - TimeThreshold1_, npe)));
0146         finalToA1 += smearing_stat_thr1;
0147         finalToA2 += smearing_stat_thr1 + smearing_stat_thr2;
0148       }
0149
0150       // --- Add in quadrature the uncertainties due to the SiPM timing resolution, the SiPM DCR,
0151       //     the electronic noise and the clock distribution:
0152       float sigma2_tot_thr1 = SPTR2_ / TimeThreshold1_ + sigma2_DCR(npe) + sigma2_electronics(npe) + SigmaClock2_;
0153       float sigma2_tot_thr2 = SPTR2_ / TimeThreshold2_ + sigma2_DCR(npe) + sigma2_electronics(npe) + SigmaClock2_;
0154
0155       // --- Smear the arrival times using the correlated uncertainties:
0156       float smearing_thr1_uncorr = CLHEP::RandGaussQ::shoot(hre, 0., sqrt(sigma2_tot_thr1));
0157       float smearing_thr2_uncorr = CLHEP::RandGaussQ::shoot(hre, 0., sqrt(sigma2_tot_thr2));
0158
0159       finalToA1 += cosPhi_ * smearing_thr1_uncorr + sinPhi_ * smearing_thr2_uncorr;
0160       finalToA2 += sinPhi_ * smearing_thr1_uncorr + cosPhi_ * smearing_thr2_uncorr;
0161
0162       //Smear the energy according to TOFHIR energy branch measured resolution
0163       float tofhir_ampnoise_relsigma = ROOT::Math::Chebyshev4(npe,
0164                                                               sigmaRelTOFHIRenergy_[0],
0165                                                               sigmaRelTOFHIRenergy_[1],
0166                                                               sigmaRelTOFHIRenergy_[2],
0167                                                               sigmaRelTOFHIRenergy_[3],
0168                                                               sigmaRelTOFHIRenergy_[4]);
0169       float smearing_tofhir = CLHEP::RandGaussQ::shoot(hre, 0., tofhir_ampnoise_relsigma);
0170       // the amplitude resolution already includes the photostatistics fluctuation, use the original average deposit
0171       chargeColl[iside] = (it->second).hit_info[2 * iside][iBX] * Npe_to_pC_ *
0172                           (1. + smearing_tofhir);  // the p.e. number is here converted to pC
0173
0174       toa1[iside] = finalToA1;
0175       toa2[iside] = finalToA2;
0176
0177     }  // iside loop
0178
0179     //run the shaper to create a new data frame
0180     BTLDataFrame rawDataFrame(it->first.detid_);
0181     runTrivialShaper(rawDataFrame, chargeColl, toa1, toa2, it->first.row_, it->first.column_);
0182     updateOutput(output, rawDataFrame);
0183
0184   }  // MTDSimHitDataAccumulator loop
0185 }
0186
0187 void BTLElectronicsSim::runTrivialShaper(BTLDataFrame& dataFrame,
0188                                          const mtd::MTDSimHitData& chargeColl,
0189                                          const mtd::MTDSimHitData& toa1,
0190                                          const mtd::MTDSimHitData& toa2,
0191                                          const uint8_t row,
0192                                          const uint8_t col) const {
0193   bool debug = debug_;
0194 #ifdef EDM_ML_DEBUG
0195   for (int it = 0; it < (int)(chargeColl.size()); it++)
0196     debug |= (chargeColl[it] > adcThreshold_MIP_);
0197 #endif
0198
0199   if (debug)
0200     edm::LogVerbatim("BTLElectronicsSim") << "[runTrivialShaper]" << std::endl;
0201
0202   //set new ADCs
0203   for (int it = 0; it < (int)(chargeColl.size()); it++) {
0204     BTLSample newSample;
0205     newSample.set(false, false, 0, 0, 0, row, col);
0206
0207     //brute force saturation, maybe could to better with an exponential like saturation
0208     const uint32_t adc = std::min((uint32_t)std::floor(chargeColl[it] / adcLSB_MIP_), adcBitSaturation_);
0209     const uint32_t tdc_time1 = std::min((uint32_t)std::floor(toa1[it] / toaLSB_ns_), tdcBitSaturation_);
0210     const uint32_t tdc_time2 = std::min((uint32_t)std::floor(toa2[it] / toaLSB_ns_), tdcBitSaturation_);
0211
0212     newSample.set(
0213         chargeColl[it] > adcThreshold_MIP_, tdc_time1 == tdcBitSaturation_, tdc_time2, tdc_time1, adc, row, col);
0214     dataFrame.setSample(it, newSample);
0215
0216     if (debug)
0217       edm::LogVerbatim("BTLElectronicsSim") << adc << " (" << chargeColl[it] << "/" << adcLSB_MIP_ << ") ";
0218   }
0219
0220   if (debug) {
0221     std::ostringstream msg;
0222     dataFrame.print(msg);
0223     edm::LogVerbatim("BTLElectronicsSim") << msg.str() << std::endl;
0224   }
0225 }
0226
0227 void BTLElectronicsSim::updateOutput(BTLDigiCollection& coll, const BTLDataFrame& rawDataFrame) const {
0228   BTLDataFrame dataFrame(rawDataFrame.id());
0229   dataFrame.resize(dfSIZE);
0230   bool putInEvent(false);
0231   for (int it = 0; it < dfSIZE; ++it) {
0232     dataFrame.setSample(it, rawDataFrame[it]);
0233     if (it == 0)
0234       putInEvent = rawDataFrame[it].threshold();
0235   }
0236
0237   if (putInEvent) {
0238     coll.push_back(dataFrame);
0239   }
0240 }
0241
0242 float BTLElectronicsSim::sigma2_pe(const float& Q, const float& R) const {
0243   float OneOverR = 1. / R;
0244   float OneOverR2 = OneOverR * OneOverR;
0245
0246   // --- This is Eq. (17) from Nucl. Instr. Meth. A 564 (2006) 185
0247   float sigma2 = Q * OneOverR2 *
0248                  (1. + 2. * (Q + 1.) * OneOverR + (Q + 1.) * (6. * Q + 11) * OneOverR2 +
0249                   (Q + 1.) * (Q + 2.) * (2. * Q + 5.) * OneOverR2 * OneOverR);
0250
0251   return sigma2;
0252 }
0253
0254 float BTLElectronicsSim::sigma_stochastic(const float& npe) const { return testBeamMIPTimeRes_ / std::sqrt(npe); }
0255
0256 float BTLElectronicsSim::sigma2_DCR(const float& npe) const {
0257   // trick to safely switch off the electronics contribution for resolution studies
0258
0259   if (DCRxRiseTime_ == 0.) {
0260     return 0.;
0261   }
0262   return DCRxRiseTime_ * ScintillatorDecayTime2_ / npe / npe;
0263 }
0264
0265 float BTLElectronicsSim::sigma2_electronics(const float npe) const {
0266   // trick to safely switch off the electronics contribution for resolution studies
0267
0268   if (SigmaElectronicNoise2_ == 0.) {
0269     return 0.;
0270   }
0271   return SigmaElectronicNoise2_ * ScintillatorDecayTime2_ / npe / npe;
0272 }