Merge pull request #13835 from lgray/reworked_rechit_scheme

HGCal Reworked rechit scheme

RecoLocalCalo/HGCalRecAlgos/interface/HGCalRecHitSimpleAlgo.h

@@ -56,6 +56,5 @@ class HGCalRecHitSimpleAlgo : public HGCalRecHitAbsAlgo {
 private:
   float adcToGeVConstant_;
   bool  adcToGeVConstantIsSet_;
-
 };
 #endif

RecoLocalCalo/HGCalRecAlgos/interface/HGCalUncalibRecHitRecWeightsAlgo.h

@@ -17,20 +17,36 @@
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
 
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+
 template<class C> class HGCalUncalibRecHitRecWeightsAlgo 
 {
  public:
   // destructor
   virtual ~HGCalUncalibRecHitRecWeightsAlgo<C>() { };
 
-  virtual void set_isSiFESim(const bool isSiFE) { isSiFESim_ = isSiFE; }
+  void set_isSiFESim(const bool isSiFE) { isSiFESim_ = isSiFE; }
+
+  void set_ADCLSB(const double adclsb) { adcLSB_ = adclsb; }
+  void set_TDCLSB(const double tdclsb) { tdcLSB_ = tdclsb; }
 
-  virtual void set_ADCLSB(const double adclsb) { adcLSB_ = adclsb; }
-  virtual void set_TDCLSB(const double tdclsb) { tdcLSB_ = tdclsb; }
+  void set_toaLSBToNS(const double lsb2ns) { toaLSBToNS_ = lsb2ns; }
 
-  virtual void set_toaLSBToNS(const double lsb2ns) { toaLSBToNS_ = lsb2ns; }
+  void set_tdcOnsetfC(const double tdcOnset) { tdcOnsetfC_ = tdcOnset; }
 
-  virtual void set_tdcOnsetfC(const double tdcOnset) { tdcOnsetfC_ = tdcOnset; }
+  void set_fCPerMIP(const std::vector<double>& fCPerMIP) { 
+    if( std::any_of(fCPerMIP.cbegin(), 
+                    fCPerMIP.cend(), 
+                    [](double conv){ return conv <= 0.0; }) ) {
+      throw cms::Exception("BadConversionFactor") << "At least one of fCPerMIP is zero!" << std::endl;
+    }
+    fCPerMIP_ = fCPerMIP; 
+  }
+
+  void setGeometry(const HGCalGeometry* geom) { 
+    if ( geom ) ddd_ = &(geom->topology().dddConstants());
+    else ddd_ = nullptr;
+  }
 
   /// Compute parameters
   virtual HGCUncalibratedRecHit makeRecHit( const C& dataFrame ) {
@@ -51,30 +67,37 @@ template<class C> class HGCalUncalibRecHitRecWeightsAlgo
         // need to increment by one so TDC doesn't overlap with ADC last bin
 	amplitude_ = ( std::floor(tdcOnsetfC_/adcLSB_) + 1.0 )* adcLSB_ + double(sample.data()) * tdcLSB_;
 	jitter_    = double(sample.toa()) * toaLSBToNS_;
-	LogDebug("HGCUncalibratedRecHit") << "TDC+: set the energy to: " << amplitude_ << ' ' << sample.data() 
+	LogDebug("HGCUncalibratedRecHit") << "TDC+: set the charge to: " << amplitude_ << ' ' << sample.data() 
                                           << ' ' << tdcLSB_ << std::endl
-                                          << "TDC+: set the jitter to: " << jitter_ << ' ' 
+                                          << "TDC+: set the ToA to: " << jitter_ << ' ' 
                                           << sample.toa() << ' ' << toaLSBToNS_ << ' '
                                           << " flag=" << flag << std::endl;
-      } 
-      else {
+      } else {
 	amplitude_ = double(sample.data()) * adcLSB_;
-	LogDebug("HGCUncalibratedRecHit") << "ADC+: set the energy to: " << amplitude_ << ' ' << sample.data() 
+	LogDebug("HGCUncalibratedRecHit") << "ADC+: set the charge to: " << amplitude_ << ' ' << sample.data() 
                                           << ' ' << adcLSB_ << ' ' << std::endl;
       }
-    }
-    else {
+    } else {
       amplitude_ = double(sample.data()) * adcLSB_;
-      LogDebug("HGCUncalibratedRecHit") << "ADC+: set the energy to: " << amplitude_ << ' ' << sample.data() 
+      LogDebug("HGCUncalibratedRecHit") << "ADC+: set the charge to: " << amplitude_ << ' ' << sample.data() 
                                         << ' ' << adcLSB_ << ' ' << std::endl;
     }
 
