Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-17 23:58:03

0001 
0002 #include "SimCalorimetry/HGCalSimProducers/interface/HGCFEElectronics.h"
0003 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0004 #include "FWCore/Utilities/interface/transform.h"
0005 
0006 #include "vdt/vdtMath.h"
0007 
0008 using namespace hgc_digi;
0009 
0010 //
0011 template <class DFr>
0012 HGCFEElectronics<DFr>::HGCFEElectronics(const edm::ParameterSet& ps)
0013     : fwVersion_{ps.getParameter<uint32_t>("fwVersion")},
0014       adcPulse_{},
0015       pulseAvgT_{},
0016       tdcForToAOnset_fC_{},
0017       adcSaturation_fC_{-1.0},
0018       adcLSB_fC_{},
0019       tdcLSB_fC_{},
0020       tdcSaturation_fC_{-1.0},
0021       adcThreshold_fC_{},
0022       tdcOnset_fC_{},
0023       toaLSB_ns_{},
0024       tdcResolutionInNs_{1e-9},  // set time resolution very small by default
0025       targetMIPvalue_ADC_{},
0026       jitterNoise_ns_{},
0027       jitterConstant_ns_{},
0028       eventTimeOffset_ns_{{0.02, 0.02, 0.02}},
0029       noise_fC_{},
0030       toaMode_(WEIGHTEDBYE) {
0031   edm::LogVerbatim("HGCFE") << "[HGCFEElectronics] running with version " << fwVersion_ << std::endl;
0032   if (ps.exists("adcPulse")) {
0033     auto temp = ps.getParameter<std::vector<double> >("adcPulse");
0034     for (unsigned i = 0; i < temp.size(); ++i) {
0035       adcPulse_[i] = (float)temp[i];
0036     }
0037     // normalize adc pulse
0038     for (unsigned i = 0; i < adcPulse_.size(); ++i) {
0039       adcPulse_[i] = adcPulse_[i] / adcPulse_[2];
0040     }
0041     temp = ps.getParameter<std::vector<double> >("pulseAvgT");
0042     for (unsigned i = 0; i < temp.size(); ++i) {
0043       pulseAvgT_[i] = (float)temp[i];
0044     }
0045   }
0046   if (ps.exists("adcNbits")) {
0047     uint32_t adcNbits = ps.getParameter<uint32_t>("adcNbits");
0048     adcSaturation_fC_ = ps.getParameter<double>("adcSaturation_fC");
0049     adcLSB_fC_ = adcSaturation_fC_ / pow(2., adcNbits);
0050     edm::LogVerbatim("HGCFE") << "[HGCFEElectronics] " << adcNbits << " bit ADC defined"
0051                               << " with LSB=" << adcLSB_fC_ << " saturation to occur @ " << adcSaturation_fC_
0052                               << std::endl;
0053   }
0054 
0055   if (ps.exists("tdcNbits")) {
0056     tdcNbits_ = ps.getParameter<uint32_t>("tdcNbits");
0057     setTDCfsc(ps.getParameter<double>("tdcSaturation_fC"));
0058     edm::LogVerbatim("HGCFE") << "[HGCFEElectronics] " << tdcNbits_ << " bit TDC defined with LSB=" << tdcLSB_fC_
0059                               << " saturation to occur @ " << tdcSaturation_fC_
0060                               << " (NB lowered by 1 part in a million)" << std::endl;
0061   }
0062   if (ps.exists("targetMIPvalue_ADC"))
0063     targetMIPvalue_ADC_ = ps.getParameter<uint32_t>("targetMIPvalue_ADC");
0064   if (ps.exists("adcThreshold_fC"))
0065     adcThreshold_fC_ = ps.getParameter<double>("adcThreshold_fC");
0066   if (ps.exists("tdcOnset_fC"))
0067     tdcOnset_fC_ = ps.getParameter<double>("tdcOnset_fC");
0068   if (ps.exists("tdcForToAOnset_fC")) {
0069     auto temp = ps.getParameter<std::vector<double> >("tdcForToAOnset_fC");
0070     if (temp.size() == tdcForToAOnset_fC_.size()) {
0071       std::copy_n(temp.begin(), temp.size(), tdcForToAOnset_fC_.begin());
0072     } else {
0073       throw cms::Exception("BadConfiguration") << " HGCFEElectronics wrong size for ToA thresholds ";
0074     }
0075   }
0076   if (ps.exists("toaLSB_ns"))
0077     toaLSB_ns_ = ps.getParameter<double>("toaLSB_ns");
0078   if (ps.exists("tdcChargeDrainParameterisation")) {
0079     for (auto val : ps.getParameter<std::vector<double> >("tdcChargeDrainParameterisation")) {
0080       tdcChargeDrainParameterisation_.push_back((float)val);
0081     }
0082   }
0083   if (ps.exists("tdcResolutionInPs"))
0084     tdcResolutionInNs_ = ps.getParameter<double>("tdcResolutionInPs") * 1e-3;  // convert to ns
0085   if (ps.exists("toaMode"))
0086     toaMode_ = ps.getParameter<uint32_t>("toaMode");
0087 
0088   if (ps.exists("jitterNoise_ns")) {
0089     auto temp = ps.getParameter<std::vector<double> >("jitterNoise_ns");
0090     if (temp.size() == jitterNoise_ns_.size()) {
0091       std::copy_n(temp.begin(), temp.size(), jitterNoise_ns_.begin());
0092     } else {
0093       throw cms::Exception("BadConfiguration") << " HGCFEElectronics wrong size for ToA jitterNoise ";
0094     }
0095   }
0096   if (ps.exists("jitterConstant_ns")) {
0097     auto temp = ps.getParameter<std::vector<double> >("jitterConstant_ns");
0098     if (temp.size() == jitterConstant_ns_.size()) {
0099       std::copy_n(temp.begin(), temp.size(), jitterConstant_ns_.begin());
0100     } else {
0101       throw cms::Exception("BadConfiguration") << " HGCFEElectronics wrong size for ToA jitterConstant ";
0102     }
0103   }
0104 }
0105 
0106 //
0107 template <class DFr>
0108 void HGCFEElectronics<DFr>::runTrivialShaper(
0109     DFr& dataFrame, HGCSimHitData& chargeColl, uint32_t thrADC, float lsbADC, uint32_t gainIdx, float maxADC) {
0110   bool debug(false);
0111 
0112 #ifdef EDM_ML_DEBUG
0113   for (int it = 0; it < (int)(chargeColl.size()); it++)
0114     debug |= (chargeColl[it] > adcThreshold_fC_);
0115 #endif
0116 
0117   if (debug)
0118     edm::LogVerbatim("HGCFE") << "[runTrivialShaper]" << std::endl;
0119 
0120   if (lsbADC < 0)
0121     lsbADC = adcLSB_fC_;
0122   if (maxADC < 0)
0123     // lower adcSaturation_fC_ by one part in a million
0124     // to ensure largest charge converted in bits is 0xfff==4095, not 0x1000
0125     // no effect on charges loewer than; no impact on cpu time, only done once
0126     maxADC = adcSaturation_fC_ * (1 - 1e-6);
0127   for (int it = 0; it < (int)(chargeColl.size()); it++) {
0128     //brute force saturation, maybe could to better with an exponential like saturation
0129     const uint32_t adc = std::floor(std::min(chargeColl[it], maxADC) / lsbADC);
0130     HGCSample newSample;
0131     newSample.set(adc > thrADC, false, gainIdx, 0, adc);
0132     dataFrame.setSample(it, newSample);
0133 
0134     if (debug)
0135       edm::LogVerbatim("HGCFE") << adc << " (" << chargeColl[it] << "/" << adcLSB_fC_ << ") ";
0136   }
0137 
0138   if (debug) {
0139     std::ostringstream msg;
0140     dataFrame.print(msg);
0141     edm::LogVerbatim("HGCFE") << msg.str() << std::endl;
0142   }
0143 }
0144 
0145 //
0146 template <class DFr>
0147 void HGCFEElectronics<DFr>::runSimpleShaper(DFr& dataFrame,
0148                                             HGCSimHitData& chargeColl,
0149                                             uint32_t thrADC,
0150                                             float lsbADC,
0151                                             uint32_t gainIdx,
0152                                             float maxADC,
0153                                             const hgc_digi::FEADCPulseShape& adcPulse) {
0154   //convolute with pulse shape to compute new ADCs
0155   newCharge.fill(0.f);
0156   bool debug(false);
0157   for (int it = 0; it < (int)(chargeColl.size()); it++) {
0158     const float charge(chargeColl[it]);
0159     if (charge == 0.f)
0160       continue;
0161 
0162 #ifdef EDM_ML_DEBUG
0163     debug |= (charge > adcThreshold_fC_);
0164 #endif
0165 
0166     if (debug)
0167       edm::LogVerbatim("HGCFE") << "\t Redistributing SARS ADC" << charge << " @ " << it;
0168 
0169     for (int ipulse = -2; ipulse < (int)(adcPulse.size()) - 2; ipulse++) {
0170       if (it + ipulse < 0)
0171         continue;
0172       if (it + ipulse >= (int)(dataFrame.size()))
0173         continue;
0174       const float chargeLeak = charge * adcPulse[(ipulse + 2)];
0175       newCharge[it + ipulse] += chargeLeak;
0176 
0177       if (debug)
0178         edm::LogVerbatim("HGCFE") << " | " << it + ipulse << " " << chargeLeak;
0179     }
0180 
0181     if (debug)
0182       edm::LogVerbatim("HGCFE") << std::endl;
0183   }
0184 
0185   for (int it = 0; it < (int)(newCharge.size()); it++) {
0186     //brute force saturation, maybe could to better with an exponential like saturation
0187     const uint32_t adc = std::floor(std::min(newCharge[it], maxADC) / lsbADC);
0188     HGCSample newSample;
0189     newSample.set(adc > thrADC, false, gainIdx, 0, adc);
0190     dataFrame.setSample(it, newSample);
0191 
0192     if (debug)
0193       edm::LogVerbatim("HGCFE") << adc << " (" << std::min(newCharge[it], maxADC) << "/" << lsbADC << " ) ";
0194   }
0195 
0196   if (debug) {
0197     std::ostringstream msg;
0198     dataFrame.print(msg);
0199     edm::LogVerbatim("HGCFE") << msg.str() << std::endl;
0200   }
0201 }
0202 
0203 //
0204 template <class DFr>
0205 void HGCFEElectronics<DFr>::runShaperWithToT(DFr& dataFrame,
0206                                              HGCSimHitData& chargeColl,
0207                                              HGCSimHitData& toaColl,
0208                                              CLHEP::HepRandomEngine* engine,
0209                                              uint32_t thrADC,
0210                                              float lsbADC,
0211                                              uint32_t gainIdx,
0212                                              float maxADC,
0213                                              int thickness,
0214                                              float tdcOnsetAuto,
0215                                              float noiseWidth,
0216                                              const hgc_digi::FEADCPulseShape& adcPulse) {
0217   busyFlags.fill(false);
0218   totFlags.fill(false);
0219   toaFlags.fill(false);
0220   newCharge.fill(0.f);
0221   toaFromToT.fill(0.f);
0222 
0223 #ifdef EDM_ML_DEBUG
0224   constexpr bool debug_state(true);
0225 #else
0226   constexpr bool debug_state(false);
0227 #endif
0228 
0229   bool debug = debug_state;
0230 
0231   float timeToA = 0.f;
0232 
0233   //configure the ADC <-> TDC transition depending on the value passed as argument
0234   double tdcOnset(maxADC);
0235   if (maxADC < 0) {
0236     maxADC = adcSaturation_fC_;
0237     tdcOnset = tdcOnset_fC_;
0238   }
0239 
0240   //configure the ADC LSB depending on the value passed as argument
0241   if (lsbADC < 0)
0242     lsbADC = adcLSB_fC_;
0243 
0244   //first look at time
0245   //for pileup look only at intime signals
0246   //ToA is in central BX if fired -- std::floor(BX/25.)+9;
0247   int fireBX = 9;
0248   //noise fluctuation on charge is added after ToA computation
0249   //do not recheck the ToA firing threshold tdcForToAOnset_fC_[thickness-1] not to bias the efficiency
0250   //to be done properly with realistic ToA shaper and jitter for the moment accounted in the smearing
0251   if (toaColl[fireBX] != 0.f) {
0252     timeToA = toaColl[fireBX];
0253     float sensor_noise = noiseWidth <= 0 ? noise_fC_[thickness - 1] : noiseWidth;
0254     float noise = jitterNoise_ns_[thickness - 1] * sensor_noise;
0255     float jitter = chargeColl[fireBX] == 0 ? 0 : (noise / chargeColl[fireBX]);
0256     if (jitter != 0)
0257       timeToA = CLHEP::RandGaussQ::shoot(engine, timeToA, jitter);
0258     else if (tdcResolutionInNs_ != 0)
0259       timeToA = CLHEP::RandGaussQ::shoot(engine, timeToA, tdcResolutionInNs_);
0260     timeToA += eventTimeOffset_ns_[thickness - 1];
0261     if (timeToA >= 0.f && timeToA <= 25.f)
0262       toaFlags[fireBX] = true;
0263   }
0264 
0265   //now look at charge
0266   //first identify bunches which will trigger ToT
0267   //if(debug_state) edm::LogVerbatim("HGCFE") << "[runShaperWithToT]" << std::endl;
0268   for (int it = 0; it < (int)(chargeColl.size()); ++it) {
0269     debug = debug_state;
0270     //if already flagged as busy it can't be re-used to trigger the ToT
0271     if (busyFlags[it])
0272       continue;
0273 
0274     if (tdcOnsetAuto < 0) {
0275       tdcOnsetAuto = tdcOnset_fC_;
0276     }
0277     //if below TDC onset will be handled by SARS ADC later
0278     float charge = chargeColl[it];
0279     if (charge < tdcOnset) {
0280       debug = false;
0281       continue;
0282     }
0283 
0284     //raise TDC mode for charge computation
0285     //ToA anyway fired independently will be sorted out with realistic ToA dedicated shaper
0286     float toa = timeToA;
0287     totFlags[it] = true;
0288 
0289     if (debug)
0290       edm::LogVerbatim("HGCFE") << "\t q=" << charge << " fC with <toa>=" << toa << " ns, triggers ToT @ " << it
0291                                 << std::endl;
0292 
0293     //compute total charge to be integrated and integration time
0294     //needs a loop as ToT will last as long as there is charge to dissipate
0295     int busyBxs(0);
0296     float totalCharge(charge), finalToA(toa), integTime(0);
0297     while (true) {
0298       //compute integration time in ns and # bunches
0299       //float newIntegTime(0);
0300       int poffset = 0;
0301       float charge_offset = 0.f;
0302       const float charge_kfC(totalCharge * 1e-3);
0303       if (charge_kfC < tdcChargeDrainParameterisation_[3]) {
0304         //newIntegTime=tdcChargeDrainParameterisation_[0]*pow(charge_kfC,2)+tdcChargeDrainParameterisation_[1]*charge_kfC+tdcChargeDrainParameterisation_[2];
0305       } else if (charge_kfC < tdcChargeDrainParameterisation_[7]) {
0306         poffset = 4;
0307         charge_offset = tdcChargeDrainParameterisation_[3];
0308         //newIntegTime=tdcChargeDrainParameterisation_[4]*pow(charge_kfC-tdcChargeDrainParameterisation_[3],2)+tdcChargeDrainParameterisation_[5]*(charge_kfC-tdcChargeDrainParameterisation_[3])+tdcChargeDrainParameterisation_[6];
0309       } else {
0310         poffset = 8;
0311         charge_offset = tdcChargeDrainParameterisation_[7];
0312         //newIntegTime=tdcChargeDrainParameterisation_[8]*pow(charge_kfC-tdcChargeDrainParameterisation_[7],2)+tdcChargeDrainParameterisation_[9]*(charge_kfC-tdcChargeDrainParameterisation_[7])+tdcChargeDrainParameterisation_[10];
0313       }
0314       const float charge_mod = charge_kfC - charge_offset;
0315       const float newIntegTime =
0316           ((tdcChargeDrainParameterisation_[poffset] * charge_mod + tdcChargeDrainParameterisation_[poffset + 1]) *
0317                charge_mod +
0318            tdcChargeDrainParameterisation_[poffset + 2]);
0319 
0320       const int newBusyBxs = std::floor(newIntegTime / 25.f) + 1;
0321 
0322       //if no update is needed regarding the number of bunches,
0323       //then the ToT integration time has converged
0324       integTime = newIntegTime;
0325       if (newBusyBxs == busyBxs)
0326         break;
0327 
0328       //update charge integrated during ToT
0329       if (debug) {
0330         if (busyBxs == 0)
0331           edm::LogVerbatim("HGCFE") << "\t Intial busy estimate=" << integTime << " ns = " << newBusyBxs << " bxs"
0332                                     << std::endl;
0333         else
0334           edm::LogVerbatim("HGCFE") << "\t ...integrated charge overflows initial busy estimate, interating again"
0335                                     << std::endl;
0336       }
0337 
0338       //update number of busy bunches
0339       busyBxs = newBusyBxs;
0340 
0341       //reset charge to be integrated
0342       totalCharge = charge;
0343       if (toaMode_ == WEIGHTEDBYE)
0344         finalToA = toa * charge;
0345 
0346       //add leakage from previous bunches in SARS ADC mode
0347       for (int jt = 0; jt < it; ++jt) {
0348         const unsigned int deltaT = (it - jt);
0349         if ((deltaT + 2) >= adcPulse.size() || chargeColl[jt] == 0.f || totFlags[jt] || busyFlags[jt])
0350           continue;
0351 
0352         const float leakCharge = chargeColl[jt] * adcPulse[deltaT + 2];
0353         totalCharge += leakCharge;
0354         if (toaMode_ == WEIGHTEDBYE)
0355           finalToA += leakCharge * pulseAvgT_[deltaT + 2];
0356 
0357         if (debug)
0358           edm::LogVerbatim("HGCFE") << "\t\t leaking " << chargeColl[jt] << " fC @ deltaT=-" << deltaT << " -> +"
0359                                     << leakCharge << " with avgT=" << pulseAvgT_[deltaT + 2] << std::endl;
0360       }
0361 
0362       //add contamination from posterior bunches
0363       for (int jt = it + 1; jt < it + busyBxs && jt < dataFrame.size(); ++jt) {
0364         //this charge will be integrated in TDC mode
0365         //disable for SARS ADC
0366         busyFlags[jt] = true;
0367 
0368         const float extraCharge = chargeColl[jt];
0369         if (extraCharge == 0.f)
0370           continue;
0371         if (debug)
0372           edm::LogVerbatim("HGCFE") << "\t\t adding " << extraCharge << " fC @ deltaT=+" << (jt - it) << std::endl;
0373 
0374         totalCharge += extraCharge;
0375         if (toaMode_ == WEIGHTEDBYE)
0376           finalToA += extraCharge * toaColl[jt];
0377       }
0378 
0379       //finalize ToA contamination
0380       if (toaMode_ == WEIGHTEDBYE)
0381         finalToA /= totalCharge;
0382     }
0383 
0384     newCharge[it] = (totalCharge - tdcOnset);
0385 
0386     if (debug)
0387       edm::LogVerbatim("HGCFE") << "\t Final busy estimate=" << integTime << " ns = " << busyBxs << " bxs" << std::endl
0388                                 << "\t Total integrated=" << totalCharge << " fC <toa>=" << toaFromToT[it]
0389                                 << " (raw=" << finalToA << ") ns " << std::endl;
0390 
0391     //last fC (tdcOnset) are dissipated trough pulse
0392     if (it + busyBxs < (int)(newCharge.size())) {
0393       const float deltaT2nextBx((busyBxs * 25 - integTime));
0394       const float tdcOnsetLeakage(tdcOnset * vdt::fast_expf(-deltaT2nextBx / tdcChargeDrainParameterisation_[11]));
0395       if (debug)
0396         edm::LogVerbatim("HGCFE") << "\t Leaking remainder of TDC onset " << tdcOnset << " fC, to be dissipated in "
0397                                   << deltaT2nextBx << " DeltaT/tau=" << deltaT2nextBx << " / "
0398                                   << tdcChargeDrainParameterisation_[11] << " ns, adds " << tdcOnsetLeakage << " fC @ "
0399                                   << it + busyBxs << " bx (first free bx)" << std::endl;
0400       newCharge[it + busyBxs] += tdcOnsetLeakage;
0401     }
0402   }
0403 
0404   //including the leakage from bunches in SARS ADC when not declared busy or in ToT
0405   auto runChargeSharing = [&]() {
0406     int ipulse = 0;
0407     for (int it = 0; it < (int)(chargeColl.size()); ++it) {
0408       //if busy, charge has been already integrated
0409       //if(debug) edm::LogVerbatim("HGCFE") << "\t SARS ADC pulse activated @ " << it << " : ";
0410       if (!totFlags[it] && !busyFlags[it]) {
0411         const int start = std::max(0, 2 - it);
0412         const int stop = std::min((int)adcPulse.size(), (int)newCharge.size() - it + 2);
0413         for (ipulse = start; ipulse < stop; ++ipulse) {
0414           const int itoffset = it + ipulse - 2;
0415           //notice that if the channel is already busy,
0416           //it has already been affected by the leakage of the SARS ADC
0417           //if(totFlags[itoffset] || busyFlags[itoffset]) continue;
0418           if (!totFlags[itoffset] && !busyFlags[itoffset]) {
0419             newCharge[itoffset] += chargeColl[it] * adcPulse[ipulse];
0420           }
0421           //if(debug) edm::LogVerbatim("HGCFE") << " | " << itoffset << " " << chargeColl[it]*adcPulse[ipulse] << "( " << chargeColl[it] << "->";
0422           //if(debug) edm::LogVerbatim("HGCFE") << newCharge[itoffset] << ") ";
0423         }
0424       }
0425 
0426       if (debug)
0427         edm::LogVerbatim("HGCFE") << std::endl;
0428     }
0429   };
0430   runChargeSharing();
0431 
0432   //For the future need to understand how to deal with toa for out of time signals
0433   //and for that should keep track of the BX firing the ToA somewhere (also to restore the use of finalToA)
0434   /*
0435   float finalToA(0.);
0436   for(int it=0; it<(int)(newCharge.size()); it++){
0437     if(toaFlags[it]){
0438       finalToA = toaFromToT[it];
0439       //to avoid +=25 for small negative time taken as 0
0440       while(finalToA < -1.e-5)  finalToA+=25.f;
0441       while(finalToA > 25.f) finalToA-=25.f;
0442       toaFromToT[it] = finalToA;
0443     }
0444   }
0445   */
0446   //timeToA is already in 0-25ns range by construction
0447 
0448   //set new ADCs and ToA
0449   if (debug)
0450     edm::LogVerbatim("HGCFE") << "\t final result : ";
0451   for (int it = 0; it < (int)(newCharge.size()); it++) {
0452     if (debug)
0453       edm::LogVerbatim("HGCFE") << chargeColl[it] << " -> " << newCharge[it] << " ";
0454 
0455     HGCSample newSample;
0456     if (totFlags[it] || busyFlags[it]) {
0457       if (totFlags[it]) {
0458         //brute force saturation, maybe could to better with an exponential like saturation
0459         const float saturatedCharge(std::min(newCharge[it], tdcSaturation_fC_));
0460         //working version for in-time PU and signal
0461         newSample.set(
0462             true, true, gainIdx, (uint16_t)(timeToA / toaLSB_ns_), (uint16_t)(std::floor(saturatedCharge / tdcLSB_fC_)));
0463         if (toaFlags[it])
0464           newSample.setToAValid(true);
0465       } else {
0466         newSample.set(false, true, gainIdx, 0, 0);
0467       }
0468     } else {
0469       //brute force saturation, maybe could to better with an exponential like saturation
0470       const uint16_t adc = std::floor(std::min(newCharge[it], maxADC) / lsbADC);
0471       //working version for in-time PU and signal
0472       newSample.set(adc > thrADC, false, gainIdx, (uint16_t)(timeToA / toaLSB_ns_), adc);
0473       if (toaFlags[it])
0474         newSample.setToAValid(true);
0475     }
0476     dataFrame.setSample(it, newSample);
0477   }
0478 
0479   if (debug) {
0480     std::ostringstream msg;
0481     dataFrame.print(msg);
0482     edm::LogVerbatim("HGCFE") << msg.str() << std::endl;
0483   }
0484 }
0485 
0486 // cause the compiler to generate the appropriate code
0487 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0488 template class HGCFEElectronics<HGCalDataFrame>;