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