Add in the digitization code for HGCal

SimCalorimetry/HGCalSimProducers/BuildFile.xml

@@ -0,0 +1,12 @@
+<use   name="FWCore/Framework"/>
+<use   name="FWCore/ParameterSet"/>
+<use   name="FWCore/Utilities"/>
+<use   name="SimGeneral/MixingModule"/>
+<use   name="DataFormats/ForwardDetId"/>
+<use   name="Geometry/Records"/>
+<use   name="Geometry/HGCalGeometry"/>
+<use   name="hepmc"/>
+<use   name="CLHEP"/>
+<export>
+  <lib   name="1"/>
+</export>

SimCalorimetry/HGCalSimProducers/interface/HGCDigitizer.h

@@ -0,0 +1,103 @@
+#ifndef HGCalSimProducers_HGCDigitizer_h
+#define HGcalSimProducers_HGCDigitizer_h
+
+#include "DataFormats/DetId/interface/DetId.h"
+#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+
+#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+
+#include "SimCalorimetry/HGCalSimProducers/interface/HGCEEDigitizer.h"
+#include "SimCalorimetry/HGCalSimProducers/interface/HGCHEfrontDigitizer.h"
+#include "SimCalorimetry/HGCalSimProducers/interface/HGCHEbackDigitizer.h"
+#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+
+#include <vector>
+#include <map>
+#include <memory>
+
+class PCaloHit;
+class PileUpEventPrincipal;
+
+namespace edm {
+  class ConsumesCollector;
+}
+
+namespace CLHEP {
+  class HepRandomEngine;
+}
+
+class HGCDigitizer {
+
+public:
+
+  explicit HGCDigitizer(const edm::ParameterSet& ps, edm::ConsumesCollector& iC);
+  ~HGCDigitizer() { }
+
+  /**
+     @short handle SimHit accumulation
+   */
+  void accumulate(edm::Event const& e, edm::EventSetup const& c, CLHEP::HepRandomEngine*);
+  void accumulate(PileUpEventPrincipal const& e, edm::EventSetup const& c, CLHEP::HepRandomEngine*);
+  void accumulate(edm::Handle<edm::PCaloHitContainer> const &hits, int bxCrossing,const edm::ESHandle<HGCalGeometry> &geom, CLHEP::HepRandomEngine*);
+
+  /**
+     @short actions at the start/end of event
+   */
+  void initializeEvent(edm::Event const& e, edm::EventSetup const& c);
+  void finalizeEvent(edm::Event& e, edm::EventSetup const& c, CLHEP::HepRandomEngine*);
+
+  /**
+   */
+  bool producesEEDigis()       { return (mySubDet_==ForwardSubdetector::HGCEE);  }
+  bool producesHEfrontDigis()  { return (mySubDet_==ForwardSubdetector::HGCHEF); }
+  bool producesHEbackDigis()   { return (mySubDet_==ForwardSubdetector::HGCHEB); }
+  std::string digiCollection() { return digiCollection_; }
+
+  /**
+      @short actions at the start/end of run
+   */
+  void beginRun(const edm::EventSetup & es);
+  void endRun();
+
+private :
+
+  //used for initialization
+  bool checkValidDetIds_;
+
+  //input/output names
+  std::string hitCollection_,digiCollection_;
+
+  //digitization type (it's up to the specializations to decide it's meaning)
+  int digitizationType_;
+
+  //handle sim hits
+  int maxSimHitsAccTime_;
+  int bxTime_;
+  std::unique_ptr<HGCSimHitDataAccumulator> simHitAccumulator_;  
+  void resetSimHitDataAccumulator();
+
+  //digitizers
+  std::unique_ptr<HGCEEDigitizer>      theHGCEEDigitizer_;
+  std::unique_ptr<HGCHEbackDigitizer>  theHGCHEbackDigitizer_;
+  std::unique_ptr<HGCHEfrontDigitizer> theHGCHEfrontDigitizer_;
+
+  //subdetector id
+  ForwardSubdetector mySubDet_;
+
+  //misc switches
+  bool useAllChannels_;
+  int  verbosity_;
+
+  //reference speed to evaluate time of arrival at the sensititive detector, assuming the center of CMS
+  float refSpeed_;
+
+  //delay to apply after evaluating time of arrival at the sensitive detector
+  float tofDelay_;
+};
+
+#endif

SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h

