Merge pull request #26811 from deguio/HGCHEbackDigitizer

improvements to the HGCHEback digitizer and implementation of a end-of-life scenario for the HGC scintillator section

L1Trigger/L1THGCal/python/hgcalVFEProducer_cfi.py

@@ -44,10 +44,11 @@
                      linnBits = cms.uint32(16),
                      siliconCellLSB_fC =  cms.double( triggerCellLsbBeforeCompression*(2**triggerCellTruncationBits) ),
                      scintillatorCellLSB_MIP = cms.double(float(adcSaturationBH_MIP.value())/(2**float(adcNbitsBH.value()))),
+                     noiseSilicon = cms.PSet(),
+                     noiseScintillator = cms.PSet(),
                      # cell thresholds before TC sums
                      # Cut at 3sigma of the noise
-                     thresholdsSilicon = cms.vdouble([3.*x for x in digiparam.HGCAL_noise_fC.values.value()]),
-                     thresholdScintillator = cms.double(3.*digiparam.HGCAL_noise_MIP.value.value()),
+                     noiseThreshold = cms.double(3), # in units of sigmas of the noise
                      # Floating point compression
                      exponentBits = cms.uint32(4),
                      mantissaBits = cms.uint32(4),
@@ -59,6 +60,13 @@
                      thickCorr = cms.double(thicknessCorrection_200)
                      )
 
+# isolate these refs in case they aren't available in some other WF
+from Configuration.Eras.Modifier_phase2_hgcal_cff import phase2_hgcal
+phase2_hgcal.toModify(vfe_proc,
+    noiseSilicon = cms.PSet(refToPSet_ = cms.string("HGCAL_noise_fC")),
+    noiseScintillator = cms.PSet(refToPSet_ = cms.string("HGCAL_noise_heback")),
+)
+
 hgcalVFEProducer = cms.EDProducer(
         "HGCalVFEProducer",
         eeDigis = cms.InputTag('simHGCalUnsuppressedDigis:EE'),

L1Trigger/L1THGCal/src/veryfrontend/HGCalVFESummationImpl.cc

@@ -3,19 +3,26 @@
 HGCalVFESummationImpl::HGCalVFESummationImpl(const edm::ParameterSet& conf)
     : thickness_corrections_(conf.getParameter<std::vector<double>>("ThicknessCorrections")),
       lsb_silicon_fC_(conf.getParameter<double>("siliconCellLSB_fC")),
-      lsb_scintillator_MIP_(conf.getParameter<double>("scintillatorCellLSB_MIP")),
-      thresholds_silicon_(conf.getParameter<std::vector<double>>("thresholdsSilicon")),
-      threshold_scintillator_(conf.getParameter<double>("thresholdScintillator")) {
+      lsb_scintillator_MIP_(conf.getParameter<double>("scintillatorCellLSB_MIP")) {
   const unsigned nThickness = 3;
   if (thickness_corrections_.size() != nThickness) {
     throw cms::Exception("Configuration")
         << thickness_corrections_.size() << " thickness corrections are given instead of " << nThickness
         << " (the number of sensor thicknesses)";
   }
+  thresholds_silicon_ =
+      conf.getParameter<edm::ParameterSet>("noiseSilicon").getParameter<std::vector<double>>("values");
   if (thresholds_silicon_.size() != nThickness) {
     throw cms::Exception("Configuration") << thresholds_silicon_.size() << " silicon thresholds are given instead of "
                                           << nThickness << " (the number of sensor thicknesses)";
   }
+  threshold_scintillator_ = conf.getParameter<edm::ParameterSet>("noiseScintillator").getParameter<double>("noise_MIP");
+  const auto threshold = conf.getParameter<double>("noiseThreshold");
+  std::transform(
+      thresholds_silicon_.begin(), thresholds_silicon_.end(), thresholds_silicon_.begin(), [threshold](auto noise) {
+        return noise * threshold;
+      });
+  threshold_scintillator_ *= threshold;
 }
 
 void HGCalVFESummationImpl::triggerCellSums(const HGCalTriggerGeometryBase& geometry,

RecoLocalCalo/HGCalRecProducers/interface/HGCalCLUEAlgo.h

@@ -42,7 +42,7 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
      fcPerMip_(ps.getParameter<std::vector<double> >("fcPerMip")),
      fcPerEle_(ps.getParameter<double>("fcPerEle")),
      nonAgedNoises_(ps.getParameter<edm::ParameterSet>("noises").getParameter<std::vector<double> >("values")),
-     noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("value")),
+     noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
      initialized_(false),
      cells_(2*(maxlayer+1)),
      numberOfClustersPerLayer_(2*(maxlayer+1),0)
@@ -99,7 +99,11 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
     descNestedNoises.add<std::vector<double> >("values", {});
     iDesc.add<edm::ParameterSetDescription>("noises", descNestedNoises);
     edm::ParameterSetDescription descNestedNoiseMIP;
-    descNestedNoiseMIP.add<double>("value", 0 );
+    descNestedNoiseMIP.add<bool>("scaleByDose", false );
+    iDesc.add<edm::ParameterSetDescription>("scaleByDose", descNestedNoiseMIP);
+    descNestedNoiseMIP.add<std::string>("doseMap", "" );
+    iDesc.add<edm::ParameterSetDescription>("doseMap", descNestedNoiseMIP);
+    descNestedNoiseMIP.add<double>("noise_MIP", 1./100. );
     iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
   }
 

RecoLocalCalo/HGCalRecProducers/interface/HGCalImagingAlgo.h

@@ -48,7 +48,7 @@ class HGCalImagingAlgo : public HGCalClusteringAlgoBase
      fcPerMip_(ps.getParameter<std::vector<double> >("fcPerMip")),
      fcPerEle_(ps.getParameter<double>("fcPerEle")),
      nonAgedNoises_(ps.getParameter<edm::ParameterSet>("noises").getParameter<std::vector<double> >("values")),
-     noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("value")),
+     noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
      initialized_(false),
      points_(2*(maxlayer+1)),
      minpos_(2*(maxlayer+1),{ {0.0f,0.0f} }),
