Merge pull request cms-sw#40634 from VinInn/SpeedupHCalPreMix4

Spedup HCAL Digitization

Author: cmsbuild

CondFormats/HcalObjects/src/HcalSiPMCharacteristics.cc

@@ -77,6 +77,7 @@ int HcalSiPMCharacteristics::getPixels(int type) const {
 std::vector<float> HcalSiPMCharacteristics::getNonLinearities(int type) const {
   const HcalSiPMCharacteristics::PrecisionItem* item = findByType(type);
   std::vector<float> pars;
+  pars.reserve(3);
   if (item) {
     pars.push_back(item->parLin1_);
     pars.push_back(item->parLin2_);

SimCalorimetry/HcalSimAlgos/interface/HcalSiPM.h

@@ -19,17 +19,17 @@ namespace CLHEP {
   class HepRandomEngine;
 }
 
-class HcalSiPM {
+class HcalSiPM final {
 public:
   HcalSiPM(int nCells = 1, double tau = 15.);
 
-  virtual ~HcalSiPM();
+  ~HcalSiPM();
 
-  void resetSiPM() { std::fill(theSiPM.begin(), theSiPM.end(), -999.); }
-  virtual double hitCells(CLHEP::HepRandomEngine*, unsigned int pes, double tempDiff = 0., double photonTime = 0.);
+  void resetSiPM() { std::fill(theSiPM.begin(), theSiPM.end(), -999.f); }
+  double hitCells(CLHEP::HepRandomEngine*, unsigned int pes, double tempDiff = 0., double photonTime = 0.);
 
-  virtual double totalCharge() const { return totalCharge(theLastHitTime); }
-  virtual double totalCharge(double time) const;
+  double totalCharge() const { return totalCharge(theLastHitTime); }
+  double totalCharge(double time) const;
 
   int getNCells() const { return theCellCount; }
   double getTau() const { return theTau; }
@@ -52,11 +52,10 @@ class HcalSiPM {
   unsigned int addCrossTalkCells(CLHEP::HepRandomEngine* engine, unsigned int in_pes);
 
   //numerical random generation from Borel-Tanner distribution
-  double Borel(unsigned int n, double lambda, unsigned int k);
   const cdfpair& BorelCDF(unsigned int k);
 
   unsigned int theCellCount;
-  std::vector<double> theSiPM;
+  std::vector<float> theSiPM;
   double theTau;
   double theTauInv;
   double theCrossTalk;

SimCalorimetry/HcalSimAlgos/interface/HcalSiPMHitResponse.h

@@ -19,7 +19,7 @@ class PCaloHitCompareTimes {
   bool operator()(const PCaloHit* a, const PCaloHit* b) const { return a->time() < b->time(); }
 };
 
-class HcalSiPMHitResponse : public CaloHitResponse {
+class HcalSiPMHitResponse final : public CaloHitResponse {
 public:
   HcalSiPMHitResponse(const CaloVSimParameterMap* parameterMap,
                       const CaloShapes* shapes,

SimCalorimetry/HcalSimAlgos/interface/HcalSiPMShape.h

@@ -2,17 +2,24 @@
 #ifndef HcalSimAlgos_HcalSiPMShape_h
 #define HcalSimAlgos_HcalSiPMShape_h
 
+#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVShape.h"
 #include <vector>
 
-class HcalSiPMShape : public CaloVShape {
+class HcalSiPMShape final : public CaloVShape {
 public:
   HcalSiPMShape(unsigned int signalShape = 206);
   HcalSiPMShape(const HcalSiPMShape& other);
 
   ~HcalSiPMShape() override {}
 
-  double operator()(double time) const override;
+  int nBins() const { return nBins_; }
+  double operator[](int i) const { return nt_[i]; }
+
+  double operator()(double time) const override {
+    int jtime(time * HcalPulseShapes::invDeltaTSiPM_ + 0.5);
+    return (jtime >= 0 && jtime < nBins_) ? nt_[jtime] : 0;
+  }
 
   double timeToRise() const override { return 0.0; }
 

SimCalorimetry/HcalSimAlgos/src/HcalSiPM.cc

@@ -12,69 +12,104 @@
 using std::vector;
 //345678911234567892123456789312345678941234567895123456789612345678971234567898
 HcalSiPM::HcalSiPM(int nCells, double tau)
-    : theCellCount(nCells), theSiPM(nCells, 1.), theCrossTalk(0.), theTempDep(0.), theLastHitTime(-1.), nonlin(nullptr) {
+    : theCellCount(nCells),
+      theSiPM(nCells, -999.f),
+      theCrossTalk(0.),
+      theTempDep(0.),
+      theLastHitTime(-1.),
+      nonlin(nullptr) {
   setTau(tau);
   assert(theCellCount > 0);
-  resetSiPM();
 }
 
 HcalSiPM::~HcalSiPM() {
   if (nonlin)
     delete nonlin;
 }
 
-//================================================================================
-//implementation of Borel-Tanner distribution
-double HcalSiPM::Borel(unsigned int n, double lambda, unsigned int k) {
-  if (n < k)
-    return 0;
-  double dn = double(n);
-  double dk = double(k);
-  double dnk = dn - dk;
-  double ldn = lambda * dn;
-  double logb = -ldn + dnk * log(ldn) - TMath::LnGamma(dnk + 1);
-  double b = 0;
-  if (logb >= -20) {  // protect against underflow
-    b = (dk / dn);
-    if ((n - k) < 100)
-      b *= (exp(-ldn) * pow(ldn, dnk)) / TMath::Factorial(n - k);
-    else
-      b *= exp(logb);
+namespace {
+
+  /*
+  //================================================================================
+  //original implementation of Borel-Tanner distribution
+  // here for reference
+  constexpr double originalBorel(unsigned int n, double lambda, unsigned int k) {
+    if (n < k)
+      return 0;
+    double dn = double(n);
+    double dk = double(k);
+    double dnk = dn - dk;
+    double ldn = lambda * dn;
+    double logb = -ldn + dnk * log(ldn) - TMath::LnGamma(dnk + 1);
+    double b = 0;
+    if (logb >= -20) {  // protect against underflow
+      b = (dk / dn);
+      if ((n - k) < 100)
+        b *= (exp(-ldn) * pow(ldn, dnk)) / TMath::Factorial(n - k);
+      else
+        b *= exp(logb);
+    }
+    return b;
+  }
+  */
+
+  using FLOAT = double;
+  //================================================================================
+  //modified implementation of Borel-Tanner distribution
+  constexpr double Borel(unsigned int i, FLOAT lambda, unsigned int k, double iFact) {
+    auto n = k + i;
+    FLOAT dn = FLOAT(n);
+    FLOAT dk = FLOAT(k);
+    FLOAT dnk = FLOAT(i);
+
+    FLOAT ldn = lambda * dn;
+    FLOAT b0 = (dk / dn);
+    FLOAT b = 0;
+    if (i < 100) {
+      b = b0 * (std::exp(-ldn) * std::pow(ldn, dnk)) / iFact;
+    } else {
+      FLOAT logb = -ldn + dnk * std::log(ldn) - std::log(iFact);
+      // protect against underflow
+      b = (logb >= -20.) ? b0 * std::exp(logb) : 0;
+    }
+    return b;
   }
-  return b;
-}
+
+}  // namespace
 
 const HcalSiPM::cdfpair& HcalSiPM::BorelCDF(unsigned int k) {
   // EPSILON determines the min and max # of xtalk cells that can be
   // simulated.
-  static const double EPSILON = 1e-6;
-  typename cdfmap::const_iterator it;
-  it = borelcdfs.find(k);
+  constexpr double EPSILON = 1e-6;
+  constexpr uint32_t maxCDFsize = 170;  // safe max to avoid factorial to be infinite
+  auto it = borelcdfs.find(k);
   if (it == borelcdfs.end()) {
     vector<double> cdf;
+    cdf.reserve(64);
 
     // Find the first n=k+i value for which cdf[i] > EPSILON
     unsigned int i;
-    double b = 0., sumb = 0.;
-    for (i = 0;; i++) {
-      b = Borel(k + i, theCrossTalk, k);
-      sumb += b;
+    double sumb = 0.;
+    double iFact = 1.;
+    for (i = 0; i < maxCDFsize; i++) {
+      if (i > 0)
+        iFact *= double(i);
+      sumb += Borel(i, theCrossTalk, k, iFact);
       if (sumb >= EPSILON)
         break;
     }
 
     cdf.push_back(sumb);
     unsigned int borelstartn = i;
 
-    // calculate cdf[i]
-    for (++i;; ++i) {
-      b = Borel(k + i, theCrossTalk, k);
-      sumb += b;
+    // calculate cdf[i]  limit to 170 to avoid iFact to become infinite
+    for (++i; i < maxCDFsize; ++i) {
+      iFact *= double(i);
+      sumb += Borel(i, theCrossTalk, k, iFact);
       cdf.push_back(sumb);
-      if (1 - sumb < EPSILON)
+      if (1. - sumb < EPSILON)
         break;
     }
-
     it = (borelcdfs.emplace(k, make_pair(borelstartn, cdf))).first;
   }
 
@@ -141,7 +176,6 @@ double HcalSiPM::totalCharge(double time) const {
 }
 
 void HcalSiPM::setNCells(int nCells) {
-  assert(nCells > 0);
   theCellCount = nCells;
   theSiPM.resize(nCells);
   resetSiPM();

SimCalorimetry/HcalSimAlgos/src/HcalSiPMHitResponse.cc

@@ -205,59 +205,69 @@ CaloSamples HcalSiPMHitResponse::makeSiPMSignal(DetId const& id,
   unsigned int pe(0);
   double hitPixels(0.), elapsedTime(0.);
 
-  auto& sipmPulseShape(shapeMap[pars.signalShape(id)]);
+  auto const& sipmPulseShape(shapeMap[pars.signalShape(id)]);
 
-  std::list<std::pair<double, double> > pulses;
-  std::list<std::pair<double, double> >::iterator pulse;
-  double timeDiff, pulseBit;
   LogDebug("HcalSiPMHitResponse") << "makeSiPMSignal for " << HcalDetId(id);
 
-  for (unsigned int tbin(0); tbin < photonTimeBins.size(); ++tbin) {
+  const int nptb = photonTimeBins.size();
+  double sum[nptb];
+  for (auto i = 0; i < nptb; ++i)
+    sum[i] = 0;
+  for (int tbin(0); tbin < nptb; ++tbin) {
     pe = photonTimeBins[tbin];
+    if (pe <= 0)
+      continue;
     preciseBin = tbin;
     sampleBin = preciseBin / nbins;
-    if (pe > 0) {
-      //skip saturation/recovery and pulse smearing for premix stage 1
-      if (PreMixDigis and HighFidelityPreMix) {
-        signal[sampleBin] += pe;
-        signal.preciseAtMod(preciseBin) += pe;
-        elapsedTime += dt;
-        continue;
-      }
-
-      hitPixels = theSiPM.hitCells(engine, pe, 0., elapsedTime);
-      LogDebug("HcalSiPMHitResponse") << " elapsedTime: " << elapsedTime << " sampleBin: " << sampleBin
-                                      << " preciseBin: " << preciseBin << " pe: " << pe << " hitPixels: " << hitPixels;
-      if (pars.doSiPMSmearing()) {
-        pulses.push_back(std::pair<double, double>(elapsedTime, hitPixels));
-      } else {
-        signal[sampleBin] += hitPixels;
-        signal.preciseAtMod(preciseBin) += 0.6 * hitPixels;
-        if (preciseBin > 0)
-          signal.preciseAtMod(preciseBin - 1) += 0.2 * hitPixels;
-        if (preciseBin < signal.preciseSize() - 1)
-          signal.preciseAtMod(preciseBin + 1) += 0.2 * hitPixels;
-      }
+    //skip saturation/recovery and pulse smearing for premix stage 1
+    if (PreMixDigis and HighFidelityPreMix) {
+      signal[sampleBin] += pe;
+      signal.preciseAtMod(preciseBin) += pe;
+      elapsedTime += dt;
+      continue;
     }
 
-    if (pars.doSiPMSmearing()) {
-      pulse = pulses.begin();
-      while (pulse != pulses.end()) {
-        timeDiff = elapsedTime - pulse->first;
-        pulseBit = sipmPulseShape(timeDiff) * pulse->second;
-        LogDebug("HcalSiPMHitResponse") << " pulse t: " << pulse->first << " pulse A: " << pulse->second
-                                        << " timeDiff: " << timeDiff << " pulseBit: " << pulseBit;
-        signal[sampleBin] += pulseBit;
-        signal.preciseAtMod(preciseBin) += pulseBit;
-
-        if (timeDiff > 1 && sipmPulseShape(timeDiff) < 1e-7)
-          pulse = pulses.erase(pulse);
-        else
-          ++pulse;
+    hitPixels = theSiPM.hitCells(engine, pe, 0., elapsedTime);
+    LogDebug("HcalSiPMHitResponse") << " elapsedTime: " << elapsedTime << " sampleBin: " << sampleBin
+                                    << " preciseBin: " << preciseBin << " pe: " << pe << " hitPixels: " << hitPixels;
+    if (!pars.doSiPMSmearing()) {
+      signal[sampleBin] += hitPixels;
+      signal.preciseAtMod(preciseBin) += 0.6 * hitPixels;
+      if (preciseBin > 0)
+        signal.preciseAtMod(preciseBin - 1) += 0.2 * hitPixels;
+      if (preciseBin < signal.preciseSize() - 1)
+        signal.preciseAtMod(preciseBin + 1) += 0.2 * hitPixels;
+    } else {
+      // add "my" smearing to future bins...
+      // this loop can vectorize....
+      for (auto i = tbin; i < nptb; ++i) {
+        auto itdiff = i - tbin;
+        if (itdiff == sipmPulseShape.nBins())
+          break;
+        auto shape = sipmPulseShape[itdiff];
+        auto pulseBit = shape * hitPixels;
+        sum[i] += pulseBit;
+        if (shape < 1.e-7 && itdiff > int(HcalPulseShapes::invDeltaTSiPM_))
+          break;
       }
     }
     elapsedTime += dt;
   }
+  if (pars.doSiPMSmearing())
+    for (auto i = 0; i < nptb; ++i) {
+      auto iSampleBin = i / nbins;
+      signal[iSampleBin] += sum[i];
+      signal.preciseAtMod(i) += sum[i];
+    }
+
+#ifdef EDM_ML_DEBUG
+  LogDebug("HcalSiPMHitResponse") << nbins << ' ' << nptb << ' ' << HcalDetId(id);
+  for (auto i = 0; i < nptb; ++i) {
+    auto iSampleBin = (nbins > 1) ? i / nbins : i;
+    LogDebug("HcalSiPMHitResponse") << i << ' ' << iSampleBin << ' ' << signal[iSampleBin] << ' '
+                                    << signal.preciseAtMod(i);
+  }
+#endif
 
   return signal;
 }

SimCalorimetry/HcalSimAlgos/src/HcalSiPMShape.cc

@@ -10,13 +10,6 @@ HcalSiPMShape::HcalSiPMShape(unsigned int signalShape)
 
 HcalSiPMShape::HcalSiPMShape(const HcalSiPMShape& other) : CaloVShape(other), nBins_(other.nBins_), nt_(other.nt_) {}
 
-double HcalSiPMShape::operator()(double time) const {
-  int jtime(time * HcalPulseShapes::invDeltaTSiPM_ + 0.5);
-  if (jtime >= 0 && jtime < nBins_)
-    return nt_[jtime];
-  return 0.;
-}
-
 void HcalSiPMShape::computeShape(unsigned int signalShape) {
   //grab correct function pointer based on shape
   double (*analyticPulseShape)(double);