@@ -0,0 +1,167 @@
+#ifndef SimCalorimetry_HGCSimProducers_hgcdigitizerbase
+#define SimCalorimetry_HGCSimProducers_hgcdigitizerbase
+
+#include <array>
+#include <iostream>
+#include <vector>
+#include <memory>
+#include <unordered_map>
+
+#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
+#include "FWCore/Utilities/interface/EDMException.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "CLHEP/Random/RandGaussQ.h"
+
+typedef float HGCSimEn_t;
+typedef std::array<HGCSimEn_t,6> HGCSimHitData;
+typedef std::unordered_map<uint32_t, HGCSimHitData> HGCSimHitDataAccumulator;
+
+template <class D>
+class HGCDigitizerBase {
+
+public:
+
+  typedef edm::SortedCollection<D> DColl;
+
+  /**
+     @short CTOR
+  */
+ HGCDigitizerBase(const edm::ParameterSet &ps) {
+    myCfg_         = ps.getParameter<edm::ParameterSet>("digiCfg"); 
+    bxTime_        = ps.getParameter<int32_t>("bxTime");
+    doTimeSamples_ = myCfg_.getParameter< bool >("doTimeSamples");
+    mipInKeV_      = myCfg_.getParameter<double>("mipInKeV");
+    lsbInMIP_      = myCfg_.getParameter<double>("lsbInMIP");
+    mip2noise_     = myCfg_.getParameter<double>("mip2noise");
+    adcThreshold_  = myCfg_.getParameter< uint32_t >("adcThreshold");
+    shaperN_       = myCfg_.getParameter< double >("shaperN");
+    shaperTau_     = myCfg_.getParameter< double >("shaperTau");
+  }
+
+  /**
+     @short steer digitization mode
+   */
+  void run(std::auto_ptr<DColl> &digiColl,HGCSimHitDataAccumulator &simData,uint32_t digitizationType, CLHEP::HepRandomEngine* engine) {
+    if(digitizationType==0) runTrivial(digiColl,simData,engine);
+    else                    runDigitizer(digiColl,simData,digitizationType,engine);
+  }
+
+
+  /**
+     @short a trivial digitization: sum energies and digitize without noise
+  */
+  void runTrivial(std::auto_ptr<DColl> &coll,HGCSimHitDataAccumulator &simData, CLHEP::HepRandomEngine* engine) {
+    for(HGCSimHitDataAccumulator::iterator it=simData.begin();
+	it!=simData.end(); it++) {
+
+      //create a new data frame
+      D rawDataFrame( it->first );
+
+      for(size_t i=0; i<it->second.size(); i++) {
+	//convert total energy GeV->keV->ADC counts
+	double totalEn=(it->second)[i]*1e6;
+
+	//add noise (in keV)
+	double noiseEn=CLHEP::RandGaussQ::shoot(engine,0,mipInKeV_/mip2noise_);
+	totalEn += noiseEn;
+	if(totalEn<0) totalEn=0;
+
+	//round to integer (sample will saturate the value according to available bits)
+	uint16_t totalEnInt = floor( (totalEn/mipInKeV_) / lsbInMIP_ );
+
+	//0 gain for the moment
+	HGCSample singleSample;
+	singleSample.set(0, totalEnInt);
+
+	rawDataFrame.setSample(i, singleSample);
+      }
+
+      //run the shaper
+      runShaper(rawDataFrame);
+
+      //update the output according to the final shape
+      updateOutput(coll,rawDataFrame);
+    }
+  }
+
+  /**
+     @short prepares the output according to the number of time samples to produce
+  */
+  void updateOutput(std::auto_ptr<DColl> &coll,D rawDataFrame) {
+    size_t itIdx(4); //index to the in-time digi - this could be configurable in a future version
+
+    //check if in-time sample is above threshold and put result into the event
+    if(doTimeSamples_) {
+      if(rawDataFrame[itIdx].adc() < adcThreshold_ ) return;
+      coll->push_back(rawDataFrame);
+    } else {
+      //create a new data frame containing only the in-time digi
+      D singleRawDataFrame( rawDataFrame.id() );
+      singleRawDataFrame.resize(1);
+
+      HGCSample singleSample;
+      singleSample.set(rawDataFrame[itIdx].gain(),rawDataFrame[itIdx].adc());
+      singleRawDataFrame.setSample(0, singleSample);
+      if(singleRawDataFrame[0].adc() < adcThreshold_ ) return;
+      coll->push_back(singleRawDataFrame);
+    }
+  }
+
+  /**
+     @short applies a shape to each time sample and propagates the tails to the subsequent time samples
+   */
+  void runShaper(D &dataFrame) {
+    std::vector<uint16_t> oldADC(dataFrame.size());
+    for(int it=0; it<dataFrame.size(); it++) {
+      uint16_t gain=dataFrame[it].gain();
+      oldADC[it]=dataFrame[it].adc();
+      uint16_t newADC(oldADC[it]);
+
+      if(shaperN_*shaperTau_>0){
+	for(int jt=0; jt<it; jt++) {
+	  float relTime(bxTime_*(it-jt)+shaperN_*shaperTau_);	
+	  newADC += uint16_t(oldADC[jt]*pow(relTime/(shaperN_*shaperTau_),shaperN_)*exp(-(relTime-shaperN_*shaperTau_)/shaperTau_));	      
+	}
+      }
+
+      HGCSample newSample;
+      newSample.set(gain,newADC);
+      dataFrame.setSample(it,newSample);
+    }
+  }
+
+  /**
+     @short to be specialized by top class
+   */
+  virtual void runDigitizer(std::auto_ptr<DColl> &coll,HGCSimHitDataAccumulator &simData,uint32_t digitizerType, CLHEP::HepRandomEngine* engine) {
+    throw cms::Exception("HGCDigitizerBaseException") << " Failed to find specialization of runDigitizer";
+  }
+
+  /**
+     @short DTOR
+  */
+  ~HGCDigitizerBase() { };
+
+  //baseline configuration
+  edm::ParameterSet myCfg_;
+
+  //minimum ADC counts to produce DIGIs
+  uint32_t adcThreshold_;
+
+  //parameters for the trivial digitization scheme
+  double mipInKeV_, lsbInMIP_, mip2noise_;
+
+  //parameters for trivial shaping
+  double shaperN_, shaperTau_;
+
+  //bunch time
+  int bxTime_;
+
+  //if true will put both in time and out-of-time samples in the event
+  bool doTimeSamples_;
+
+private:
+
+};
+
+#endif

