Merge pull request #34252 from TaeunKwon/FixingSaturation

Updates of the TP reconstruction algorithm to fix the TP energy saturation effect

Author: cmsbuild

SimCalorimetry/HcalTrigPrimAlgos/interface/HcalTriggerPrimitiveAlgo.h

@@ -89,6 +89,7 @@ class HcalTriggerPrimitiveAlgo {
   void setNumFilterPresamplesHEQIE11(int presamples) { numberOfFilterPresamplesHEQIE11_ = presamples; }
 
   void setUpgradeFlags(bool hb, bool he, bool hf);
+  void setFixSaturationFlag(bool fix_saturation);
   void overrideParameters(const edm::ParameterSet& ps);
 
 private:
@@ -110,7 +111,10 @@ class HcalTriggerPrimitiveAlgo {
   /// adds the actual digis
   void analyze(IntegerCaloSamples& samples, HcalTriggerPrimitiveDigi& result);
   // 2017 and later: QIE11
-  void analyzeQIE11(IntegerCaloSamples& samples, HcalTriggerPrimitiveDigi& result, const HcalFinegrainBit& fg_algo);
+  void analyzeQIE11(IntegerCaloSamples& samples,
+                    std::vector<bool> sample_saturation,
+                    HcalTriggerPrimitiveDigi& result,
+                    const HcalFinegrainBit& fg_algo);
   // Version 0: RCT
   void analyzeHF(IntegerCaloSamples& samples, HcalTriggerPrimitiveDigi& result, const int hf_lumi_shift);
   // Version 1: 1x1
@@ -163,6 +167,9 @@ class HcalTriggerPrimitiveAlgo {
   typedef std::map<HcalTrigTowerDetId, IntegerCaloSamples> SumMap;
   SumMap theSumMap;
 
+  typedef std::map<HcalTrigTowerDetId, std::vector<bool>> SatMap;
+  SatMap theSatMap;
+
   struct HFDetails {
     IntegerCaloSamples long_fiber;
     IntegerCaloSamples short_fiber;
@@ -210,6 +217,8 @@ class HcalTriggerPrimitiveAlgo {
   bool upgrade_he_ = false;
   bool upgrade_hf_ = false;
 
+  bool fix_saturation_ = false;
+
   edm::ParameterSet override_parameters_;
 
   bool override_adc_hf_ = false;
@@ -250,6 +259,7 @@ void HcalTriggerPrimitiveAlgo::run(const HcalTPGCoder* incoder,
   conditions_ = conditions;
 
   theSumMap.clear();
+  theSatMap.clear();
   theTowerMapFGSum.clear();
   HF_Veto.clear();
   fgMap_.clear();
@@ -293,13 +303,18 @@ void HcalTriggerPrimitiveAlgo::run(const HcalTPGCoder* incoder,
       if (fgMap_.find(item.first) != fgMap_.end()) {
         analyze(item.second, result.back());
       } else if (fgUpgradeMap_.find(item.first) != fgUpgradeMap_.end()) {
-        analyzeQIE11(item.second, result.back(), fg_algo);
+        SatMap::iterator item_sat = theSatMap.find(detId);
+        if (item_sat == theSatMap.end())
+          analyzeQIE11(item.second, std::vector<bool>(), result.back(), fg_algo);
+        else
+          analyzeQIE11(item.second, item_sat->second, result.back(), fg_algo);
       }
     }
   }
 
   // Free up some memory
   theSumMap.clear();
+  theSatMap.clear();
   theTowerMapFGSum.clear();
   HF_Veto.clear();
   fgMap_.clear();

SimCalorimetry/HcalTrigPrimAlgos/src/HcalTriggerPrimitiveAlgo.cc

@@ -76,6 +76,8 @@ void HcalTriggerPrimitiveAlgo::setUpgradeFlags(bool hb, bool he, bool hf) {
   upgrade_hf_ = hf;
 }
 
+void HcalTriggerPrimitiveAlgo::setFixSaturationFlag(bool fix_saturation) { fix_saturation_ = fix_saturation; }
+
 void HcalTriggerPrimitiveAlgo::overrideParameters(const edm::ParameterSet& ps) {
   override_parameters_ = ps;
 
@@ -278,6 +280,7 @@ void HcalTriggerPrimitiveAlgo::addSignal(const QIE11DataFrame& frame) {
 void HcalTriggerPrimitiveAlgo::addSignal(const IntegerCaloSamples& samples) {
   HcalTrigTowerDetId id(samples.id());
   SumMap::iterator itr = theSumMap.find(id);
+
   if (itr == theSumMap.end()) {
     theSumMap.insert(std::make_pair(id, samples));
   } else {
@@ -286,6 +289,29 @@ void HcalTriggerPrimitiveAlgo::addSignal(const IntegerCaloSamples& samples) {
       (itr->second)[i] += samples[i];
     }
   }
+
+  // if fix_saturation == true, keep track of tower with saturated input LUT
+  if (fix_saturation_) {
+    SatMap::iterator itr_sat = theSatMap.find(id);
+
+    assert((itr == theSumMap.end()) == (itr_sat == theSatMap.end()));
+
+    if (itr_sat == theSatMap.end()) {
+      vector<bool> check_sat;
+      for (int i = 0; i < samples.size(); ++i) {
+        if (!(samples[i] < QIE11_LINEARIZATION_ET)) {
+          check_sat.push_back(true);
+        } else
+          check_sat.push_back(false);
+      }
+      theSatMap.insert(std::make_pair(id, check_sat));
+    } else {
+      for (int i = 0; i < samples.size(); ++i) {
+        if (!(samples[i] < QIE11_LINEARIZATION_ET))
+          (itr_sat->second)[i] = true;
+      }
+    }
+  }
 }
 
 void HcalTriggerPrimitiveAlgo::analyze(IntegerCaloSamples& samples, HcalTriggerPrimitiveDigi& result) {
@@ -371,6 +397,7 @@ void HcalTriggerPrimitiveAlgo::analyze(IntegerCaloSamples& samples, HcalTriggerP
 }
 
 void HcalTriggerPrimitiveAlgo::analyzeQIE11(IntegerCaloSamples& samples,
+                                            vector<bool> sample_saturation,
                                             HcalTriggerPrimitiveDigi& result,
                                             const HcalFinegrainBit& fg_algo) {
   HcalDetId detId(samples.id());
@@ -398,15 +425,26 @@ void HcalTriggerPrimitiveAlgo::analyzeQIE11(IntegerCaloSamples& samples,
   IntegerCaloSamples sum(samples.id(), samples.size());
 
   std::vector<HcalTrigTowerDetId> ids = theTrigTowerGeometry->towerIds(detId);
+
+  // keep track of tower with saturated energy and force the total TP saturated
+  bool force_saturation[samples.size()];
+  for (int i = 0; i < samples.size(); i++) {
+    force_saturation[i] = false;
+  }
+
   //slide algo window
   for (unsigned int ibin = 0; ibin < dgSamples - shrink; ++ibin) {
     int algosumvalue = 0;
+    bool check_sat = false;
     for (unsigned int i = 0; i < filterSamples; i++) {
       //add up value * scale factor
       // In addition, divide by two in the 10 degree phi segmentation region
       // to mimic 5 degree segmentation for the trigger
       unsigned int sample = samples[ibin + i];
-      if (sample > QIE11_MAX_LINEARIZATION_ET)
+
+      if (fix_saturation_ && (sample_saturation.size() > ibin + i))
+        check_sat = (sample_saturation[ibin + i] | (sample > QIE11_MAX_LINEARIZATION_ET));
+      else if (sample > QIE11_MAX_LINEARIZATION_ET)
         sample = QIE11_MAX_LINEARIZATION_ET;
 
       // Usually use a segmentation factor of 1.0 but for ieta >= 21 use 0.5
@@ -426,6 +464,9 @@ void HcalTriggerPrimitiveAlgo::analyzeQIE11(IntegerCaloSamples& samples,
     //else if (algosumvalue>QIE11_LINEARIZATION_ET) sum[ibin]=QIE11_LINEARIZATION_ET;
     else
       sum[ibin] = algosumvalue;  //assign value to sum[]
+
+    if (check_sat)
+      force_saturation[ibin] = true;
   }
 
   std::vector<int> finegrain(tpSamples, false);
@@ -448,6 +489,8 @@ void HcalTriggerPrimitiveAlgo::analyzeQIE11(IntegerCaloSamples& samples,
 
     if (isPeak) {
       output[ibin] = std::min<unsigned int>(sum[idx], QIE11_MAX_LINEARIZATION_ET);
+      if (fix_saturation_ && force_saturation[idx])
+        output[ibin] = QIE11_MAX_LINEARIZATION_ET;
     } else {
       // Not a peak
       output[ibin] = 0;

SimCalorimetry/HcalTrigPrimProducers/python/hcaltpdigi_cfi.py

@@ -4,6 +4,7 @@
 from Configuration.Eras.Modifier_run2_HE_2017_cff import run2_HE_2017
 from Configuration.Eras.Modifier_run2_HF_2017_cff import run2_HF_2017
 from Configuration.Eras.Modifier_run3_HB_cff import run3_HB
+from Configuration.Eras.Modifier_run3_common_cff import run3_common
 
 LSParameter =cms.untracked.PSet(
 HcalFeatureHFEMBit= cms.bool(False),
@@ -66,6 +67,8 @@
     upgradeHB = cms.bool(False),
     upgradeHE = cms.bool(False),
 
+    applySaturationFix = cms.bool(False), # Apply the TP energy saturation fix for Peak Finder Algorithm only for Run3 
+
     # parameters = cms.untracked.PSet(
     #     FGVersionHBHE=cms.uint32(0),
     #     TDCMask=cms.uint64(0xFFFFFFFFFFFFFFFF),
@@ -98,4 +101,5 @@
 )
 run2_HF_2017.toModify(tpScales.HF, NCTShift=cms.int32(2))
 run3_HB.toModify(simHcalTriggerPrimitiveDigis, upgradeHB=cms.bool(True))
+run3_common.toModify(simHcalTriggerPrimitiveDigis, applySaturationFix=cms.bool(True))
 run3_HB.toModify(tpScales.HBHE, LSBQIE11Overlap=cms.double(1/16.))

SimCalorimetry/HcalTrigPrimProducers/src/HcalTrigPrimDigiProducer.cc

@@ -50,6 +50,7 @@ HcalTrigPrimDigiProducer::HcalTrigPrimDigiProducer(const edm::ParameterSet& ps)
     theAlgo_.overrideParameters(pset);
   }
   theAlgo_.setUpgradeFlags(upgrades[0], upgrades[1], upgrades[2]);
+  theAlgo_.setFixSaturationFlag(ps.getParameter<bool>("applySaturationFix"));
 
   HFEMB_ = false;
   if (ps.exists("LSConfig")) {