Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-12 23:30:40

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