File indexing completed on 2024-04-06 12:29:32
0001 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSiPM.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "CLHEP/Random/RandGaussQ.h"
0005 #include "CLHEP/Random/RandPoissonQ.h"
0006 #include "CLHEP/Random/RandFlat.h"
0007 #include "TMath.h"
0008 #include <cmath>
0009 #include <cassert>
0010 #include <utility>
0011
0012 using std::vector;
0013
0014 HcalSiPM::HcalSiPM(int nCells, double tau)
0015 : theCellCount(nCells),
0016 theSiPM(nCells, -999.f),
0017 theCrossTalk(0.),
0018 theTempDep(0.),
0019 theLastHitTime(-1.),
0020 nonlin(nullptr) {
0021 setTau(tau);
0022 assert(theCellCount > 0);
0023 }
0024
0025 HcalSiPM::~HcalSiPM() {
0026 if (nonlin)
0027 delete nonlin;
0028 }
0029
0030 namespace {
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056 using FLOAT = double;
0057
0058
0059 constexpr double Borel(unsigned int i, FLOAT lambda, unsigned int k, double iFact) {
0060 auto n = k + i;
0061 FLOAT dn = FLOAT(n);
0062 FLOAT dk = FLOAT(k);
0063 FLOAT dnk = FLOAT(i);
0064
0065 FLOAT ldn = lambda * dn;
0066 FLOAT b0 = (dk / dn);
0067 FLOAT b = 0;
0068 if (i < 100) {
0069 b = b0 * (std::exp(-ldn) * std::pow(ldn, dnk)) / iFact;
0070 } else {
0071 FLOAT logb = -ldn + dnk * std::log(ldn) - std::log(iFact);
0072
0073 b = (logb >= -20.) ? b0 * std::exp(logb) : 0;
0074 }
0075 return b;
0076 }
0077
0078 }
0079
0080 const HcalSiPM::cdfpair& HcalSiPM::BorelCDF(unsigned int k) {
0081
0082
0083 constexpr double EPSILON = 1e-6;
0084 constexpr uint32_t maxCDFsize = 170;
0085 auto it = borelcdfs.find(k);
0086 if (it == borelcdfs.end()) {
0087 vector<double> cdf;
0088 cdf.reserve(64);
0089
0090
0091 unsigned int i;
0092 double sumb = 0.;
0093 double iFact = 1.;
0094 for (i = 0; i < maxCDFsize; i++) {
0095 if (i > 0)
0096 iFact *= double(i);
0097 sumb += Borel(i, theCrossTalk, k, iFact);
0098 if (sumb >= EPSILON)
0099 break;
0100 }
0101
0102 cdf.push_back(sumb);
0103 unsigned int borelstartn = i;
0104
0105
0106 for (++i; i < maxCDFsize; ++i) {
0107 iFact *= double(i);
0108 sumb += Borel(i, theCrossTalk, k, iFact);
0109 cdf.push_back(sumb);
0110 if (1. - sumb < EPSILON)
0111 break;
0112 }
0113 it = (borelcdfs.emplace(k, make_pair(borelstartn, cdf))).first;
0114 }
0115
0116 return it->second;
0117 }
0118
0119 unsigned int HcalSiPM::addCrossTalkCells(CLHEP::HepRandomEngine* engine, unsigned int in_pes) {
0120 const cdfpair& cdf = BorelCDF(in_pes);
0121
0122 double U = CLHEP::RandFlat::shoot(engine);
0123 std::vector<double>::const_iterator up;
0124 up = std::lower_bound(cdf.second.cbegin(), cdf.second.cend(), U);
0125
0126 LogDebug("HcalSiPM") << "cdf size = " << cdf.second.size() << ", U = " << U << ", in_pes = " << in_pes
0127 << ", 2ndary_pes = " << (up - cdf.second.cbegin() + cdf.first);
0128
0129
0130 return (up - cdf.second.cbegin() + cdf.first);
0131 }
0132
0133
0134
0135 double HcalSiPM::hitCells(CLHEP::HepRandomEngine* engine, unsigned int pes, double tempDiff, double photonTime) {
0136
0137
0138
0139
0140
0141
0142
0143 if ((theCrossTalk > 0.) && (theCrossTalk < 1.))
0144 pes += addCrossTalkCells(engine, pes);
0145
0146
0147
0148
0149
0150 if (theTau <= 0)
0151 return pes;
0152
0153 unsigned int pixel;
0154 double sum(0.), hit(0.);
0155 for (unsigned int pe(0); pe < pes; ++pe) {
0156 pixel = CLHEP::RandFlat::shootInt(engine, theCellCount);
0157 hit = (theSiPM[pixel] < 0.) ? 1.0 : (cellCharge(photonTime - theSiPM[pixel]));
0158 sum += hit * (1 + (tempDiff * theTempDep));
0159 theSiPM[pixel] = photonTime;
0160 }
0161
0162 theLastHitTime = photonTime;
0163
0164 return sum;
0165 }
0166
0167 double HcalSiPM::totalCharge(double time) const {
0168
0169
0170 double tot(0.), hit(0.);
0171 for (unsigned int i = 0; i < theCellCount; ++i) {
0172 hit = (theSiPM[i] < 0.) ? 1. : cellCharge(time - theSiPM[i]);
0173 tot += hit;
0174 }
0175 return tot;
0176 }
0177
0178 void HcalSiPM::setNCells(int nCells) {
0179 theCellCount = nCells;
0180 theSiPM.resize(nCells);
0181 resetSiPM();
0182 }
0183
0184 void HcalSiPM::setTau(double tau) {
0185 theTau = tau;
0186 if (theTau > 0)
0187 theTauInv = 1. / theTau;
0188 else
0189 theTauInv = 0;
0190 }
0191
0192 void HcalSiPM::setCrossTalk(double xTalk) {
0193
0194
0195 double oldCrossTalk = theCrossTalk;
0196
0197 if ((xTalk < 0) || (xTalk >= 1)) {
0198 theCrossTalk = 0.;
0199 } else {
0200 theCrossTalk = xTalk;
0201 }
0202
0203
0204 if (theCrossTalk != oldCrossTalk) {
0205 borelcdfs.clear();
0206 if (theCrossTalk > 0)
0207 for (int k = 1; k <= 100; k++)
0208 BorelCDF(k);
0209 }
0210 }
0211
0212 void HcalSiPM::setTemperatureDependence(double dTemp) {
0213
0214 theTempDep = dTemp;
0215 }
0216
0217 double HcalSiPM::cellCharge(double deltaTime) const {
0218 if (deltaTime <= 0.)
0219 return 0.;
0220 if (deltaTime * theTauInv > 10.)
0221 return 1.;
0222 double result(1. - std::exp(-deltaTime * theTauInv));
0223 return (result > 0.99) ? 1.0 : result;
0224 }
0225
0226 void HcalSiPM::setSaturationPars(const std::vector<float>& pars) {
0227 if (nonlin)
0228 delete nonlin;
0229
0230 nonlin = new HcalSiPMnonlinearity(pars);
0231 }