@@ -114,7 +114,11 @@ static void fillPSetDescription(edm::ParameterSetDescription& iDesc) {
   descNestedNoises.add<std::vector<double> >("values", {});
   iDesc.add<edm::ParameterSetDescription>("noises", descNestedNoises);
   edm::ParameterSetDescription descNestedNoiseMIP;
-  descNestedNoiseMIP.add<double>("value", 0 );
+  descNestedNoiseMIP.add<bool>("scaleByDose", false );
+  iDesc.add<edm::ParameterSetDescription>("scaleByDose", descNestedNoiseMIP);
+  descNestedNoiseMIP.add<std::string>("doseMap", "" );
+  iDesc.add<edm::ParameterSetDescription>("doseMap", descNestedNoiseMIP);
+  descNestedNoiseMIP.add<double>("noise_MIP", 1./100. );
   iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
 }
 

RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.cc

@@ -52,7 +52,7 @@ HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet&ps) :
     hgcEE_cce_ = ps.getParameter<edm::ParameterSet>("HGCEE_cce").getParameter< std::vector<double> > ("values");
     hgcHEF_noise_fC_ = ps.getParameter<edm::ParameterSet>("HGCHEF_noise_fC").getParameter < std::vector<double> > ("values");
     hgcHEF_cce_ = ps.getParameter<edm::ParameterSet>("HGCHEF_cce").getParameter< std::vector<double> > ("values");
-    hgcHEB_noise_MIP_ = ps.getParameter<edm::ParameterSet>("HGCHEB_noise_MIP").getParameter<double>("value");
+    hgcHEB_noise_MIP_ = ps.getParameter<edm::ParameterSet>("HGCHEB_noise_MIP").getParameter<double>("noise_MIP");
 
     // don't produce rechit if detid is a ghost one
     rangeMatch_ = ps.getParameter<uint32_t>("rangeMatch");

RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py