SimCalorimetry/HGCalSimProducers/interface/HGCEEDigitizer.h

@@ -0,0 +1,17 @@
+#ifndef SimCalorimetry_HGCSimProducers_hgceedigitizer
+#define SimCalorimetry_HGCSimProducers_hgceedigitizer
+
+#include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h"
+#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
+
+class HGCEEDigitizer : public HGCDigitizerBase<HGCEEDataFrame> {
+
+public:
+  HGCEEDigitizer(const edm::ParameterSet& ps);
+  void runDigitizer(std::auto_ptr<HGCEEDigiCollection> &digiColl,HGCSimHitDataAccumulator &simData,uint32_t digitizationType, CLHEP::HepRandomEngine* engine);
+  ~HGCEEDigitizer();
+private:
+
+};
+
+#endif 

SimCalorimetry/HGCalSimProducers/interface/HGCHEbackDigitizer.h

@@ -0,0 +1,23 @@
+#ifndef SimCalorimetry_HGCSimProducers_hgchebackdigitizer
+#define SimCalorimetry_HGCSimProducers_hgchebackdigitizer
+
+#include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h"
+#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
+
+
+class HGCHEbackDigitizer : public HGCDigitizerBase<HGCHEDataFrame> {
+
+public:
+
+  HGCHEbackDigitizer(const edm::ParameterSet& ps);
+  void runDigitizer(std::auto_ptr<HGCHEDigiCollection> &digiColl,HGCSimHitDataAccumulator &simData,uint32_t digitizationType, CLHEP::HepRandomEngine* engine);
+  ~HGCHEbackDigitizer();
+
+private:
+
+  //calice-like digitization parameters
+  float nPEperMIP_, nTotalPE_, xTalk_, sdPixels_;
+  void runCaliceLikeDigitizer(std::auto_ptr<HGCHEDigiCollection> &digiColl,HGCSimHitDataAccumulator &simData, CLHEP::HepRandomEngine* engine);
+};
+
+#endif 

SimCalorimetry/HGCalSimProducers/interface/HGCHEfrontDigitizer.h

@@ -0,0 +1,17 @@
+#ifndef SimCalorimetry_HGCSimProducers_hgchefrontdigitizer
+#define SimCalorimetry_HGCSimProducers_hgchefrontdigitizer
+
+#include "SimCalorimetry/HGCalSimProducers/interface/HGCDigitizerBase.h"
+#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
+
+class HGCHEfrontDigitizer : public HGCDigitizerBase<HGCHEDataFrame> {
+
+public:
+  HGCHEfrontDigitizer(const edm::ParameterSet& ps);
+  void runDigitizer(std::auto_ptr<HGCHEDigiCollection> &digiColl,HGCSimHitDataAccumulator &simData,uint32_t digitizationType, CLHEP::HepRandomEngine* engine);
+  ~HGCHEfrontDigitizer();
+private:
+
+};
+
+#endif 

SimCalorimetry/HGCalSimProducers/plugins/BuildFile.xml

@@ -0,0 +1,18 @@
+<use name="FWCore/Framework"/>
+<use name="FWCore/Utilities"/>
+<use name="Geometry/Records"/>
+<use name="MagneticField/Records"/>
+<use name="DetectorDescription/Core"/>
+<use name="DetectorDescription/Parser"/>
+<use name="DetectorDescription/OfflineDBLoader"/>
+<use name="SimG4CMS/Calo"/>
+<use name="CondFormats/GeometryObjects"/>
+<use name="DataFormats/ForwardDetId"/>
+<use name="DataFormats/Candidate"/>
+<use name="DataFormats/HepMCCandidate"/>
+<use name="SimGeneral/MixingModule"/>
+<use name="CommonTools/UtilAlgos"/>
+<use name="SimCalorimetry/HGCalSimProducers"/>
+<use name="geant4core"/>
+<use name="root"/>
+<flags EDM_PLUGIN="1"/>