HGCalSimAlgos/src/HGCalSciNoiseMap.cc

0001 #include "SimCalorimetry/HGCalSimAlgos/interface/HGCalSciNoiseMap.h"
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003 #include <fstream>
0004
0005 //
0006 HGCalSciNoiseMap::HGCalSciNoiseMap()
0007     : refEdge_(3.),
0008       ignoreSiPMarea_(false),
0009       overrideSiPMarea_(false),
0010       ignoreTileArea_(false),
0011       ignoreDoseScale_(false),
0012       ignoreFluenceScale_(false),
0013       ignoreNoise_(false),
0014       refDarkCurrent_(0.5),
0015       aimMIPtoADC_(15),
0016       maxSiPMPE_(8888) {
0017   //number of photo electrons per MIP per scintillator type (irradiated, based on testbeam results)
0018   //reference is a 30*30 mm^2 tile and 2 mm^2 SiPM (with 15um pixels), at the 2 V over-voltage
0019   //based on https://indico.cern.ch/event/927798/contributions/3900921/attachments/2054679/3444966/2020Jun10_sn_scenes.pdf
0020   nPEperMIP_[CAST] = 80.31 / 2;     //cast
0021   nPEperMIP_[MOULDED] = 57.35 / 2;  //moulded
0022
0023   //full scale charge per gain in nPE
0024   //this is chosen for now such that the ref. MIP peak is at N ADC counts
0025   //to be changed once the specs are fully defined such that the algorithm chooses the best gain
0026   //for the MIP peak to be close to N ADC counts (similar to Si) or applying other specific criteria
0027   fscADCPerGain_[GAIN_2] = nPEperMIP_[CAST] * 1024. / aimMIPtoADC_;      //2 mm^2  SiPM
0028   fscADCPerGain_[GAIN_4] = 2 * nPEperMIP_[CAST] * 1024. / aimMIPtoADC_;  //4mm^2   SiPM
0029
0030   //lsb: adc has 10 bits -> 1024 counts at max
0031   for (size_t i = 0; i < GAINRANGE_N; i++)
0032     lsbPerGain_[i] = fscADCPerGain_[i] / 1024.f;
0033 }
0034
0035 //
0036 void HGCalSciNoiseMap::setDoseMap(const std::string& fullpath, const unsigned int algo) {
0037   //decode bits of the algo word
0038   ignoreSiPMarea_ = ((algo >> IGNORE_SIPMAREA) & 0x1);
0039   overrideSiPMarea_ = ((algo >> OVERRIDE_SIPMAREA) & 0x1);
0040   ignoreTileArea_ = ((algo >> IGNORE_TILEAREA) & 0x1);
0041   ignoreDoseScale_ = ((algo >> IGNORE_DOSESCALE) & 0x1);
0042   ignoreFluenceScale_ = ((algo >> IGNORE_FLUENCESCALE) & 0x1);
0043   ignoreNoise_ = ((algo >> IGNORE_NOISE) & 0x1);
0044   ignoreTileType_ = ((algo >> IGNORE_TILETYPE) & 0x1);
0045
0046   //call base class method
0047   HGCalRadiationMap::setDoseMap(fullpath, algo);
0048 }
0049
0050 //
0051 void HGCalSciNoiseMap::setSipmMap(const std::string& fullpath) { sipmMap_ = readSipmPars(fullpath); }
0052
0053 //
0054 void HGCalSciNoiseMap::setNpePerMIP(float npePerMIP) {
0055   nPEperMIP_[CAST] = npePerMIP;
0056   nPEperMIP_[MOULDED] = npePerMIP;
0057 }
0058
0059 //
0060 std::unordered_map<int, float> HGCalSciNoiseMap::readSipmPars(const std::string& fullpath) {
0061   std::unordered_map<int, float> result;
0062   //no file means default sipm size
0063   if (fullpath.empty())
0064     return result;
0065
0066   edm::FileInPath fp(fullpath);
0067   std::ifstream infile(fp.fullPath());
0068   if (!infile.is_open()) {
0069     throw cms::Exception("FileNotFound") << "Unable to open '" << fullpath << "'" << std::endl;
0070   }
0071   std::string line;
0072   while (getline(infile, line)) {
0073     int layer;
0074     float boundary;
0075
0076     //space-separated
0077     std::stringstream linestream(line);
0078     linestream >> layer >> boundary;
0079
0080     result[layer] = boundary;
0081   }
0082   return result;
0083 }
0084
0085 //
0086 void HGCalSciNoiseMap::setReferenceDarkCurrent(double idark) { refDarkCurrent_ = idark; }
0087
0088 //
0089 HGCalSciNoiseMap::SiPMonTileCharacteristics HGCalSciNoiseMap::scaleByDose(const HGCScintillatorDetId& cellId,
0090                                                                           const double radius,
0091                                                                           int aimMIPtoADC,
0092                                                                           GainRange_t gainPreChoice) {
0093   int layer = cellId.layer();
0094   bool hasDoseMap(!(getDoseMap().empty()));
0095
0096   //LIGHT YIELD
0097   double lyScaleFactor(1.f);
0098   //formula is: A = A0 * exp( -D^0.65 / 199.6)
0099   //where A0 is the response of the undamaged detector, D is the dose
0100   if (!ignoreDoseScale_ && hasDoseMap) {
0101     double cellDose = getDoseValue(DetId::HGCalHSc, layer, radius);  //in kRad
0102     constexpr double expofactor = 1. / 199.6;
0103     const double dosespower = 0.65;
0104     lyScaleFactor = std::exp(-std::pow(cellDose, dosespower) * expofactor);
0105   }
0106
0107   //NOISE
0108   double noise(0.f);
0109   if (!ignoreNoise_) {
0110     double cellFluence = getFluenceValue(DetId::HGCalHSc, layer, radius);  //in 1-Mev-equivalent neutrons per cm2
0111
0112     //MODEL 1 : formula is N = 2.18 * sqrt(F * A / 2e13)
0113     //where F is the fluence and A is the SiPM area (scaling with the latter is done below)
0114     if (refDarkCurrent_ < 0) {
0115       noise = 2.18;
0116       if (!ignoreFluenceScale_ && hasDoseMap) {
0117         constexpr double fluencefactor = 2. / (2 * 1e13);  //reference SiPM area = 2mm^2
0118         noise *= sqrt(cellFluence * fluencefactor);
0119       }
0120     }
0121
0122     //MODEL 2 : formula is  3.16 *  sqrt( (Idark * 1e-12) / (qe * gain) * (F / F0) )
0123     //where F is the fluence (neq/cm2), gain is the SiPM gain, qe is the electron charge (C), Idark is dark current (mA)
0124     else {
0125       constexpr double refFluence(2.0E+13);
0126       constexpr double refGain(235000.);
0127       double Rdark = (refDarkCurrent_ * 1E-12) / (CLHEP::e_SI * refGain);
0128       if (!ignoreFluenceScale_ && hasDoseMap)
0129         Rdark *= (cellFluence / refFluence);
0130       noise = 3.16 * sqrt(Rdark);
0131     }
0132   }
0133
0134   //ADDITIONAL SCALING FACTORS
0135   double tileAreaSF = scaleByTileArea(cellId, radius);
0136   std::pair<double, HGCalSciNoiseMap::GainRange_t> sipm = scaleBySipmArea(cellId, radius, gainPreChoice);
0137   double sipmAreaSF = sipm.first;
0138   HGCalSciNoiseMap::GainRange_t gain = sipm.second;
0139
0140   lyScaleFactor *= tileAreaSF * sipmAreaSF;
0141   noise *= sqrt(sipmAreaSF);
0142
0143   //final signal depending on scintillator type
0144   double S(nPEperMIP_[CAST]);
0145   if (!ignoreTileType_ && cellId.type() == 2)
0146     S = nPEperMIP_[MOULDED];
0147   S *= lyScaleFactor;
0148
0149   HGCalSciNoiseMap::SiPMonTileCharacteristics sipmChar;
0150   sipmChar.s = S;
0151   sipmChar.lySF = lyScaleFactor;
0152   sipmChar.n = noise;
0153   sipmChar.gain = gain;
0154   sipmChar.thrADC = std::floor(0.5 * S / lsbPerGain_[gain]);
0155   sipmChar.ntotalPE = maxSiPMPE_ * sipmAreaSF;
0156   sipmChar.xtalk = refXtalk_;
0157   return sipmChar;
0158 }
0159
0160 //
0161 double HGCalSciNoiseMap::scaleByTileArea(const HGCScintillatorDetId& cellId, const double radius) {
0162   double scaleFactor(1.f);
0163
0164   if (ignoreTileArea_)
0165     return scaleFactor;
0166
0167   [[clang::suppress]]
0168   double edge(refEdge_);  //start with reference 3cm of edge
0169   if (cellId.type() == 0) {
0170     constexpr double factor = 2 * M_PI * 1. / 360.;
0171     edge = radius * factor;  //1 degree
0172   } else {
0173     constexpr double factor = 2 * M_PI * 1. / 288.;
0174     edge = radius * factor;  //1.25 degrees
0175   }
0176   scaleFactor = refEdge_ / edge;
0177   return scaleFactor;
0178 }
0179
0180 //
0181 std::pair<double, HGCalSciNoiseMap::GainRange_t> HGCalSciNoiseMap::scaleBySipmArea(
0182     const HGCScintillatorDetId& cellId, const double radius, const HGCalSciNoiseMap::GainRange_t& gainPreChoice) {
0183   //start with the prechosen gain
0184   //if auto then override it according to the SiPM area
0185   HGCalSciNoiseMap::GainRange_t gain(gainPreChoice);
0186   if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
0187     gain = GainRange_t::GAIN_2;
0188
0189   double scaleFactor(1.f);
0190
0191   if (ignoreSiPMarea_)
0192     return std::pair<double, HGCalSciNoiseMap::GainRange_t>(scaleFactor, gain);
0193
0194   //use sipm area boundary map
0195   if (overrideSiPMarea_) {
0196     int layer = cellId.layer();
0197     if (sipmMap_.count(layer) > 0 && radius < sipmMap_[layer]) {
0198       scaleFactor = 2.f;
0199       if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
0200         gain = GainRange_t::GAIN_4;
0201     }
0202   }
0203   //read from DetId
0204   else {
0205     int sipm = cellId.sipm();
0206     if (sipm == 0) {
0207       scaleFactor = 2.f;
0208       if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
0209         gain = GainRange_t::GAIN_4;
0210     }
0211   }
0212
0213   return std::pair<double, HGCalSciNoiseMap::GainRange_t>(scaleFactor, gain);
0214 }