EcalRecProducers/plugins/EcalDetailedTimeRecHitProducer.cc

0001 /** \class  EcalDetailedTimeRecHitProducer
0002  *   produce ECAL rechits associating them with improved ecalTimeDigis
0003  *  \author Paolo Meridiani, INFN Roma
0004  **/
0005
0006 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"
0007 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
0008 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
0009 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0010 #include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
0011 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0014 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0015 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0016 #include "DataFormats/EcalDigi/interface/EcalTimeDigi.h"
0017 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0018 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0019 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0020 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0021 #include "DataFormats/VertexReco/interface/Vertex.h"
0022 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0023 #include "FWCore/Framework/interface/ESHandle.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/Utilities/interface/ESGetToken.h"
0030 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0031 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0032 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0033 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0034 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0035 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0036 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalRecHitAbsAlgo.h"
0037 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalRecHitSimpleAlgo.h"
0038 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0039 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0040
0041 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0042 #include <CLHEP/Units/SystemOfUnits.h>
0043
0044 #include <cmath>
0045 #include <iostream>
0046 #include <memory>
0047 #include <vector>
0048
0049 class EcalDetailedTimeRecHitProducer : public edm::stream::EDProducer<> {
0050 public:
0051   explicit EcalDetailedTimeRecHitProducer(const edm::ParameterSet& ps);
0052   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0053
0054 private:
0055   //Functions to correct the TOF from the EcalDigi which is not corrected for the vertex position
0056   double deltaTimeOfFlight(GlobalPoint& vertex, const DetId& detId, int layer) const;
0057
0058   const CaloGeometry* m_geometry;
0059
0060   edm::EDGetTokenT<EBRecHitCollection> EBRecHitCollection_;  // secondary name given to collection of EBrechits
0061   edm::EDGetTokenT<EERecHitCollection> EERecHitCollection_;  // secondary name given to collection of EErechits
0062
0063   edm::EDGetTokenT<reco::VertexCollection> recoVertex_;
0064   edm::EDGetTokenT<edm::SimVertexContainer> simVertex_;
0065   bool correctForVertexZPosition_;
0066   bool useMCTruthVertex_;
0067
0068   int ebTimeLayer_;
0069   int eeTimeLayer_;
0070
0071   edm::EDGetTokenT<EcalTimeDigiCollection>
0072       ebTimeDigiCollection_;  // secondary name given to collection of EB uncalib rechits
0073   edm::EDGetTokenT<EcalTimeDigiCollection>
0074       eeTimeDigiCollection_;  // secondary name given to collection of EE uncalib rechits
0075
0076   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometry_;
0077
0078   std::string EBDetailedTimeRecHitCollection_;  // secondary name to be given to EB collection of hits
0079   std::string EEDetailedTimeRecHitCollection_;  // secondary name to be given to EE collection of hits
0080 };
0081
0082 EcalDetailedTimeRecHitProducer::EcalDetailedTimeRecHitProducer(const edm::ParameterSet& ps) : m_geometry(nullptr) {
0083   EBRecHitCollection_ = consumes<EBRecHitCollection>(ps.getParameter<edm::InputTag>("EBRecHitCollection"));
0084   EERecHitCollection_ = consumes<EERecHitCollection>(ps.getParameter<edm::InputTag>("EERecHitCollection"));
0085
0086   ebTimeDigiCollection_ = consumes<EcalTimeDigiCollection>(ps.getParameter<edm::InputTag>("EBTimeDigiCollection"));
0087   eeTimeDigiCollection_ = consumes<EcalTimeDigiCollection>(ps.getParameter<edm::InputTag>("EETimeDigiCollection"));
0088
0089   caloGeometry_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0090
0091   EBDetailedTimeRecHitCollection_ = ps.getParameter<std::string>("EBDetailedTimeRecHitCollection");
0092   EEDetailedTimeRecHitCollection_ = ps.getParameter<std::string>("EEDetailedTimeRecHitCollection");
0093
0094   correctForVertexZPosition_ = ps.getParameter<bool>("correctForVertexZPosition");
0095   useMCTruthVertex_ = ps.getParameter<bool>("useMCTruthVertex");
0096   if (correctForVertexZPosition_) {
0097     if (not useMCTruthVertex_) {
0098       recoVertex_ = consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("recoVertex"));
0099     } else {
0100       simVertex_ = consumes<edm::SimVertexContainer>(ps.getParameter<edm::InputTag>("simVertex"));
0101     }
0102   }
0103
0104   ebTimeLayer_ = ps.getParameter<int>("EBTimeLayer");
0105   eeTimeLayer_ = ps.getParameter<int>("EETimeLayer");
0106
0107   produces<EBRecHitCollection>(EBDetailedTimeRecHitCollection_);
0108   produces<EERecHitCollection>(EEDetailedTimeRecHitCollection_);
0109 }
0110
0111 void EcalDetailedTimeRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0112   using namespace edm;
0113   using namespace reco;
0114
0115   edm::ESHandle<CaloGeometry> hGeometry = es.getHandle(caloGeometry_);
0116
0117   m_geometry = hGeometry.product();
0118
0119   Handle<EBRecHitCollection> pEBRecHits;
0120   Handle<EERecHitCollection> pEERecHits;
0121
0122   const EBRecHitCollection* EBRecHits = nullptr;
0123   const EERecHitCollection* EERecHits = nullptr;
0124
0125   evt.getByToken(EBRecHitCollection_, pEBRecHits);
0126   if (pEBRecHits.isValid()) {
0127     EBRecHits = pEBRecHits.product();  // get a ptr to the product
0128 #ifdef DEBUG
0129     LogDebug("EcalRecHitDebug") << "total # EB rechits to be re-calibrated: " << EBRecHits->size();
0130 #endif
0131   }
0132
0133   evt.getByToken(EERecHitCollection_, pEERecHits);
0134   if (pEERecHits.isValid()) {
0135     EERecHits = pEERecHits.product();  // get a ptr to the product
0136 #ifdef DEBUG
0137     LogDebug("EcalRecHitDebug") << "total # EE uncalibrated rechits to be re-calibrated: " << EERecHits->size();
0138 #endif
0139   }
0140
0141   Handle<EcalTimeDigiCollection> pEBTimeDigis;
0142   Handle<EcalTimeDigiCollection> pEETimeDigis;
0143
0144   const EcalTimeDigiCollection* ebTimeDigis = nullptr;
0145   const EcalTimeDigiCollection* eeTimeDigis = nullptr;
0146
0147   evt.getByToken(ebTimeDigiCollection_, pEBTimeDigis);
0148   //evt.getByToken( digiProducer_, pEBTimeDigis);
0149   if (pEBTimeDigis.isValid()) {
0150     ebTimeDigis = pEBTimeDigis.product();  // get a ptr to the produc
0151     edm::LogInfo("EcalDetailedTimeRecHitInfo") << "total # ebTimeDigis: " << ebTimeDigis->size();
0152   }
0153
0154   evt.getByToken(eeTimeDigiCollection_, pEETimeDigis);
0155   //evt.getByToken( digiProducer_, pEETimeDigis);
0156   if (pEETimeDigis.isValid()) {
0157     eeTimeDigis = pEETimeDigis.product();  // get a ptr to the product
0158     edm::LogInfo("EcalDetailedTimeRecHitInfo") << "total # eeTimeDigis: " << eeTimeDigis->size();
0159   }
0160   // collection of rechits to put in the event
0161   std::unique_ptr<EBRecHitCollection> EBDetailedTimeRecHits(new EBRecHitCollection);
0162   std::unique_ptr<EERecHitCollection> EEDetailedTimeRecHits(new EERecHitCollection);
0163
0164   std::unique_ptr<GlobalPoint> vertex;
0165
0166   if (correctForVertexZPosition_) {
0167     if (!useMCTruthVertex_) {
0168       //Get the first reco vertex
0169       // get primary vertices
0170
0171       edm::Handle<VertexCollection> VertexHandle;
0172       evt.getByToken(recoVertex_, VertexHandle);
0173
0174       if (VertexHandle.isValid()) {
0175         if (!(*VertexHandle).empty())  //at least 1 vertex
0176         {
0177           const reco::Vertex* myVertex = &(*VertexHandle)[0];
0178           vertex = std::make_unique<GlobalPoint>(myVertex->x(), myVertex->y(), myVertex->z());
0179         }
0180       }
0181
0182     } else {
0183       edm::Handle<SimVertexContainer> VertexHandle;
0184       evt.getByToken(simVertex_, VertexHandle);
0185
0186       if (VertexHandle.isValid()) {
0187         if (!(*VertexHandle).empty())  //at least 1 vertex
0188         {
0189           assert((*VertexHandle)[0].vertexId() == 0);
0190           const SimVertex* myVertex = &(*VertexHandle)[0];
0191           vertex = std::make_unique<GlobalPoint>(
0192               myVertex->position().x(), myVertex->position().y(), myVertex->position().z());
0193         }
0194       }
0195     }
0196   }
0197
0198   if (EBRecHits && ebTimeDigis) {
0199     // loop over uncalibrated rechits to make calibrated ones
0200     for (EBRecHitCollection::const_iterator it = EBRecHits->begin(); it != EBRecHits->end(); ++it) {
0201       EcalRecHit aHit((*it));
0202       EcalTimeDigiCollection::const_iterator timeDigi = ebTimeDigis->find((*it).id());
0203       if (timeDigi != ebTimeDigis->end()) {
0204         if (timeDigi->sampleOfInterest() >= 0) {
0205           float myTime = (*timeDigi)[timeDigi->sampleOfInterest()];
0206           //Vertex corrected ToF
0207           if (vertex) {
0208             aHit.setTime(myTime + deltaTimeOfFlight(*vertex, (*it).id(), ebTimeLayer_));
0209           } else
0210             //Uncorrected ToF
0211             aHit.setTime(myTime);
0212         }
0213       }
0214       // leave standard time if no timeDigi is associated (e.g. noise recHits)
0215       EBDetailedTimeRecHits->push_back(aHit);
0216     }
0217   }
0218
0219   if (EERecHits && eeTimeDigis) {
0220     // loop over uncalibrated rechits to make calibrated ones
0221     for (EERecHitCollection::const_iterator it = EERecHits->begin(); it != EERecHits->end(); ++it) {
0222       EcalRecHit aHit(*it);
0223       EcalTimeDigiCollection::const_iterator timeDigi = eeTimeDigis->find((*it).id());
0224       if (timeDigi != eeTimeDigis->end()) {
0225         if (timeDigi->sampleOfInterest() >= 0) {
0226           float myTime = (*timeDigi)[timeDigi->sampleOfInterest()];
0227           //Vertex corrected ToF
0228           if (vertex) {
0229             aHit.setTime(myTime + deltaTimeOfFlight(*vertex, (*it).id(), eeTimeLayer_));
0230           } else
0231             //Uncorrected ToF
0232             aHit.setTime(myTime);
0233         }
0234       }
0235       EEDetailedTimeRecHits->push_back(aHit);
0236     }
0237   }
0238   // put the collection of recunstructed hits in the event
0239   LogInfo("EcalDetailedTimeRecHitInfo") << "total # EB rechits: " << EBDetailedTimeRecHits->size();
0240   LogInfo("EcalDetailedTimeRecHitInfo") << "total # EE rechits: " << EEDetailedTimeRecHits->size();
0241
0242   evt.put(std::move(EBDetailedTimeRecHits), EBDetailedTimeRecHitCollection_);
0243   evt.put(std::move(EEDetailedTimeRecHits), EEDetailedTimeRecHitCollection_);
0244 }
0245
0246 double EcalDetailedTimeRecHitProducer::deltaTimeOfFlight(GlobalPoint& vertex, const DetId& detId, int layer) const {
0247   auto cellGeometry(m_geometry->getGeometry(detId));
0248   assert(nullptr != cellGeometry);
0249   GlobalPoint layerPos =
0250       cellGeometry->getPosition(double(layer) + 0.5);  //depth in mm in the middle of the layer position
0251   GlobalVector tofVector = layerPos - vertex;
0252   return (layerPos.mag() * CLHEP::cm - tofVector.mag() * CLHEP::cm) / (float)c_light;
0253 }
0254
0255 #include "FWCore/Framework/interface/MakerMacros.h"
0256 DEFINE_FWK_MODULE(EcalDetailedTimeRecHitProducer);