Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:28:39

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), theSiPM(nCells, 1.), theCrossTalk(0.), theTempDep(0.), theLastHitTime(-1.), nonlin(nullptr) {
0016   setTau(tau);
0017   assert(theCellCount > 0);
0018   resetSiPM();
0019 }
0020 
0021 HcalSiPM::~HcalSiPM() {
0022   if (nonlin)
0023     delete nonlin;
0024 }
0025 
0026 //================================================================================
0027 //implementation of Borel-Tanner distribution
0028 double HcalSiPM::Borel(unsigned int n, double lambda, unsigned int k) {
0029   if (n < k)
0030     return 0;
0031   double dn = double(n);
0032   double dk = double(k);
0033   double dnk = dn - dk;
0034   double ldn = lambda * dn;
0035   double logb = -ldn + dnk * log(ldn) - TMath::LnGamma(dnk + 1);
0036   double b = 0;
0037   if (logb >= -20) {  // protect against underflow
0038     b = (dk / dn);
0039     if ((n - k) < 100)
0040       b *= (exp(-ldn) * pow(ldn, dnk)) / TMath::Factorial(n - k);
0041     else
0042       b *= exp(logb);
0043   }
0044   return b;
0045 }
0046 
0047 const HcalSiPM::cdfpair& HcalSiPM::BorelCDF(unsigned int k) {
0048   // EPSILON determines the min and max # of xtalk cells that can be
0049   // simulated.
0050   static const double EPSILON = 1e-6;
0051   typename cdfmap::const_iterator it;
0052   it = borelcdfs.find(k);
0053   if (it == borelcdfs.end()) {
0054     vector<double> cdf;
0055 
0056     // Find the first n=k+i value for which cdf[i] > EPSILON
0057     unsigned int i;
0058     double b = 0., sumb = 0.;
0059     for (i = 0;; i++) {
0060       b = Borel(k + i, theCrossTalk, k);
0061       sumb += b;
0062       if (sumb >= EPSILON)
0063         break;
0064     }
0065 
0066     cdf.push_back(sumb);
0067     unsigned int borelstartn = i;
0068 
0069     // calculate cdf[i]
0070     for (++i;; ++i) {
0071       b = Borel(k + i, theCrossTalk, k);
0072       sumb += b;
0073       cdf.push_back(sumb);
0074       if (1 - sumb < EPSILON)
0075         break;
0076     }
0077 
0078     it = (borelcdfs.emplace(k, make_pair(borelstartn, cdf))).first;
0079   }
0080 
0081   return it->second;
0082 }
0083 
0084 unsigned int HcalSiPM::addCrossTalkCells(CLHEP::HepRandomEngine* engine, unsigned int in_pes) {
0085   const cdfpair& cdf = BorelCDF(in_pes);
0086 
0087   double U = CLHEP::RandFlat::shoot(engine);
0088   std::vector<double>::const_iterator up;
0089   up = std::lower_bound(cdf.second.cbegin(), cdf.second.cend(), U);
0090 
0091   LogDebug("HcalSiPM") << "cdf size = " << cdf.second.size() << ", U = " << U << ", in_pes = " << in_pes
0092                        << ", 2ndary_pes = " << (up - cdf.second.cbegin() + cdf.first);
0093 
0094   // returns the number of secondary pes produced
0095   return (up - cdf.second.cbegin() + cdf.first);
0096 }
0097 
0098 //================================================================================
0099 
0100 double HcalSiPM::hitCells(CLHEP::HepRandomEngine* engine, unsigned int pes, double tempDiff, double photonTime) {
0101   // response to light impulse with pes input photons.  The return is the number
0102   // of micro-pixels hit.  If a fraction other than 0. is supplied then the
0103   // micro-pixel doesn't fully discharge.  The tempDiff is the temperature
0104   // difference from nominal and is used to modify the relative strength of a
0105   // hit pixel.  Pixels which are fractionally charged return a fractional
0106   // number of hit pixels.
0107 
0108   if ((theCrossTalk > 0.) && (theCrossTalk < 1.))
0109     pes += addCrossTalkCells(engine, pes);
0110 
0111   // Account for saturation - disabled in lieu of recovery model below
0112   //pes = nonlin->getPixelsFired(pes);
0113 
0114   //disable saturation/recovery model for bad tau values
0115   if (theTau <= 0)
0116     return pes;
0117 
0118   unsigned int pixel;
0119   double sum(0.), hit(0.);
0120   for (unsigned int pe(0); pe < pes; ++pe) {
0121     pixel = CLHEP::RandFlat::shootInt(engine, theCellCount);
0122     hit = (theSiPM[pixel] < 0.) ? 1.0 : (cellCharge(photonTime - theSiPM[pixel]));
0123     sum += hit * (1 + (tempDiff * theTempDep));
0124     theSiPM[pixel] = photonTime;
0125   }
0126 
0127   theLastHitTime = photonTime;
0128 
0129   return sum;
0130 }
0131 
0132 double HcalSiPM::totalCharge(double time) const {
0133   // sum of the micro-pixels.  NP is a fully charged device.
0134   // 0 is a fullly depleted device.
0135   double tot(0.), hit(0.);
0136   for (unsigned int i = 0; i < theCellCount; ++i) {
0137     hit = (theSiPM[i] < 0.) ? 1. : cellCharge(time - theSiPM[i]);
0138     tot += hit;
0139   }
0140   return tot;
0141 }
0142 
0143 void HcalSiPM::setNCells(int nCells) {
0144   assert(nCells > 0);
0145   theCellCount = nCells;
0146   theSiPM.resize(nCells);
0147   resetSiPM();
0148 }
0149 
0150 void HcalSiPM::setTau(double tau) {
0151   theTau = tau;
0152   if (theTau > 0)
0153     theTauInv = 1. / theTau;
0154   else
0155     theTauInv = 0;
0156 }
0157 
0158 void HcalSiPM::setCrossTalk(double xTalk) {
0159   // set the cross-talk probability
0160 
0161   double oldCrossTalk = theCrossTalk;
0162 
0163   if ((xTalk < 0) || (xTalk >= 1)) {
0164     theCrossTalk = 0.;
0165   } else {
0166     theCrossTalk = xTalk;
0167   }
0168 
0169   // Recalculate the crosstalk CDFs
0170   if (theCrossTalk != oldCrossTalk) {
0171     borelcdfs.clear();
0172     if (theCrossTalk > 0)
0173       for (int k = 1; k <= 100; k++)
0174         BorelCDF(k);
0175   }
0176 }
0177 
0178 void HcalSiPM::setTemperatureDependence(double dTemp) {
0179   // set the temperature dependence
0180   theTempDep = dTemp;
0181 }
0182 
0183 double HcalSiPM::cellCharge(double deltaTime) const {
0184   if (deltaTime <= 0.)
0185     return 0.;
0186   if (deltaTime * theTauInv > 10.)
0187     return 1.;
0188   double result(1. - std::exp(-deltaTime * theTauInv));
0189   return (result > 0.99) ? 1.0 : result;
0190 }
0191 
0192 void HcalSiPM::setSaturationPars(const std::vector<float>& pars) {
0193   if (nonlin)
0194     delete nonlin;
0195 
0196   nonlin = new HcalSiPMnonlinearity(pars);
0197 }