CMS 3D CMS Logo

HcalSiPM.cc
Go to the documentation of this file.
3 
4 #include "CLHEP/Random/RandGaussQ.h"
5 #include "CLHEP/Random/RandPoissonQ.h"
6 #include "CLHEP/Random/RandFlat.h"
7 #include "TMath.h"
8 #include <cmath>
9 #include <cassert>
10 #include <utility>
11 
12 using std::vector;
13 //345678911234567892123456789312345678941234567895123456789612345678971234567898
15  : theCellCount(nCells), theSiPM(nCells, 1.), theCrossTalk(0.), theTempDep(0.), theLastHitTime(-1.), nonlin(nullptr) {
16  setTau(tau);
17  assert(theCellCount > 0);
18  resetSiPM();
19 }
20 
22  if (nonlin)
23  delete nonlin;
24 }
25 
26 //================================================================================
27 //implementation of Borel-Tanner distribution
28 double HcalSiPM::Borel(unsigned int n, double lambda, unsigned int k) {
29  if (n < k)
30  return 0;
31  double dn = double(n);
32  double dk = double(k);
33  double dnk = dn - dk;
34  double ldn = lambda * dn;
35  double logb = -ldn + dnk * log(ldn) - TMath::LnGamma(dnk + 1);
36  double b = 0;
37  if (logb >= -20) { // protect against underflow
38  b = (dk / dn);
39  if ((n - k) < 100)
40  b *= (exp(-ldn) * pow(ldn, dnk)) / TMath::Factorial(n - k);
41  else
42  b *= exp(logb);
43  }
44  return b;
45 }
46 
47 const HcalSiPM::cdfpair& HcalSiPM::BorelCDF(unsigned int k) {
48  // EPSILON determines the min and max # of xtalk cells that can be
49  // simulated.
50  static const double EPSILON = 1e-6;
51  typename cdfmap::const_iterator it;
52  it = borelcdfs.find(k);
53  if (it == borelcdfs.end()) {
54  vector<double> cdf;
55 
56  // Find the first n=k+i value for which cdf[i] > EPSILON
57  unsigned int i;
58  double b = 0., sumb = 0.;
59  for (i = 0;; i++) {
60  b = Borel(k + i, theCrossTalk, k);
61  sumb += b;
62  if (sumb >= EPSILON)
63  break;
64  }
65 
66  cdf.push_back(sumb);
67  unsigned int borelstartn = i;
68 
69  // calculate cdf[i]
70  for (++i;; ++i) {
71  b = Borel(k + i, theCrossTalk, k);
72  sumb += b;
73  cdf.push_back(sumb);
74  if (1 - sumb < EPSILON)
75  break;
76  }
77 
78  it = (borelcdfs.emplace(k, make_pair(borelstartn, cdf))).first;
79  }
80 
81  return it->second;
82 }
83 
84 unsigned int HcalSiPM::addCrossTalkCells(CLHEP::HepRandomEngine* engine, unsigned int in_pes) {
85  const cdfpair& cdf = BorelCDF(in_pes);
86 
87  double U = CLHEP::RandFlat::shoot(engine);
88  std::vector<double>::const_iterator up;
89  up = std::lower_bound(cdf.second.cbegin(), cdf.second.cend(), U);
90 
91  LogDebug("HcalSiPM") << "cdf size = " << cdf.second.size() << ", U = " << U << ", in_pes = " << in_pes
92  << ", 2ndary_pes = " << (up - cdf.second.cbegin() + cdf.first);
93 
94  // returns the number of secondary pes produced
95  return (up - cdf.second.cbegin() + cdf.first);
96 }
97 
98 //================================================================================
99 
100 double HcalSiPM::hitCells(CLHEP::HepRandomEngine* engine, unsigned int pes, double tempDiff, double photonTime) {
101  // response to light impulse with pes input photons. The return is the number
102  // of micro-pixels hit. If a fraction other than 0. is supplied then the
103  // micro-pixel doesn't fully discharge. The tempDiff is the temperature
104  // difference from nominal and is used to modify the relative strength of a
105  // hit pixel. Pixels which are fractionally charged return a fractional
106  // number of hit pixels.
107 
108  if ((theCrossTalk > 0.) && (theCrossTalk < 1.))
109  pes += addCrossTalkCells(engine, pes);
110 
111  // Account for saturation - disabled in lieu of recovery model below
112  //pes = nonlin->getPixelsFired(pes);
113 
114  //disable saturation/recovery model for bad tau values
115  if (theTau <= 0)
116  return pes;
117 
118  unsigned int pixel;
119  double sum(0.), hit(0.);
120  for (unsigned int pe(0); pe < pes; ++pe) {
121  pixel = CLHEP::RandFlat::shootInt(engine, theCellCount);
122  hit = (theSiPM[pixel] < 0.) ? 1.0 : (cellCharge(photonTime - theSiPM[pixel]));
123  sum += hit * (1 + (tempDiff * theTempDep));
124  theSiPM[pixel] = photonTime;
125  }
126 
127  theLastHitTime = photonTime;
128 
129  return sum;
130 }
131 
132 double HcalSiPM::totalCharge(double time) const {
133  // sum of the micro-pixels. NP is a fully charged device.
134  // 0 is a fullly depleted device.
135  double tot(0.), hit(0.);
136  for (unsigned int i = 0; i < theCellCount; ++i) {
137  hit = (theSiPM[i] < 0.) ? 1. : cellCharge(time - theSiPM[i]);
138  tot += hit;
139  }
140  return tot;
141 }
142 
144  assert(nCells > 0);
146  theSiPM.resize(nCells);
147  resetSiPM();
148 }
149 
150 void HcalSiPM::setTau(double tau) {
151  theTau = tau;
152  if (theTau > 0)
153  theTauInv = 1. / theTau;
154  else
155  theTauInv = 0;
156 }
157 
159  // set the cross-talk probability
160 
161  double oldCrossTalk = theCrossTalk;
162 
163  if ((xTalk < 0) || (xTalk >= 1)) {
164  theCrossTalk = 0.;
165  } else {
167  }
168 
169  // Recalculate the crosstalk CDFs
170  if (theCrossTalk != oldCrossTalk) {
171  borelcdfs.clear();
172  if (theCrossTalk > 0)
173  for (int k = 1; k <= 100; k++)
174  BorelCDF(k);
175  }
176 }
177 
179  // set the temperature dependence
180  theTempDep = dTemp;
181 }
182 
183 double HcalSiPM::cellCharge(double deltaTime) const {
184  if (deltaTime <= 0.)
185  return 0.;
186  if (deltaTime * theTauInv > 10.)
187  return 1.;
188  double result(1. - std::exp(-deltaTime * theTauInv));
189  return (result > 0.99) ? 1.0 : result;
190 }
191 
192 void HcalSiPM::setSaturationPars(const std::vector<float>& pars) {
193  if (nonlin)
194  delete nonlin;
195 
196  nonlin = new HcalSiPMnonlinearity(pars);
197 }
Definition: BitonicSort.h:7
double theTempDep
Definition: HcalSiPM.h:63
virtual ~HcalSiPM()
Definition: HcalSiPM.cc:21
void setSaturationPars(const std::vector< float > &pars)
Definition: HcalSiPM.cc:192
double theTauInv
Definition: HcalSiPM.h:61
void resetSiPM()
Definition: HcalSiPM.h:28
caConstants::TupleMultiplicity const CAHitNtupletGeneratorKernelsGPU::HitToTuple const cms::cuda::AtomicPairCounter GPUCACell const *__restrict__ uint32_t const *__restrict__ nCells
assert(be >=bs)
double theCrossTalk
Definition: HcalSiPM.h:62
std::vector< double > theSiPM
Definition: HcalSiPM.h:59
const cdfpair & BorelCDF(unsigned int k)
Definition: HcalSiPM.cc:47
double Borel(unsigned int n, double lambda, unsigned int k)
Definition: HcalSiPM.cc:28
virtual double hitCells(CLHEP::HepRandomEngine *, unsigned int pes, double tempDiff=0., double photonTime=0.)
Definition: HcalSiPM.cc:100
void setNCells(int nCells)
Definition: HcalSiPM.cc:143
void setTau(double tau)
Definition: HcalSiPM.cc:150
double theLastHitTime
Definition: HcalSiPM.h:64
void setCrossTalk(double xtalk)
Definition: HcalSiPM.cc:158
cdfmap borelcdfs
Definition: HcalSiPM.h:68
const float EPSILON
Definition: AnglesUtil.h:22
unsigned int addCrossTalkCells(CLHEP::HepRandomEngine *engine, unsigned int in_pes)
Definition: HcalSiPM.cc:84
HcalSiPMnonlinearity * nonlin
Definition: HcalSiPM.h:66
double b
Definition: hdecay.h:118
virtual double totalCharge() const
Definition: HcalSiPM.h:31
double theTau
Definition: HcalSiPM.h:60
unsigned int theCellCount
Definition: HcalSiPM.h:58
HcalSiPM(int nCells=1, double tau=15.)
Definition: HcalSiPM.cc:14
void setTemperatureDependence(double tempDep)
Definition: HcalSiPM.cc:178
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
std::pair< unsigned int, std::vector< double > > cdfpair
Definition: HcalSiPM.h:46
double cellCharge(double deltaTime) const
Definition: HcalSiPM.cc:183
#define LogDebug(id)