+    int thickness = 1;
+    if( ddd_ != nullptr ) {
+      HGCalDetId hid(dataFrame.id());
+      thickness = ddd_->waferTypeL(hid.wafer());
+    }
+    amplitude_ = amplitude_/fCPerMIP_[thickness-1];
+
     LogDebug("HGCUncalibratedRecHit") << "Final uncalibrated amplitude : " << amplitude_ << std::endl;
     return HGCUncalibratedRecHit( dataFrame.id(), amplitude_, pedestal_, jitter_, chi2_, flag);
    }
 
  private:
-   double adcLSB_, tdcLSB_, fCToMIP_, toaLSBToNS_, tdcOnsetfC_;
-   bool   isSiFESim_;
+   double adcLSB_, tdcLSB_, toaLSBToNS_, tdcOnsetfC_;
+   bool   isSiFESim_;  
+   std::vector<double> fCPerMIP_;
+   const HGCalDDDConstants* ddd_;
 };
 #endif

RecoLocalCalo/HGCalRecProducers/BuildFile.xml

@@ -1,5 +1,6 @@
 <use   name="FWCore/MessageLogger"/>
 <use   name="FWCore/Framework"/>
+<use   name="Geometry/HGCalGeometry"/>
 <use   name="clhep"/>
 <export>
   <lib   name="1"/>

RecoLocalCalo/HGCalRecProducers/plugins/BuildFile.xml

@@ -8,6 +8,7 @@
 <use   name="CondFormats/DataRecord"/>
 <use   name="FWCore/MessageLogger"/>
 <use   name="FWCore/MessageService"/>
+<use   name="Geometry/HGCalGeometry"/>
 <library   file="*.cc" name="RecoLocalCaloHGCalRecProducersPlugins">
   <flags   EDM_PLUGIN="1"/>
 </library>

RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.cc

@@ -13,10 +13,12 @@ HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet&ps) :
   constexpr float keV2GeV = 1e-6;
   // HGCee constants 
   HGCEE_keV2DIGI_   =  ps.getParameter<double>("HGCEE_keV2DIGI");
+  HGCEE_fCPerMIP_   =  ps.getParameter<std::vector<double> >("HGCEE_fCPerMIP");
   hgceeUncalib2GeV_ = keV2GeV/HGCEE_keV2DIGI_;
 
   // HGChef constants
   HGCHEF_keV2DIGI_   =  ps.getParameter<double>("HGCHEF_keV2DIGI");
+  HGCHEF_fCPerMIP_    =  ps.getParameter<std::vector<double> >("HGCHEF_fCPerMIP");
   hgchefUncalib2GeV_ = keV2GeV/HGCHEF_keV2DIGI_;
 
   // HGCheb constants
@@ -25,6 +27,13 @@ HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet&ps) :
 }
 
 void HGCalRecHitWorkerSimple::set(const edm::EventSetup& es) {
+  edm::ESHandle<HGCalGeometry> hgceeGeoHandle; 
+  edm::ESHandle<HGCalGeometry> hgchefGeoHandle; 
+  es.get<IdealGeometryRecord>().get("HGCalEESensitive",hgceeGeoHandle); 
+  es.get<IdealGeometryRecord>().get("HGCalHESiliconSensitive",hgchefGeoHandle); 
+  ddds_[0] = &(hgceeGeoHandle->topology().dddConstants());
+  ddds_[1] = &(hgchefGeoHandle->topology().dddConstants());
+  ddds_[2] = nullptr;  
 }
 
 