@@ -221,7 +221,7 @@
     HGCEE_cce = hgceeDigitizer.digiCfg.chargeCollectionEfficiencies,
     HGCHEF_noise_fC = hgchefrontDigitizer.digiCfg.noise_fC,
     HGCHEF_cce = hgchefrontDigitizer.digiCfg.chargeCollectionEfficiencies,
-    HGCHEB_noise_MIP = hgchebackDigitizer.digiCfg.noise_MIP,
+    HGCHEB_noise_MIP = hgchebackDigitizer.digiCfg.noise,
     HGCHFNose_noise_fC = hfnoseDigitizer.digiCfg.noise_fC,
     HGCHFNose_cce = hfnoseDigitizer.digiCfg.chargeCollectionEfficiencies,
 

RecoLocalCalo/HGCalRecProducers/python/hgcalLayerClusters_cff.py

@@ -16,4 +16,4 @@
 hgcalLayerClusters.plugin.thicknessCorrection = cms.vdouble(HGCalRecHit.thicknessCorrection)
 hgcalLayerClusters.plugin.fcPerEle = cms.double(fC_per_ele)
 hgcalLayerClusters.plugin.noises = cms.PSet(refToPSet_ = cms.string('HGCAL_noises'))
-hgcalLayerClusters.plugin.noiseMip = hgchebackDigitizer.digiCfg.noise_MIP
+hgcalLayerClusters.plugin.noiseMip = hgchebackDigitizer.digiCfg.noise

SimCalorimetry/HGCalSimProducers/data/doseParams_3000fb_fluka-3.5.15.9.txt

@@ -0,0 +1,14 @@
+9  3.985908 -0.0233164 0.000479769 -8.90993e-06 5.74269e-08	  16.460299 -0.0282926 0.000176 -1.24226e-06 3.33661e-09
+10 3.865145 -0.0157737 1.4102e-05 -2.59944e-08 5.18596e-09	  16.074041 -0.0142934 -2.3325e-05 -9.91417e-08 1.05114e-09
+11 3.797408 -0.0110597 -0.000150187 1.70552e-06 -9.14724e-10	15.856707 -0.00698478 -0.000129121 5.11442e-07 -1.75663e-10
+12 3.851879 -0.00268387 -0.000502753 6.91762e-06 -2.62058e-08	15.514695 0.00534509 -0.000304267 1.51454e-06 -2.17969e-09
+13 3.637403 -0.0102584 -0.000280915 4.39397e-06 -1.67187e-08	15.235543 0.0140316 -0.000430826 2.25819e-06 -3.69355e-09
+14 3.502875 -0.00930028 -0.000349823 5.39975e-06 -2.10104e-08	15.364650 0.00885852 -0.000386321 2.09248e-06 -3.47966e-09
+15 3.408253 -0.00713129 -0.000358892 4.95141e-06 -1.79598e-08	15.591748 0.00056066 -0.000302279 1.70913e-06 -2.84797e-09
+16 3.308403 -0.010831 -0.000302115 4.5388e-06 -1.67407e-08	  15.353311 0.00786759 -0.000395923 2.15057e-06 -3.55091e-09
+17 3.274672 -0.00658889 -0.000350272 4.3462e-06 -1.44797e-08	15.555453 0.00110278 -0.000332765 1.8694e-06 -3.08188e-09
+18 3.083276 -0.0102271 -0.00028511 3.5795e-06 -1.08933e-08	  15.822784 -0.00815451 -0.000231267 1.35734e-06 -2.15243e-09
+19 3.125504 -0.0073409 -0.000282136 3.16049e-06 -9.34773e-09	15.909287 -0.0112805 -0.000205724 1.24919e-06 -1.99359e-09
+20 2.995494 -0.0193198 -6.97825e-06 6.81663e-07 -1.76918e-09	16.115456 -0.0185292 -0.000131819 9.26615e-07 -1.49273e-09
+21 2.973583 -0.00764613 -0.000184954 1.73145e-06 -4.03636e-09	16.356612 -0.0288606 -1.22903e-05 3.72635e-07 -6.06651e-10
+22 2.821734 -0.0158869 1.38538e-06 3.83847e-07 -1.34869e-09	  16.758791 -0.0481384 0.000220859 -6.4673e-07 8.8165e-10

