Add in CCE degradation and add HGCal to aging customizations

DataFormats/HGCDigi/interface/HGCSample.h

@@ -33,6 +33,9 @@ class HGCSample {
   void setData(uint16_t data)           { setWord(data, kDataMask,   kDataShift);   }
   void set(bool thr, bool mode,uint16_t toa, uint16_t data) 
   { 
+    toa = ( toa > kToAMask ? kToAMask : toa );
+    data = ( data > kDataMask ? kDataMask : data);
+
     value_ = ( ( (uint32_t)thr  & kThreshMask ) << kThreshShift | 
                ( (uint32_t)mode & kModeMask   ) << kModeShift   |
                ( (uint32_t)toa  & kToAMask    ) << kToAShift    | 
@@ -64,6 +67,7 @@ class HGCSample {
    */
   void setWord(uint32_t word, uint32_t mask, uint32_t pos)
   {
+    if( word > mask ) word = mask; // deal with saturation
     //clear required bits
     const uint32_t masked_word = (word & mask) << pos;
     value_ &= ~(masked_word); 

RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.cc

@@ -51,7 +51,9 @@ HGCalRecHitWorkerSimple::HGCalRecHitWorkerSimple(const edm::ParameterSet&ps) :
 
 
     hgcEE_noise_fC_ = ps.getParameter < std::vector<double> > ("HGCEE_noise_fC");
+    hgcEE_cce_ = ps.getParameter< std::vector<double> > ("HGCEE_cce");
     hgcHEF_noise_fC_ = ps.getParameter < std::vector<double> > ("HGCHEF_noise_fC");
+    hgcHEF_cce_ = ps.getParameter< std::vector<double> > ("HGCHEF_cce");
     hgcHEB_noise_MIP_ = ps.getParameter<double>("HGCHEB_noise_MIP");
 
     // don't produce rechit if detid is a ghost one
@@ -98,20 +100,23 @@ bool HGCalRecHitWorkerSimple::run(const edm::Event & evt, const HGCUncalibratedR
     float sigmaNoiseGeV = 0.f;
     unsigned int layer = tools_->getLayerWithOffset(detid);
     HGCalDetId hid;
+    float cce_correction = 1.0;
 
     switch (detid.subdetId())
     {
     case HGCEE:
         rechitMaker_->setADCToGeVConstant(float(hgceeUncalib2GeV_));
         hid = detid;
         thickness = ddds_[hid.subdetId() - 3]->waferTypeL(hid.wafer());
+        cce_correction = hgcEE_cce_[thickness - 1];
         sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness]
                     * hgcEE_noise_fC_[thickness - 1] / hgcEE_fCPerMIP_[thickness - 1];
         break;
     case HGCHEF:
         rechitMaker_->setADCToGeVConstant(float(hgchefUncalib2GeV_));
         hid = detid;
         thickness = ddds_[hid.subdetId() - 3]->waferTypeL(hid.wafer());
+        cce_correction = hgcHEF_cce_[thickness - 1];
         sigmaNoiseGeV = 1e-3 * weights_[layer] * rcorr_[thickness]
                     * hgcHEF_noise_fC_[thickness - 1] / hgcHEF_fCPerMIP_[thickness - 1];
         break;
@@ -129,10 +134,10 @@ bool HGCalRecHitWorkerSimple::run(const edm::Event & evt, const HGCUncalibratedR
     // 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]);
+    const double new_E = myrechit.energy() * (thickness == -1 ? 1.0 : rcorr_[thickness])/cce_correction;
 
 
-    myrechit.setEnergy(new_E);
+    myrechit.setEnergy(new_E); 
     myrechit.setSignalOverSigmaNoise(new_E/sigmaNoiseGeV);
     result.push_back(myrechit);
 

RecoLocalCalo/HGCalRecProducers/plugins/HGCalRecHitWorkerSimple.h

@@ -19,17 +19,19 @@
 class HGCalRecHitWorkerSimple : public HGCalRecHitWorkerBaseClass {
  public:
   HGCalRecHitWorkerSimple(const edm::ParameterSet&);
-  virtual ~HGCalRecHitWorkerSimple();                       
+  ~HGCalRecHitWorkerSimple() override;                       
 
-  void set(const edm::EventSetup& es);
-  bool run(const edm::Event& evt, const HGCUncalibratedRecHit& uncalibRH, HGCRecHitCollection & result);
+  void set(const edm::EventSetup& es) override;
+  bool run(const edm::Event& evt, const HGCUncalibratedRecHit& uncalibRH, HGCRecHitCollection & result) override;
 
  protected:
 
   double hgcEE_keV2DIGI_,  hgceeUncalib2GeV_;
   std::vector<double> hgcEE_fCPerMIP_;
+  std::vector<double> hgcEE_cce_;
   double hgcHEF_keV2DIGI_, hgchefUncalib2GeV_;
   std::vector<double> hgcHEF_fCPerMIP_;
+  std::vector<double> hgcHEF_cce_;
   double hgcHEB_keV2DIGI_, hgchebUncalib2GeV_;
   bool hgcEE_isSiFE_, hgcHEF_isSiFE_, hgcHEB_isSiFE_;
 

RecoLocalCalo/HGCalRecProducers/python/HGCalRecHit_cfi.py

@@ -85,7 +85,9 @@
 
     thicknessCorrection = cms.vdouble(1.132,1.092,1.084), # 100, 200, 300 um 
     HGCEE_noise_fC = hgceeDigitizer.digiCfg.noise_fC,
+    HGCEE_cce = hgceeDigitizer.digiCfg.chargeCollectionEfficiency,
     HGCHEF_noise_fC = hgchefrontDigitizer.digiCfg.noise_fC,
+    HGCHEF_cce = hgchefrontDigitizer.digiCfg.chargeCollectionEfficiency,
     HGCHEB_noise_MIP = hgchebackDigitizer.digiCfg.noise_MIP,
     # algo
     algo = cms.string("HGCalRecHitWorkerSimple")

SLHCUpgradeSimulations/Configuration/python/aging.py

@@ -8,6 +8,16 @@ def getHcalDigitizer(process):
         return process.mix.digitizers.hcal
     return None
 
+def getHGCalDigitizer(process,section):
+    if hasattr(process,'mix') and hasattr(process.mix,'digitizers'):
+        if section == 'EE' and hasattr(process.mix.digitizers,'hgceeDigitizer'):
+            return process.mix.digitizers.hgceeDigitizer
+        elif section == 'FH' and hasattr(process.mix.digitizers,'hgchefrontDigitizer'):
+            return process.mix.digitizers.hgchefrontDigitizer
+        elif section == 'BH' and hasattr(process.mix.digitizers,'hgchebackDigitizer'):
+            return process.mix.digitizers.hgchebackDigitizer
+    return None
+
 # change assumptions about lumi rate
 def setScenarioHLLHC(module,scenarioHLLHC):
     if scenarioHLLHC=="nominal":
@@ -49,6 +59,13 @@ def ageHF(process,turnon):
         process.es_hardcode.HFRecalibration = cms.bool(turnon)
     return process
 
+def agedHGCal(process):
+    from SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi import HGCal_setEndOfLifeNoise
+    HGCal_setEndOfLifeNoise(getHGCalDigitizer(process,'EE'))
+    HGCal_setEndOfLifeNoise(getHGCalDigitizer(process,'FH'))
+    HGCal_setEndOfLifeNoise(getHGCalDigitizer(process,'BH'))
+    return process
+
 # needs lumi to set proper ZS thresholds (tbd)
 def ageSiPM(process,turnon,lumi):
     process.es_hardcode.hbUpgrade.doRadiationDamage = turnon
@@ -182,14 +199,17 @@ def customise_aging_1000(process):
 def customise_aging_3000(process):
     process=ageHcal(process,3000,5.0e34,"nominal")
     process=ageEcal(process,3000,5.0e34)
+    process=agedHGCal(process)
     return process
 
 def customise_aging_3000_ultimate(process):
     process=ageHcal(process,3000,7.5e34,"ultimate")
     process=ageEcal(process,3000,7.5e34)
+    process=agedHGCal(process)
     return process
 
 def customise_aging_4500_ultimate(process):
     process=ageHcal(process,4500,7.5e34,"ultimate")
     process=ageEcal(process,4500,7.5e34)
+    process=agedHGCal(process)
     return process

SimCalorimetry/HGCalSimProducers/interface/HGCDigitizer.h

@@ -114,6 +114,8 @@ private :
   //average occupancies
   std::array<double,3> averageOccupancies_;
   uint32_t nEvents_;
+
+  std::vector<float> cce_;
 };
 
 

SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h

@@ -87,6 +87,9 @@ class HGCDigitizerBase {
 
   //noise level
   std::vector<float> noise_fC_;
+
+  //charge collection efficiency
+  std::vector<double> cce_;
 
   //front-end electronics model
   std::unique_ptr<HGCFEElectronics<DFr> > myFEelectronics_;

SimCalorimetry/HGCalSimProducers/interface/HGCEEDigitizer.h

@@ -11,7 +11,7 @@ class HGCEEDigitizer : public HGCDigitizerBase<HGCEEDataFrame> {
   void runDigitizer(std::unique_ptr<HGCEEDigiCollection> &digiColl,hgc::HGCSimHitDataAccumulator &simData,
 		    const CaloSubdetectorGeometry* theGeom, const std::unordered_set<DetId>& validIds,
 		    uint32_t digitizationType, CLHEP::HepRandomEngine* engine) override; 
-  ~HGCEEDigitizer();
+  ~HGCEEDigitizer() override;
 private:
 
 };

SimCalorimetry/HGCalSimProducers/interface/HGCFEElectronics.h

@@ -33,13 +33,13 @@ class HGCFEElectronics
      @short switches according to the firmware version
    */
   inline void runShaper(DFr &dataFrame, hgc::HGCSimHitData& chargeColl, 
-                        hgc::HGCSimHitData& toa, int thickness, CLHEP::HepRandomEngine* engine)
+                        hgc::HGCSimHitData& toa, int thickness, CLHEP::HepRandomEngine* engine, float cce = 1.0)
   {    
     switch(fwVersion_)
       {
-      case SIMPLE :  { runSimpleShaper(dataFrame,chargeColl, thickness);      break; }
-      case WITHTOT : { runShaperWithToT(dataFrame,chargeColl,toa, thickness, engine); break; }
-      default :      { runTrivialShaper(dataFrame,chargeColl, thickness);     break; }
+      case SIMPLE :  { runSimpleShaper(dataFrame,chargeColl, thickness, cce);      break; }
+      case WITHTOT : { runShaperWithToT(dataFrame,chargeColl,toa, thickness, engine, cce); break; }
+      default :      { runTrivialShaper(dataFrame,chargeColl, thickness, cce);     break; }
       }
   }
 
@@ -55,18 +55,18 @@ class HGCFEElectronics
   /**
      @short converts charge to digis without pulse shape
    */
-  void runTrivialShaper(DFr &dataFrame, hgc::HGCSimHitData& chargeColl, int thickness);
+  void runTrivialShaper(DFr &dataFrame, hgc::HGCSimHitData& chargeColl, int thickness, float cce = 1.0);
 
   /**
      @short applies a shape to each time sample and propagates the tails to the subsequent time samples
    */
-  void runSimpleShaper(DFr &dataFrame, hgc::HGCSimHitData& chargeColl, int thickness);
+  void runSimpleShaper(DFr &dataFrame, hgc::HGCSimHitData& chargeColl, int thickness, float cce = 1.0);
 
   /**
      @short implements pulse shape and switch to time over threshold including deadtime
    */
   void runShaperWithToT(DFr &dataFrame, hgc::HGCSimHitData& chargeColl, 
-                        hgc::HGCSimHitData& toa, int thickness, CLHEP::HepRandomEngine* engine);
+                        hgc::HGCSimHitData& toa, int thickness, CLHEP::HepRandomEngine* engine, float cce = 1.0);
 
   /**
      @short returns how ToT will be computed
@@ -87,6 +87,7 @@ class HGCFEElectronics
   float adcSaturation_fC_, adcLSB_fC_, tdcLSB_fC_, tdcSaturation_fC_,
     adcThreshold_fC_, tdcOnset_fC_, toaLSB_ns_, tdcResolutionInNs_; 
   uint32_t toaMode_;
+  bool thresholdFollowsMIP_;
   //caches
   std::array<bool,hgc::nSamples>  busyFlags, totFlags;
   hgc::HGCSimHitData newCharge, toaFromToT;

SimCalorimetry/HGCalSimProducers/interface/HGCHEbackDigitizer.h

@@ -12,7 +12,7 @@ class HGCHEbackDigitizer : public HGCDigitizerBase<HGCBHDataFrame>
   void runDigitizer(std::unique_ptr<HGCBHDigiCollection> &digiColl,hgc::HGCSimHitDataAccumulator &simData,
 		    const CaloSubdetectorGeometry* theGeom, const std::unordered_set<DetId>& validIds,
 		    uint32_t digitizationType, CLHEP::HepRandomEngine* engine) override;
-  ~HGCHEbackDigitizer();
+  ~HGCHEbackDigitizer() override;
 
  private:
 

SimCalorimetry/HGCalSimProducers/interface/HGCHEfrontDigitizer.h

@@ -11,7 +11,7 @@ class HGCHEfrontDigitizer : public HGCDigitizerBase<HGCHEDataFrame> {
   void runDigitizer(std::unique_ptr<HGCHEDigiCollection> &digiColl, hgc::HGCSimHitDataAccumulator &simData,
 		    const CaloSubdetectorGeometry* theGeom, const std::unordered_set<DetId>& validIds,
 		    uint32_t digitizationType, CLHEP::HepRandomEngine* engine) override;
-  ~HGCHEfrontDigitizer();
+  ~HGCHEfrontDigitizer() override;
 private:
 
 };

SimCalorimetry/HGCalSimProducers/python/hgcalDigitizer_cfi.py

@@ -3,7 +3,9 @@
 # Base configurations for HGCal digitizers
 eV_per_eh_pair = 3.62
 fC_per_ele     = 1.6020506e-4
-nonAgedNoises = [2100.0,2100.0,1600.0] #100,200,300 um (in electrons)
+nonAgedCCEs    = [1.0, 1.0, 1.0]
+nonAgedNoises  = [2100.0,2100.0,1600.0] #100,200,300 um (in electrons)
+thresholdTracksMIP = False
 
 # ECAL
 hgceeDigitizer = cms.PSet( 
@@ -20,6 +22,8 @@
     verbosity         = cms.untracked.uint32(0),
     digiCfg = cms.PSet( 
         keV2fC           = cms.double(0.044259), #1000 eV/3.62 (eV per e) / 6.24150934e3 (e per fC)
+
+        chargeCollectionEfficiency = cms.vdouble( nonAgedCCEs ),
         noise_fC         = cms.vdouble( [x*fC_per_ele for x in nonAgedNoises] ), #100,200,300 um
         doTimeSamples    = cms.bool(False),                                         
         feCfg   = cms.PSet( 
@@ -42,6 +46,7 @@
             # raise threshold flag (~MIP/2) this is scaled 
             # for different thickness
             adcThreshold_fC   = cms.double(0.672),
+            thresholdFollowsMIP        = cms.bool(thresholdTracksMIP),
             # raise usage of TDC and mode flag (from J. Kaplon)
             tdcOnset_fC       = cms.double(60) ,
             # LSB for time of arrival estimate from TDC in ns
@@ -72,6 +77,7 @@
     verbosity         = cms.untracked.uint32(0),
     digiCfg = cms.PSet(        
         keV2fC           = cms.double(0.044259), #1000 eV / 3.62 (eV per e) / 6.24150934e3 (e per fC)
+        chargeCollectionEfficiency = cms.vdouble( nonAgedCCEs ),
         noise_fC         = cms.vdouble( [x*fC_per_ele for x in nonAgedNoises] ), #100,200,300 um
         doTimeSamples    = cms.bool(False),                                         
         feCfg   = cms.PSet( 
@@ -93,6 +99,7 @@
             # raise threshold flag (~MIP/2) this is scaled 
             # for different thickness
             adcThreshold_fC   = cms.double(0.672),
+            thresholdFollowsMIP        = cms.bool(thresholdTracksMIP),
             # raise usage of TDC and mode flag (from J. Kaplon)
             tdcOnset_fC       = cms.double(60) ,
             # LSB for time of arrival estimate from TDC in ns
@@ -138,15 +145,18 @@
             # ADC saturation : in this case we use the same variable but fC=MIP
             adcSaturation_fC = cms.double(1024.0),
             # 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.50),
+            thresholdFollowsMIP = cms.bool(False)
             )
         )                              
     )
 
 #function to set noise to aged HGCal
+endOfLifeCCEs = [0.5, 0.5, 0.7]
 endOfLifeNoises = [2400.0,2250.0,1750.0]
 def HGCal_setEndOfLifeNoise(digitizer):
     if( digitizer.digiCollection != "HGCDigisHEback" ):
         digitizer.digiCfg.noise_fC = cms.vdouble( [x*fC_per_ele for x in endOfLifeNoises] )
+        digitizer.digiCfg.chargeCollectionEfficiency = cms.double(endOfLifeCCEs)
     else: #use S/N of 7 for SiPM readout
         digitizer.digiCfg.noise_MIP = cms.vdouble( 1.0/5.0 )

SimCalorimetry/HGCalSimProducers/src/HGCDigitizer.cc

@@ -20,6 +20,7 @@
 
 #include <algorithm>
 #include <boost/foreach.hpp>
+#include "FWCore/Utilities/interface/transform.h"
 
 using namespace hgc_digi;
 
@@ -69,7 +70,7 @@ namespace {
     const auto& topo     = geom->topology();
     const auto& dddConst = topo.dddConstants();
 
-    int subdet, layer, cell, sec, subsec, zp;
+    int subdet(DetId(simId).subdetId()), layer, cell, sec, subsec, zp;
 
     const bool isSqr = (dddConst.geomMode() == HGCalGeometryMode::Square);
     if (isSqr) {
@@ -91,6 +92,26 @@ namespace {
 
     return result;
   }  
+
+  float getCCE(const HGCalGeometry* geom,
+	       const DetId& detid,
+	       const std::vector<float>&cces) {
+    if( cces.empty() ) return 1.f;
+    const auto& topo     = geom->topology();
+    const auto& dddConst = topo.dddConstants();
+    uint32_t id(detid.rawId());
+    HGCalDetId hid(id);
+    int wafer = HGCalDetId(id).wafer();
+    int waferTypeL = dddConst.waferTypeL(wafer);  
+    return cces[waferTypeL-1];
+  }
+
+  float getCCE(const HcalGeometry* geom,
+	       const DetId& id,
+	       const std::vector<float>&cces) {
+    return 1.f;
+  }
+
 }
 
 //
@@ -110,10 +131,21 @@ HGCDigitizer::HGCDigitizer(const edm::ParameterSet& ps,
   digitizationType_  = ps.getParameter< uint32_t >("digitizationType");
   verbosity_         = ps.getUntrackedParameter< uint32_t >("verbosity",0);
   tofDelay_          = ps.getParameter< double >("tofDelay");  
-  
+
   std::unordered_set<DetId>().swap(validIds_);
 
   iC.consumes<std::vector<PCaloHit> >(edm::InputTag("g4SimHits",hitCollection_));
+  const auto& myCfg_ = ps.getParameter<edm::ParameterSet>("digiCfg");
+
+  if( myCfg_.existsAs<std::vector<double> >( "chargeCollectionEfficiencies" ) ) {
+    cce_.clear();
+    const auto& temp = myCfg_.getParameter<std::vector<double> >("chargeCollectionEfficiencies");
+    for( double cce : temp ) {
+      cce_.push_back(cce);
+    }
+  } else {
+    std::vector<float>().swap(cce_);
+  }
 
   if(hitCollection_.find("HitsEE")!=std::string::npos) { 
     mySubDet_=ForwardSubdetector::HGCEE;  
@@ -296,9 +328,9 @@ void HGCDigitizer::accumulate(edm::Handle<edm::PCaloHitContainer> const &hits,
 
     if (verbosity_>0) {
       if (producesEEDigis())
-	edm::LogInfo("HGCDigitizer") << " i/p " << std::hex << the_hit.id() << std::dec << " o/p " << id.rawId() << std::endl;
+	edm::LogInfo("HGCDigitizer") << " i/p " << std::hex << the_hit.id() << " o/p " << id.rawId() << std::dec << std::endl;
       else
-	edm::LogInfo("HGCDigitizer") << " i/p " << std::hex << the_hit.id() << std::dec << " o/p " << id.rawId() << std::endl;
+	edm::LogInfo("HGCDigitizer") << " i/p " << std::hex << the_hit.id() << " o/p " << id.rawId() << std::dec << std::endl;
     }
 
     if( 0 != id.rawId() ) {      
@@ -322,8 +354,8 @@ void HGCDigitizer::accumulate(edm::Handle<edm::PCaloHitContainer> const &hits,
 
     const float toa    = std::get<2>(hitRefs[i]);
     const PCaloHit &hit=hits->at( hitidx );     
-    const float charge = hit.energy()*1e6*keV2fC;
-      
+    const float charge = hit.energy()*1e6*keV2fC*getCCE(geom,id,cce_);
+
     //distance to the center of the detector
     const float dist2center( getPositionDistance(geom,id) );