@@ -34,13 +43,16 @@ HGCalRecHitWorkerSimple::run( const edm::Event & evt,
                               HGCRecHitCollection & result ) {
   DetId detid=uncalibRH.id();  
   uint32_t recoFlag = 0;
+  const std::vector<double>* fCPerMIP = nullptr;
 
   switch( detid.subdetId() ) {
   case HGCEE:
     rechitMaker_->setADCToGeVConstant(float(hgceeUncalib2GeV_) );
+    fCPerMIP = &HGCEE_fCPerMIP_;
     break;
   case HGCHEF:
     rechitMaker_->setADCToGeVConstant(float(hgchefUncalib2GeV_) );
+    fCPerMIP = &HGCHEF_fCPerMIP_;
     break;
   case HGCHEB:
     rechitMaker_->setADCToGeVConstant(float(hgchebUncalib2GeV_) );
@@ -51,8 +63,16 @@ HGCalRecHitWorkerSimple::run( const edm::Event & evt,
   }
 
   // make the rechit and put in the output collection
-  if (recoFlag == 0) {
+  if (recoFlag == 0) {    
     HGCRecHit myrechit( rechitMaker_->makeRecHit(uncalibRH, 0) );    
+    HGCalDetId hid(detid);
+    if( fCPerMIP != nullptr ) {
+      const int thk = ddds_[hid.subdetId()-3]->waferTypeL(hid.wafer());
+      // units out of rechit maker are MIP * (GeV/fC)
+      // so multiple
+      const double new_E = myrechit.energy()*(*fCPerMIP)[thk-1];
+      myrechit.setEnergy(new_E);
+    }    
     result.push_back(myrechit);
   }
 

RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.h

@@ -13,6 +13,8 @@
 
 #include "FWCore/Framework/interface/ESHandle.h"
 
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+
 class HGCalRecHitWorkerSimple : public HGCalRecHitWorkerBaseClass {
  public:
   HGCalRecHitWorkerSimple(const edm::ParameterSet&);
@@ -24,9 +26,13 @@ class HGCalRecHitWorkerSimple : public HGCalRecHitWorkerBaseClass {
  protected:
 
   double HGCEE_keV2DIGI_,  hgceeUncalib2GeV_;
+  std::vector<double> HGCEE_fCPerMIP_;
   double HGCHEF_keV2DIGI_, hgchefUncalib2GeV_;
+  std::vector<double> HGCHEF_fCPerMIP_;
   double HGCHEB_keV2DIGI_, hgchebUncalib2GeV_;
 
+  std::array<const HGCalDDDConstants*, 3> ddds_;
+
   std::vector<int> v_chstatus_;
 
   std::vector<int> v_DB_reco_flags_;

RecoLocalCalo/HGCalRecProducers/plugins/HGCalUncalibRecHitWorkerWeights.cc

@@ -16,6 +16,7 @@ void configureIt(const edm::ParameterSet& conf,
   constexpr char tdcSaturation[]    = "tdcSaturation";
   constexpr char tdcOnset[]         = "tdcOnset";
   constexpr char toaLSB_ns[]        = "toaLSB_ns";
+  constexpr char fCPerMIP[]          = "fCPerMIP";
 
   if( conf.exists(isSiFE) ) {
     maker.set_isSiFESim(conf.getParameter<bool>(isSiFE));
@@ -49,6 +50,13 @@ void configureIt(const edm::ParameterSet& conf,
   } else {
     maker.set_toaLSBToNS(-1.);
   }
+
+  if( conf.exists(fCPerMIP) ) {
+    maker.set_fCPerMIP(conf.getParameter<std::vector<double> >(fCPerMIP));
+  } else {
+    maker.set_fCPerMIP(std::vector<double>({1.0}));
+  }
+
 }
 
 HGCalUncalibRecHitWorkerWeights::HGCalUncalibRecHitWorkerWeights(const edm::ParameterSet&ps) :
@@ -65,7 +73,13 @@ HGCalUncalibRecHitWorkerWeights::HGCalUncalibRecHitWorkerWeights(const edm::Para
 void
 HGCalUncalibRecHitWorkerWeights::set(const edm::EventSetup& es)
 {
-
+  edm::ESHandle<HGCalGeometry> hgceeGeoHandle; 
+  edm::ESHandle<HGCalGeometry> hgchefGeoHandle; 
+  es.get<IdealGeometryRecord>().get("HGCalEESensitive",hgceeGeoHandle) ; 
+  es.get<IdealGeometryRecord>().get("HGCalHESiliconSensitive",hgchefGeoHandle) ; 
+  uncalibMaker_ee_.setGeometry(hgceeGeoHandle.product());
+  uncalibMaker_hef_.setGeometry(hgchefGeoHandle.product());
+  uncalibMaker_heb_.setGeometry(nullptr);
 }
 
 

RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py

@@ -1,6 +1,8 @@
 import FWCore.ParameterSet.Config as cms
 
 from SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi import *
+from RecoLocalCalo.HGCalRecProducers.HGCalUncalibRecHit_cfi import *
+
 # HGCAL rechit producer
 HGCalRecHit = cms.EDProducer(
     "HGCalRecHitProducer",
@@ -13,7 +15,9 @@
 
     # digi constants
     HGCEE_keV2DIGI  = hgceeDigitizer.digiCfg.keV2fC,
+    HGCEE_fCPerMIP  = HGCalUncalibRecHit.HGCEEConfig.fCPerMIP,
     HGCHEF_keV2DIGI = hgchefrontDigitizer.digiCfg.keV2fC,
+    HGCHEF_fCPerMIP = HGCalUncalibRecHit.HGCHEFConfig.fCPerMIP,
     HGCHEB_keV2DIGI = hgchebackDigitizer.digiCfg.keV2MIP,
 
     # algo

RecoLocalCalo/HGCalRecProducers/python/HGCalUncalibRecHit_cfi.py

@@ -21,7 +21,8 @@
         tdcNbits      = hgceeDigitizer.digiCfg.feCfg.tdcNbits,
         tdcSaturation = hgceeDigitizer.digiCfg.feCfg.tdcSaturation_fC,
         tdcOnset      = hgceeDigitizer.digiCfg.feCfg.tdcOnset_fC,
-        toaLSB_ns     = hgceeDigitizer.digiCfg.feCfg.toaLSB_ns
+        toaLSB_ns     = hgceeDigitizer.digiCfg.feCfg.toaLSB_ns,
+        fCPerMIP      = cms.vdouble(1.25,2.57,3.88) #100um, 200um, 300um
         ),
 
     HGCHEFConfig = cms.PSet(
@@ -33,7 +34,8 @@
         tdcNbits      = hgchefrontDigitizer.digiCfg.feCfg.tdcNbits,
         tdcSaturation = hgchefrontDigitizer.digiCfg.feCfg.tdcSaturation_fC,
         tdcOnset      = hgchefrontDigitizer.digiCfg.feCfg.tdcOnset_fC,
-        toaLSB_ns     = hgchefrontDigitizer.digiCfg.feCfg.toaLSB_ns
+        toaLSB_ns     = hgchefrontDigitizer.digiCfg.feCfg.toaLSB_ns,
+        fCPerMIP      = cms.vdouble(1.25,2.57,3.88) #100um, 200um, 300um
         ),
 
     HGCHEBConfig = cms.PSet(

RecoLocalCalo/HGCalRecProducers/test/testHGCalRecoLocal_cfg.py

@@ -11,8 +11,8 @@
 process.load('FWCore.MessageService.MessageLogger_cfi')
 process.load('Configuration.EventContent.EventContent_cff')
 process.load('SimGeneral.MixingModule.mixNoPU_cfi')
-process.load('Configuration.Geometry.GeometryExtended2023Reco_cff')
-process.load('Configuration.Geometry.GeometryExtended2023_cff')
+process.load('Configuration.Geometry.GeometryExtended2023DevReco_cff')
+process.load('Configuration.Geometry.GeometryExtended2023Dev_cff')
 process.load('Configuration.StandardSequences.MagneticField_38T_PostLS1_cff')
 process.load('Configuration.StandardSequences.Generator_cff')
 process.load('IOMC.EventVertexGenerators.VtxSmearedGauss_cfi')

RecoParticleFlow/PFClusterProducer/interface/PFRecHitQTests.h

@@ -638,4 +638,80 @@ class PFRecHitQTestThresholdInMIPs : public PFRecHitQTestBase {
   }
 };
 
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+class PFRecHitQTestThresholdInThicknessNormalizedMIPs : public PFRecHitQTestBase {
+ public:
+  PFRecHitQTestThresholdInThicknessNormalizedMIPs() :
+    geometryInstance_(""),
+    recHitEnergy_keV_(0.),
+    threshold_(0.),
+    mip_(0.),
+    recHitEnergyMultiplier_(0.) {
+    }
+
+  PFRecHitQTestThresholdInThicknessNormalizedMIPs(const edm::ParameterSet& iConfig) :
+    PFRecHitQTestBase(iConfig),
+    geometryInstance_(iConfig.getParameter<std::string>("geometryInstance")),
+    recHitEnergy_keV_(iConfig.getParameter<bool>("recHitEnergyIs_keV")),
+    threshold_(iConfig.getParameter<double>("thresholdInMIPs")),
+    mip_(iConfig.getParameter<double>("mipValueInkeV")),
+    recHitEnergyMultiplier_(iConfig.getParameter<double>("recHitEnergyMultiplier")) { 
+    }
+
+    void beginEvent(const edm::Event& event,const edm::EventSetup& iSetup) {
+      edm::ESHandle<HGCalGeometry> geoHandle;
+      iSetup.get<IdealGeometryRecord>().get(geometryInstance_,geoHandle);
+      ddd_ = &(geoHandle->topology().dddConstants());
+    }
+
+    bool test(reco::PFRecHit& hit,const EcalRecHit& rh,bool& clean) {
+      throw cms::Exception("WrongDetector")
+	<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
+      return false;
+    }
+    bool test(reco::PFRecHit& hit,const HBHERecHit& rh,bool& clean) {
+      throw cms::Exception("WrongDetector")
+	<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
+      return false;
+    }
+
+    bool test(reco::PFRecHit& hit,const HFRecHit& rh,bool& clean) {
+      throw cms::Exception("WrongDetector")
+	<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
+      return false;
+    }
+    bool test(reco::PFRecHit& hit,const HORecHit& rh,bool& clean) {
+      throw cms::Exception("WrongDetector")
+	<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
+      return false;
+    }
+
+    bool test(reco::PFRecHit& hit,const CaloTower& rh,bool& clean) {
+      throw cms::Exception("WrongDetector")
+	<< "PFRecHitQTestThresholdInMIPs only works for HGCAL!";
+      return false;
+    }
+
+    bool test(reco::PFRecHit& hit,const HGCRecHit& rh,bool& clean) {
+      const double newE = ( recHitEnergy_keV_ ? 
+			    1.0e-6*rh.energy()*recHitEnergyMultiplier_ :
+			    rh.energy()*recHitEnergyMultiplier_ );
+      const int wafer = HGCalDetId(rh.detid()).wafer();
+      const float mult  = (float) ddd_->waferTypeL(wafer); // 1 for 100um, 2 for 200um, 3 for 300um
+      hit.setEnergy(newE);      
+      return pass(hit,mult);
+    }
+
+ protected:
+    const std::string geometryInstance_;
+    const bool recHitEnergy_keV_;
+    const double threshold_,mip_,recHitEnergyMultiplier_;    
+    const HGCalDDDConstants*  ddd_;
+
+    bool pass(const reco::PFRecHit& hit, const float mult) {
+      const double hitValueInMIPs = 1e6*hit.energy()/(mult*mip_);
+      return hitValueInMIPs > threshold_;
+    }
+};
+
 #endif