Merge pull request #30422 from franzoni/hgcal_eol_pulse_update_112X_bis

HGCAL end of life DIGITISATION update

SimCalorimetry/HGCalSimAlgos/interface/HGCalRadiationMap.h

@@ -7,36 +7,38 @@
 #include "vdt/vdtMath.h"
 #include <string>
 
-typedef std::array<double, 8> radiiVec;
-
 /**
    @class HGCalRadiationMap
    @short parses a txt file with dose/fluence parameters and provides functions for noise, etc.
  */
 class HGCalRadiationMap {
 public:
   struct DoseParameters {
-    DoseParameters() : a_(0.), b_(0.), c_(0.), d_(0.), e_(0.), f_(0.), g_(0.), h_(0.), i_(0.), j_(0.) {}
-    double a_, b_, c_, d_, e_, f_, g_, h_, i_, j_;
+    DoseParameters()
+        : a_(0.), b_(0.), c_(0.), d_(0.), e_(0.), doff_(0.), f_(0.), g_(0.), h_(0.), i_(0.), j_(0.), foff_(0.) {}
+    double a_, b_, c_, d_, e_, doff_, f_, g_, h_, i_, j_, foff_;
   };
 
-  HGCalRadiationMap(){};
+  HGCalRadiationMap();
   ~HGCalRadiationMap(){};
 
   typedef std::map<std::pair<int, int>, DoseParameters> doseParametersMap;
 
   void setGeometry(const CaloSubdetectorGeometry *);
   void setDoseMap(const std::string &, const unsigned int);
 
-  double getDoseValue(const int, const int, const radiiVec &, bool logVal = false);
-  double getFluenceValue(const int, const int, const radiiVec &, bool logVal = false);
+  double computeRadius(const HGCScintillatorDetId &);
+
+  double getDoseValue(const int, const int, const double, bool logVal = false);
+  double getFluenceValue(const int, const int, const double, bool logVal = false);
 
   const unsigned int &algo() { return algo_; }
   const HGCalGeometry *geom() { return hgcalGeom_; }
   const HGCalTopology *topo() { return hgcalTopology_; }
   const HGCalDDDConstants *ddd() { return hgcalDDD_; }
 
   inline const doseParametersMap &getDoseMap() { return doseMap_; }
+  inline void setFluenceScaleFactor(double val) { fluenceSFlog10_ = log10(val); }
 
 private:
   doseParametersMap readDosePars(const std::string &);
@@ -48,6 +50,7 @@ class HGCalRadiationMap {
   doseParametersMap doseMap_;
   //conversion from grey to krad
   const double grayToKrad_ = 0.1;
+  double fluenceSFlog10_;
 };
 
 #endif

SimCalorimetry/HGCalSimAlgos/interface/HGCalSciNoiseMap.h

@@ -18,18 +18,18 @@ class HGCalSciNoiseMap : public HGCalRadiationMap {
   /**
      @short returns the signal scaling and the noise
   */
-  double scaleByTileArea(const HGCScintillatorDetId&, const radiiVec&);
-  double scaleBySipmArea(const HGCScintillatorDetId&, const double&);
-  std::pair<double, double> scaleByDose(const HGCScintillatorDetId&, const radiiVec&);
+  double scaleByTileArea(const HGCScintillatorDetId &, const double);
+  double scaleBySipmArea(const HGCScintillatorDetId &, const double);
+  std::pair<double, double> scaleByDose(const HGCScintillatorDetId &, const double);
 
-  radiiVec computeRadius(const HGCScintillatorDetId&);
-  void setSipmMap(const std::string&);
+  void setSipmMap(const std::string &);
 
 private:
-  std::unordered_map<int, float> readSipmPars(const std::string&);
+  std::unordered_map<int, float> readSipmPars(const std::string &);
 
   //size of the reference scintillator tile
   const double refEdge_;
+
   //sipm size boundaries
   std::unordered_map<int, float> sipmMap_;
 };

SimCalorimetry/HGCalSimAlgos/interface/HGCalSiNoiseMap.h

@@ -2,10 +2,12 @@
 #define simcalorimetry_hgcalsimalgos_hgcalsinoisemap
 
 #include "SimCalorimetry/HGCalSimAlgos/interface/HGCalRadiationMap.h"
+#include "SimCalorimetry/HGCalSimProducers/interface/HGCFEElectronics.h"
 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
 #include <string>
 #include <array>
+#include <unordered_map>
 
 /**
    @class HGCalSiNoiseMap
@@ -14,13 +16,19 @@
 class HGCalSiNoiseMap : public HGCalRadiationMap {
 public:
   enum GainRange_t { q80fC, q160fC, q320fC, AUTO };
-  enum NoiseMapAlgoBits_t { FLUENCE, CCE, NOISE };
+  enum NoiseMapAlgoBits_t { FLUENCE, CCE, NOISE, PULSEPERGAIN, CACHEDOP };
+
+  struct SiCellOpCharacteristicsCore {
+    SiCellOpCharacteristicsCore() : cce(0.), noise(0.), gain(0), thrADC(0) {}
+    float cce, noise;
+    unsigned short gain, thrADC;
+  };
 
   struct SiCellOpCharacteristics {
-    SiCellOpCharacteristics()
-        : lnfluence(0.), fluence(0.), ileak(0.), cce(1.), noise(0.), mipfC(0), gain(0), mipADC(0), thrADC(0) {}
-    double lnfluence, fluence, ileak, cce, noise, mipfC;
-    unsigned int gain, mipADC, thrADC;
+    SiCellOpCharacteristics() : lnfluence(0.), fluence(0.), ileak(0.), enc_s(0.), enc_p(0.), mipfC(0), mipADC(0) {}
+    SiCellOpCharacteristicsCore core;
+    double lnfluence, fluence, ileak, enc_s, enc_p, mipfC;
+    unsigned int mipADC;
   };
 
   HGCalSiNoiseMap();
@@ -47,13 +55,35 @@ class HGCalSiNoiseMap : public HGCalRadiationMap {
   */
   void setDoseMap(const std::string &, const unsigned int &);
 
+  /**
+     @short specialization of the base class method which sets the geometry so that it can instantiate an operation
+     cache the first time it is called - intrinsically related to the valid detIds in the geometry
+     the filling of the cache is ignored by configuration or if it has already been filled
+   */
+  void setGeometry(const CaloSubdetectorGeometry *, GainRange_t gain = GainRange_t::AUTO, int aimMIPtoADC = 10);
+
   /**
      @short returns the charge collection efficiency and noise
      if gain range is set to auto, it will find the most appropriate gain to put the mip peak close to 10 ADC counts
   */
+  const SiCellOpCharacteristicsCore getSiCellOpCharacteristicsCore(const HGCSiliconDetId &did,
+                                                                   GainRange_t gain,
+                                                                   int aimMIPtoADC);
+  const SiCellOpCharacteristicsCore getSiCellOpCharacteristicsCore(const HGCSiliconDetId &did) {
+    return getSiCellOpCharacteristicsCore(did, defaultGain_, defaultAimMIPtoADC_);
+  }
   SiCellOpCharacteristics getSiCellOpCharacteristics(const HGCSiliconDetId &did,
                                                      GainRange_t gain = GainRange_t::AUTO,
                                                      int aimMIPtoADC = 10);
+  SiCellOpCharacteristics getSiCellOpCharacteristics(double &cellCap,
+                                                     double &cellVol,
+                                                     double &mipEqfC,
+                                                     std::vector<double> &cceParam,
+                                                     int &subdet,
+                                                     int &layer,
+                                                     double &radius,
+                                                     GainRange_t &gain,
+                                                     int &aimMIPtoADC);
 
   std::array<double, 3> &getMipEqfC() { return mipEqfC_; }
   std::array<double, 3> &getCellCapacitance() { return cellCapacitance_; }
@@ -62,36 +92,52 @@ class HGCalSiNoiseMap : public HGCalRadiationMap {
   std::vector<double> &getIleakParam() { return ileakParam_; }
   std::vector<std::vector<double> > &getENCsParam() { return encsParam_; }
   std::vector<double> &getLSBPerGain() { return lsbPerGain_; }
+  void setDefaultADCPulseShape(const hgc_digi::FEADCPulseShape &adcPulse) { defaultADCPulse_ = adcPulse; };
+  const hgc_digi::FEADCPulseShape &adcPulseForGain(GainRange_t gain) {
+    if (ignoreGainDependentPulse_)
+      return defaultADCPulse_;
+    return adcPulses_[gain];
+  };
   std::vector<double> &getMaxADCPerGain() { return chargeAtFullScaleADCPerGain_; }
+  double getENCpad(double ileak);
+  void setCachedOp(bool flag) { activateCachedOp_ = flag; }
+
+  inline void setENCCommonNoiseSubScale(double val) { encCommonNoiseSub_ = val; }
 
 private:
+  GainRange_t defaultGain_;
+  int defaultAimMIPtoADC_;
+
+  //cache of SiCellOpCharacteristics
+  std::map<uint32_t, SiCellOpCharacteristicsCore> siopCache_;
+
   //vector of three params, per sensor type: 0:120 [mum], 1:200, 2:300
   std::array<double, 3> mipEqfC_, cellCapacitance_, cellVolume_;
   std::vector<std::vector<double> > cceParam_;
 
   //leakage current/volume vs fluence
   std::vector<double> ileakParam_;
 
-  //shaper noise param
-  const double encpScale_;
-
   //common noise subtraction noise (final scaling value)
-  const double encCommonNoiseSub_;
+  double encCommonNoiseSub_;
 
   //electron charge in fC
   const double qe2fc_;
 
-  //electronics noise (series+parallel) polynomial coeffs;
+  //electronics noise (series+parallel) polynomial coeffs and ADC pulses;
   std::vector<std::vector<double> > encsParam_;
+  hgc_digi::FEADCPulseShape defaultADCPulse_;
+  std::vector<hgc_digi::FEADCPulseShape> adcPulses_;
 
   //lsb
   std::vector<double> lsbPerGain_, chargeAtFullScaleADCPerGain_;
 
   //conversions
   const double unitToMicro_ = 1.e6;
+  const double unitToMicroLog_ = log(unitToMicro_);
 
-  //flags used to disable specific components of the Si operation parameters
-  bool ignoreFluence_, ignoreCCE_, ignoreNoise_;
+  //flags used to disable specific components of the Si operation parameters or usage of operation cache
+  bool ignoreFluence_, ignoreCCE_, ignoreNoise_, ignoreGainDependentPulse_, activateCachedOp_;
 };
 
 #endif

SimCalorimetry/HGCalSimAlgos/src/HGCalRadiationMap.cc

@@ -2,6 +2,9 @@
 #include "FWCore/ParameterSet/interface/FileInPath.h"
 #include <fstream>
 
+//
+HGCalRadiationMap::HGCalRadiationMap() : fluenceSFlog10_(0.) {}
+
 //
 void HGCalRadiationMap::setDoseMap(const std::string& fullpath, const unsigned int algo) {
   doseMap_ = readDosePars(fullpath);
@@ -16,18 +19,39 @@ void HGCalRadiationMap::setGeometry(const CaloSubdetectorGeometry* geom) {
 }
 
 //
-double HGCalRadiationMap::getDoseValue(const int subdet, const int layer, const radiiVec& radius, bool logVal) {
+double HGCalRadiationMap::computeRadius(const HGCScintillatorDetId& cellId) {
+  GlobalPoint global = geom()->getPosition(cellId);
+  return std::sqrt(std::pow(global.x(), 2) + std::pow(global.y(), 2));
+}
+
+//
+double HGCalRadiationMap::getDoseValue(const int subdet, const int layer, const double radius, bool logVal) {
   std::pair<int, int> key(subdet, layer);
-  double cellDoseLog10 = doseMap_[key].a_ + doseMap_[key].b_ * radius[4] + doseMap_[key].c_ * radius[5] +
-                         doseMap_[key].d_ * radius[6] + doseMap_[key].e_ * radius[7];
+
+  double r(radius - doseMap_[key].doff_);
+  double r2(r * r);
+  double r3(r2 * r);
+  double r4(r3 * r);
+
+  double cellDoseLog10 =
+      doseMap_[key].a_ + doseMap_[key].b_ * r + doseMap_[key].c_ * r2 + doseMap_[key].d_ * r3 + doseMap_[key].e_ * r4;
+
   return logVal ? cellDoseLog10 * M_LN10 + log(grayToKrad_) : std::pow(10, cellDoseLog10) * grayToKrad_;
 }
 
 //
-double HGCalRadiationMap::getFluenceValue(const int subdet, const int layer, const radiiVec& radius, bool logVal) {
+double HGCalRadiationMap::getFluenceValue(const int subdet, const int layer, const double radius, bool logVal) {
   std::pair<int, int> key(subdet, layer);
-  double cellFluenceLog10 = doseMap_[key].f_ + doseMap_[key].g_ * radius[0] + doseMap_[key].h_ * radius[1] +
-                            doseMap_[key].i_ * radius[2] + doseMap_[key].j_ * radius[3];
+
+  double r(radius - doseMap_[key].foff_);
+  double r2(r * r);
+  double r3(r2 * r);
+  double r4(r3 * r);
+
+  double cellFluenceLog10 =
+      doseMap_[key].f_ + doseMap_[key].g_ * r + doseMap_[key].h_ * r2 + doseMap_[key].i_ * r3 + doseMap_[key].j_ * r4;
+  cellFluenceLog10 += fluenceSFlog10_;
+
   return logVal ? cellFluenceLog10 * M_LN10 : std::pow(10, cellFluenceLog10);
 }
 
@@ -54,7 +78,7 @@ std::map<std::pair<int, int>, HGCalRadiationMap::DoseParameters> HGCalRadiationM
     //space-separated
     std::stringstream linestream(line);
     linestream >> subdet >> layer >> dosePars.a_ >> dosePars.b_ >> dosePars.c_ >> dosePars.d_ >> dosePars.e_ >>
-        dosePars.f_ >> dosePars.g_ >> dosePars.h_ >> dosePars.i_ >> dosePars.j_;
+        dosePars.doff_ >> dosePars.f_ >> dosePars.g_ >> dosePars.h_ >> dosePars.i_ >> dosePars.j_ >> dosePars.foff_;
 
     std::pair<int, int> key(subdet, layer);
     result[key] = dosePars;

SimCalorimetry/HGCalSimAlgos/src/HGCalSciNoiseMap.cc

@@ -35,7 +35,7 @@ std::unordered_map<int, float> HGCalSciNoiseMap::readSipmPars(const std::string&
 }
 
 //
-std::pair<double, double> HGCalSciNoiseMap::scaleByDose(const HGCScintillatorDetId& cellId, const radiiVec& radius) {
+std::pair<double, double> HGCalSciNoiseMap::scaleByDose(const HGCScintillatorDetId& cellId, const double radius) {
   if (getDoseMap().empty())
     return std::make_pair(1., 0.);
 
@@ -58,22 +58,22 @@ std::pair<double, double> HGCalSciNoiseMap::scaleByDose(const HGCScintillatorDet
   return std::make_pair(scaleFactor, noise);
 }
 
-double HGCalSciNoiseMap::scaleByTileArea(const HGCScintillatorDetId& cellId, const radiiVec& radius) {
+double HGCalSciNoiseMap::scaleByTileArea(const HGCScintillatorDetId& cellId, const double radius) {
   double edge;
   if (cellId.type() == 0) {
     constexpr double factor = 2 * M_PI * 1. / 360.;
-    edge = radius[0] * factor;  //1 degree
+    edge = radius * factor;  //1 degree
   } else {
     constexpr double factor = 2 * M_PI * 1. / 288.;
-    edge = radius[0] * factor;  //1.25 degrees
+    edge = radius * factor;  //1.25 degrees
   }
 
   double scaleFactor = refEdge_ / edge;  //assume reference 3cm of edge
 
   return scaleFactor;
 }
 
-double HGCalSciNoiseMap::scaleBySipmArea(const HGCScintillatorDetId& cellId, const double& radius) {
+double HGCalSciNoiseMap::scaleBySipmArea(const HGCScintillatorDetId& cellId, const double radius) {
   if (sipmMap_.empty())
     return 1.;
 
@@ -83,20 +83,3 @@ double HGCalSciNoiseMap::scaleBySipmArea(const HGCScintillatorDetId& cellId, con
   else
     return 1.;
 }
-
-radiiVec HGCalSciNoiseMap::computeRadius(const HGCScintillatorDetId& cellId) {
-  GlobalPoint global = geom()->getPosition(cellId);
-
-  double radius2 = std::pow(global.x(), 2) + std::pow(global.y(), 2);  //in cm
-  double radius4 = std::pow(radius2, 2);
-  double radius = sqrt(radius2);
-  double radius3 = radius2 * radius;
-
-  double radius_m100 = radius - 100;
-  double radius_m100_2 = std::pow(radius_m100, 2);
-  double radius_m100_3 = radius_m100_2 * radius_m100;
-  double radius_m100_4 = std::pow(radius_m100_2, 2);
-
-  radiiVec radii{{radius, radius2, radius3, radius4, radius_m100, radius_m100_2, radius_m100_3, radius_m100_4}};
-  return radii;
-}