EgammaHLTProducers/plugins/HLTEcalPFClusterIsolationProducer.cc

0001 #include <iostream>
0002 #include <vector>
0003 #include <memory>
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0013
0014 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0015 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0016
0017 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0019
0020 #include "DataFormats/Math/interface/deltaR.h"
0021
0022 #include "RecoEgamma/EgammaIsolationAlgos/interface/EcalPFClusterIsolation.h"
0023 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0024
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032
0033 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0034 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0035
0036 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0037 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateIsolation.h"
0038
0039 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0040 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0041
0042 template <typename T1>
0043 class HLTEcalPFClusterIsolationProducer : public edm::stream::EDProducer<> {
0044   typedef std::vector<T1> T1Collection;
0045   typedef edm::Ref<T1Collection> T1Ref;
0046   typedef edm::AssociationMap<edm::OneToValue<std::vector<T1>, float>> T1IsolationMap;
0047
0048 public:
0049   explicit HLTEcalPFClusterIsolationProducer(const edm::ParameterSet&);
0050   ~HLTEcalPFClusterIsolationProducer() override;
0051
0052   void produce(edm::Event&, const edm::EventSetup&) override;
0053   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0054
0055 private:
0056   bool computedRVeto(T1Ref candRef, reco::PFClusterRef pfclu);
0057
0058   edm::EDGetTokenT<T1Collection> recoCandidateProducer_;
0059   const edm::EDGetTokenT<reco::PFClusterCollection> pfClusterProducer_;
0060   const edm::EDGetTokenT<double> rhoProducer_;
0061
0062   const double drMax_;
0063   const double drVetoBarrel_;
0064   const double drVetoEndcap_;
0065   const double etaStripBarrel_;
0066   const double etaStripEndcap_;
0067   const double energyBarrel_;
0068   const double energyEndcap_;
0069
0070   const bool doRhoCorrection_;
0071   const double rhoMax_;
0072   const double rhoScale_;
0073   const std::vector<double> effectiveAreas_;
0074   const std::vector<double> absEtaLowEdges_;
0075 };
0076
0077 template <typename T1>
0078 HLTEcalPFClusterIsolationProducer<T1>::HLTEcalPFClusterIsolationProducer(const edm::ParameterSet& config)
0079     : pfClusterProducer_(consumes<reco::PFClusterCollection>(config.getParameter<edm::InputTag>("pfClusterProducer"))),
0080       rhoProducer_(consumes<double>(config.getParameter<edm::InputTag>("rhoProducer"))),
0081       drMax_(config.getParameter<double>("drMax")),
0082       drVetoBarrel_(config.getParameter<double>("drVetoBarrel")),
0083       drVetoEndcap_(config.getParameter<double>("drVetoEndcap")),
0084       etaStripBarrel_(config.getParameter<double>("etaStripBarrel")),
0085       etaStripEndcap_(config.getParameter<double>("etaStripEndcap")),
0086       energyBarrel_(config.getParameter<double>("energyBarrel")),
0087       energyEndcap_(config.getParameter<double>("energyEndcap")),
0088       doRhoCorrection_(config.getParameter<bool>("doRhoCorrection")),
0089       rhoMax_(config.getParameter<double>("rhoMax")),
0090       rhoScale_(config.getParameter<double>("rhoScale")),
0091       effectiveAreas_(config.getParameter<std::vector<double>>("effectiveAreas")),
0092       absEtaLowEdges_(config.getParameter<std::vector<double>>("absEtaLowEdges")) {
0093   if (doRhoCorrection_) {
0094     if (absEtaLowEdges_.size() != effectiveAreas_.size())
0095       throw cms::Exception("IncompatibleVects") << "absEtaLowEdges and effectiveAreas should be of the same size. \n";
0096
0097     if (absEtaLowEdges_.at(0) != 0.0)
0098       throw cms::Exception("IncompleteCoverage") << "absEtaLowEdges should start from 0. \n";
0099
0100     for (unsigned int aIt = 0; aIt < absEtaLowEdges_.size() - 1; aIt++) {
0101       if (!(absEtaLowEdges_.at(aIt) < absEtaLowEdges_.at(aIt + 1)))
0102         throw cms::Exception("ImproperBinning") << "absEtaLowEdges entries should be in increasing order. \n";
0103     }
0104   }
0105
0106   std::string recoCandidateProducerName = "recoCandidateProducer";
0107   if ((typeid(HLTEcalPFClusterIsolationProducer<T1>) ==
0108        typeid(HLTEcalPFClusterIsolationProducer<reco::RecoEcalCandidate>)))
0109     recoCandidateProducerName = "recoEcalCandidateProducer";
0110
0111   recoCandidateProducer_ = consumes<T1Collection>(config.getParameter<edm::InputTag>(recoCandidateProducerName));
0112   produces<T1IsolationMap>();
0113 }
0114
0115 template <typename T1>
0116 HLTEcalPFClusterIsolationProducer<T1>::~HLTEcalPFClusterIsolationProducer() {}
0117
0118 template <typename T1>
0119 void HLTEcalPFClusterIsolationProducer<T1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0120   std::string recoCandidateProducerName = "recoCandidateProducer";
0121   if ((typeid(HLTEcalPFClusterIsolationProducer<T1>) ==
0122        typeid(HLTEcalPFClusterIsolationProducer<reco::RecoEcalCandidate>)))
0123     recoCandidateProducerName = "recoEcalCandidateProducer";
0124
0125   edm::ParameterSetDescription desc;
0126   desc.add<edm::InputTag>(recoCandidateProducerName, edm::InputTag("hltL1SeededRecoEcalCandidatePF"));
0127   desc.add<edm::InputTag>("pfClusterProducer", edm::InputTag("hltParticleFlowClusterECAL"));
0128   desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAllCalo"));
0129   desc.add<bool>("doRhoCorrection", false);
0130   desc.add<double>("rhoMax", 9.9999999E7);
0131   desc.add<double>("rhoScale", 1.0);
0132   desc.add<double>("drMax", 0.3);
0133   desc.add<double>("drVetoBarrel", 0.0);
0134   desc.add<double>("drVetoEndcap", 0.0);
0135   desc.add<double>("etaStripBarrel", 0.0);
0136   desc.add<double>("etaStripEndcap", 0.0);
0137   desc.add<double>("energyBarrel", 0.0);
0138   desc.add<double>("energyEndcap", 0.0);
0139   desc.add<std::vector<double>>("effectiveAreas", {0.29, 0.21});  // 2016 post-ichep sinEle default

