HcalSimAlgos/src/HcalSiPM.cc

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 //345678911234567892123456789312345678941234567895123456789612345678971234567898
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   //original implementation of Borel-Tanner distribution
0035   // here for reference
0036   constexpr double originalBorel(unsigned int n, double lambda, unsigned int k) {
0037     if (n < k)
0038       return 0;
0039     double dn = double(n);
0040     double dk = double(k);
0041     double dnk = dn - dk;
0042     double ldn = lambda * dn;
0043     double logb = -ldn + dnk * log(ldn) - TMath::LnGamma(dnk + 1);
0044     double b = 0;
0045     if (logb >= -20) {  // protect against underflow
0046       b = (dk / dn);
0047       if ((n - k) < 100)
0048         b *= (exp(-ldn) * pow(ldn, dnk)) / TMath::Factorial(n - k);
0049       else
0050         b *= exp(logb);
0051     }
0052     return b;
0053   }
0054   */
0055
0056   using FLOAT = double;
0057   //================================================================================
0058   //modified implementation of Borel-Tanner distribution
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       // protect against underflow
0073       b = (logb >= -20.) ? b0 * std::exp(logb) : 0;
0074     }
0075     return b;
0076   }
0077
0078 }  // namespace
0079
0080 const HcalSiPM::cdfpair& HcalSiPM::BorelCDF(unsigned int k) {
0081   // EPSILON determines the min and max # of xtalk cells that can be
0082   // simulated.
0083   constexpr double EPSILON = 1e-6;
0084   constexpr uint32_t maxCDFsize = 170;  // safe max to avoid factorial to be infinite
0085   auto it = borelcdfs.find(k);
0086   if (it == borelcdfs.end()) {
0087     vector<double> cdf;
0088     cdf.reserve(64);
0089
0090     // Find the first n=k+i value for which cdf[i] > EPSILON
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     // calculate cdf[i]  limit to 170 to avoid iFact to become infinite
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   // returns the number of secondary pes produced
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   // response to light impulse with pes input photons.  The return is the number
0137   // of micro-pixels hit.  If a fraction other than 0. is supplied then the
0138   // micro-pixel doesn't fully discharge.  The tempDiff is the temperature
0139   // difference from nominal and is used to modify the relative strength of a
0140   // hit pixel.  Pixels which are fractionally charged return a fractional
0141   // number of hit pixels.
0142
0143   if ((theCrossTalk > 0.) && (theCrossTalk < 1.))
0144     pes += addCrossTalkCells(engine, pes);
0145
0146   // Account for saturation - disabled in lieu of recovery model below
0147   //pes = nonlin->getPixelsFired(pes);
0148
0149   //disable saturation/recovery model for bad tau values
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   // sum of the micro-pixels.  NP is a fully charged device.
0169   // 0 is a fullly depleted device.
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   // set the cross-talk probability
0194
0195   double oldCrossTalk = theCrossTalk;
0196
0197   if ((xTalk < 0) || (xTalk >= 1)) {
0198     theCrossTalk = 0.;
0199   } else {
0200     theCrossTalk = xTalk;
0201   }
0202
0203   // Recalculate the crosstalk CDFs
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   // set the temperature dependence
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 }