diff --git a/Configuration/Eras/python/Era_Phase2C1_timing_cff.py b/Configuration/Eras/python/Era_Phase2C1_timing_cff.py new file mode 100644 index 0000000000000..a567047a044d6 --- /dev/null +++ b/Configuration/Eras/python/Era_Phase2C1_timing_cff.py @@ -0,0 +1,7 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.Eras.Era_Phase2C1_cff import Phase2C1 +from Configuration.Eras.Modifier_phase2_timing_cff import phase2_timing + +Phase2C1_timing = cms.ModifierChain(Phase2C1, phase2_timing) + diff --git a/Configuration/Eras/python/Era_Phase2C2_timing_cff.py b/Configuration/Eras/python/Era_Phase2C2_timing_cff.py new file mode 100644 index 0000000000000..abcc898538be4 --- /dev/null +++ b/Configuration/Eras/python/Era_Phase2C2_timing_cff.py @@ -0,0 +1,7 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.Eras.Era_Phase2C2_cff import Phase2C2 +from Configuration.Eras.Modifier_phase2_timing_cff import phase2_timing + +Phase2C2_timing = cms.ModifierChain(Phase2C2, phase2_timing) + diff --git a/Configuration/Eras/python/Modifier_phase2_timing_cff.py b/Configuration/Eras/python/Modifier_phase2_timing_cff.py new file mode 100644 index 0000000000000..3b0c24c641185 --- /dev/null +++ b/Configuration/Eras/python/Modifier_phase2_timing_cff.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms + +phase2_timing = cms.Modifier() + diff --git a/Configuration/PyReleaseValidation/python/relval_steps.py b/Configuration/PyReleaseValidation/python/relval_steps.py index f1403532d935c..dc56ab0e2b839 100644 --- a/Configuration/PyReleaseValidation/python/relval_steps.py +++ b/Configuration/PyReleaseValidation/python/relval_steps.py @@ -1603,6 +1603,9 @@ def lhegensim(fragment,howMuch): defaultDataSets['2023D1']='' defaultDataSets['2023D2']='' defaultDataSets['2023D3']='' +defaultDataSets['2023D1Timing']='' +defaultDataSets['2023D2Timing']='' +defaultDataSets['2023D3Timing']='' keys=defaultDataSets.keys() for key in keys: diff --git a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py index 895cafb4d48f1..8abd401c087ee 100644 --- a/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py +++ b/Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py @@ -13,9 +13,15 @@ '2023D1', '2023D1PU', '2023D2', - '2023D2PU', - '2023D3', + '2023D2PU', + '2023D3', '2023D3PU', + '2023D1Timing', + '2023D1TimingPU', + '2023D2Timing', + '2023D2TimingPU', + '2023D3Timing', + '2023D3TimingPU' ] upgradeSteps=[ @@ -62,23 +68,35 @@ 'Custom' : 'SLHCUpgradeSimulations/Configuration/combinedCustoms.cust_2023tilted', 'Era' : 'Phase2C1', 'ScenToRun' : ['GenSimFull','DigiFull','RecoFullGlobal','HARVESTFullGlobal'], - }, + }, '2023D2' : { 'Geom' : 'Extended2023D2', 'GT' : 'auto:run2_mc', 'Custom' : 'SLHCUpgradeSimulations/Configuration/combinedCustoms.cust_2023flat', 'Era' : 'Phase2C1', 'ScenToRun' : ['GenSimFull','DigiFull','RecoFullGlobal','HARVESTFullGlobal'], - }, + }, '2023D3' : { 'Geom' : 'Extended2023D3', 'GT' : 'auto:run2_mc', 'Custom' : 'SLHCUpgradeSimulations/Configuration/combinedCustoms.cust_2023tilted', 'Era' : 'Phase2C2', 'ScenToRun' : ['GenSimFull','DigiFull','RecoFullGlobal', 'HARVESTFullGlobal'], - }, + } } + +#Timing (later we can alter geometry, etc, if need be) +upgradeProperties[2023]['2023D1Timing'] = deepcopy(upgradeProperties[2023]['2023D1']) +upgradeProperties[2023]['2023D1Timing']['Era'] = 'Phase2C1_timing' +upgradeProperties[2023]['2023D2Timing'] = deepcopy(upgradeProperties[2023]['2023D2']) +upgradeProperties[2023]['2023D2Timing']['Era'] = 'Phase2C1_timing' +upgradeProperties[2023]['2023D3Timing'] = deepcopy(upgradeProperties[2023]['2023D3']) +upgradeProperties[2023]['2023D3Timing']['Era'] = 'Phase2C2_timing' + +#standard PU sequences +upgradeProperties[2023]['2017PU'] = deepcopy(upgradeProperties[2017]['2017']) +upgradeProperties[2023]['2017PU']['ScenToRun'] = ['GenSimFull','DigiFullPU','RecoFullPU','HARVESTFullPU'] upgradeProperties[2023]['2023D1PU'] = deepcopy(upgradeProperties[2023]['2023D1']) upgradeProperties[2023]['2023D1PU']['ScenToRun'] = ['GenSimFull','DigiFullPU','RecoFullGlobalPU', 'HARVESTFullGlobalPU'] upgradeProperties[2023]['2023D2PU'] = deepcopy(upgradeProperties[2023]['2023D2']) @@ -86,6 +104,15 @@ upgradeProperties[2023]['2023D3PU'] = deepcopy(upgradeProperties[2023]['2023D3']) upgradeProperties[2023]['2023D3PU']['ScenToRun'] = ['GenSimFull','DigiFullPU','RecoFullGlobalPU', 'HARVESTFullGlobalPU'] +#Timing PU (for now copy ScenToRun of standard PU) +upgradeProperties[2023]['2023D1TimingPU'] = deepcopy(upgradeProperties[2023]['2023D1Timing']) +upgradeProperties[2023]['2023D1TimingPU']['ScenToRun'] = deepcopy(upgradeProperties[2023]['2023D1PU']['ScenToRun']) +upgradeProperties[2023]['2023D2TimingPU'] = deepcopy(upgradeProperties[2023]['2023D2Timing']) +upgradeProperties[2023]['2023D2TimingPU']['ScenToRun'] = deepcopy(upgradeProperties[2023]['2023D2PU']['ScenToRun']) +upgradeProperties[2023]['2023D3TimingPU'] = deepcopy(upgradeProperties[2023]['2023D3Timing']) +upgradeProperties[2023]['2023D3TimingPU']['ScenToRun'] = deepcopy(upgradeProperties[2023]['2023D3PU']['ScenToRun']) + + from Configuration.PyReleaseValidation.relval_steps import Kby upgradeFragments=['FourMuPt_1_200_pythia8_cfi', diff --git a/Configuration/StandardSequences/python/Eras.py b/Configuration/StandardSequences/python/Eras.py index dacd32adfee76..db3af16ddff9d 100644 --- a/Configuration/StandardSequences/python/Eras.py +++ b/Configuration/StandardSequences/python/Eras.py @@ -20,7 +20,9 @@ def __init__(self): 'Run2_2017_trackingPhase1PU70', 'Run3', 'Phase2C1', - 'Phase2C2'] + 'Phase2C2', + 'Phase2C1_timing', + 'Phase2C2_timing'] internalUseMods = ['run2_common', 'run2_25ns_specific', 'run2_50ns_specific', 'run2_HI_specific', @@ -28,12 +30,11 @@ def __init__(self): 'run2_HE_2017', 'stage2L1Trigger', 'phase1Pixel', 'run3_GEM', 'phase2_common', 'phase2_tracker', - 'phase2_hgcal', 'phase2_muon', + 'phase2_hgcal', 'phase2_muon', 'phase2_timing', 'phase2_hcal', 'trackingLowPU', 'trackingPhase1', 'trackingPhase1PU70', 'ctpps_2016', 'trackingPhase2PU140', 'tracker_apv_vfp30_2016'] - for e in allEras: eObj=getattr(__import__('Configuration.Eras.Era_'+e+'_cff',globals(),locals(),[e],0),e) self.addEra(e,eObj) @@ -78,8 +79,7 @@ def inspect(self,name=None,onlyChosen=False,details=True): if name is None: if not onlyChosen or getattr(self,e).isChosen(): self.inspectEra(e,details) - - + eras=Eras() diff --git a/DataFormats/CaloRecHit/interface/CaloRecHit.h b/DataFormats/CaloRecHit/interface/CaloRecHit.h index 6cd565a3c1ee0..82e06705822f0 100644 --- a/DataFormats/CaloRecHit/interface/CaloRecHit.h +++ b/DataFormats/CaloRecHit/interface/CaloRecHit.h @@ -17,6 +17,7 @@ class CaloRecHit { float energy() const { return energy_; } void setEnergy(float energy) { energy_=energy; } float time() const { return time_; } + void setTime(float time) { time_=time; } const DetId& detid() const { return id_; } uint32_t flags() const { return flags_; } void setFlags(uint32_t flags) { flags_=flags; } diff --git a/DataFormats/EcalDigi/interface/EcalDigiCollections.h b/DataFormats/EcalDigi/interface/EcalDigiCollections.h index a1fd2d33ec976..34071b3fa36b4 100644 --- a/DataFormats/EcalDigi/interface/EcalDigiCollections.h +++ b/DataFormats/EcalDigi/interface/EcalDigiCollections.h @@ -4,6 +4,7 @@ #include "DataFormats/EcalDigi/interface/EBDataFrame.h" #include "DataFormats/EcalDigi/interface/EEDataFrame.h" #include "DataFormats/EcalDigi/interface/ESDataFrame.h" +#include "DataFormats/EcalDigi/interface/EcalTimeDigi.h" #include "DataFormats/EcalDigi/interface/EcalTriggerPrimitiveDigi.h" #include "DataFormats/EcalDigi/interface/EcalTrigPrimCompactColl.h" #include "DataFormats/EcalDigi/interface/EcalPseudoStripInputDigi.h" @@ -104,6 +105,7 @@ void swap(ESDigiCollection& lhs, ESDigiCollection& rhs) { lhs.swap(rhs); } +typedef edm::SortedCollection EcalTimeDigiCollection; typedef edm::SortedCollection EcalTrigPrimDigiCollection; typedef edm::SortedCollection EcalPSInputDigiCollection; diff --git a/DataFormats/EcalDigi/interface/EcalTimeDigi.h b/DataFormats/EcalDigi/interface/EcalTimeDigi.h new file mode 100644 index 0000000000000..eb649f886a73d --- /dev/null +++ b/DataFormats/EcalDigi/interface/EcalTimeDigi.h @@ -0,0 +1,50 @@ +#ifndef _DataFormats_EcalDigi_ECALTIMEDIGI_H_ +#define _DataFormats_EcalDigi_ECALTIMEDIGI_H_ + +#include +#include +#include "DataFormats/DetId/interface/DetId.h" + +class EcalTimeDigi { + public: + typedef DetId key_type; ///< For the sorted collection + + EcalTimeDigi(); // for persistence + explicit EcalTimeDigi(const DetId& id); + + void swap(EcalTimeDigi& rh) { + std::swap(id_,rh.id_); + std::swap(size_,rh.size_); + std::swap(data_,rh.data_); + } + + const DetId& id() const { return id_; } + int size() const { return size_; } + + const float& operator[](unsigned int i) const { return data_[i]; } + const float& sample(unsigned int i) const { return data_[i]; } + + void setSize(unsigned int size); + void setSample(unsigned int i, const float sam) { data_[i]=sam; } + void setSampleOfInterest(int i) { sampleOfInterest_=i; } + + /// Gets the BX==0 sample. If =-1 then it means that only OOT hits are present + int sampleOfInterest() const { return sampleOfInterest_; } + +private: + DetId id_; + unsigned int size_; + int sampleOfInterest_; + std::vector data_; +}; + + +inline void swap(EcalTimeDigi& lh, EcalTimeDigi& rh) { + lh.swap(rh); +} + +std::ostream& operator<<(std::ostream& s, const EcalTimeDigi& digi); + + + +#endif diff --git a/DataFormats/EcalDigi/src/EcalTimeDigi.cc b/DataFormats/EcalDigi/src/EcalTimeDigi.cc new file mode 100644 index 0000000000000..82a71fe6e4ad9 --- /dev/null +++ b/DataFormats/EcalDigi/src/EcalTimeDigi.cc @@ -0,0 +1,22 @@ +#include "DataFormats/EcalDigi/interface/EcalTimeDigi.h" + +namespace { + constexpr unsigned int MAXSAMPLES = 10; +} + +EcalTimeDigi::EcalTimeDigi() : id_(0), size_(0), sampleOfInterest_(-1), data_(MAXSAMPLES) { +} + +EcalTimeDigi::EcalTimeDigi(const DetId& id) : id_(id), + size_(0), sampleOfInterest_(-1), data_(MAXSAMPLES) { +} + +void EcalTimeDigi::setSize(unsigned int size) { + if (size>MAXSAMPLES) size_=MAXSAMPLES; + else size_=size; + data_.resize(size_); +} + + + + diff --git a/DataFormats/EcalDigi/src/classes.h b/DataFormats/EcalDigi/src/classes.h index 907731a3ed953..0f52217473aa9 100644 --- a/DataFormats/EcalDigi/src/classes.h +++ b/DataFormats/EcalDigi/src/classes.h @@ -6,11 +6,14 @@ namespace DataFormats_EcalDigi { std::vector vMGPA_; std::vector vFEM_; std::vector vESSample_; + std::vector vETS_; std::vector vETPS_; std::vector vEPSIS_; std::vector vMD_; + std::vector vTD_; edm::SortedCollection vES_; + edm::SortedCollection vETDP_; edm::SortedCollection vETP_; edm::SortedCollection vEPSI_; edm::SortedCollection vEBSRF_; @@ -18,10 +21,12 @@ namespace DataFormats_EcalDigi { edm::SortedCollection vEPN_; edm::SortedCollection vMDS_; EcalMatacqDigi Matacq_; + EcalTimeDigi Time_; EBDigiCollection theEB_; EEDigiCollection theEE_; ESDigiCollection theES_; + EcalTimeDigiCollection theEBTime_; EcalTrigPrimDigiCollection theETP_; EcalTrigPrimCompactColl theETP2_; @@ -34,6 +39,7 @@ namespace DataFormats_EcalDigi { edm::Wrapper anotherEBw_; edm::Wrapper anotherEEw_; edm::Wrapper anotherESw_; + edm::Wrapper anotherETDw_; edm::Wrapper anotherETPw_; edm::Wrapper anotherETP2w_; edm::Wrapper anotherEBSRFw_; @@ -42,6 +48,7 @@ namespace DataFormats_EcalDigi { edm::Wrapper anotherMDw_; edm::Wrapper< edm::SortedCollection > theESw_; + edm::Wrapper< edm::SortedCollection > theETDw_; edm::Wrapper< edm::SortedCollection > theETPw_; edm::Wrapper< edm::SortedCollection > theEPSIw_; edm::Wrapper< edm::SortedCollection > theEBSRFw_; diff --git a/DataFormats/EcalDigi/src/classes_def.xml b/DataFormats/EcalDigi/src/classes_def.xml index bc104d4307322..718e90bf61503 100644 --- a/DataFormats/EcalDigi/src/classes_def.xml +++ b/DataFormats/EcalDigi/src/classes_def.xml @@ -44,9 +44,14 @@ + + + + + @@ -56,6 +61,7 @@ + @@ -83,6 +89,7 @@ + @@ -97,6 +104,7 @@ + diff --git a/DataFormats/EcalRecHit/interface/EcalRecHit.h b/DataFormats/EcalRecHit/interface/EcalRecHit.h index 5882dac8dabb6..27d149fbe1a07 100644 --- a/DataFormats/EcalRecHit/interface/EcalRecHit.h +++ b/DataFormats/EcalRecHit/interface/EcalRecHit.h @@ -68,6 +68,7 @@ class EcalRecHit { float energy() const { return energy_; } void setEnergy(float energy) { energy_=energy; } float time() const { return time_; } + void setTime(float time) { time_ = time; } const DetId& detid() const { return id_; } /// get the id diff --git a/RecoLocalCalo/Configuration/python/ecalLocalRecoSequence_cff.py b/RecoLocalCalo/Configuration/python/ecalLocalRecoSequence_cff.py index cb0ebd7026905..0a8d6440ab982 100644 --- a/RecoLocalCalo/Configuration/python/ecalLocalRecoSequence_cff.py +++ b/RecoLocalCalo/Configuration/python/ecalLocalRecoSequence_cff.py @@ -18,6 +18,7 @@ from RecoLocalCalo.EcalRecProducers.ecalCompactTrigPrim_cfi import * from RecoLocalCalo.EcalRecProducers.ecalTPSkim_cfi import * +from RecoLocalCalo.EcalRecProducers.ecalDetailedTimeRecHit_cfi import * #ecalUncalibRecHitSequence = cms.Sequence(ecalGlobalUncalibRecHit* # ecalDetIdToBeRecovered) @@ -33,3 +34,7 @@ ecalLocalRecoSequence = cms.Sequence(ecalUncalibRecHitSequence* ecalRecHitSequence) +from Configuration.StandardSequences.Eras import eras +from RecoLocalCalo.EcalRecProducers.ecalDetailedTimeRecHit_cfi import * +_phase2_timing_ecalRecHitSequence = cms.Sequence( ecalRecHitSequence.copy() + ecalDetailedTimeRecHit ) +eras.phase2_timing.toReplaceWith( ecalRecHitSequence, _phase2_timing_ecalRecHitSequence ) diff --git a/RecoLocalCalo/Configuration/python/ecalLocalReco_EventContent_cff.py b/RecoLocalCalo/Configuration/python/ecalLocalReco_EventContent_cff.py index afa1e18acb6b4..9c335c7f14749 100644 --- a/RecoLocalCalo/Configuration/python/ecalLocalReco_EventContent_cff.py +++ b/RecoLocalCalo/Configuration/python/ecalLocalReco_EventContent_cff.py @@ -30,3 +30,12 @@ ) ) +from Configuration.StandardSequences.Eras import eras +#mods for timing +_phase2_timing_EcalOutputCommands = ['keep *_mix_EBTimeDigi_*', + 'keep *_mix_EETimeDigi_*', + 'keep *_ecalDetailedTimeRecHit_*_*'] + +eras.phase2_timing.toModify( ecalLocalRecoFEVT, outputCommands = ecalLocalRecoFEVT.outputCommands + _phase2_timing_EcalOutputCommands ) +eras.phase2_timing.toModify( ecalLocalRecoRECO, outputCommands = ecalLocalRecoRECO.outputCommands + _phase2_timing_EcalOutputCommands ) + diff --git a/RecoLocalCalo/EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.cc b/RecoLocalCalo/EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.cc new file mode 100644 index 0000000000000..a916f6177517f --- /dev/null +++ b/RecoLocalCalo/EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.cc @@ -0,0 +1,238 @@ +/** \class EcalDetailedTimeRecHitProducer + * produce ECAL detailed time Rechits + * \author Paolo Meridiani + * + **/ + +#include "RecoLocalCalo/EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.h" +#include "RecoLocalCalo/EcalRecAlgos/interface/EcalRecHitSimpleAlgo.h" +#include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h" +#include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h" +#include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h" +#include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h" +#include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h" +#include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h" + +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include +#include +#include + +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/CaloGeometry/interface/CaloGenericDetId.h" +#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h" +#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h" + +#include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHit.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h" +#include "DataFormats/EcalDigi/interface/EcalTimeDigi.h" + + + +#include "CLHEP/Units/GlobalPhysicalConstants.h" +#include "CLHEP/Units/GlobalSystemOfUnits.h" + +EcalDetailedTimeRecHitProducer::EcalDetailedTimeRecHitProducer(const edm::ParameterSet& ps) : + m_geometry(0) +{ + EBRecHitCollection_ = consumes( ps.getParameter("EBRecHitCollection") ); + EERecHitCollection_ = consumes( ps.getParameter("EERecHitCollection") ); + + ebTimeDigiCollection_ = consumes( ps.getParameter("EBTimeDigiCollection") ); + eeTimeDigiCollection_ = consumes( ps.getParameter("EETimeDigiCollection") ); + + EBDetailedTimeRecHitCollection_ = ps.getParameter("EBDetailedTimeRecHitCollection"); + EEDetailedTimeRecHitCollection_ = ps.getParameter("EEDetailedTimeRecHitCollection"); + + correctForVertexZPosition_ = ps.getParameter("correctForVertexZPosition"); + useMCTruthVertex_ = ps.getParameter("useMCTruthVertex"); + recoVertex_ = consumes( ps.getParameter("recoVertex") ); + simVertex_ = consumes( ps.getParameter("simVertex") ); + + ebTimeLayer_ = ps.getParameter("EBTimeLayer"); + eeTimeLayer_ = ps.getParameter("EETimeLayer"); + + produces< EBRecHitCollection >(EBDetailedTimeRecHitCollection_); + produces< EERecHitCollection >(EEDetailedTimeRecHitCollection_); +} + +EcalDetailedTimeRecHitProducer::~EcalDetailedTimeRecHitProducer() { +} + + +void EcalDetailedTimeRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) +{ + using namespace edm; + using namespace reco; + + edm::ESHandle hGeometry ; + es.get().get( hGeometry ) ; + + m_geometry = &*hGeometry; + + Handle< EBRecHitCollection > pEBRecHits; + Handle< EERecHitCollection > pEERecHits; + + const EBRecHitCollection* EBRecHits = 0; + const EERecHitCollection* EERecHits = 0; + + evt.getByToken( EBRecHitCollection_, pEBRecHits); + if ( pEBRecHits.isValid() ) { + EBRecHits = pEBRecHits.product(); // get a ptr to the product +#ifdef DEBUG + LogDebug("EcalRecHitDebug") << "total # EB rechits to be re-calibrated: " << EBRecHits->size(); +#endif + } + + evt.getByToken( EERecHitCollection_, pEERecHits); + if ( pEERecHits.isValid() ) { + EERecHits = pEERecHits.product(); // get a ptr to the product +#ifdef DEBUG + LogDebug("EcalRecHitDebug") << "total # EE uncalibrated rechits to be re-calibrated: " << EERecHits->size(); +#endif + } + + Handle< EcalTimeDigiCollection > pEBTimeDigis; + Handle< EcalTimeDigiCollection > pEETimeDigis; + + const EcalTimeDigiCollection* ebTimeDigis =0; + const EcalTimeDigiCollection* eeTimeDigis =0; + + evt.getByToken( ebTimeDigiCollection_, pEBTimeDigis); + //evt.getByToken( digiProducer_, pEBTimeDigis); + if ( pEBTimeDigis.isValid() ) { + ebTimeDigis = pEBTimeDigis.product(); // get a ptr to the produc + edm::LogInfo("EcalDetailedTimeRecHitInfo") << "total # ebTimeDigis: " << ebTimeDigis->size() ; + } + + evt.getByToken( eeTimeDigiCollection_, pEETimeDigis); + //evt.getByToken( digiProducer_, pEETimeDigis); + if ( pEETimeDigis.isValid() ) { + eeTimeDigis = pEETimeDigis.product(); // get a ptr to the product + edm::LogInfo("EcalDetailedTimeRecHitInfo") << "total # eeTimeDigis: " << eeTimeDigis->size() ; + } + // collection of rechits to put in the event + std::unique_ptr< EBRecHitCollection > EBDetailedTimeRecHits( new EBRecHitCollection ); + std::unique_ptr< EERecHitCollection > EEDetailedTimeRecHits( new EERecHitCollection ); + + std::unique_ptr vertex; + + if (correctForVertexZPosition_) + { + if (!useMCTruthVertex_) + { + //Get the first reco vertex + // get primary vertices + + edm::Handle VertexHandle; + evt.getByToken(recoVertex_, VertexHandle); + + if ( VertexHandle.isValid() ) + { + if ((*VertexHandle).size()>0) //at least 1 vertex + { + const reco::Vertex* myVertex= &(*VertexHandle)[0]; + vertex.reset( new GlobalPoint(myVertex->x(),myVertex->y(),myVertex->z()) ); + } + } + + } + else + { + edm::Handle VertexHandle; + evt.getByToken(simVertex_, VertexHandle); + + if ( VertexHandle.isValid() ) + { + if ((*VertexHandle).size()>0) //at least 1 vertex + { + assert ((*VertexHandle)[0].vertexId() == 0); + const SimVertex* myVertex= &(*VertexHandle)[0]; + vertex.reset( new GlobalPoint(myVertex->position().x(),myVertex->position().y(),myVertex->position().z()) ); + } + } + + } + } + + if (EBRecHits && ebTimeDigis) { + // loop over uncalibrated rechits to make calibrated ones + for(EBRecHitCollection::const_iterator it = EBRecHits->begin(); it != EBRecHits->end(); ++it) { + EcalRecHit aHit( (*it) ); + EcalTimeDigiCollection::const_iterator timeDigi=ebTimeDigis->find((*it).id()); + if (timeDigi!=ebTimeDigis->end()) + { + if (timeDigi->sampleOfInterest()>=0) + { + float myTime=(*timeDigi)[timeDigi->sampleOfInterest()]; + //Vertex corrected ToF + if (vertex) + { + aHit.setTime(myTime+deltaTimeOfFlight(*vertex,(*it).id(),ebTimeLayer_)); + } + else + //Uncorrected ToF + aHit.setTime(myTime); + } + } + // leave standard time if no timeDigi is associated (e.g. noise recHits) + EBDetailedTimeRecHits->push_back( aHit ); + } + } + + if (EERecHits && eeTimeDigis) + { + // loop over uncalibrated rechits to make calibrated ones + for(EERecHitCollection::const_iterator it = EERecHits->begin(); + it != EERecHits->end(); ++it) { + + EcalRecHit aHit( *it ); + EcalTimeDigiCollection::const_iterator timeDigi=eeTimeDigis->find((*it).id()); + if (timeDigi!=eeTimeDigis->end()) + { + if (timeDigi->sampleOfInterest()>=0) + { + float myTime=(*timeDigi)[timeDigi->sampleOfInterest()]; + //Vertex corrected ToF + if (vertex) + { + aHit.setTime(myTime+deltaTimeOfFlight(*vertex,(*it).id(),eeTimeLayer_)); + } + else + //Uncorrected ToF + aHit.setTime(myTime); + } + } + EEDetailedTimeRecHits->push_back( aHit ); + } + } + // put the collection of recunstructed hits in the event + LogInfo("EcalDetailedTimeRecHitInfo") << "total # EB rechits: " << EBDetailedTimeRecHits->size(); + LogInfo("EcalDetailedTimeRecHitInfo") << "total # EE rechits: " << EEDetailedTimeRecHits->size(); + + evt.put( std::move(EBDetailedTimeRecHits), EBDetailedTimeRecHitCollection_ ); + evt.put( std::move(EEDetailedTimeRecHits), EEDetailedTimeRecHitCollection_ ); + +} + +double EcalDetailedTimeRecHitProducer::deltaTimeOfFlight( GlobalPoint& vertex, const DetId& detId , int layer) const +{ + const CaloCellGeometry* cellGeometry ( m_geometry->getGeometry( detId ) ) ; + assert( 0 != cellGeometry ) ; + GlobalPoint layerPos = (dynamic_cast(cellGeometry))->getPosition( double(layer)+0.5 ); //depth in mm in the middle of the layer position + GlobalVector tofVector = layerPos-vertex; + return (layerPos.mag()*cm-tofVector.mag()*cm)/(float)c_light ; +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE( EcalDetailedTimeRecHitProducer ); diff --git a/RecoLocalCalo/EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.h b/RecoLocalCalo/EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.h new file mode 100644 index 0000000000000..9e4735f284151 --- /dev/null +++ b/RecoLocalCalo/EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.h @@ -0,0 +1,59 @@ +#ifndef RecoLocalCalo_EcalRecProducers_EcalDetailedTimeRecHitProducer_HH +#define RecoLocalCalo_EcalRecProducers_EcalDetailedTimeRecHitProducer_HH +/** \class EcalDetailedTimeRecHitProducer + * produce ECAL rechits associating them with improved ecalTimeDigis + * \author Paolo Meridiani, INFN Roma + **/ + +#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "SimDataFormats/Vertex/interface/SimVertex.h" +#include "SimDataFormats/Vertex/interface/SimVertexContainer.h" + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "RecoLocalCalo/EcalRecAlgos/interface/EcalRecHitAbsAlgo.h" + + +// forward declaration + +class CaloGeometry; + +class EcalDetailedTimeRecHitProducer : public edm::stream::EDProducer<> { + + public: + explicit EcalDetailedTimeRecHitProducer(const edm::ParameterSet& ps); + ~ EcalDetailedTimeRecHitProducer(); + virtual void produce(edm::Event& evt, const edm::EventSetup& es); + + private: + + //Functions to correct the TOF from the EcalDigi which is not corrected for the vertex position + double deltaTimeOfFlight( GlobalPoint& vertex, const DetId& detId , int layer) const ; + + const CaloGeometry* m_geometry; + + edm::EDGetTokenT EBRecHitCollection_; // secondary name given to collection of EBrechits + edm::EDGetTokenT EERecHitCollection_; // secondary name given to collection of EErechits + + edm::EDGetTokenT recoVertex_; + edm::EDGetTokenT simVertex_; + bool correctForVertexZPosition_; + bool useMCTruthVertex_; + + int ebTimeLayer_; + int eeTimeLayer_; + + edm::EDGetTokenT ebTimeDigiCollection_; // secondary name given to collection of EB uncalib rechits + edm::EDGetTokenT eeTimeDigiCollection_; // secondary name given to collection of EE uncalib rechits + + std::string EBDetailedTimeRecHitCollection_; // secondary name to be given to EB collection of hits + std::string EEDetailedTimeRecHitCollection_; // secondary name to be given to EE collection of hits + +}; +#endif diff --git a/RecoLocalCalo/EcalRecProducers/python/ecalDetailedTimeRecHit_cfi.py b/RecoLocalCalo/EcalRecProducers/python/ecalDetailedTimeRecHit_cfi.py new file mode 100644 index 0000000000000..ed13d3a44bdc4 --- /dev/null +++ b/RecoLocalCalo/EcalRecProducers/python/ecalDetailedTimeRecHit_cfi.py @@ -0,0 +1,18 @@ +import FWCore.ParameterSet.Config as cms + +from SimCalorimetry.EcalSimProducers.ecalTimeDigiParameters_cff import * + +ecalDetailedTimeRecHit = cms.EDProducer("EcalDetailedTimeRecHitProducer", + EERecHitCollection = cms.InputTag("ecalRecHit","EcalRecHitsEE"), + EBRecHitCollection = cms.InputTag("ecalRecHit","EcalRecHitsEB"), + EETimeDigiCollection = cms.InputTag("mix","EETimeDigi"), + EBTimeDigiCollection = cms.InputTag("mix","EBTimeDigi"), + EBDetailedTimeRecHitCollection = cms.string('EcalRecHitsEB'), + EEDetailedTimeRecHitCollection = cms.string('EcalRecHitsEE'), + EBTimeLayer = ecal_time_digi_parameters.timeLayerBarrel, + EETimeLayer = ecal_time_digi_parameters.timeLayerEndcap, + correctForVertexZPosition=cms.bool(True), + useMCTruthVertex=cms.bool(True), + recoVertex = cms.InputTag("offlinePrimaryVerticesWithBS"), + simVertex = cms.InputTag("g4SimHits") + ) diff --git a/RecoParticleFlow/PFSimProducer/plugins/BuildFile.xml b/RecoParticleFlow/PFSimProducer/plugins/BuildFile.xml index efe3bb3873038..a81770913d136 100644 --- a/RecoParticleFlow/PFSimProducer/plugins/BuildFile.xml +++ b/RecoParticleFlow/PFSimProducer/plugins/BuildFile.xml @@ -1,4 +1,4 @@ - + diff --git a/RecoParticleFlow/PFSimProducer/plugins/ConfigurableFlatResolutionModel.cc b/RecoParticleFlow/PFSimProducer/plugins/ConfigurableFlatResolutionModel.cc new file mode 100644 index 0000000000000..ab0a8f9c42328 --- /dev/null +++ b/RecoParticleFlow/PFSimProducer/plugins/ConfigurableFlatResolutionModel.cc @@ -0,0 +1,19 @@ +#include "ResolutionModel.h" + +class ConfigurableFlatResolutionModel : public ResolutionModel { +public: + ConfigurableFlatResolutionModel( const edm::ParameterSet& conf ) : + ResolutionModel( conf ), + reso_( conf.getParameter("resolutionInNs") ) { + } + + virtual float getTimeResolution(const reco::Track&) const override { return reso_; } + virtual float getTimeResolution(const reco::PFCluster&) const override { return reso_; } + +private: + const float reso_; +}; + +DEFINE_EDM_PLUGIN(ResolutionModelFactory, + ConfigurableFlatResolutionModel, + "ConfigurableFlatResolutionModel"); diff --git a/RecoParticleFlow/PFSimProducer/plugins/EcalBarrelClusterFastTimer.cc b/RecoParticleFlow/PFSimProducer/plugins/EcalBarrelClusterFastTimer.cc new file mode 100644 index 0000000000000..c3faa5d6bc898 --- /dev/null +++ b/RecoParticleFlow/PFSimProducer/plugins/EcalBarrelClusterFastTimer.cc @@ -0,0 +1,221 @@ +// This producer eats standard PF ECAL clusters +// finds the corresponding fast-timing det IDs and attempts to +// assign a reasonable time guess + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "DataFormats/Common/interface/ValueMap.h" +#include "DataFormats/Common/interface/View.h" +#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h" +#include "DataFormats/ParticleFlowReco/interface/PFCluster.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" + +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" + +#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" + +#include + +#include "ResolutionModel.h" +#include "CLHEP/Units/SystemOfUnits.h" +#include "FWCore/Utilities/interface/isFinite.h" +#include "CLHEP/Random/RandGauss.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/RandomNumberGenerator.h" + +class EcalBarrelClusterFastTimer : public edm::global::EDProducer<> { +public: + EcalBarrelClusterFastTimer(const edm::ParameterSet&); + ~EcalBarrelClusterFastTimer() { } + + virtual void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; + +private: + // inputs + edm::EDGetTokenT ebTimeHitsToken_; + edm::EDGetTokenT > ebClustersToken_; + edm::EDGetTokenT timedVertexToken_; + // options + std::vector > _resolutions; + const float minFraction_, minEnergy_; + const unsigned ecalDepth_; + // functions + std::pair getTimeForECALPFCluster(const EcalRecHitCollection&,const reco::PFCluster&) const; + float correctTimeToVertex(const float intime, const DetId& timeDet, const reco::Vertex& vtx, + const CaloSubdetectorGeometry* ecalGeom) const; +}; + +DEFINE_FWK_MODULE(EcalBarrelClusterFastTimer); + +namespace { + static const std::string resolution("Resolution"); + + constexpr float cm_per_ns = 29.9792458; + + template + void writeValueMap(edm::Event &iEvent, + const edm::Handle > & handle, + const std::vector & values, + const std::string & label) { + std::auto_ptr > valMap(new edm::ValueMap()); + typename edm::ValueMap::Filler filler(*valMap); + filler.insert(handle, values.begin(), values.end()); + filler.fill(); + iEvent.put(valMap, label); + } +} + +EcalBarrelClusterFastTimer::EcalBarrelClusterFastTimer(const edm::ParameterSet& conf) : + ebTimeHitsToken_(consumes( conf.getParameter("ebTimeHits") ) ), + ebClustersToken_(consumes >( conf.getParameter("ebClusters") ) ), + timedVertexToken_(consumes( conf.getParameter("timedVertices") ) ), + minFraction_( conf.getParameter("minFractionToConsider") ), + minEnergy_(conf.getParameter("minEnergyToConsider") ), + ecalDepth_(conf.getParameter("ecalDepth") ) +{ + // setup resolution models + const std::vector& resos = conf.getParameterSetVector("resolutionModels"); + for( const auto& reso : resos ) { + const std::string& name = reso.getParameter("modelName"); + ResolutionModel* resomod = ResolutionModelFactory::get()->create(name,reso); + _resolutions.emplace_back( resomod ); + + // times and time resolutions for general tracks + produces > >(name); // indexed by primary vertex then PF cluster + produces >(name+resolution); + } + // get RNG engine + edm::Service rng; + if (!rng.isAvailable()){ + throw cms::Exception("Configuration") + << "EcalBarrelClusterFastTimer::EcalBarrelClusterFastTimer() - RandomNumberGeneratorService is not present in configuration file.\n" + << "Add the service in the configuration file or remove the modules that require it."; + } +} + +void EcalBarrelClusterFastTimer::produce(edm::StreamID sid, edm::Event& evt, const edm::EventSetup& es) const { + // get RNG engine + edm::Service rng; + auto rng_engine = &(rng->getEngine(sid)); + + edm::Handle > clustersH; + edm::Handle timehitsH; + edm::Handle verticesH; + + evt.getByToken(ebClustersToken_,clustersH); + evt.getByToken(ebTimeHitsToken_,timehitsH); + evt.getByToken(timedVertexToken_,verticesH); + + edm::ESHandle geoHandle; + es.get().get(geoHandle); + const auto* barrelGeom = + geoHandle->getSubdetectorGeometry(DetId::Ecal,EcalBarrel); + + const auto& clusters = *clustersH; + const auto& timehits = *timehitsH; + const auto& vertices = *verticesH; + + std::vector > times; // perfect times keyed to cluster index + times.reserve(clusters.size()); + + for( const auto& cluster : clusters ) { + times.push_back( getTimeForECALPFCluster( timehits, cluster ) ); + } + + for( const auto& reso : _resolutions ) { + const std::string& name = reso->name(); + std::vector resolutions; + std::vector > smeared_times; + resolutions.reserve(clusters.size()); + smeared_times.reserve(clusters.size()); + std::auto_ptr > > outP( new std::vector< std::vector > ); + auto& out = *outP; + + // smear once then correct to multiple vertices + for( unsigned i = 0 ; i < clusters.size(); ++i ) { + const float theresolution = reso->getTimeResolution(clusters[i]); + + smeared_times.emplace_back( CLHEP::RandGauss::shoot(rng_engine, times[i].first, theresolution), times[i].second ); + resolutions.push_back( theresolution ); + } + + // now loop over vertices and back-out the correction to (0,0,0) + out.resize(vertices.size()); + for( unsigned i = 0; i < vertices.size(); ++i ) { + const reco::Vertex& vtx = vertices[i]; + out[i].reserve(clusters.size()); + for( unsigned j = 0; j < clusters.size(); ++j ) { + const auto& basetime = smeared_times[j]; + const float corrtime = correctTimeToVertex(basetime.first,basetime.second,vtx,barrelGeom); + out[i].push_back(corrtime); + } + } + + evt.put(outP,name); + writeValueMap(evt,clustersH,resolutions,name+resolution); + } + +} + +std::pair EcalBarrelClusterFastTimer::getTimeForECALPFCluster(const EcalRecHitCollection& timehits, const reco::PFCluster& cluster) const { + + const auto& rhfs = cluster.recHitFractions(); + unsigned best_hit = 0; + float best_energy = -1.f; + for( const auto& rhf : rhfs ) { + const auto& hitref = rhf.recHitRef(); + const unsigned detid = hitref->detId(); + const auto fraction = rhf.fraction(); + const auto energy = hitref->energy(); + if( fraction < minFraction_ || energy < minEnergy_ ) continue; + auto timehit = timehits.find(detid); + float e_hit = rhf.fraction() * hitref->energy(); + if( e_hit > best_energy && timehit->isTimeValid() ) { + best_energy = e_hit; + best_hit = detid; + } + } + + float best_time_guess = -1.0; + if( best_energy > 0. ) { + best_time_guess = timehits.find(best_hit)->time(); + } + + //std::cout << "EcalBarrelFastTimer: " << best_time_guess << ' ' << best_energy << ' ' << best_hit << std::endl; + + return std::make_pair(best_time_guess,DetId(best_hit)); +} + +float EcalBarrelClusterFastTimer::correctTimeToVertex(const float intime, const DetId& timeDet, const reco::Vertex& vtx, + const CaloSubdetectorGeometry* ecalGeom) const { + if( timeDet.rawId() == 0 ) return -999.; + // correct the cluster time from 0,0,0 to the primary vertex given + const CaloCellGeometry* cellGeometry ( ecalGeom->getGeometry( timeDet ) ) ; + if( nullptr == cellGeometry ) { + throw cms::Exception("BadECALBarrelCell") + << timeDet << " is not a valid ECAL Barrel DetId!"; + } + //depth in mm in the middle of the layer position; + GlobalPoint layerPos = (dynamic_cast(cellGeometry))->getPosition( ecalDepth_+0.5 ); + const math::XYZPoint layerPos_cm( layerPos.x(), layerPos.y(), layerPos.z() ); + const math::XYZVector to_center = layerPos_cm - math::XYZPoint(0.,0.,0.); + const math::XYZVector to_vtx = layerPos_cm - vtx.position(); + + /* + std::cout << intime << ' ' << to_center.r()/cm_per_ns << ' ' << to_vtx.r()/cm_per_ns + << ' ' << intime + to_center.r()/cm_per_ns - to_vtx.r()/cm_per_ns << std::endl; + */ + + return intime + to_center.r()/cm_per_ns - to_vtx.r()/cm_per_ns; +} diff --git a/RecoParticleFlow/PFSimProducer/plugins/PerfectResolutionModel.cc b/RecoParticleFlow/PFSimProducer/plugins/PerfectResolutionModel.cc new file mode 100644 index 0000000000000..26829f734c91b --- /dev/null +++ b/RecoParticleFlow/PFSimProducer/plugins/PerfectResolutionModel.cc @@ -0,0 +1,14 @@ +#include "ResolutionModel.h" + +class PerfectResolutionModel : public ResolutionModel { +public: + PerfectResolutionModel( const edm::ParameterSet& conf ) : ResolutionModel( conf ) {} + + virtual float getTimeResolution(const reco::Track&) const override { return 1e-6; } + virtual float getTimeResolution(const reco::PFCluster&) const override { return 1e-6; } + +}; + +DEFINE_EDM_PLUGIN(ResolutionModelFactory, + PerfectResolutionModel, + "PerfectResolutionModel"); diff --git a/RecoParticleFlow/PFSimProducer/plugins/ResolutionModel.cc b/RecoParticleFlow/PFSimProducer/plugins/ResolutionModel.cc new file mode 100644 index 0000000000000..7b47fd80b0de7 --- /dev/null +++ b/RecoParticleFlow/PFSimProducer/plugins/ResolutionModel.cc @@ -0,0 +1,4 @@ +#include "ResolutionModel.h" + +EDM_REGISTER_PLUGINFACTORY(ResolutionModelFactory, + "ResolutionModelFactory"); diff --git a/RecoParticleFlow/PFSimProducer/plugins/ResolutionModel.h b/RecoParticleFlow/PFSimProducer/plugins/ResolutionModel.h new file mode 100644 index 0000000000000..101507a844f68 --- /dev/null +++ b/RecoParticleFlow/PFSimProducer/plugins/ResolutionModel.h @@ -0,0 +1,33 @@ +#ifndef __RecoFTL_FastTimingKludge_ResolutionModel_h__ +#define __RecoFTL_FastTimingKludge_ResolutionModel_h__ + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/ParticleFlowReco/interface/PFCluster.h" + +#include +#include + +class ResolutionModel { + public: + ResolutionModel(const edm::ParameterSet& conf): + _modelName(conf.getParameter("modelName")) { + } + virtual ~ResolutionModel() { } + // get rid of things we should never use... + ResolutionModel(const ResolutionModel&) = delete; + ResolutionModel& operator=(const ResolutionModel&) = delete; + + virtual float getTimeResolution(const reco::Track&) const { return -1.f; } + virtual float getTimeResolution(const reco::PFCluster&) const { return -1.f; } + + const std::string& name() const { return _modelName; } + + private: + const std::string _modelName; +}; + +#include "FWCore/PluginManager/interface/PluginFactory.h" +typedef edmplugin::PluginFactory< ResolutionModel* (const edm::ParameterSet&) > ResolutionModelFactory; + +#endif diff --git a/RecoParticleFlow/PFSimProducer/python/ecalBarrelClusterFastTimer_cfi.py b/RecoParticleFlow/PFSimProducer/python/ecalBarrelClusterFastTimer_cfi.py new file mode 100644 index 0000000000000..fe4f8d09bf791 --- /dev/null +++ b/RecoParticleFlow/PFSimProducer/python/ecalBarrelClusterFastTimer_cfi.py @@ -0,0 +1,12 @@ +import FWCore.ParameterSet.Config as cms + +ecalBarrelClusterFastTimer = cms.EDProducer( + 'EcalBarrelClusterFastTimer', + ebTimeHits = cms.InputTag('ecalDetailedTimeRecHit:EcalRecHitsEB'), + ebClusters = cms.InputTag('particleFlowClusterECAL'), + timedVertices = cms.InputTag('offlinePrimaryVertices4D'), + minFractionToConsider = cms.double(0.1), + minEnergyToConsider = cms.double(0.0), + ecalDepth = cms.double(7.0), + resolutionModels = cms.VPSet( cms.PSet( modelName = cms.string('PerfectResolutionModel') ) ) + ) diff --git a/SLHCUpgradeSimulations/Configuration/python/combinedCustoms.py b/SLHCUpgradeSimulations/Configuration/python/combinedCustoms.py index d8d49532bafd8..081ce61d1b2b5 100644 --- a/SLHCUpgradeSimulations/Configuration/python/combinedCustoms.py +++ b/SLHCUpgradeSimulations/Configuration/python/combinedCustoms.py @@ -3,7 +3,6 @@ from SLHCUpgradeSimulations.Configuration.customise_mixing import customise_NoCrossing from SLHCUpgradeSimulations.Configuration.phase1TkCustoms import customise as customisePhase1Tk from SLHCUpgradeSimulations.Configuration.HCalCustoms import customise_HcalPhase1, customise_HcalPhase0 - from SLHCUpgradeSimulations.Configuration.fixMissingUpgradeGTPayloads import fixRPCConditions from SLHCUpgradeSimulations.Configuration.phase2TkTilted import customise as customiseTiltedTK @@ -21,8 +20,6 @@ def cust_2023flat(process): process=customiseFlatTK(process) return process - -######Below are the customized used for SLHC releases def cust_2019(process): process=customisePostLS1(process,displayDeprecationWarning=False) process=customisePhase1Tk(process) diff --git a/SimCalorimetry/Configuration/python/SimCalorimetry_EventContent_cff.py b/SimCalorimetry/Configuration/python/SimCalorimetry_EventContent_cff.py index 36c64be1a771c..58ef42496f8ce 100644 --- a/SimCalorimetry/Configuration/python/SimCalorimetry_EventContent_cff.py +++ b/SimCalorimetry/Configuration/python/SimCalorimetry_EventContent_cff.py @@ -28,3 +28,9 @@ # eras.run2_common.toModify( SimCalorimetryFEVTDEBUG.outputCommands, func=lambda outputCommands: outputCommands.append('keep *_simHcalUnsuppressedDigis_*_*') ) eras.run2_common.toModify( SimCalorimetryRAW.outputCommands, func=lambda outputCommands: outputCommands.append('keep *_simHcalUnsuppressedDigis_*_*') ) + +eras.phase2_timing.toModify(SimCalorimetryFEVTDEBUG.outputCommands, func=lambda outputCommands: outputCommands.append('keep *_mix_EETimeDigi_*') ) +eras.phase2_timing.toModify(SimCalorimetryFEVTDEBUG.outputCommands, func=lambda outputCommands: outputCommands.append('keep *_mix_EBTimeDigi_*') ) + +eras.phase2_timing.toModify(SimCalorimetryRAW.outputCommands, func=lambda outputCommands: outputCommands.append('keep *_mix_EETimeDigi_*') ) +eras.phase2_timing.toModify(SimCalorimetryRAW.outputCommands, func=lambda outputCommands: outputCommands.append('keep *_mix_EBTimeDigi_*') ) diff --git a/SimCalorimetry/EcalSimAlgos/interface/EcalTimeMapDigitizer.h b/SimCalorimetry/EcalSimAlgos/interface/EcalTimeMapDigitizer.h new file mode 100644 index 0000000000000..add1b2eedf13d --- /dev/null +++ b/SimCalorimetry/EcalSimAlgos/interface/EcalTimeMapDigitizer.h @@ -0,0 +1,159 @@ +#ifndef EcalSimAlgos_EcalTimeMapDigitizer_h +#define EcalSimAlgos_EcalTimeMapDigitizer_h + +/** Turns hits into digis. Assumes that + there's an ElectroncsSim class with the + interface analogToDigital(const CaloSamples &, Digi &); +*/ + +#include "CalibFormats/CaloObjects/interface/CaloTSamples.h" +#include "SimDataFormats/CrossingFrame/interface/MixCollection.h" +#include "SimDataFormats/CaloHit/interface/PCaloHit.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "CalibFormats/CaloObjects/interface/CaloTSamplesBase.h" +#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h" + +class CaloSubdetectorGeometry ; + +class EcalTimeMapDigitizer +{ + public: + + struct time_average{ + + static const unsigned short time_average_capacity=10; + + DetId id; + float average_time[time_average_capacity]; + unsigned int nhits[time_average_capacity]; + float tot_energy[time_average_capacity]; + + time_average(const DetId& myId): + id(myId) + { + for (unsigned int i(0);i0) + average_time[i]=average_time[i]/tot_energy[i]; + else + average_time[i]=0; + } + }; + + void + setZero() + { + for (unsigned int i(0);i0) + return false; + } + return true; + }; + + }; + + typedef time_average TimeSamples; + + typedef EcalTimeDigi Digi; + + typedef std::vector< unsigned int > VecInd ; + + explicit EcalTimeMapDigitizer(EcalSubdetector myDet); + + virtual ~EcalTimeMapDigitizer(); + + void add(const std::vector & hits, int bunchCrossing); + + void setGeometry( const CaloSubdetectorGeometry* geometry ) ; + + void initializeMap(); + + void run(EcalTimeDigiCollection& output ); + + void finalizeHits(); + + inline void setTimeLayerId( const int& layerId ) { m_timeLayerId = layerId; }; + + int getTimeLayerId() { return m_timeLayerId; }; + + EcalSubdetector subdetector() { return m_subDet; }; + + int minBunch() const ; + + int maxBunch() const ; + + void blankOutUsedSamples() ; // blank out previously used elements + + +/* const CaloVHitFilter* hitFilter() const ; */ + + private: + + TimeSamples* findSignal( const DetId& detId ); + + VecInd& index() ; + + const VecInd& index() const ; + + unsigned int samplesSize() const; + + unsigned int samplesSizeAll() const; + + const TimeSamples* operator[]( unsigned int i ) const; + + TimeSamples* operator[]( unsigned int i ); + + TimeSamples* vSam( unsigned int i ); + + TimeSamples* vSamAll( unsigned int i ); + + const TimeSamples* vSamAll( unsigned int i ) const; + + EcalSubdetector m_subDet; + //time difference between bunches + + static const int BUNCHSPACE=25; + + static const float MIN_ENERGY_THRESHOLD; //50 KeV threshold to consider a valid hit in the timing detector + + static const int m_minBunch=-4; + static const int m_maxBunch=5; + + int m_timeLayerId; + + const CaloSubdetectorGeometry* m_geometry ; + + double timeOfFlight( const DetId& detId , int layer) const; + + std::vector m_vSam; + + VecInd m_index; + +}; + +#endif diff --git a/SimCalorimetry/EcalSimAlgos/src/EBHitResponse.cc b/SimCalorimetry/EcalSimAlgos/src/EBHitResponse.cc index 135720c898717..19050ea891482 100644 --- a/SimCalorimetry/EcalSimAlgos/src/EBHitResponse.cc +++ b/SimCalorimetry/EcalSimAlgos/src/EBHitResponse.cc @@ -124,6 +124,7 @@ EBHitResponse::putAPDSignal( const DetId& detId , double EBHitResponse::apdSignalAmplitude( const PCaloHit& hit, CLHEP::HepRandomEngine* engine ) const { + // std::cout << "*** " << hit.depth() << std::endl; assert( 1 == hit.depth() || 2 == hit.depth() ) ; @@ -228,7 +229,7 @@ void EBHitResponse::add( const PCaloHit& hit, CLHEP::HepRandomEngine* engine ) { if (!edm::isNotFinite( hit.time() ) && ( 0 == hitFilter() || hitFilter()->accepts( hit ) ) ) { - if( 0 == hit.depth() ) // for now take only nonAPD hits + if( 0 == hit.depth() || hit.depth() >=100 ) // for now take only nonAPD hits { if( !m_apdOnly ) putAnalogSignal( hit, engine ) ; } @@ -269,40 +270,40 @@ EBHitResponse::run( MixCollection& hits, CLHEP::HepRandomEngine* engin ( 0 == hitFilter() || hitFilter()->accepts( hit ) ) ) { - if( 0 == hit.depth() ) // for now take only nonAPD hits - { - if( !m_apdOnly ) putAnalogSignal( hit, engine ) ; - } - else // APD hits here - { + if( 0 == hit.depth() || hit.depth() >=100 ) // for now take only nonAPD hits + { + if( !m_apdOnly ) putAnalogSignal( hit, engine ) ; + } + else // APD hits here + { if( apdParameters()->addToBarrel() || m_apdOnly ) - { - const unsigned int icell ( EBDetId( hit.id() ).denseIndex() ) ; - m_apdNpeVec[ icell ] += apdSignalAmplitude( hit, engine ) ; - if( 0 == m_apdTimeVec[ icell ] ) m_apdTimeVec[ icell ] = hit.time() ; - } - } + { + const unsigned int icell ( EBDetId( hit.id() ).denseIndex() ) ; + m_apdNpeVec[ icell ] += apdSignalAmplitude( hit, engine ) ; + if( 0 == m_apdTimeVec[ icell ] ) m_apdTimeVec[ icell ] = hit.time() ; + } + } } } - + if( apdParameters()->addToBarrel() || m_apdOnly ) - { - for( unsigned int i ( 0 ) ; i != bSize ; ++i ) - { - if( 0 < m_apdNpeVec[i] ) + { + for( unsigned int i ( 0 ) ; i != bSize ; ++i ) { - putAPDSignal( EBDetId::detIdFromDenseIndex( i ), - m_apdNpeVec[i] , - m_apdTimeVec[i] ) ; - - // now zero out for next time - m_apdNpeVec[i] = 0. ; - m_apdTimeVec[i] = 0. ; + if( 0 < m_apdNpeVec[i] ) + { + putAPDSignal( EBDetId::detIdFromDenseIndex( i ), + m_apdNpeVec[i] , + m_apdTimeVec[i] ) ; + + // now zero out for next time + m_apdNpeVec[i] = 0. ; + m_apdTimeVec[i] = 0. ; + } } - } - } + } } unsigned int diff --git a/SimCalorimetry/EcalSimAlgos/src/EcalTimeMapDigitizer.cc b/SimCalorimetry/EcalSimAlgos/src/EcalTimeMapDigitizer.cc new file mode 100644 index 0000000000000..c939b402e9754 --- /dev/null +++ b/SimCalorimetry/EcalSimAlgos/src/EcalTimeMapDigitizer.cc @@ -0,0 +1,342 @@ +#include "SimCalorimetry/EcalSimAlgos/interface/EcalTimeMapDigitizer.h" + +#include "Geometry/CaloGeometry/interface/CaloGenericDetId.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h" +#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h" + +//#include "FWCore/ServiceRegistry/interface/Service.h" +//#include "FWCore/Utilities/interface/RandomNumberGenerator.h" +//#include "CLHEP/Random/RandPoissonQ.h" +//#include "CLHEP/Random/RandGaussQ.h" + +#include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/Utilities/interface/isFinite.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "CLHEP/Units/GlobalPhysicalConstants.h" +#include "CLHEP/Units/GlobalSystemOfUnits.h" + +#include + +//#define ecal_time_debug 1 + +const float EcalTimeMapDigitizer::MIN_ENERGY_THRESHOLD=5e-5; //50 KeV threshold to consider a valid hit in the timing detector + +EcalTimeMapDigitizer::EcalTimeMapDigitizer(EcalSubdetector myDet): + m_subDet(myDet), + m_geometry(0) +{ +// edm::Service rng ; +// if ( !rng.isAvailable() ) +// { +// throw cms::Exception("Configuration") +// << "EcalTimeMapDigitizer requires the RandomNumberGeneratorService\n" +// "which is not present in the configuration file. You must add the service\n" +// "in the configuration file or remove the modules that require it."; +// } +// m_RandPoisson = new CLHEP::RandPoissonQ( rng->getEngine() ) ; +// m_RandGauss = new CLHEP::RandGaussQ( rng->getEngine() ) ; + + unsigned int size=0; + DetId detId(0); + + //Initialize the map + if (myDet==EcalBarrel) + { + size=EBDetId::kSizeForDenseIndexing; + detId=EBDetId::detIdFromDenseIndex( 0 ) ; + } + else if (myDet==EcalEndcap) + { + size=EEDetId::kSizeForDenseIndexing; + detId=EEDetId::detIdFromDenseIndex( 0 ) ; + } + else + edm::LogError("TimeDigiError") << "[EcalTimeMapDigitizer]::ERROR::This subdetector " << myDet << " is not implemented"; + + + assert(m_maxBunch-m_minBunch+1<=10); + assert(m_minBunch<=0); + + m_vSam.reserve( size ) ; + + for( unsigned int i ( 0 ) ; i != size ; ++i ) + { + // m_vSam.emplace_back(CaloGenericDetId( detId.det(), detId.subdetId(), i ) , + // m_maxBunch-m_minBunch+1, abs(m_minBunch) ); + m_vSam.emplace_back(TimeSamples(CaloGenericDetId( detId.det(), detId.subdetId(), i ))); + } + + + edm::LogInfo("TimeDigiInfo") << "[EcalTimeDigitizer]::Subdetector " << m_subDet << "::Reserved size for time digis " << m_vSam.size(); + +#ifdef ecal_time_debug + std::cout << "[EcalTimeDigitizer]::Subdetector " << m_subDet << "::Reserved size for time digis " << m_vSam.size() << std::endl; +#endif + +} + +EcalTimeMapDigitizer::~EcalTimeMapDigitizer() +{ +} + +void +EcalTimeMapDigitizer::add(const std::vector & hits, int bunchCrossing) { + if(bunchCrossing>=m_minBunch && bunchCrossing<=m_maxBunch ) { + + for(std::vector::const_iterator it = hits.begin(), itEnd = hits.end(); it != itEnd; ++it) { + //here goes the map logic + + if (edm::isNotFinite( (*it).time() ) ) + continue; + + //Just consider only the hits belonging to the specified time layer + if ((*it).depth()!=100+m_timeLayerId) //modified layerId start from 100 + continue; + + if ((*it).energy()zero() ) m_index.push_back( di ) ; + return result ; +} + + + +void +EcalTimeMapDigitizer::setGeometry( const CaloSubdetectorGeometry* geometry ) +{ + m_geometry = geometry ; +} + + +void +EcalTimeMapDigitizer::blankOutUsedSamples() // blank out previously used elements +{ + const unsigned int size ( m_index.size() ) ; + + + for( unsigned int i ( 0 ) ; i != size ; ++i ) + { +#ifdef ecal_time_debug + std::cout << "Zeroing " << m_index[i] << std::endl; +#endif + vSamAll( m_index[i] )->setZero() ; + } + + m_index.erase( m_index.begin() , // done and make ready to start over + m_index.end() ) ; +} + + +void +EcalTimeMapDigitizer::finalizeHits() +{ + //Getting the size from the actual hits + const unsigned int size ( m_index.size() ) ; + + //Here just averaging the MC truth + for( unsigned int i ( 0 ) ; i != size ; ++i ) + { +#ifdef ecal_time_debug + std::cout << "Averaging " << m_index[i] << std::endl; +#endif + vSamAll( m_index[i] )->calculateAverage() ; +#ifdef ecal_time_debug + for ( unsigned int j ( 0 ) ; j != vSamAll( m_index[i] )->time_average_capacity; ++j ) + std::cout << j << "\t" << vSamAll( m_index[i] )->average_time[j] << "\t" << vSamAll( m_index[i] )->nhits[j] << "\t" << vSamAll( m_index[i] )->tot_energy[j] << std::endl; +#endif + } + + //Noise can be added here + +} + +void +EcalTimeMapDigitizer::initializeMap() +{ +#ifdef ecal_time_debug + std::cout << "[EcalTimeMapDigitizer]::Zeroing the used samples" << std::endl; +#endif + blankOutUsedSamples(); +} + + +void +EcalTimeMapDigitizer::run( EcalTimeDigiCollection& output ) +{ +#ifdef ecal_time_debug + std::cout << "[EcalTimeMapDigitizer]::Finalizing hits and fill output collection" << std::endl; +#endif + + //Do the time averages and add noise if simulated + finalizeHits(); + + //Until we do now simulated noise we can get the size from the index vector + const unsigned int ssize ( m_index.size() ) ; + + output.reserve( ssize ) ; + + for( unsigned int i ( 0 ) ; i != ssize ; ++i ) + { + + #ifdef ecal_time_debug + std::cout << "----- in digi loop " << i << std::endl; + #endif + + output.push_back( Digi(vSamAll( m_index[i] )->id) ); + + unsigned int nTimeHits=0; + float timeHits[vSamAll( m_index[i] )->time_average_capacity]; + unsigned int timeBX[vSamAll( m_index[i] )->time_average_capacity]; + + for ( unsigned int j ( 0 ) ; j != vSamAll( m_index[i] )->time_average_capacity; ++j ) //here sampling on the OOTPU + { + if (vSamAll( m_index[i] )->nhits[j]>0) + { + timeHits[nTimeHits]= vSamAll( m_index[i] )->average_time[j]; + timeBX[nTimeHits++]= m_minBunch+j; + } + } + + output.back().setSize(nTimeHits); + + for (unsigned int j ( 0 ) ; j != nTimeHits; ++j ) //filling the !zero hits + { + output.back().setSample(j,timeHits[j]); + if (timeBX[j]==0) + { +#ifdef ecal_time_debug + std::cout << "setting interesting sample " << j << std::endl; +#endif + output.back().setSampleOfInterest(j); + } + } + + #ifdef ecal_time_debug + std::cout << "digi " << output.back().id().rawId() << "\t" << output.back().size(); + if (output.back().sampleOfInterest()>0) + std::cout << "\tBX0 time " << output.back().sample(output.back().sampleOfInterest()) << std::endl; + else + std::cout << "\tNo in time hits" << std::endl; + #endif + } +#ifdef ecal_time_debug + std::cout << "[EcalTimeMapDigitizer]::Output collection size " << output.size() << std::endl; +#endif + +} + + +double +EcalTimeMapDigitizer::timeOfFlight( const DetId& detId , int layer) const +{ + //not using the layer yet + const CaloCellGeometry* cellGeometry ( m_geometry->getGeometry( detId ) ) ; + assert( 0 != cellGeometry ) ; + GlobalPoint layerPos = (dynamic_cast(cellGeometry))->getPosition( double(layer)+0.5 ); //depth in mm in the middle of the layer position + return layerPos.mag()*cm/c_light ; +} + + +unsigned int +EcalTimeMapDigitizer::samplesSize() const +{ + return m_vSam.size() ; +} + +unsigned int +EcalTimeMapDigitizer::samplesSizeAll() const +{ + return m_vSam.size() ; +} + +const EcalTimeMapDigitizer::TimeSamples* +EcalTimeMapDigitizer::operator[]( unsigned int i ) const +{ + return &m_vSam[ i ] ; +} + +EcalTimeMapDigitizer::TimeSamples* +EcalTimeMapDigitizer::operator[]( unsigned int i ) +{ + return &m_vSam[ i ] ; +} + +EcalTimeMapDigitizer::TimeSamples* +EcalTimeMapDigitizer::vSam( unsigned int i ) +{ + return &m_vSam[ i ] ; +} + +EcalTimeMapDigitizer::TimeSamples* +EcalTimeMapDigitizer::vSamAll( unsigned int i ) +{ + return &m_vSam[ i ] ; +} + +const EcalTimeMapDigitizer::TimeSamples* +EcalTimeMapDigitizer::vSamAll( unsigned int i ) const +{ + return &m_vSam[ i ] ; +} + +int +EcalTimeMapDigitizer::minBunch() const +{ + return m_minBunch ; +} + +int +EcalTimeMapDigitizer::maxBunch() const +{ + return m_maxBunch ; +} + +EcalTimeMapDigitizer::VecInd& +EcalTimeMapDigitizer::index() +{ + return m_index ; +} + +const EcalTimeMapDigitizer::VecInd& +EcalTimeMapDigitizer::index() const +{ + return m_index ; +} + +// const EcalTimeMapDigitizer::TimeSamples* +// EcalTimeMapDigitizer::findDetId( const DetId& detId ) const +// { +// const unsigned int di ( CaloGenericDetId( detId ).denseIndex() ) ; +// return vSamAll( di ) ; +// } diff --git a/SimCalorimetry/EcalSimProducers/interface/EcalTimeDigiProducer.h b/SimCalorimetry/EcalSimProducers/interface/EcalTimeDigiProducer.h new file mode 100644 index 0000000000000..422cb8073e272 --- /dev/null +++ b/SimCalorimetry/EcalSimProducers/interface/EcalTimeDigiProducer.h @@ -0,0 +1,68 @@ +#ifndef SimCalorimetry_EcalSimProducers_EcalTimeDigiProducer_h +#define SimCalorimetry_EcalSimProducers_EcalTimeDigiProducer_h + + +#include "DataFormats/Math/interface/Error.h" + +#include "SimCalorimetry/EcalSimAlgos/interface/EcalDigitizerTraits.h" +#include "SimGeneral/MixingModule/interface/DigiAccumulatorMixMod.h" +#include "SimDataFormats/CaloHit/interface/PCaloHit.h" +#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h" +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" + +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include + +class ESDigitizer ; + +class CaloGeometry ; +class EcalSimParameterMap ; +class PileUpEventPrincipal ; +class EcalTimeMapDigitizer; + +namespace edm { + class Event; + class EventSetup; + template class Handle; + class ParameterSet; +} + +class EcalTimeDigiProducer : public DigiAccumulatorMixMod { + public: + + EcalTimeDigiProducer( const edm::ParameterSet& params , edm::stream::EDProducerBase& mixMod, edm::ConsumesCollector&); + virtual ~EcalTimeDigiProducer(); + + virtual void initializeEvent(edm::Event const& e, edm::EventSetup const& c); + virtual void accumulate(edm::Event const& e, edm::EventSetup const& c); + virtual void accumulate(PileUpEventPrincipal const& e, edm::EventSetup const& c, edm::StreamID const&); + virtual void finalizeEvent(edm::Event& e, edm::EventSetup const& c); + + private: + + typedef edm::Handle > HitsHandle; + void accumulateCaloHits(HitsHandle const& ebHandle, HitsHandle const& eeHandle, int bunchCrossing); + + void checkGeometry(const edm::EventSetup& eventSetup) ; + + void updateGeometry() ; + + const std::string m_EBdigiCollection ; + const std::string m_EEdigiCollection ; + const edm::InputTag m_hitsProducerTagEB ; + const edm::InputTag m_hitsProducerTagEE ; + const edm::EDGetTokenT > m_hitsProducerTokenEB ; + const edm::EDGetTokenT > m_hitsProducerTokenEE ; + + private: + int m_timeLayerEB; + int m_timeLayerEE; + const CaloGeometry* m_Geometry ; + + EcalTimeMapDigitizer* m_BarrelDigitizer; + EcalTimeMapDigitizer* m_EndcapDigitizer; +}; + +#endif diff --git a/SimCalorimetry/EcalSimProducers/plugins/SealModule.cc b/SimCalorimetry/EcalSimProducers/plugins/SealModule.cc index 4b7641c709e74..17db6a577e626 100644 --- a/SimCalorimetry/EcalSimProducers/plugins/SealModule.cc +++ b/SimCalorimetry/EcalSimProducers/plugins/SealModule.cc @@ -1,6 +1,8 @@ #include "FWCore/Framework/interface/MakerMacros.h" #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixModFactory.h" #include "SimCalorimetry/EcalSimProducers/interface/EcalDigiProducer.h" +#include "SimCalorimetry/EcalSimProducers/interface/EcalTimeDigiProducer.h" DEFINE_DIGI_ACCUMULATOR(EcalDigiProducer); +DEFINE_DIGI_ACCUMULATOR(EcalTimeDigiProducer); diff --git a/SimCalorimetry/EcalSimProducers/python/ecalTimeDigiParameters_cff.py b/SimCalorimetry/EcalSimProducers/python/ecalTimeDigiParameters_cff.py new file mode 100644 index 0000000000000..538128ce0096b --- /dev/null +++ b/SimCalorimetry/EcalSimProducers/python/ecalTimeDigiParameters_cff.py @@ -0,0 +1,10 @@ +import FWCore.ParameterSet.Config as cms + +ecal_time_digi_parameters = cms.PSet( + hitsProducerEB = cms.InputTag('g4SimHits:EcalHitsEB'), + hitsProducerEE = cms.InputTag('g4SimHits:EcalHitsEE'), + EBtimeDigiCollection = cms.string('EBTimeDigi'), + EEtimeDigiCollection = cms.string('EETimeDigi'), + timeLayerBarrel = cms.int32(7), + timeLayerEndcap = cms.int32(3) + ) diff --git a/SimCalorimetry/EcalSimProducers/src/EcalTimeDigiProducer.cc b/SimCalorimetry/EcalSimProducers/src/EcalTimeDigiProducer.cc new file mode 100644 index 0000000000000..43338e8fdabd5 --- /dev/null +++ b/SimCalorimetry/EcalSimProducers/src/EcalTimeDigiProducer.cc @@ -0,0 +1,163 @@ +#include "FWCore/Framework/interface/Event.h" +#include "SimCalorimetry/EcalSimProducers/interface/EcalTimeDigiProducer.h" +#include "SimCalorimetry/EcalSimAlgos/interface/EcalTimeMapDigitizer.h" + +#include "CalibFormats/CaloObjects/interface/CaloSamples.h" + +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h" + +#include "FWCore/Framework/interface/ESHandle.h" +#include "DataFormats/Common/interface/Handle.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h" + + +//#define ecal_time_debug 1 + +EcalTimeDigiProducer::EcalTimeDigiProducer( const edm::ParameterSet& params, edm::stream::EDProducerBase& mixMod,edm::ConsumesCollector& sumes) : + DigiAccumulatorMixMod(), + m_EBdigiCollection ( params.getParameter("EBtimeDigiCollection") ) , + m_EEdigiCollection ( params.getParameter("EEtimeDigiCollection") ) , + m_hitsProducerTagEB ( params.getParameter("hitsProducerEB" ) ) , + m_hitsProducerTagEE ( params.getParameter("hitsProducerEE" ) ) , + m_hitsProducerTokenEB ( sumes.consumes >( m_hitsProducerTagEB) ) , + m_hitsProducerTokenEE ( sumes.consumes >( m_hitsProducerTagEE) ) , + m_timeLayerEB ( params.getParameter ("timeLayerBarrel") ), + m_timeLayerEE ( params.getParameter ("timeLayerEndcap") ), + m_Geometry ( 0 ) +{ + mixMod.produces(m_EBdigiCollection); + mixMod.produces(m_EEdigiCollection); + + m_BarrelDigitizer = new EcalTimeMapDigitizer(EcalBarrel); + m_EndcapDigitizer = new EcalTimeMapDigitizer(EcalEndcap); + +#ifdef ecal_time_debug + std::cout << "[EcalTimeDigiProducer]::Create EB " << m_EBdigiCollection << " and EE " << m_EEdigiCollection << " collections and digitizers" << std::endl; +#endif + + m_BarrelDigitizer->setTimeLayerId(m_timeLayerEB); + m_EndcapDigitizer->setTimeLayerId(m_timeLayerEE); +} + +EcalTimeDigiProducer::~EcalTimeDigiProducer() +{ +} + +void +EcalTimeDigiProducer::initializeEvent(edm::Event const& event, edm::EventSetup const& eventSetup) { + checkGeometry( eventSetup ); + // checkCalibrations( event, eventSetup ); + // here the methods to clean the maps + m_BarrelDigitizer->initializeMap(); + m_EndcapDigitizer->initializeMap(); +} + +void +EcalTimeDigiProducer::accumulateCaloHits(HitsHandle const& ebHandle, HitsHandle const& eeHandle, int bunchCrossing) { + // accumulate the simHits and do the averages in a given layer per bunch crossing + if(ebHandle.isValid()) { + m_BarrelDigitizer->add(*ebHandle.product(), bunchCrossing); + } + + if(eeHandle.isValid()) { + m_EndcapDigitizer->add(*eeHandle.product(), bunchCrossing); + } + +} + +void +EcalTimeDigiProducer::accumulate(edm::Event const& e, edm::EventSetup const& eventSetup) { + // Step A: Get Inputs + edm::Handle > ebHandle; + e.getByToken(m_hitsProducerTokenEB, ebHandle); + + edm::Handle > eeHandle; + e.getByToken(m_hitsProducerTokenEE, eeHandle); + +#ifdef ecal_time_debug + std::cout << "[EcalTimeDigiProducer]::Accumulate Hits HS event" << std::endl; +#endif + + accumulateCaloHits(ebHandle, eeHandle, 0); +} + +void +EcalTimeDigiProducer::accumulate(PileUpEventPrincipal const& e, edm::EventSetup const& eventSetup, edm::StreamID const&) { + edm::Handle > ebHandle; + e.getByLabel(m_hitsProducerTagEB, ebHandle); + + edm::Handle > eeHandle; + e.getByLabel(m_hitsProducerTagEE, eeHandle); + +#ifdef ecal_time_debug + std::cout << "[EcalTimeDigiProducer]::Accumulate Hits for BC " << e.bunchCrossing() << std::endl; +#endif + accumulateCaloHits(ebHandle, eeHandle, e.bunchCrossing()); +} + +void +EcalTimeDigiProducer::finalizeEvent(edm::Event& event, edm::EventSetup const& eventSetup) { + std::unique_ptr barrelResult ( new EcalTimeDigiCollection() ) ; + std::unique_ptr endcapResult ( new EcalTimeDigiCollection() ) ; + +#ifdef ecal_time_debug + std::cout << "[EcalTimeDigiProducer]::finalizeEvent" << std::endl; +#endif + + // here basically just put everything in the final collections + m_BarrelDigitizer->run( *barrelResult ) ; + +#ifdef ecal_time_debug + std::cout << "[EcalTimeDigiProducer]::EB Digi size " << barrelResult->size() << std::endl; +#endif + + edm::LogInfo("TimeDigiInfo") << "EB time Digis: " << barrelResult->size() ; + + m_EndcapDigitizer->run( *endcapResult ) ; + +#ifdef ecal_time_debug + std::cout << "[EcalTimeDigiProducer]::EE Digi size " << endcapResult->size() << std::endl; +#endif + + edm::LogInfo("TimeDigiInfo") << "EE Digis: " << endcapResult->size() ; + +#ifdef ecal_time_debug + std::cout << "[EcalTimeDigiProducer]::putting collections into the event " << std::endl; +#endif + + event.put( std::move(barrelResult), m_EBdigiCollection ) ; + event.put( std::move(endcapResult), m_EEdigiCollection ) ; +} + +void +EcalTimeDigiProducer::checkGeometry( const edm::EventSetup & eventSetup ) +{ + // TODO find a way to avoid doing this every event + edm::ESHandle hGeometry ; + eventSetup.get().get( hGeometry ) ; + + const CaloGeometry* pGeometry = &*hGeometry; + + if( pGeometry != m_Geometry ) + { + m_Geometry = pGeometry; + updateGeometry(); + } +} + +void +EcalTimeDigiProducer::updateGeometry() +{ + m_BarrelDigitizer->setGeometry( + m_Geometry->getSubdetectorGeometry( DetId::Ecal, EcalBarrel ) ) ; + m_EndcapDigitizer->setGeometry( + m_Geometry->getSubdetectorGeometry( DetId::Ecal, EcalEndcap ) ) ; +} diff --git a/SimG4CMS/Calo/interface/ECalSD.h b/SimG4CMS/Calo/interface/ECalSD.h index 7d212e31e14d2..aa55721d44180 100644 --- a/SimG4CMS/Calo/interface/ECalSD.h +++ b/SimG4CMS/Calo/interface/ECalSD.h @@ -32,13 +32,12 @@ class ECalSD : public CaloSD { virtual ~ECalSD(); virtual double getEnergyDeposit(G4Step*); virtual uint16_t getRadiationLength(G4Step *); + virtual uint16_t getLayerIDForTimeSim(G4Step *); virtual uint32_t setDetUnitId(G4Step*); void setNumberingScheme(EcalNumberingScheme*); virtual int getTrackID(G4Track*); virtual uint16_t getDepth(G4Step*); - private: - void initMap(G4String, const DDCompactView &); double curve_LY(G4Step*); double crystalLength(G4LogicalVolume*); @@ -48,9 +47,10 @@ class ECalSD : public CaloSD { const DDsvalues_type&); std::vector getStringArray(const std::string&, const DDsvalues_type&); - + bool isEB; + bool isEE; EcalNumberingScheme * numberingScheme; - bool useWeight, storeTrack, storeRL; + bool useWeight, storeTrack, storeRL, storeLayerTimeSim; bool useBirk, useBirkL3; double birk1, birk2, birk3, birkSlope, birkCut; double slopeLY; @@ -60,6 +60,7 @@ class ECalSD : public CaloSD { EcalBaseNumber theBaseNumber; EnergyResolutionVsLumi ageing; bool ageingWithSlopeLY; + }; #endif // ECalSD_h diff --git a/SimG4CMS/Calo/src/ECalSD.cc b/SimG4CMS/Calo/src/ECalSD.cc index 93a0dd91a8b45..22cd609d7f211 100644 --- a/SimG4CMS/Calo/src/ECalSD.cc +++ b/SimG4CMS/Calo/src/ECalSD.cc @@ -66,6 +66,9 @@ ECalSD::ECalSD(G4String name, const DDCompactView & cpv, bool isItTB = m_EC.getUntrackedParameter("TestBeam", false); bool nullNS = m_EC.getUntrackedParameter("NullNumbering", false); storeRL = m_EC.getUntrackedParameter("StoreRadLength", false); + + //Changes for improved timing simulation + storeLayerTimeSim = m_EC.getUntrackedParameter("StoreLayerTimeSim", false); ageingWithSlopeLY = m_EC.getUntrackedParameter("AgeingWithSlopeLY", false); if(ageingWithSlopeLY) ageing.setLumies(p.getParameter("ECalSD").getParameter("DelivLuminosity"), @@ -93,13 +96,24 @@ ECalSD::ECalSD(G4String name, const DDCompactView & cpv, EcalNumberingScheme* scheme=0; if (nullNS) scheme = 0; - else if (name == "EcalHitsEB") scheme = dynamic_cast(new EcalBarrelNumberingScheme()); - else if (name == "EcalHitsEE") scheme = dynamic_cast(new EcalEndcapNumberingScheme()); - else if (name == "EcalHitsES") { - if (isItTB) scheme = dynamic_cast(new ESTBNumberingScheme()); - else scheme = dynamic_cast(new EcalPreshowerNumberingScheme()); - useWeight = false; - } else {edm::LogWarning("EcalSim") << "ECalSD: ReadoutName not supported\n";} + else if (name == "EcalHitsEB") + { + scheme = dynamic_cast(new EcalBarrelNumberingScheme()); + isEB=1; + } + else if (name == "EcalHitsEE") + { + scheme = dynamic_cast(new EcalEndcapNumberingScheme()); + isEE=1; + } + else if (name == "EcalHitsES") + { + if (isItTB) scheme = dynamic_cast(new ESTBNumberingScheme()); + else scheme = dynamic_cast(new EcalPreshowerNumberingScheme()); + useWeight = false; + } + else {edm::LogWarning("EcalSim") << "ECalSD: ReadoutName not supported\n";} + if (scheme) setNumberingScheme(scheme); #ifdef DebugLog @@ -120,14 +134,15 @@ ECalSD::ECalSD(G4String name, const DDCompactView & cpv, } edm::LogInfo("EcalSim") << "ECalSD:: Suppression Flag " << suppressHeavy - << " protons below " << kmaxProton << " MeV," - << " neutrons below " << kmaxNeutron << " MeV and" - << " ions below " << kmaxIon << " MeV\n" - << " Depth1 Name = " << depth1Name - << " and Depth2 Name = " << depth2Name; - + << "\tprotons below " << kmaxProton << " MeV," + << "\tneutrons below " << kmaxNeutron << " MeV" + << "\tions below " << kmaxIon << " MeV" + << "\n\tDepth1 Name = " << depth1Name + << "\tDepth2 Name = " << depth2Name + << "\n\tstoreRL" << storeRL + << "\tstoreLayerTimeSim " << storeLayerTimeSim + << "\n\ttime Granularity " << p.getParameter("ECalSD").getParameter("TimeSliceUnit") << " ns"; if (useWeight) initMap(name,cpv); - } ECalSD::~ECalSD() { @@ -224,14 +239,11 @@ int ECalSD::getTrackID(G4Track* aTrack) { uint16_t ECalSD::getDepth(G4Step * aStep) { G4LogicalVolume* lv = aStep->GetPreStepPoint()->GetTouchable()->GetVolume(0)->GetLogicalVolume(); - uint16_t ret = 0; - if (any(useDepth1,lv)) ret = 1; - else if (any(useDepth2,lv)) ret = 2; - else if (storeRL) ret = getRadiationLength(aStep); -#ifdef DebugLog - LogDebug("EcalSim") << "Volume " << lv->GetName() << " Depth " << ret; -#endif - return ret; + if (any(useDepth1,lv)) return 1; + else if (any(useDepth2,lv)) return 2; + else if (storeRL) return getRadiationLength(aStep); + else if (storeLayerTimeSim) return getLayerIDForTimeSim(aStep); + return 0; } uint16_t ECalSD::getRadiationLength(G4Step * aStep) { @@ -253,6 +265,44 @@ uint16_t ECalSD::getRadiationLength(G4Step * aStep) { return thisX0; } +uint16_t ECalSD::getLayerIDForTimeSim(G4Step * aStep) +{ + constexpr char refl[] = "refl"; + float layerSize = 1*cm; //layer size in cm + if (!isEB && !isEE) + return 0; + + if (aStep != NULL ) { + G4StepPoint* hitPoint = aStep->GetPreStepPoint(); + G4LogicalVolume* lv = hitPoint->GetTouchable()->GetVolume(0)->GetLogicalVolume(); + G4ThreeVector localPoint = setToLocal(hitPoint->GetPosition(), + hitPoint->GetTouchable()); + double crlength = crystalLength(lv); + double detz; + + const auto& name = lv->GetName(); + + if( name.size() > 4 && name.compare(name.size()-4,4,refl) == 0 ) + { + if (isEB) + detz = (float)(0.5*crlength + localPoint.z()); + else + detz = (float)(0.5*crlength - localPoint.z()); + } + else + { + if (isEB) + detz = (float)(0.5*crlength - localPoint.z()); + else + detz = (float)(0.5*crlength + localPoint.z()); + } + if (detz<0) + detz=0; + return 100+(int)detz/layerSize; + } + return 0; +} + uint32_t ECalSD::setDetUnitId(G4Step * aStep) { if (numberingScheme == 0) { return EBDetId(1,1)(); @@ -417,13 +467,6 @@ double ECalSD::curve_LY(G4Step* aStep) { << " take weight = " << weight; } } -#ifdef DebugLog - LogDebug("EcalSim") << "ECalSD, light coll curve : " << dapd - << " crlength = " << crlength - << " crystal name = " << lv->GetName() - << " z of localPoint = " << localPoint.z() - << " take weight = " << weight; -#endif return weight; } diff --git a/SimG4Core/Configuration/python/SimG4Core_cff.py b/SimG4Core/Configuration/python/SimG4Core_cff.py index 8dce116832452..4517d913aefe1 100644 --- a/SimG4Core/Configuration/python/SimG4Core_cff.py +++ b/SimG4Core/Configuration/python/SimG4Core_cff.py @@ -8,3 +8,6 @@ from Configuration.StandardSequences.Eras import eras eras.phase2_hcal.toModify( g4SimHits, HCalSD = dict( TestNumberingScheme = True ) ) +eras.phase2_timing.toModify( g4SimHits.ECalSD, + StoreLayerTimeSim = cms.untracked.bool(True), + TimeSliceUnit = cms.double(0.001) ) diff --git a/SimGeneral/MixingModule/python/digitizers_cfi.py b/SimGeneral/MixingModule/python/digitizers_cfi.py index 030d3aa999e21..71172a0739e16 100644 --- a/SimGeneral/MixingModule/python/digitizers_cfi.py +++ b/SimGeneral/MixingModule/python/digitizers_cfi.py @@ -52,6 +52,10 @@ ) eras.phase2_common.toModify( theDigitizers, castor = None ) + +from SimGeneral.MixingModule.ecalTimeDigitizer_cfi import ecalTimeDigitizer +eras.phase2_timing.toModify( theDigitizers, + ecalTime = ecalTimeDigitizer.clone() ) theDigitizersValid = cms.PSet( theDigitizers, @@ -63,3 +67,5 @@ eras.phase2_hgcal.toModify( theDigitizersValid, calotruth = cms.PSet( caloParticles ) ) + + diff --git a/SimGeneral/MixingModule/python/ecalTimeDigitizer_cfi.py b/SimGeneral/MixingModule/python/ecalTimeDigitizer_cfi.py new file mode 100644 index 0000000000000..6c786830386ec --- /dev/null +++ b/SimGeneral/MixingModule/python/ecalTimeDigitizer_cfi.py @@ -0,0 +1,9 @@ +import FWCore.ParameterSet.Config as cms + +from SimCalorimetry.EcalSimProducers.ecalTimeDigiParameters_cff import * + +ecalTimeDigitizer = cms.PSet( + ecal_time_digi_parameters, + accumulatorType = cms.string("EcalTimeDigiProducer"), +) +