HcalIsolatedTrackReco/plugins/EcalIsolatedParticleCandidateProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EcalIsolatedParticleCandidateProducer
0004 // Class:      EcalIsolatedParticleCandidateProducer
0005 //
0006 /**\class EcalIsolatedParticleCandidateProducer EcalIsolatedParticleCandidateProducer.cc Calibration/EcalIsolatedParticleCandidateProducer/src/EcalIsolatedParticleCandidateProducer.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Grigory Safronov
0015 //         Created:  Thu Jun  7 17:21:58 MSD 2007
0016 //
0017 //
0018
0019 // system include files
0020 #include <cmath>
0021 #include <memory>
0022
0023 // user include files
0024
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/global/EDProducer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0032
0033 #include "DataFormats/Common/interface/Ref.h"
0034 #include "DataFormats/DetId/interface/DetId.h"
0035 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0036 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0037 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0038 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0039
0040 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0041 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0042
0043 //#define EDM_ML_DEBUG
0044
0045 //
0046 // class decleration
0047 //
0048
0049 class EcalIsolatedParticleCandidateProducer : public edm::global::EDProducer<> {
0050 public:
0051   explicit EcalIsolatedParticleCandidateProducer(const edm::ParameterSet&);
0052   ~EcalIsolatedParticleCandidateProducer() override;
0053
0054   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0055
0056 private:
0057   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0058   void beginJob() override {}
0059   void endJob() override {}
0060
0061   const double InConeSize_;
0062   const double OutConeSize_;
0063   const double hitCountEthr_;
0064   const double hitEthr_;
0065
0066   const edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_l1tau_;
0067   const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_hlt_;
0068   const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0069   const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0070
0071   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0072
0073   // ----------member data ---------------------------
0074 };
0075
0076 EcalIsolatedParticleCandidateProducer::EcalIsolatedParticleCandidateProducer(const edm::ParameterSet& conf)
0077     : InConeSize_(conf.getParameter<double>("EcalInnerConeSize")),
0078       OutConeSize_(conf.getParameter<double>("EcalOuterConeSize")),
0079       hitCountEthr_(conf.getParameter<double>("ECHitCountEnergyThreshold")),
0080       hitEthr_(conf.getParameter<double>("ECHitEnergyThreshold")),
0081       tok_l1tau_(consumes<l1extra::L1JetParticleCollection>(conf.getParameter<edm::InputTag>("L1eTauJetsSource"))),
0082       tok_hlt_(consumes<trigger::TriggerFilterObjectWithRefs>(conf.getParameter<edm::InputTag>("L1GTSeedLabel"))),
0083       tok_EB_(consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("EBrecHitCollectionLabel"))),
0084       tok_EE_(consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("EErecHitCollectionLabel"))),
0085       tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()) {
0086   //register your products
0087   produces<reco::IsolatedPixelTrackCandidateCollection>();
0088 }
0089
0090 EcalIsolatedParticleCandidateProducer::~EcalIsolatedParticleCandidateProducer() {
0091   // do anything here that needs to be done at desctruction time
0092   // (e.g. close files, deallocate resources etc.)
0093 }
0094
0095 //
0096 // member functions
0097 //
0098
0099 // ------------ method called to produce the data  ------------
0100 void EcalIsolatedParticleCandidateProducer::produce(edm::StreamID,
0101                                                     edm::Event& iEvent,
0102                                                     const edm::EventSetup& iSetup) const {
0103 #ifdef EDM_ML_DEBUG
0104   edm::LogVerbatim("HcalIsoTrack") << "get tau";
0105 #endif
0106   const edm::Handle<l1extra::L1JetParticleCollection>& l1Taus = iEvent.getHandle(tok_l1tau_);
0107
0108 #ifdef EDM_ML_DEBUG
0109   edm::LogVerbatim("HcalIsoTrack") << "get geom";
0110 #endif
0111   const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0112
0113 #ifdef EDM_ML_DEBUG
0114   edm::LogVerbatim("HcalIsoTrack") << "get ec rechit";
0115 #endif
0116   const edm::Handle<EcalRecHitCollection>& ecalEB = iEvent.getHandle(tok_EB_);
0117   const edm::Handle<EcalRecHitCollection>& ecalEE = iEvent.getHandle(tok_EE_);
0118
0119 #ifdef EDM_ML_DEBUG
0120   edm::LogVerbatim("HcalIsoTrack") << "get l1 trig obj";
0121 #endif
0122
0123   const edm::Handle<trigger::TriggerFilterObjectWithRefs>& l1trigobj = iEvent.getHandle(tok_hlt_);
0124
0125   std::vector<edm::Ref<l1t::TauBxCollection> > l1tauobjref;
0126   std::vector<edm::Ref<l1t::JetBxCollection> > l1jetobjref;
0127
0128   l1trigobj->getObjects(trigger::TriggerL1Tau, l1tauobjref);
0129   l1trigobj->getObjects(trigger::TriggerL1Jet, l1jetobjref);
0130
0131   double ptTriggered = -10;
0132   double etaTriggered = -100;
0133   double phiTriggered = -100;
0134
0135 #ifdef EDM_ML_DEBUG
0136   edm::LogVerbatim("HcalIsoTrack") << "find highest pT triggered obj";
0137 #endif
0138   for (unsigned int p = 0; p < l1tauobjref.size(); p++) {
0139     if (l1tauobjref[p]->pt() > ptTriggered) {
0140       ptTriggered = l1tauobjref[p]->pt();
0141       phiTriggered = l1tauobjref[p]->phi();
0142       etaTriggered = l1tauobjref[p]->eta();
0143     }
0144   }
0145   for (unsigned int p = 0; p < l1jetobjref.size(); p++) {
0146     if (l1jetobjref[p]->pt() > ptTriggered) {
0147       ptTriggered = l1jetobjref[p]->pt();
0148       phiTriggered = l1jetobjref[p]->phi();
0149       etaTriggered = l1jetobjref[p]->eta();
0150     }
0151   }
0152
0153   auto iptcCollection = std::make_unique<reco::IsolatedPixelTrackCandidateCollection>();
0154
0155 #ifdef EDM_ML_DEBUG
0156   edm::LogVerbatim("HcalIsoTrack") << "loop over l1taus";
0157 #endif
0158   for (l1extra::L1JetParticleCollection::const_iterator tit = l1Taus->begin(); tit != l1Taus->end(); tit++) {
0159     double dphi = fabs(tit->phi() - phiTriggered);
0160     if (dphi > M_PI)
0161       dphi = 2 * M_PI - dphi;
0162     double Rseed = sqrt(pow(etaTriggered - tit->eta(), 2) + dphi * dphi);
0163     if (Rseed < 1.2)
0164       continue;
0165     int nhitOut = 0;
0166     int nhitIn = 0;
0167     double OutEnergy = 0;
0168     double InEnergy = 0;
0169 #ifdef EDM_ML_DEBUG
0170     edm::LogVerbatim("HcalIsoTrack") << "loops over rechits";
0171 #endif
0172     for (EcalRecHitCollection::const_iterator eItr = ecalEB->begin(); eItr != ecalEB->end(); eItr++) {
0173       double phiD, R;
0174       const GlobalPoint& pos = geo->getPosition(eItr->detid());
0175       double phihit = pos.phi();
0176       double etahit = pos.eta();
0177       phiD = fabs(phihit - tit->phi());
0178       if (phiD > M_PI)
0179         phiD = 2 * M_PI - phiD;
0180       R = sqrt(pow(etahit - tit->eta(), 2) + phiD * phiD);
0181
0182       if (R < OutConeSize_ && R > InConeSize_ && eItr->energy() > hitCountEthr_) {
0183         nhitOut++;
0184       }
0185       if (R < InConeSize_ && eItr->energy() > hitCountEthr_) {
0186         nhitIn++;
0187       }
0188
0189       if (R < OutConeSize_ && R > InConeSize_ && eItr->energy() > hitEthr_) {
0190         OutEnergy += eItr->energy();
0191       }
0192       if (R < InConeSize_ && eItr->energy() > hitEthr_) {
0193         InEnergy += eItr->energy();
0194       }
0195     }
0196
0197     for (EcalRecHitCollection::const_iterator eItr = ecalEE->begin(); eItr != ecalEE->end(); eItr++) {
0198       double phiD, R;
0199       const GlobalPoint& pos = geo->getPosition(eItr->detid());
0200       double phihit = pos.phi();
0201       double etahit = pos.eta();
0202       phiD = fabs(phihit - tit->phi());
0203       if (phiD > M_PI)
0204         phiD = 2 * M_PI - phiD;
0205       R = sqrt(pow(etahit - tit->eta(), 2) + phiD * phiD);
0206       if (R < OutConeSize_ && R > InConeSize_ && eItr->energy() > hitCountEthr_) {
0207         nhitOut++;
0208       }
0209       if (R < InConeSize_ && eItr->energy() > hitCountEthr_) {
0210         nhitIn++;
0211       }
0212       if (R < OutConeSize_ && R > InConeSize_ && eItr->energy() > hitEthr_) {
0213         OutEnergy += eItr->energy();
0214       }
0215       if (R < InConeSize_ && eItr->energy() > hitEthr_) {
0216         InEnergy += eItr->energy();
0217       }
0218     }
0219 #ifdef EDM_ML_DEBUG
0220     edm::LogVerbatim("HcalIsoTrack") << "create and push_back candidate";
0221 #endif
0222     reco::IsolatedPixelTrackCandidate newca(
0223         l1extra::L1JetParticleRef(l1Taus, tit - l1Taus->begin()), InEnergy, OutEnergy, nhitIn, nhitOut);
0224     iptcCollection->push_back(newca);
0225   }
0226
0227   //Use the ExampleData to create an ExampleData2 which
0228   // is put into the Event
0229
0230 #ifdef EDM_ML_DEBUG
0231   edm::LogVerbatim("HcalIsoTrack") << "put cand into event";
0232 #endif
0233   iEvent.put(std::move(iptcCollection));
0234 }
0235
0236 void EcalIsolatedParticleCandidateProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0237   edm::ParameterSetDescription desc;
0238   desc.add<double>("ECHitEnergyThreshold", 0.05);
0239   desc.add<edm::InputTag>("L1eTauJetsSource", edm::InputTag("l1extraParticles", "Tau"));
0240   desc.add<edm::InputTag>("L1GTSeedLabel", edm::InputTag("l1sIsolTrack"));
0241   desc.add<edm::InputTag>("EBrecHitCollectionLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0242   desc.add<edm::InputTag>("EErecHitCollectionLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0243   desc.add<double>("ECHitCountEnergyThreshold", 0.5);
0244   desc.add<double>("EcalInnerConeSize", 0.3);
0245   desc.add<double>("EcalOuterConeSize", 0.7);
0246   descriptions.add("ecalIsolPartProd", desc);
0247 }
0248
0249 #include "FWCore/PluginManager/interface/ModuleDef.h"
0250 #include "FWCore/Framework/interface/MakerMacros.h"
0251
0252 DEFINE_FWK_MODULE(EcalIsolatedParticleCandidateProducer);