HGCalSimAlgos/interface/HGCalSiNoiseMap.icc

0001 #include "SimCalorimetry/HGCalSimAlgos/interface/HGCalSiNoiseMap.h"
0002 #include "TMath.h"
0003
0004 //
0005 template <typename T>
0006 HGCalSiNoiseMap<T>::HGCalSiNoiseMap()
0007     : defaultGain_(GainRange_t::AUTO),
0008       defaultAimMIPtoADC_(10),
0009       encCommonNoiseSub_(sqrt(1.0)),
0010       qe2fc_(1.60217646E-4),
0011       ignoreFluence_(false),
0012       ignoreCCE_(false),
0013       ignoreNoise_(false),
0014       ignoreGainDependentPulse_(false),
0015       activateCachedOp_(false) {
0016   //q80fC
0017   encsParam_.push_back({636., 15.6, 0.0328});   //2nd order polynomial coefficients as function of capacitance
0018   chargeAtFullScaleADCPerGain_.push_back(80.);  //the num of fC (charge) which corresponds to the max ADC value
0019   adcPulses_.push_back({{0., 0., 1.0, 0.066 / 0.934, 0., 0.}});  //in-time bunch is the 3rd entry in the array
0020   //q160fC
0021   encsParam_.push_back({1045., 8.74, 0.0685});
0022   chargeAtFullScaleADCPerGain_.push_back(160.);
0023   adcPulses_.push_back({{0., 0., 1.0, 0.153 / 0.847, 0., 0.}});
0024   // q320fC
0025   encsParam_.push_back({1915., 2.79, 0.0878});
0026   chargeAtFullScaleADCPerGain_.push_back(320.);
0027   adcPulses_.push_back({{0., 0., 1.0, 0.0963 / 0.9037, 0., 0.}});
0028
0029   //start with a default value
0030   defaultADCPulse_ = adcPulses_[(GainRange_t)q160fC];
0031
0032   // adc has 10 bits -> 1024 counts at max ( >0 baseline to be handled)
0033   for (auto i : chargeAtFullScaleADCPerGain_)
0034     lsbPerGain_.push_back(i / 1024.);
0035
0036   //fine sensors: 120 mum -  67: MPV of charge[number of e-]/mum for a mip in silicon; srouce PDG
0037   const double mipEqfC_120 = 120. * 67. * qe2fc_;
0038   mipEqfC_[0] = mipEqfC_120;
0039   const double cellCapacitance_120 = 50;
0040   cellCapacitance_[0] = cellCapacitance_120;
0041   const double cellVolume_120 = 0.56 * (120.e-4);
0042   cellVolume_[0] = cellVolume_120;
0043
0044   //thin sensors: 200 mum
0045   const double mipEqfC_200 = 200. * 70. * qe2fc_;
0046   mipEqfC_[1] = mipEqfC_200;
0047   const double cellCapacitance_200 = 65;
0048   cellCapacitance_[1] = cellCapacitance_200;
0049   const double cellVolume_200 = 1.26 * (200.e-4);
0050   cellVolume_[1] = cellVolume_200;
0051
0052   //thick sensors: 300 mum
0053   const double mipEqfC_300 = 300. * 73. * qe2fc_;
0054   mipEqfC_[2] = mipEqfC_300;
0055   const double cellCapacitance_300 = 45;
0056   cellCapacitance_[2] = cellCapacitance_300;
0057   const double cellVolume_300 = 1.26 * (300.e-4);
0058   cellVolume_[2] = cellVolume_300;
0059 }
0060
0061 //
0062 template <typename T>
0063 void HGCalSiNoiseMap<T>::setDoseMap(const std::string &fullpath, const unsigned int &algo) {
0064   //decode bits in the algo word
0065   ignoreFluence_ = ((algo >> FLUENCE) & 0x1);
0066   ignoreCCE_ = ((algo >> CCE) & 0x1);
0067   ignoreNoise_ = ((algo >> NOISE) & 0x1);
0068   ignoreGainDependentPulse_ = ((algo >> PULSEPERGAIN) & 0x1);
0069   activateCachedOp_ = ((algo >> CACHEDOP) & 0x1);
0070
0071   //call base class method
0072   HGCalRadiationMap::setDoseMap(fullpath, algo);
0073 }
0074
0075 //
0076 template <typename T>
0077 double HGCalSiNoiseMap<T>::getENCpad(double ileak) {
0078   if (ileak > 45.40)
0079     return 23.30 * ileak + 1410.04;
0080   else if (ileak > 38.95)
0081     return 30.07 * ileak + 1156.76;
0082   else if (ileak > 32.50)
0083     return 38.58 * ileak + 897.94;
0084   else if (ileak > 26.01)
0085     return 193.67 * pow(ileak, 0.70) + 21.12;
0086   else if (ileak > 19.59)
0087     return 167.60 * pow(ileak, 0.77);
0088   else if (ileak > 13.06)
0089     return 162.35 * pow(ileak, 0.82);
0090   else if (ileak > 6.53)
0091     return 202.73 * pow(ileak, 0.81);
0092   else
0093     return 457.15 * pow(ileak, 0.57);
0094 }
0095
0096 //
0097 template <typename T>
0098 double HGCalSiNoiseMap<T>::getTDCOnsetAuto(uint32_t gainIdx) {
0099   if (gainIdx < 3) {
0100     return chargeAtFullScaleADCPerGain_[gainIdx] - 1e-6;
0101   }
0102
0103   return -1;
0104 }
0105
0106 //
0107 template <typename T>
0108 void HGCalSiNoiseMap<T>::setGeometry(const CaloSubdetectorGeometry *hgcGeom, GainRange_t gain, int aimMIPtoADC) {
0109   //call base class method
0110   HGCalRadiationMap::setGeometry(hgcGeom);
0111
0112   defaultGain_ = gain;
0113   defaultAimMIPtoADC_ = aimMIPtoADC;
0114
0115   //exit if cache is to be ignored
0116   if (!activateCachedOp_)
0117     return;
0118
0119   //fill cache if it's not filled
0120   if (!siopCache_.empty())
0121     return;
0122
0123   const auto &validDetIds = geom()->getValidDetIds();
0124   for (const auto &did : validDetIds) {
0125     unsigned int rawId(did.rawId());
0126     T hgcDetId(rawId);
0127
0128     //compute and store in cache
0129     siopCache_.emplace(rawId, getSiCellOpCharacteristicsCore(hgcDetId));
0130   }
0131 }
0132
0133 //
0134 template <typename T>
0135 const typename HGCalSiNoiseMap<T>::SiCellOpCharacteristicsCore HGCalSiNoiseMap<T>::getSiCellOpCharacteristicsCore(
0136     const T &cellId, GainRange_t gain, int aimMIPtoADC) {
0137   //re-compute
0138   if (!activateCachedOp_)
0139     return getSiCellOpCharacteristics(cellId, gain, aimMIPtoADC).core;
0140
0141   uint32_t key = cellId.rawId();
0142
0143   return siopCache_[key];
0144 }
0145
0146 //
0147 template <typename T>
0148 typename HGCalSiNoiseMap<T>::SiCellOpCharacteristics HGCalSiNoiseMap<T>::getSiCellOpCharacteristics(const T &cellId,
0149                                                                                                     GainRange_t gain,
0150                                                                                                     int aimMIPtoADC) {
0151   //decode cell properties
0152   int layer = cellId.layer();
0153   unsigned int cellThick = cellId.type();
0154   double cellCap(cellCapacitance_[cellThick]);
0155   double cellVol(cellVolume_[cellThick]);
0156   double mipEqfC(mipEqfC_[cellThick]);
0157   std::vector<double> &cceParam = cceParam_[cellThick];
0158
0159   //location of the cell
0160   int subdet = cellId.subdet();
0161   const auto &xy(ddd()->locateCell(cellId.zside(),
0162                                    cellId.layer(),
0163                                    cellId.waferU(),
0164                                    cellId.waferV(),
0165                                    cellId.cellU(),
0166                                    cellId.cellV(),
0167                                    true,
0168                                    true,
0169                                    false,
0170                                    false,
0171                                    false));
0172   double radius = sqrt(std::pow(xy.first, 2) + std::pow(xy.second, 2));  //in cm
0173
0174   //call baseline method and add to cache
0175   return getSiCellOpCharacteristics(cellCap, cellVol, mipEqfC, cceParam, subdet, layer, radius, gain, aimMIPtoADC);
0176 }
0177
0178 //
0179 template <typename T>
0180 typename HGCalSiNoiseMap<T>::SiCellOpCharacteristics HGCalSiNoiseMap<T>::getSiCellOpCharacteristics(
0181     double &cellCap,
0182     double &cellVol,
0183     double &mipEqfC,
0184     std::vector<double> &cceParam,
0185     int &subdet,
0186     int &layer,
0187     double &radius,
0188     GainRange_t &gain,
0189     int &aimMIPtoADC) {
0190   SiCellOpCharacteristics siop;
0191
0192   //leakage current and CCE [muA]
0193   if (ignoreFluence_) {
0194     siop.fluence = 0;
0195     siop.lnfluence = -1;
0196     siop.ileak = exp(ileakParam_[1]) * cellVol * unitToMicro_;
0197     siop.core.cce = 1;
0198   } else {
0199     if (getDoseMap().empty()) {
0200       throw cms::Exception("BadConfiguration")
0201           << " Fluence is required but no DoseMap has been passed to HGCalSiNoiseMap";
0202       return siop;
0203     }
0204
0205     siop.lnfluence = getFluenceValue(subdet, layer, radius, true);
0206     siop.fluence = exp(siop.lnfluence);
0207
0208     double conv(log(cellVol) + unitToMicroLog_);
0209     siop.ileak = exp(ileakParam_[0] * siop.lnfluence + ileakParam_[1] + conv);
0210
0211     //charge collection efficiency
0212     if (ignoreCCE_) {
0213       siop.core.cce = 1.0;
0214     } else {
0215       //in the region where measurements are available we have used a linear approximation in log (f)
0216       //the extrapolation to the lower fluence regions is done with an error function matched at the boundary
0217       //this is a bit more expensive than the previous approximation (lin-lin approx) but it fits
0218       //better the TDR, TTU and the new 2021 data
0219       if (siop.fluence > cceParam[0]) {
0220         siop.core.cce = cceParam[2] * siop.lnfluence + cceParam[1];
0221       } else {
0222         double at = cceParam[2] * log(cceParam[0]) + cceParam[1];
0223         double bt = -cceParam[0] / TMath::ErfcInverse(1. / at);
0224         siop.core.cce = at * TMath::Erfc((siop.fluence - cceParam[0]) / bt);
0225       }
0226
0227       siop.core.cce = std::max((float)0., siop.core.cce);
0228     }
0229   }
0230
0231   //determine the gain to apply accounting for cce
0232   //algo:  start with the most favored = lowest gain possible (=highest range)
0233   //       test for the other gains in the preferred order
0234   //       the first to yield <=15 ADC counts is taken
0235   //       this relies on the fact that these gains shift the mip peak by factors of 2
0236   //       in the presence of more gains 15 should be updated accordingly
0237   //note:  move computation to higher granularity level (ROC, trigger tower, once decided)
0238   double S(siop.core.cce * mipEqfC);
0239   if (gain == GainRange_t::AUTO) {
0240     gain = GainRange_t::q320fC;
0241     std::vector<GainRange_t> orderedGainChoice = {GainRange_t::q160fC, GainRange_t::q80fC};
0242     for (const auto &igain : orderedGainChoice) {
0243       double mipPeakADC(S / lsbPerGain_[igain]);
0244       if (mipPeakADC > 16)
0245         break;
0246       gain = igain;
0247     }
0248   }
0249
0250   //fill in the parameters of the struct
0251   siop.core.gain = gain;
0252   siop.mipfC = S;
0253   siop.mipADC = std::floor(S / lsbPerGain_[gain]);
0254   siop.core.thrADC = std::floor(S / 2. / lsbPerGain_[gain]);
0255
0256   //build noise estimate
0257   if (ignoreNoise_) {
0258     siop.core.noise = 0.0;
0259     siop.enc_s = 0.0;
0260     siop.enc_p = 0.0;
0261   } else {
0262     siop.enc_s = encsParam_[gain][0] + encsParam_[gain][1] * cellCap + encsParam_[gain][2] * pow(cellCap, 2);
0263     siop.enc_p = getENCpad(siop.ileak);
0264     siop.core.noise = hypot(siop.enc_p * encCommonNoiseSub_, siop.enc_s) * qe2fc_;
0265   }
0266
0267   return siop;
0268 }