Phase 2 - HGCal - Adding Energy above noise in RecHit

DataFormats/HGCRecHit/interface/HGCRecHit.h

@@ -74,6 +74,8 @@ class HGCRecHit : public CaloRecHit {
 
   float chi2() const;
   float outOfTimeChi2() const;
+  float signalOverSigmaNoise() const;
+
 
   // set the energy for out of time events
   // (only energy >= 0 will be stored)
@@ -83,6 +85,8 @@ class HGCRecHit : public CaloRecHit {
   void setChi2( float chi2 );
   void setOutOfTimeChi2( float chi2 );
   void setOutOfTimeEnergy( float energy );
+  void setSignalOverSigmaNoise( float sOverNoise );
+
 
   void setTimeError( uint8_t timeErrBits );
 
@@ -102,6 +106,7 @@ class HGCRecHit : public CaloRecHit {
 
   /// store rechit condition (see Flags enum) in a bit-wise way 
   uint32_t flagBits_;
+  uint8_t signalOverSigmaNoise_;
 };
 
 std::ostream& operator<<(std::ostream& s, const HGCRecHit& hit);

DataFormats/HGCRecHit/src/HGCRecHit.cc

@@ -14,12 +14,12 @@ HGCRecHit::HGCRecHit(const DetId& id, float energy, float time, uint32_t flags,
 
 float HGCRecHit::chi2() const {
   uint32_t rawChi2 = 0x7F & (flags()>>4);
-  return (float)rawChi2 / (float)((1<<7)-1) * 64.;
+  return (float)rawChi2 / (float)((1<<7)-1) * 64.f;
 }
 
 float HGCRecHit::outOfTimeChi2() const {
   uint32_t rawChi2Prob = 0x7F & (flags()>>24);
-  return (float)rawChi2Prob / (float)((1<<7)-1) * 64.;
+  return (float)rawChi2Prob / (float)((1<<7)-1) * 64.f;
 }
 
 float HGCRecHit::outOfTimeEnergy() const {
@@ -33,13 +33,13 @@ void HGCRecHit::setChi2( float chi2 ) {
   // bound the max value of the chi2
   if ( chi2 > 64 ) chi2 = 64;
   // use 7 bits
-  uint32_t rawChi2 = lround( chi2 / 64. * ((1<<7)-1) );
+  uint32_t rawChi2 = lround( chi2 / 64.f * ((1<<7)-1) );
   // shift by 4 bits (recoFlag)
   setFlags( (~(0x7F<<4) & flags()) | ((rawChi2 & 0x7F)<<4) );
 }
 
 void HGCRecHit::setOutOfTimeEnergy( float energy ) {
-  if ( energy > 0.001 ) {
+  if ( energy > 0.001f ) {
     uint16_t exponent = lround(floor(log10(energy)))+3;
     uint16_t significand = lround(energy/pow(10,exponent-5));
     // use 13 bits (3 exponent, 10 significand)
@@ -53,11 +53,22 @@ void HGCRecHit::setOutOfTimeChi2( float chi2 ) {
   // bound the max value of chi2
   if ( chi2 > 64 ) chi2 = 64;
   // use 7 bits
-  uint32_t rawChi2 = lround( chi2 / 64. * ((1<<7)-1) );
+  uint32_t rawChi2 = lround( chi2 / 64.f * ((1<<7)-1) );
   // shift by 24 bits (recoFlag + chi2 + outOfTimeEnergy)
   setFlags( (~(0x7F<<24) & flags()) | ((rawChi2 & 0x7F)<<24) );
 }
 
+void HGCRecHit::setSignalOverSigmaNoise( float sOverNoise ) {
+    // bound the max value of sOverNoise
+    if(sOverNoise > 32.f) sOverNoise = 32.f;
+    //use 8 bits
+    signalOverSigmaNoise_ = lround( sOverNoise / 32.f * ((1<<8)-1) );
+}
+
+float HGCRecHit::signalOverSigmaNoise() const {
+    return (float)signalOverSigmaNoise_ * 0.125f;
+}
+
 void HGCRecHit::setTimeError( uint8_t timeErrBits ) {
   // take the bits and put them in the right spot
   setAux( (~0xFF & aux()) | timeErrBits );
@@ -75,7 +86,7 @@ float HGCRecHit::timeError() const {
   float LSB = 1.26008;
   uint8_t exponent = timeErrorBits>>5;
   uint8_t significand = timeErrorBits & ~(0x7<<5);
-  return pow(2.,exponent)*significand*LSB/1000.;
+  return pow(2.,exponent)*significand*LSB/1000.f;
 }
 
 bool HGCRecHit::isTimeValid() const {

DataFormats/HGCRecHit/src/classes_def.xml

@@ -9,7 +9,8 @@
   <class name="edm::RefVector<edm::SortedCollection<HGCUncalibratedRecHit,edm::StrictWeakOrdering<HGCUncalibratedRecHit> >,HGCUncalibratedRecHit,edm::refhelper::FindUsingAdvance<edm::SortedCollection<HGCUncalibratedRecHit,edm::StrictWeakOrdering<HGCUncalibratedRecHit> >,HGCUncalibratedRecHit> >"/>
   <class name="edm::RefProd<edm::SortedCollection<HGCUncalibratedRecHit,edm::StrictWeakOrdering<HGCUncalibratedRecHit> > >"/>
 
-  <class name="HGCRecHit" ClassVersion="12">
+  <class name="HGCRecHit" ClassVersion="13">
+   <version ClassVersion="13" checksum="311955435"/>
    <version ClassVersion="12" checksum="3568613480"/>
    <version ClassVersion="11" checksum="3873633606"/>
   </class>

RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.cc

@@ -9,111 +9,133 @@
 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
 
 HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet&ps) :
-  HGCalRecHitWorkerBaseClass(ps) {
-  rechitMaker_.reset( new HGCalRecHitSimpleAlgo() );
-  constexpr float keV2GeV = 1e-6;
-  // HGCee constants 
-  HGCEE_keV2DIGI_   =  ps.getParameter<double>("HGCEE_keV2DIGI");
-  HGCEE_fCPerMIP_   =  ps.getParameter<std::vector<double> >("HGCEE_fCPerMIP");
-  HGCEE_isSiFE_     =  ps.getParameter<bool>("HGCEE_isSiFE");
-  hgceeUncalib2GeV_ = keV2GeV/HGCEE_keV2DIGI_;
-
-  // HGChef constants
-  HGCHEF_keV2DIGI_   =  ps.getParameter<double>("HGCHEF_keV2DIGI");
-  HGCHEF_fCPerMIP_   =  ps.getParameter<std::vector<double> >("HGCHEF_fCPerMIP");
-  HGCHEF_isSiFE_     =  ps.getParameter<bool>("HGCHEF_isSiFE");
-  hgchefUncalib2GeV_ = keV2GeV/HGCHEF_keV2DIGI_;
-
-  // HGCheb constants
-  HGCHEB_keV2DIGI_   =  ps.getParameter<double>("HGCHEB_keV2DIGI");
-  HGCHEB_isSiFE_     =  ps.getParameter<bool>("HGCHEB_isSiFE");
-  hgchebUncalib2GeV_ = keV2GeV/HGCHEB_keV2DIGI_;
-
-  // layer weights (from Valeri/Arabella)
-  std::vector<float> weights;
-  const auto& dweights = ps.getParameter<std::vector<double> >("layerWeights");
-  for( auto weight : dweights ) {
-    weights.push_back(weight);
-  }
-  rechitMaker_->setLayerWeights(weights);
-
-  // residual correction for cell thickness
-  const auto& rcorr = ps.getParameter<std::vector<double> >("thicknessCorrection");
-  rcorr_.clear();
-  rcorr_.push_back(1.f);
-  for( auto corr : rcorr ) {
-    rcorr_.push_back(1.0/corr);
-  }
-
+        HGCalRecHitWorkerBaseClass(ps)
+{
+    rechitMaker_.reset(new HGCalRecHitSimpleAlgo());
+    tools_.reset(new hgcal::RecHitTools());
+    constexpr float keV2GeV = 1e-6;
+    // HGCee constants
+    HGCEE_keV2DIGI_ = ps.getParameter<double>("HGCEE_keV2DIGI");
+    HGCEE_fCPerMIP_ = ps.getParameter < std::vector<double> > ("HGCEE_fCPerMIP");
+    HGCEE_isSiFE_ = ps.getParameter<bool>("HGCEE_isSiFE");
+    hgceeUncalib2GeV_ = keV2GeV / HGCEE_keV2DIGI_;
+
+    // HGChef constants
+    HGCHEF_keV2DIGI_ = ps.getParameter<double>("HGCHEF_keV2DIGI");
+    HGCHEF_fCPerMIP_ = ps.getParameter < std::vector<double> > ("HGCHEF_fCPerMIP");
+    HGCHEF_isSiFE_ = ps.getParameter<bool>("HGCHEF_isSiFE");
+    hgchefUncalib2GeV_ = keV2GeV / HGCHEF_keV2DIGI_;
+
+    // HGCheb constants
+    HGCHEB_keV2DIGI_ = ps.getParameter<double>("HGCHEB_keV2DIGI");
+    HGCHEB_isSiFE_ = ps.getParameter<bool>("HGCHEB_isSiFE");
+    hgchebUncalib2GeV_ = keV2GeV / HGCHEB_keV2DIGI_;
+
+    // layer weights (from Valeri/Arabella)
+    const auto& dweights = ps.getParameter < std::vector<double> > ("layerWeights");
+    for (auto weight : dweights)
+    {
+        weights_.push_back(weight);
+    }
+
+    rechitMaker_->setLayerWeights(weights_);
+
+    // residual correction for cell thickness
+    const auto& rcorr = ps.getParameter < std::vector<double> > ("thicknessCorrection");
+    rcorr_.clear();
+    rcorr_.push_back(1.f);
+    for (auto corr : rcorr)
+    {
+        rcorr_.push_back(1.0 / corr);
+    }
+
+
+    HGCEE_noise_fC_ = ps.getParameter < std::vector<double> > ("HGCEE_noise_fC");
+    HGCHEF_noise_fC_ = ps.getParameter < std::vector<double> > ("HGCHEF_noise_fC");
+    HGCHEB_noise_MIP_ = ps.getParameter<double>("HGCHEB_noise_MIP");
+
+
 }
 
-void HGCalRecHitWorkerSimple::set(const edm::EventSetup& es) {
-  if (HGCEE_isSiFE_) {
-    edm::ESHandle<HGCalGeometry> hgceeGeoHandle; 
-    es.get<IdealGeometryRecord>().get("HGCalEESensitive",hgceeGeoHandle); 
-    ddds_[0] = &(hgceeGeoHandle->topology().dddConstants());
-  } else {
-    ddds_[0] = nullptr;
-  }
-  if (HGCHEF_isSiFE_) {
-    edm::ESHandle<HGCalGeometry> hgchefGeoHandle; 
-    es.get<IdealGeometryRecord>().get("HGCalHESiliconSensitive",hgchefGeoHandle); 
-    ddds_[1] = &(hgchefGeoHandle->topology().dddConstants());
-  } else {
-    ddds_[1] = nullptr;
-  }
-  ddds_[2] = nullptr;  
+void HGCalRecHitWorkerSimple::set(const edm::EventSetup& es)
+{
+    tools_->getEventSetup(es);
+    if (HGCEE_isSiFE_)
+    {
+        edm::ESHandle < HGCalGeometry > hgceeGeoHandle;
+        es.get<IdealGeometryRecord>().get("HGCalEESensitive", hgceeGeoHandle);
+        ddds_[0] = &(hgceeGeoHandle->topology().dddConstants());
+    }
+    else
+    {
+        ddds_[0] = nullptr;
+    }
+    if (HGCHEF_isSiFE_)
+    {
+        edm::ESHandle < HGCalGeometry > hgchefGeoHandle;
+        es.get<IdealGeometryRecord>().get("HGCalHESiliconSensitive", hgchefGeoHandle);
+        ddds_[1] = &(hgchefGeoHandle->topology().dddConstants());
+    }
+    else
+    {
+        ddds_[1] = nullptr;
+    }
+    ddds_[2] = nullptr;
 }
 
+bool HGCalRecHitWorkerSimple::run(const edm::Event & evt, const HGCUncalibratedRecHit& uncalibRH,
+        HGCRecHitCollection & result)
+{
+    DetId detid = uncalibRH.id();
+    int thickness = -1;
+    float sigmaNoiseGeV = 0.f;
+    unsigned int layer = tools_->getLayerWithOffset(detid);
+    HGCalDetId hid;
 
-bool
-HGCalRecHitWorkerSimple::run( const edm::Event & evt,
-                              const HGCUncalibratedRecHit& uncalibRH,
-                              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 HcalEndcap:
-  case HGCHEB:
-    rechitMaker_->setADCToGeVConstant(float(hgchebUncalib2GeV_) );
-    break;  
-  default:
-    throw cms::Exception("NonHGCRecHit")
-      << "Rechit with detid = " << detid.rawId() << " is not HGC!";
-  }
-
-  // make the rechit and put in the output collection
-  if (recoFlag == 0) {    
-    HGCRecHit myrechit( rechitMaker_->makeRecHit(uncalibRH, 0) );    
-    int thk = -1;
-    if( detid.subdetId() != HcalEndcap ) {
-      HGCalDetId hid(detid);
-      thk = ddds_[hid.subdetId()-3]->waferTypeL(hid.wafer());
-      // units out of rechit maker are MIP * (GeV/fC)
-      // so multiple
+    switch (detid.subdetId())
+    {
+    case HGCEE:
+        rechitMaker_->setADCToGeVConstant(float(hgceeUncalib2GeV_));
+        hid = detid;
+        thickness = ddds_[hid.subdetId() - 3]->waferTypeL(hid.wafer());
+        sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness]
+                    * HGCHEF_noise_fC_[thickness - 1] / HGCHEF_fCPerMIP_[thickness - 1];
+        break;
+    case HGCHEF:
+        rechitMaker_->setADCToGeVConstant(float(hgchefUncalib2GeV_));
+        hid = detid;
+        thickness = ddds_[hid.subdetId() - 3]->waferTypeL(hid.wafer());
+        sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness]
+                    * HGCHEF_noise_fC_[thickness - 1] / HGCHEF_fCPerMIP_[thickness - 1];
+        break;
+    case HcalEndcap:
+    case HGCHEB:
+        rechitMaker_->setADCToGeVConstant(float(hgchebUncalib2GeV_));
+        hid = detid;
+        sigmaNoiseGeV = 1e-3 * HGCHEB_noise_MIP_ * weights_[layer];
+        break;
+    default:
+        throw cms::Exception("NonHGCRecHit") << "Rechit with detid = " << detid.rawId()
+                << " is not HGC!";
     }
-    const double new_E = myrechit.energy()*(thk == -1 ? 1.0 : rcorr_[thk]);
+
+    // make the rechit and put in the output collection
+
+    HGCRecHit myrechit(rechitMaker_->makeRecHit(uncalibRH, 0));
+    const double new_E = myrechit.energy() * (thickness == -1 ? 1.0 : rcorr_[thickness]);
+
+
     myrechit.setEnergy(new_E);
+    myrechit.setSignalOverSigmaNoise(new_E/sigmaNoiseGeV);
     result.push_back(myrechit);
-  }
 
-  return true;
+    return true;
 }
 
-HGCalRecHitWorkerSimple::~HGCalRecHitWorkerSimple(){
+HGCalRecHitWorkerSimple::~HGCalRecHitWorkerSimple()
+{
 }
 
-
 #include "FWCore/Framework/interface/MakerMacros.h"
 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalRecHitWorkerFactory.h"
 DEFINE_EDM_PLUGIN( HGCalRecHitWorkerFactory, HGCalRecHitWorkerSimple, "HGCalRecHitWorkerSimple" );

RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.h

@@ -14,6 +14,7 @@
 #include "FWCore/Framework/interface/ESHandle.h"
 
 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
 
 class HGCalRecHitWorkerSimple : public HGCalRecHitWorkerBaseClass {
  public:
@@ -32,16 +33,25 @@ class HGCalRecHitWorkerSimple : public HGCalRecHitWorkerBaseClass {
   double HGCHEB_keV2DIGI_, hgchebUncalib2GeV_;
   bool HGCEE_isSiFE_, HGCHEF_isSiFE_, HGCHEB_isSiFE_;
 
+
+
+  std::vector<double> HGCEE_noise_fC_;
+  std::vector<double> HGCHEF_noise_fC_;
+  double HGCHEB_noise_MIP_;
+
+
   std::array<const HGCalDDDConstants*, 3> ddds_;
 
   std::vector<int> v_chstatus_;
 
   std::vector<int> v_DB_reco_flags_;
   bool killDeadChannels_;
 
-  std::vector<float> rcorr_;
-
+  std::vector<double> rcorr_;
+  std::vector<float> weights_;
   std::unique_ptr<HGCalRecHitSimpleAlgo> rechitMaker_;
+  std::unique_ptr<hgcal::RecHitTools> tools_;
+
 };
 
 #endif

RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py

@@ -1,5 +1,4 @@
 import FWCore.ParameterSet.Config as cms
-
 from SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi import *
 from RecoLocalCalo.HGCalRecProducers.HGCalUncalibRecHit_cfi import *
 
@@ -79,9 +78,14 @@
 
     # EM Scale calibrations
     layerWeights = dEdX_weights,
-    thicknessCorrection = cms.vdouble(1.132,1.092,1.084), # 100, 200, 300 um 
 
+    thicknessCorrection = cms.vdouble(1.132,1.092,1.084), # 100, 200, 300 um 
+    HGCEE_noise_fC = hgceeDigitizer.digiCfg.noise_fC,
+    HGCHEF_noise_fC = hgchefrontDigitizer.digiCfg.noise_fC,
+    HGCHEB_noise_MIP = hgchebackDigitizer.digiCfg.noise_MIP,
     # algo
     algo = cms.string("HGCalRecHitWorkerSimple")
 
     )
+
+