SimCalorimetry/HGCalSimProducers/interface/HGCHEbackDigitizer.h

@@ -3,6 +3,44 @@
 
 #include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h"
 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
+#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
+
+
+class HGCHEbackSignalScaler
+{
+  public:
+
+    struct DoseParameters {
+      DoseParameters(): a_(0.), b_(0.), c_(0.), d_(0.), e_(0.),
+                        f_(0.), g_(0.), h_(0.), i_(0.), j_(0.) {}
+      float a_, b_, c_, d_, e_, f_, g_, h_, i_, j_;
+    };
+
+    HGCHEbackSignalScaler() {};
+    ~HGCHEbackSignalScaler() {};
+
+    void setGeometry(const CaloSubdetectorGeometry*);
+    void setDoseMap(const std::string&);
+
+    float scaleByArea(const HGCScintillatorDetId&, const std::array<double, 8>&);
+    std::pair<float, float> scaleByDose(const HGCScintillatorDetId&, const std::array<double, 8>&);
+    double getDoseValue(const int, const std::array<double, 8>&);
+    double getFluenceValue(const int, const std::array<double, 8>&);
+    std::array<double, 8> computeRadius(const HGCScintillatorDetId&);
+
+  private:
+    std::map<int, DoseParameters> readDosePars(const std::string&);
+
+    const HGCalGeometry* hgcalGeom_;
+    std::map<int, DoseParameters> doseMap_;
+    static constexpr double greyToKrad_ = 0.1;
+    const float refEdge_ = 3; //3 cm
+
+    bool verbose_ = false;
+
+};
+
+
 
 class HGCHEbackDigitizer : public HGCDigitizerBase<HGCalDataFrame>
 {
@@ -17,11 +55,24 @@ class HGCHEbackDigitizer : public HGCDigitizerBase<HGCalDataFrame>
  private:
 
   //calice-like digitization parameters
+  uint32_t algo_;
+  bool scaleByArea_, scaleByDose_, calibDigis_;
   float keV2MIP_, noise_MIP_;
   float nPEperMIP_, nTotalPE_, xTalk_, sdPixels_;
-  void runCaliceLikeDigitizer(std::unique_ptr<HGCalDigiCollection> &digiColl,hgc::HGCSimHitDataAccumulator &simData, 
+  std::string doseMapFile_;
+  HGCHEbackSignalScaler scal_;
+
+  void runEmptyDigitizer(std::unique_ptr<HGCalDigiCollection> &digiColl,hgc::HGCSimHitDataAccumulator &simData,
+                         const CaloSubdetectorGeometry* theGeom, const std::unordered_set<DetId>& validIds,
+                         CLHEP::HepRandomEngine* engine);
+
+  void runRealisticDigitizer(std::unique_ptr<HGCalDigiCollection> &digiColl,hgc::HGCSimHitDataAccumulator &simData,
+                             const CaloSubdetectorGeometry* theGeom, const std::unordered_set<DetId>& validIds,
+                             CLHEP::HepRandomEngine* engine);
+
+  void runCaliceLikeDigitizer(std::unique_ptr<HGCalDigiCollection> &digiColl,hgc::HGCSimHitDataAccumulator &simData,
 			      const CaloSubdetectorGeometry* theGeom, const std::unordered_set<DetId>& validIds,
 			      CLHEP::HepRandomEngine* engine);
 };
 
-#endif 
+#endif

SimCalorimetry/HGCalSimProducers/python/hgcalDigitizer_cfi.py

@@ -12,8 +12,10 @@
     values = cms.vdouble( [x*fC_per_ele for x in nonAgedNoises] ), #100,200,300 um
     )
 
