CMS 3D CMS Logo

HGCalSciNoiseMap.cc
Go to the documentation of this file.
3 #include <fstream>
4 
5 //
7  : refEdge_(3.),
8  ignoreSiPMarea_(false),
9  overrideSiPMarea_(false),
10  ignoreTileArea_(false),
11  ignoreDoseScale_(false),
12  ignoreFluenceScale_(false),
13  ignoreNoise_(false),
14  refDarkCurrent_(0.5),
15  aimMIPtoADC_(15),
16  maxSiPMPE_(8888) {
17  //number of photo electrons per MIP per scintillator type (irradiated, based on testbeam results)
18  //reference is a 30*30 mm^2 tile and 2 mm^2 SiPM (with 15um pixels), at the 2 V over-voltage
19  //based on https://indico.cern.ch/event/927798/contributions/3900921/attachments/2054679/3444966/2020Jun10_sn_scenes.pdf
20  nPEperMIP_[CAST] = 80.31 / 2; //cast
21  nPEperMIP_[MOULDED] = 57.35 / 2; //moulded
22 
23  //full scale charge per gain in nPE
24  //this is chosen for now such that the ref. MIP peak is at N ADC counts
25  //to be changed once the specs are fully defined such that the algorithm chooses the best gain
26  //for the MIP peak to be close to N ADC counts (similar to Si) or applying other specific criteria
27  fscADCPerGain_[GAIN_2] = nPEperMIP_[CAST] * 1024. / aimMIPtoADC_; //2 mm^2 SiPM
28  fscADCPerGain_[GAIN_4] = 2 * nPEperMIP_[CAST] * 1024. / aimMIPtoADC_; //4mm^2 SiPM
29 
30  //lsb: adc has 10 bits -> 1024 counts at max
31  for (size_t i = 0; i < GAINRANGE_N; i++)
32  lsbPerGain_[i] = fscADCPerGain_[i] / 1024.f;
33 }
34 
35 //
36 void HGCalSciNoiseMap::setDoseMap(const std::string& fullpath, const unsigned int algo) {
37  //decode bits of the algo word
43  ignoreNoise_ = ((algo >> IGNORE_NOISE) & 0x1);
45 
46  //call base class method
48 }
49 
50 //
52 
53 //
54 void HGCalSciNoiseMap::setNpePerMIP(float npePerMIP) {
55  nPEperMIP_[CAST] = npePerMIP;
56  nPEperMIP_[MOULDED] = npePerMIP;
57 }
58 
59 //
60 std::unordered_map<int, float> HGCalSciNoiseMap::readSipmPars(const std::string& fullpath) {
61  std::unordered_map<int, float> result;
62  //no file means default sipm size
63  if (fullpath.empty())
64  return result;
65 
67  std::ifstream infile(fp.fullPath());
68  if (!infile.is_open()) {
69  throw cms::Exception("FileNotFound") << "Unable to open '" << fullpath << "'" << std::endl;
70  }
72  while (getline(infile, line)) {
73  int layer;
74  float boundary;
75 
76  //space-separated
77  std::stringstream linestream(line);
78  linestream >> layer >> boundary;
79 
80  result[layer] = boundary;
81  }
82  return result;
83 }
84 
85 //
87 
88 //
90  const double radius,
91  int aimMIPtoADC,
92  GainRange_t gainPreChoice) {
93  int layer = cellId.layer();
94  bool hasDoseMap(!(getDoseMap().empty()));
95 
96  //LIGHT YIELD
97  double lyScaleFactor(1.f);
98  //formula is: A = A0 * exp( -D^0.65 / 199.6)
99  //where A0 is the response of the undamaged detector, D is the dose
100  if (!ignoreDoseScale_ && hasDoseMap) {
101  double cellDose = getDoseValue(DetId::HGCalHSc, layer, radius); //in kRad
102  constexpr double expofactor = 1. / 199.6;
103  const double dosespower = 0.65;
104  lyScaleFactor = std::exp(-std::pow(cellDose, dosespower) * expofactor);
105  }
106 
107  //NOISE
108  double noise(0.f);
109  if (!ignoreNoise_) {
110  double cellFluence = getFluenceValue(DetId::HGCalHSc, layer, radius); //in 1-Mev-equivalent neutrons per cm2
111 
112  //MODEL 1 : formula is N = 2.18 * sqrt(F * A / 2e13)
113  //where F is the fluence and A is the SiPM area (scaling with the latter is done below)
114  if (refDarkCurrent_ < 0) {
115  noise = 2.18;
116  if (!ignoreFluenceScale_ && hasDoseMap) {
117  constexpr double fluencefactor = 2. / (2 * 1e13); //reference SiPM area = 2mm^2
118  noise *= sqrt(cellFluence * fluencefactor);
119  }
120  }
121 
122  //MODEL 2 : formula is 3.16 * sqrt( (Idark * 1e-12) / (qe * gain) * (F / F0) )
123  //where F is the fluence (neq/cm2), gain is the SiPM gain, qe is the electron charge (C), Idark is dark current (mA)
124  else {
125  constexpr double refFluence(2.0E+13);
126  constexpr double refGain(235000.);
127  double Rdark = (refDarkCurrent_ * 1E-12) / (CLHEP::e_SI * refGain);
128  if (!ignoreFluenceScale_ && hasDoseMap)
129  Rdark *= (cellFluence / refFluence);
130  noise = 3.16 * sqrt(Rdark);
131  }
132  }
133 
134  //ADDITIONAL SCALING FACTORS
135  double tileAreaSF = scaleByTileArea(cellId, radius);
136  std::pair<double, HGCalSciNoiseMap::GainRange_t> sipm = scaleBySipmArea(cellId, radius, gainPreChoice);
137  double sipmAreaSF = sipm.first;
138  HGCalSciNoiseMap::GainRange_t gain = sipm.second;
139 
140  lyScaleFactor *= tileAreaSF * sipmAreaSF;
141  noise *= sqrt(sipmAreaSF);
142 
143  //final signal depending on scintillator type
144  double S(nPEperMIP_[CAST]);
145  if (!ignoreTileType_ && cellId.type() == 2)
146  S = nPEperMIP_[MOULDED];
147  S *= lyScaleFactor;
148 
150  sipmChar.s = S;
151  sipmChar.lySF = lyScaleFactor;
152  sipmChar.n = noise;
153  sipmChar.gain = gain;
154  sipmChar.thrADC = std::floor(0.5 * S / lsbPerGain_[gain]);
155  sipmChar.ntotalPE = maxSiPMPE_ * sipmAreaSF;
156  sipmChar.xtalk = refXtalk_;
157  return sipmChar;
158 }
159 
160 //
161 double HGCalSciNoiseMap::scaleByTileArea(const HGCScintillatorDetId& cellId, const double radius) {
162  double scaleFactor(1.f);
163 
164  if (ignoreTileArea_)
165  return scaleFactor;
166 
167  double edge(refEdge_); //start with reference 3cm of edge
168  if (cellId.type() == 0) {
169  constexpr double factor = 2 * M_PI * 1. / 360.;
170  edge = radius * factor; //1 degree
171  } else {
172  constexpr double factor = 2 * M_PI * 1. / 288.;
173  edge = radius * factor; //1.25 degrees
174  }
175  scaleFactor = refEdge_ / edge;
176  return scaleFactor;
177 }
178 
179 //
180 std::pair<double, HGCalSciNoiseMap::GainRange_t> HGCalSciNoiseMap::scaleBySipmArea(
181  const HGCScintillatorDetId& cellId, const double radius, const HGCalSciNoiseMap::GainRange_t& gainPreChoice) {
182  //start with the prechosen gain
183  //if auto then override it according to the SiPM area
184  HGCalSciNoiseMap::GainRange_t gain(gainPreChoice);
185  if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
186  gain = GainRange_t::GAIN_2;
187 
188  double scaleFactor(1.f);
189 
190  if (ignoreSiPMarea_)
191  return std::pair<double, HGCalSciNoiseMap::GainRange_t>(scaleFactor, gain);
192 
193  //use sipm area boundary map
194  if (overrideSiPMarea_) {
195  int layer = cellId.layer();
196  if (sipmMap_.count(layer) > 0 && radius < sipmMap_[layer]) {
197  scaleFactor = 2.f;
198  if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
199  gain = GainRange_t::GAIN_4;
200  }
201  }
202  //read from DetId
203  else {
204  int sipm = cellId.sipm();
205  if (sipm == 0) {
206  scaleFactor = 2.f;
207  if (gainPreChoice == HGCalSciNoiseMap::GainRange_t::AUTO)
208  gain = GainRange_t::GAIN_4;
209  }
210  }
211 
212  return std::pair<double, HGCalSciNoiseMap::GainRange_t>(scaleFactor, gain);
213 }
std::unordered_map< int, float > readSipmPars(const std::string &)
parses the radius boundaries for the SiPM area assignment from a custom file
const unsigned int & algo()
double getDoseValue(const int, const int, const double, bool logVal=false)
double scaleByTileArea(const HGCScintillatorDetId &, const double)
returns the signal scaling and the noise
void setSipmMap(const std::string &)
std::array< double, GAINRANGE_N > fscADCPerGain_
constexpr int pow(int x)
Definition: conifer.h:24
constexpr int sipm() const
get/set the sipm size
std::unordered_map< int, float > sipmMap_
std::array< double, GAINRANGE_N > lsbPerGain_
T sqrt(T t)
Definition: SSEVec.h:19
constexpr G4double scaleFactor
double f[11][100]
void setReferenceDarkCurrent(double idark)
constexpr int type() const
get/set the type
#define M_PI
double getFluenceValue(const int, const int, const double, bool logVal=false)
std::pair< double, GainRange_t > scaleBySipmArea(const HGCScintillatorDetId &, const double, const GainRange_t &)
void setDoseMap(const std::string &, const unsigned int)
std::array< double, TILETYPE_N > nPEperMIP_
#define e_SI
const double refEdge_
void setDoseMap(const std::string &, const unsigned int)
void setNpePerMIP(float npePerMIP)
constexpr int layer() const
get the layer #
SiPMonTileCharacteristics scaleByDose(const HGCScintillatorDetId &, const double, const int aimMIPtoADC=15, const GainRange_t gainPreChoice=GainRange_t::AUTO)
const doseParametersMap & getDoseMap()