Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-11-07 23:58:26 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!