-HGCAL_noise_MIP = cms.PSet(
-    value = cms.double(1.0/7.0)
+HGCAL_noise_heback = cms.PSet(
+    scaleByDose = cms.bool(False),
+    doseMap = cms.string(""), #empty dose map at begin-of-life
+    noise_MIP = cms.double(1./100.)
     )
 
 HGCAL_chargeCollectionEfficiencies = cms.PSet(
@@ -154,7 +156,7 @@
     )
 
 
-# HCAL back (CALICE-like version, no pulse shape)
+# HCAL back
 hgchebackDigitizer = cms.PSet(
     accumulatorType   = cms.string("HGCDigiProducer"),
     hitCollection = cms.string("HcalHits"),
@@ -171,22 +173,26 @@
     useAllChannels    = cms.bool(True),
     verbosity         = cms.untracked.uint32(0),
     digiCfg = cms.PSet(
-        keV2MIP           = cms.double(1./616.0),
-        noise_MIP         = cms.PSet(refToPSet_ = cms.string("HGCAL_noise_MIP")),
+        #0 empty digitizer, 1 calice digitizer, 2 realistic digitizer
+        algo          = cms.uint32(2),
+        scaleByArea   = cms.bool(True),
+        noise         = cms.PSet(refToPSet_ = cms.string("HGCAL_noise_heback")), #scales both for scint raddam and sipm dark current
+        calibDigis    = cms.bool(True),
+        keV2MIP       = cms.double(1./500.0),
         doTimeSamples = cms.bool(False),
-        nPEperMIP = cms.double(11.0),
-        nTotalPE  = cms.double(1156), #1156 pixels => saturation ~600MIP
-        xTalk     = cms.double(0.25),
+        nPEperMIP = cms.double(21.0),
+        nTotalPE  = cms.double(7500),
+        xTalk     = cms.double(0.01),
         sdPixels  = cms.double(1e-6), # this is additional photostatistics noise (as implemented), not sure why it's here...
         feCfg   = cms.PSet(
             # 0 only ADC, 1 ADC with pulse shape, 2 ADC+TDC with pulse shape
             fwVersion       = cms.uint32(0),
-            # n bits for the ADC
-            adcNbits        = cms.uint32(12),
+            # n bits for the ADC (same as the silicon ROC)
+            adcNbits        = cms.uint32(10),
             # ADC saturation : in this case we use the same variable but fC=MIP
-            adcSaturation_fC = cms.double(1024.0),
+            adcSaturation_fC = cms.double(68.5), #value chosen to have 1MIP at 15ADC
             # threshold for digi production : in this case we use the same variable but fC=MIP
-            adcThreshold_fC = cms.double(0.50),
+            adcThreshold_fC = cms.double(0.5),
             thresholdFollowsMIP = cms.bool(False)
             )
         )
@@ -273,8 +279,10 @@ def HGCal_setEndOfLifeNoise(process):
     process.HGCAL_chargeCollectionEfficiencies = cms.PSet(
         values = cms.vdouble(endOfLifeCCEs)
         )
-    process.HGCAL_noise_MIP = cms.PSet(
-        value = cms.double( 1.0/5.0 )
+    process.HGCAL_noise_heback = cms.PSet(
+        scaleByDose = cms.bool(True),
+        doseMap = cms.string("SimCalorimetry/HGCalSimProducers/data/doseParams_3000fb_fluka-3.5.15.9.txt"),
+        noise_MIP = cms.double(1./5.) #uses noise map
         )
     process.HGCAL_noises = cms.PSet(
         values = cms.vdouble([x for x in endOfLifeNoises])
@@ -285,9 +293,11 @@ def HGCal_disableNoise(process):
     process.HGCAL_noise_fC = cms.PSet(
         values = cms.vdouble(0,0,0), #100,200,300 um
     )
-    process.HGCAL_noise_MIP = cms.PSet(
-        value = cms.double(0)
-    )
+    process.HGCAL_noise_heback = cms.PSet(
+        scaleByDose = cms.bool(False),
+        doseMap = cms.string(""),
+        noise_MIP = cms.double(0.) #zero noise
+        )
     process.HGCAL_noises = cms.PSet(
         values = cms.vdouble(0,0,0)
     )