Merge pull request cms-sw#31730 from apsallid/cluetweaksforregemfactors

[HGCAL] CLUE and regional factors

Author: cmsbuild

RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.cc

@@ -44,8 +44,12 @@ void HGCalCLUEAlgoT<T>::populate(const HGCRecHitCollection& hits) {
     if (dependSensor_) {
       int thickness_index = rhtools_.getSiThickIndex(detid);
       if (thickness_index == -1)
-        thickness_index = 3;
+        thickness_index = maxNumberOfThickIndices_;
+
       double storedThreshold = thresholds_[layerOnSide][thickness_index];
+      if (detid.det() == DetId::HGCalHSi || detid.subdetId() == HGCHEF) {
+        storedThreshold = thresholds_[layerOnSide][thickness_index + deltasi_index_regemfac_];
+      }
       sigmaNoise = v_sigmaNoise_[layerOnSide][thickness_index];
 
       if (hgrh.energy() < storedThreshold)
@@ -507,8 +511,9 @@ void HGCalCLUEAlgoT<T>::computeThreshold() {
   // To support the TDR geometry and also the post-TDR one (v9 onwards), we
   // need to change the logic of the vectors containing signal to noise and
   // thresholds. The first 3 indices will keep on addressing the different
-  // thicknesses of the Silicon detectors, while the last one, number 3 (the
-  // fourth) will address the Scintillators. This change will support both
+  // thicknesses of the Silicon detectors in CE_E , the next 3 indices will address
+  // the thicknesses of the Silicon detectors in CE_H, while the last one, number 6 (the
+  // seventh) will address the Scintillators. This change will support both
   // geometries at the same time.
 
   if (initialized_)
@@ -517,13 +522,13 @@ void HGCalCLUEAlgoT<T>::computeThreshold() {
   initialized_ = true;
 
   std::vector<double> dummy;
-  const unsigned maxNumberOfThickIndices = 3;
-  dummy.resize(maxNumberOfThickIndices + !isNose_, 0);  // +1 to accomodate for the Scintillators
+
+  dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0);  // +1 to accomodate for the Scintillators
   thresholds_.resize(maxlayer_, dummy);
   v_sigmaNoise_.resize(maxlayer_, dummy);
 
   for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
-    for (unsigned ithick = 0; ithick < maxNumberOfThickIndices; ++ithick) {
+    for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) {
       float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
                          (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
       thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
@@ -535,9 +540,9 @@ void HGCalCLUEAlgoT<T>::computeThreshold() {
     }
 
     if (!isNose_) {
-      float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer];
-      thresholds_[ilayer - 1][maxNumberOfThickIndices] = ecut_ * scintillators_sigmaNoise;
-      v_sigmaNoise_[ilayer - 1][maxNumberOfThickIndices] = scintillators_sigmaNoise;
+      float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer] / sciThicknessCorrection_;
+      thresholds_[ilayer - 1][maxNumberOfThickIndices_] = ecut_ * scintillators_sigmaNoise;
+      v_sigmaNoise_[ilayer - 1][maxNumberOfThickIndices_] = scintillators_sigmaNoise;
       LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " noiseMip: " << noiseMip_
                                 << " scintillators_sigmaNoise: " << scintillators_sigmaNoise << "\n";
     }

RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.h

@@ -40,12 +40,20 @@ class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
         dependSensor_(ps.getParameter<bool>("dependSensor")),
         dEdXweights_(ps.getParameter<std::vector<double>>("dEdXweights")),
         thicknessCorrection_(ps.getParameter<std::vector<double>>("thicknessCorrection")),
+        sciThicknessCorrection_(ps.getParameter<double>("sciThicknessCorrection")),
+        deltasi_index_regemfac_(ps.getParameter<int>("deltasi_index_regemfac")),
+        maxNumberOfThickIndices_(ps.getParameter<unsigned>("maxNumberOfThickIndices")),
         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>("noise_MIP")),
         use2x2_(ps.getParameter<bool>("use2x2")),
-        initialized_(false) {}
+        initialized_(false) {
+    // repeat same noises for CE-H as well
+    if (!isNose_) {
+      nonAgedNoises_.insert(std::end(nonAgedNoises_), std::begin(nonAgedNoises_), std::end(nonAgedNoises_));
+    }
+  }
 
   ~HGCalCLUEAlgoT() override {}
 
@@ -94,6 +102,9 @@ class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
     iDesc.addUntracked<unsigned int>("verbosity", 3);
     iDesc.add<std::vector<double>>("dEdXweights", {});
     iDesc.add<std::vector<double>>("thicknessCorrection", {});
+    iDesc.add<double>("sciThicknessCorrection", 0.9);
+    iDesc.add<int>("deltasi_index_regemfac", 3);
+    iDesc.add<unsigned>("maxNumberOfThickIndices", 6);
     iDesc.add<std::vector<double>>("fcPerMip", {});
     iDesc.add<double>("fcPerEle", 0.0);
     edm::ParameterSetDescription descNestedNoises;
@@ -131,6 +142,9 @@ class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
   bool dependSensor_;
   std::vector<double> dEdXweights_;
   std::vector<double> thicknessCorrection_;
+  double sciThicknessCorrection_;
+  int deltasi_index_regemfac_;
+  unsigned maxNumberOfThickIndices_;
   std::vector<double> fcPerMip_;
   double fcPerEle_;
   std::vector<double> nonAgedNoises_;

RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.cc

@@ -59,7 +59,8 @@ HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet& ps) :
     rcorr_.push_back(1.0 / corr);
   }
   // here for scintillator
