Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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   double edge(refEdge_);  //start with reference 3cm of edge
0168   if (cellId.type() == 0) {
0169     constexpr double factor = 2 * M_PI * 1. / 360.;
0170     edge = radius * factor;  //1 degree
0171   } else {
0172     constexpr double factor = 2 * M_PI * 1. / 288.;
0173     edge = radius * factor;  //1.25 degrees
0174   }
0175   scaleFactor = refEdge_ / edge;
0176   return scaleFactor;
0177 }
0178 
0179 //
0180 std::pair<double, HGCalSciNoiseMap::GainRange_t> HGCalSciNoiseMap::scaleBySipmArea(
0181     const HGCScintillatorDetId& cellId, const double radius, const HGCalSciNoiseMap::GainRange_t& gainPreChoice) {
0182   //start with the prechosen gain
0183   //if auto then override it according to the SiPM area
0184   HGCalSciNoiseMap::GainRange_t gain(gainPreChoice);
0185   if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
0186     gain = GainRange_t::GAIN_2;
0187 
0188   double scaleFactor(1.f);
0189 
0190   if (ignoreSiPMarea_)
0191     return std::pair<double, HGCalSciNoiseMap::GainRange_t>(scaleFactor, gain);
0192 
0193   //use sipm area boundary map
0194   if (overrideSiPMarea_) {
0195     int layer = cellId.layer();
0196     if (sipmMap_.count(layer) > 0 && radius < sipmMap_[layer]) {
0197       scaleFactor = 2.f;
0198       if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
0199         gain = GainRange_t::GAIN_4;
0200     }
0201   }
0202   //read from DetId
0203   else {
0204     int sipm = cellId.sipm();
0205     if (sipm == 0) {
0206       scaleFactor = 2.f;
0207       if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
0208         gain = GainRange_t::GAIN_4;
0209     }
0210   }
0211 
0212   return std::pair<double, HGCalSciNoiseMap::GainRange_t>(scaleFactor, gain);
0213 }