Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:03:44

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   double radius = sqrt(std::pow(xy.first, 2) + std::pow(xy.second, 2));  //in cm
0172 
0173   //call baseline method and add to cache
0174   return getSiCellOpCharacteristics(cellCap, cellVol, mipEqfC, cceParam, subdet, layer, radius, gain, aimMIPtoADC);
0175 }
0176 
0177 //
0178 template <typename T>
0179 typename HGCalSiNoiseMap<T>::SiCellOpCharacteristics HGCalSiNoiseMap<T>::getSiCellOpCharacteristics(
0180     double &cellCap,
0181     double &cellVol,
0182     double &mipEqfC,
0183     std::vector<double> &cceParam,
0184     int &subdet,
0185     int &layer,
0186     double &radius,
0187     GainRange_t &gain,
0188     int &aimMIPtoADC) {
0189   SiCellOpCharacteristics siop;
0190 
0191   //leakage current and CCE [muA]
0192   if (ignoreFluence_) {
0193     siop.fluence = 0;
0194     siop.lnfluence = -1;
0195     siop.ileak = exp(ileakParam_[1]) * cellVol * unitToMicro_;
0196     siop.core.cce = 1;
0197   } else {
0198     if (getDoseMap().empty()) {
0199       throw cms::Exception("BadConfiguration")
0200           << " Fluence is required but no DoseMap has been passed to HGCalSiNoiseMap";
0201       return siop;
0202     }
0203 
0204     siop.lnfluence = getFluenceValue(subdet, layer, radius, true);
0205     siop.fluence = exp(siop.lnfluence);
0206 
0207     double conv(log(cellVol) + unitToMicroLog_);
0208     siop.ileak = exp(ileakParam_[0] * siop.lnfluence + ileakParam_[1] + conv);
0209 
0210     //charge collection efficiency
0211     if (ignoreCCE_) {
0212       siop.core.cce = 1.0;
0213     } else {
0214       //in the region where measurements are available we have used a linear approximation in log (f)
0215       //the extrapolation to the lower fluence regions is done with an error function matched at the boundary
0216       //this is a bit more expensive than the previous approximation (lin-lin approx) but it fits
0217       //better the TDR, TTU and the new 2021 data
0218       if (siop.fluence > cceParam[0]) {
0219         siop.core.cce = cceParam[2] * siop.lnfluence + cceParam[1];
0220       } else {
0221         double at = cceParam[2] * log(cceParam[0]) + cceParam[1];
0222         double bt = -cceParam[0] / TMath::ErfcInverse(1. / at);
0223         siop.core.cce = at * TMath::Erfc((siop.fluence - cceParam[0]) / bt);
0224       }
0225 
0226       siop.core.cce = std::max((float)0., siop.core.cce);
0227     }
0228   }
0229 
0230   //determine the gain to apply accounting for cce
0231   //algo:  start with the most favored = lowest gain possible (=highest range)
0232   //       test for the other gains in the preferred order
0233   //       the first to yield <=15 ADC counts is taken
0234   //       this relies on the fact that these gains shift the mip peak by factors of 2
0235   //       in the presence of more gains 15 should be updated accordingly
0236   //note:  move computation to higher granularity level (ROC, trigger tower, once decided)
0237   double S(siop.core.cce * mipEqfC);
0238   if (gain == GainRange_t::AUTO) {
0239     gain = GainRange_t::q320fC;
0240     std::vector<GainRange_t> orderedGainChoice = {GainRange_t::q160fC, GainRange_t::q80fC};
0241     for (const auto &igain : orderedGainChoice) {
0242       double mipPeakADC(S / lsbPerGain_[igain]);
0243       if (mipPeakADC > 16)
0244         break;
0245       gain = igain;
0246     }
0247   }
0248 
0249   //fill in the parameters of the struct
0250   siop.core.gain = gain;
0251   siop.mipfC = S;
0252   siop.mipADC = std::floor(S / lsbPerGain_[gain]);
0253   siop.core.thrADC = std::floor(S / 2. / lsbPerGain_[gain]);
0254 
0255   //build noise estimate
0256   if (ignoreNoise_) {
0257     siop.core.noise = 0.0;
0258     siop.enc_s = 0.0;
0259     siop.enc_p = 0.0;
0260   } else {
0261     siop.enc_s = encsParam_[gain][0] + encsParam_[gain][1] * cellCap + encsParam_[gain][2] * pow(cellCap, 2);
0262     siop.enc_p = getENCpad(siop.ileak);
0263     siop.core.noise = hypot(siop.enc_p * encCommonNoiseSub_, siop.enc_s) * qe2fc_;
0264   }
0265 
0266   return siop;
0267 }