0140   desc.add<std::vector<double>>("absEtaLowEdges", {0.0, 1.479});  // Barrel, Endcap

0141   descriptions.add(defaultModuleLabel<HLTEcalPFClusterIsolationProducer<T1>>(), desc);
0142 }
0143
0144 template <typename T1>
0145 void HLTEcalPFClusterIsolationProducer<T1>::produce(edm::Event& iEvent, const edm::EventSetup&) {
0146   edm::Handle<double> rhoHandle;
0147   double rho = 0.0;
0148   if (doRhoCorrection_) {
0149     iEvent.getByToken(rhoProducer_, rhoHandle);
0150     rho = *(rhoHandle.product());
0151   }
0152
0153   if (rho > rhoMax_)
0154     rho = rhoMax_;
0155
0156   rho = rho * rhoScale_;
0157
0158   edm::Handle<T1Collection> recoCandHandle;
0159   edm::Handle<reco::PFClusterCollection> clusterHandle;
0160
0161   iEvent.getByToken(recoCandidateProducer_, recoCandHandle);
0162   iEvent.getByToken(pfClusterProducer_, clusterHandle);
0163
0164   EcalPFClusterIsolation<T1> isoAlgo(
0165       drMax_, drVetoBarrel_, drVetoEndcap_, etaStripBarrel_, etaStripEndcap_, energyBarrel_, energyEndcap_);
0166   T1IsolationMap recoCandMap(recoCandHandle);
0167
0168   for (unsigned int iReco = 0; iReco < recoCandHandle->size(); iReco++) {
0169     T1Ref candRef(recoCandHandle, iReco);
0170
0171     float sum = isoAlgo.getSum(candRef, clusterHandle);
0172
0173     if (doRhoCorrection_) {
0174       int iEA = -1;
0175       auto cEta = std::abs(candRef->eta());
0176       for (int bIt = absEtaLowEdges_.size() - 1; bIt > -1; bIt--) {
0177         if (cEta >= absEtaLowEdges_[bIt]) {
0178           iEA = bIt;
0179           break;
0180         }
0181       }
0182       sum = sum - rho * effectiveAreas_[iEA];
0183     }
0184
0185     recoCandMap.insert(candRef, sum);
0186   }
0187
0188   iEvent.put(std::make_unique<T1IsolationMap>(recoCandMap));
0189 }
0190
0191 typedef HLTEcalPFClusterIsolationProducer<reco::RecoEcalCandidate> EgammaHLTEcalPFClusterIsolationProducer;
0192 typedef HLTEcalPFClusterIsolationProducer<reco::RecoChargedCandidate> MuonHLTEcalPFClusterIsolationProducer;
0193
0194 DEFINE_FWK_MODULE(EgammaHLTEcalPFClusterIsolationProducer);
0195 DEFINE_FWK_MODULE(MuonHLTEcalPFClusterIsolationProducer);