Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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