-  rcorrscint_ = ps.getParameter<double>("sciThicknessCorrection");
+  rcorrscint_ = 1.0 / ps.getParameter<double>("sciThicknessCorrection");
+
   //This is for the index position in CE_H silicon thickness cases
   deltasi_index_regemfac_ = ps.getParameter<int>("deltasi_index_regemfac");
   const auto& rcorrnose = ps.getParameter<std::vector<double> >("thicknessNoseCorrection");
@@ -185,7 +186,7 @@ bool HGCalRecHitWorkerSimple::run(const edm::Event& evt,
       break;
     case hgcbh:
       rechitMaker_->setADCToGeVConstant(float(hgchebUncalib2GeV_));
-      sigmaNoiseGeV = 1e-3 * hgcHEB_noise_MIP_ * weights_[layer];
+      sigmaNoiseGeV = 1e-3 * hgcHEB_noise_MIP_ * weights_[layer] * rcorrscint_;
       break;
     case hgchfnose:
       rechitMaker_->setADCToGeVConstant(float(hgchfnoseUncalib2GeV_));

RecoLocalCalo/HGCalRecProducers/python/hgcalLayerClusters_cff.py

@@ -12,8 +12,10 @@
 
 hgcalLayerClusters.timeOffset = hgceeDigitizer.tofDelay
 hgcalLayerClusters.plugin.dEdXweights = cms.vdouble(dEdX.weights)
-hgcalLayerClusters.plugin.fcPerMip = cms.vdouble(HGCalUncalibRecHit.HGCEEConfig.fCPerMIP)
+hgcalLayerClusters.plugin.fcPerMip = cms.vdouble(HGCalUncalibRecHit.HGCEEConfig.fCPerMIP + HGCalUncalibRecHit.HGCHEFConfig.fCPerMIP)
 hgcalLayerClusters.plugin.thicknessCorrection = cms.vdouble(HGCalRecHit.thicknessCorrection)
+hgcalLayerClusters.plugin.sciThicknessCorrection = HGCalRecHit.sciThicknessCorrection
+hgcalLayerClusters.plugin.deltasi_index_regemfac = HGCalRecHit.deltasi_index_regemfac
 hgcalLayerClusters.plugin.fcPerEle = cms.double(fC_per_ele)
 hgcalLayerClusters.plugin.noises = cms.PSet(refToPSet_ = cms.string('HGCAL_noises'))
 hgcalLayerClusters.plugin.noiseMip = hgchebackDigitizer.digiCfg.noise