Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:09

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 "RecoEgamma/EgammaIsolationAlgos/interface/HcalPFClusterIsolation.h"
0015 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0016 
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 #include "FWCore/Framework/interface/global/EDProducer.h"
0019 
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 
0025 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0026 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0027 
0028 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0029 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateIsolation.h"
0030 
0031 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0032 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0033 
0034 template <typename T1>
0035 class HLTHcalPFClusterIsolationProducer : public edm::global::EDProducer<> {
0036   typedef std::vector<T1> T1Collection;
0037   typedef edm::Ref<T1Collection> T1Ref;
0038   typedef edm::AssociationMap<edm::OneToValue<std::vector<T1>, float>> T1IsolationMap;
0039 
0040 public:
0041   explicit HLTHcalPFClusterIsolationProducer(const edm::ParameterSet&);
0042   ~HLTHcalPFClusterIsolationProducer() override;
0043 
0044   void produce(edm::StreamID sid, edm::Event&, const edm::EventSetup&) const override;
0045   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0046 
0047 private:
0048   edm::EDGetTokenT<T1Collection> recoCandidateProducer_;
0049   const edm::EDGetTokenT<reco::PFClusterCollection> pfClusterProducerHCAL_;
0050   const edm::EDGetTokenT<double> rhoProducer_;
0051   const edm::EDGetTokenT<reco::PFClusterCollection> pfClusterProducerHFEM_;
0052   const edm::EDGetTokenT<reco::PFClusterCollection> pfClusterProducerHFHAD_;
0053 
0054   const bool useHF_;
0055 
0056   const double drMax_;
0057   const double drVetoBarrel_;
0058   const double drVetoEndcap_;
0059   const double etaStripBarrel_;
0060   const double etaStripEndcap_;
0061   const double energyBarrel_;
0062   const double energyEndcap_;
0063   const bool useEt_;
0064 
0065   const bool doRhoCorrection_;
0066   const double rhoMax_;
0067   const double rhoScale_;
0068   const std::vector<double> effectiveAreas_;
0069   const std::vector<double> absEtaLowEdges_;
0070 };
0071 
0072 template <typename T1>
0073 HLTHcalPFClusterIsolationProducer<T1>::HLTHcalPFClusterIsolationProducer(const edm::ParameterSet& config)
0074     : pfClusterProducerHCAL_(
0075           consumes<reco::PFClusterCollection>(config.getParameter<edm::InputTag>("pfClusterProducerHCAL"))),
0076       rhoProducer_(consumes<double>(config.getParameter<edm::InputTag>("rhoProducer"))),
0077       pfClusterProducerHFEM_(
0078           consumes<reco::PFClusterCollection>(config.getParameter<edm::InputTag>("pfClusterProducerHFEM"))),
0079       pfClusterProducerHFHAD_(
0080           consumes<reco::PFClusterCollection>(config.getParameter<edm::InputTag>("pfClusterProducerHFHAD"))),
0081       useHF_(config.getParameter<bool>("useHF")),
0082       drMax_(config.getParameter<double>("drMax")),
0083       drVetoBarrel_(config.getParameter<double>("drVetoBarrel")),
0084       drVetoEndcap_(config.getParameter<double>("drVetoEndcap")),
0085       etaStripBarrel_(config.getParameter<double>("etaStripBarrel")),
0086       etaStripEndcap_(config.getParameter<double>("etaStripEndcap")),
0087       energyBarrel_(config.getParameter<double>("energyBarrel")),
0088       energyEndcap_(config.getParameter<double>("energyEndcap")),
0089       useEt_(config.getParameter<bool>("useEt")),
0090       doRhoCorrection_(config.getParameter<bool>("doRhoCorrection")),
0091       rhoMax_(config.getParameter<double>("rhoMax")),
0092       rhoScale_(config.getParameter<double>("rhoScale")),
0093       effectiveAreas_(config.getParameter<std::vector<double>>("effectiveAreas")),
0094       absEtaLowEdges_(config.getParameter<std::vector<double>>("absEtaLowEdges")) {
0095   if (doRhoCorrection_) {
0096     if (absEtaLowEdges_.size() != effectiveAreas_.size())
0097       throw cms::Exception("IncompatibleVects") << "absEtaLowEdges and effectiveAreas should be of the same size. \n";
0098 
0099     if (absEtaLowEdges_.at(0) != 0.0)
0100       throw cms::Exception("IncompleteCoverage") << "absEtaLowEdges should start from 0. \n";
0101 
0102     for (unsigned int aIt = 0; aIt < absEtaLowEdges_.size() - 1; aIt++) {
0103       if (!(absEtaLowEdges_.at(aIt) < absEtaLowEdges_.at(aIt + 1)))
0104         throw cms::Exception("ImproperBinning") << "absEtaLowEdges entries should be in increasing order. \n";
0105     }
0106   }
0107 
0108   std::string recoCandidateProducerName = "recoCandidateProducer";
0109   if ((typeid(HLTHcalPFClusterIsolationProducer<T1>) ==
0110        typeid(HLTHcalPFClusterIsolationProducer<reco::RecoEcalCandidate>)))
0111     recoCandidateProducerName = "recoEcalCandidateProducer";
0112   recoCandidateProducer_ = consumes<T1Collection>(config.getParameter<edm::InputTag>(recoCandidateProducerName));
0113 
0114   produces<T1IsolationMap>();
0115 }
0116 
0117 template <typename T1>
0118 HLTHcalPFClusterIsolationProducer<T1>::~HLTHcalPFClusterIsolationProducer() {}
0119 
0120 template <typename T1>
0121 void HLTHcalPFClusterIsolationProducer<T1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0122   std::string recoCandidateProducerName = "recoCandidateProducer";
0123   if ((typeid(HLTHcalPFClusterIsolationProducer<T1>) ==
0124        typeid(HLTHcalPFClusterIsolationProducer<reco::RecoEcalCandidate>)))
0125     recoCandidateProducerName = "recoEcalCandidateProducer";
0126 
0127   edm::ParameterSetDescription desc;
0128   desc.add<edm::InputTag>(recoCandidateProducerName, edm::InputTag("hltL1SeededRecoEcalCandidatePF"));
0129   desc.add<edm::InputTag>("pfClusterProducerHCAL", edm::InputTag("hltParticleFlowClusterHCAL"));
0130   desc.ifValue(
0131       edm::ParameterDescription<bool>("useHF", false, true),
0132       true >> (edm::ParameterDescription<edm::InputTag>(
0133                    "pfClusterProducerHFEM", edm::InputTag("hltParticleFlowClusterHFEM"), true) and
0134                edm::ParameterDescription<edm::InputTag>(
0135                    "pfClusterProducerHFHAD", edm::InputTag("hltParticleFlowClusterHFHAD"), true)) or
0136           false >> (edm::ParameterDescription<edm::InputTag>("pfClusterProducerHFEM", edm::InputTag(""), true) and
0137                     edm::ParameterDescription<edm::InputTag>("pfClusterProducerHFHAD", edm::InputTag(""), true)));
0138   desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAllCalo"));
0139   desc.add<bool>("doRhoCorrection", false);
0140   desc.add<double>("rhoMax", 9.9999999E7);
0141   desc.add<double>("rhoScale", 1.0);
0142   desc.add<double>("drMax", 0.3);
0143   desc.add<double>("drVetoBarrel", 0.0);
0144   desc.add<double>("drVetoEndcap", 0.0);
0145   desc.add<double>("etaStripBarrel", 0.0);
0146   desc.add<double>("etaStripEndcap", 0.0);
0147   desc.add<double>("energyBarrel", 0.0);
0148   desc.add<double>("energyEndcap", 0.0);
0149   desc.add<bool>("useEt", true);
0150   desc.add<std::vector<double>>("effectiveAreas", {0.2, 0.25});   // 2016 post-ichep sinEle default
0151   desc.add<std::vector<double>>("absEtaLowEdges", {0.0, 1.479});  // Barrel, Endcap
0152   descriptions.add(defaultModuleLabel<HLTHcalPFClusterIsolationProducer<T1>>(), desc);
0153 }
0154 
0155 template <typename T1>
0156 void HLTHcalPFClusterIsolationProducer<T1>::produce(edm::StreamID sid,
0157                                                     edm::Event& iEvent,
0158                                                     const edm::EventSetup&) const {
0159   edm::Handle<double> rhoHandle;
0160   double rho = 0.0;
0161   if (doRhoCorrection_) {
0162     iEvent.getByToken(rhoProducer_, rhoHandle);
0163     rho = *(rhoHandle.product());
0164   }
0165 
0166   if (rho > rhoMax_)
0167     rho = rhoMax_;
0168 
0169   rho = rho * rhoScale_;
0170 
0171   edm::Handle<T1Collection> recoCandHandle;
0172 
0173   std::vector<edm::Handle<reco::PFClusterCollection>> clusterHandles;
0174   edm::Handle<reco::PFClusterCollection> clusterHcalHandle;
0175   edm::Handle<reco::PFClusterCollection> clusterHfemHandle;
0176   edm::Handle<reco::PFClusterCollection> clusterHfhadHandle;
0177 
0178   iEvent.getByToken(recoCandidateProducer_, recoCandHandle);
0179   iEvent.getByToken(pfClusterProducerHCAL_, clusterHcalHandle);
0180   //const reco::PFClusterCollection* forIsolationHcal = clusterHcalHandle.product();
0181   clusterHandles.push_back(clusterHcalHandle);
0182 
0183   if (useHF_) {
0184     iEvent.getByToken(pfClusterProducerHFEM_, clusterHfemHandle);
0185     clusterHandles.push_back(clusterHfemHandle);
0186     iEvent.getByToken(pfClusterProducerHFHAD_, clusterHfhadHandle);
0187     clusterHandles.push_back(clusterHfhadHandle);
0188   }
0189 
0190   T1IsolationMap recoCandMap(recoCandHandle);
0191   HcalPFClusterIsolation<T1> isoAlgo(
0192       drMax_, drVetoBarrel_, drVetoEndcap_, etaStripBarrel_, etaStripEndcap_, energyBarrel_, energyEndcap_, useEt_);
0193 
0194   for (unsigned int iReco = 0; iReco < recoCandHandle->size(); iReco++) {
0195     T1Ref candRef(recoCandHandle, iReco);
0196 
0197     float sum = isoAlgo.getSum(candRef, clusterHandles);
0198 
0199     if (doRhoCorrection_) {
0200       int iEA = -1;
0201       auto cEta = std::abs(candRef->eta());
0202       for (int bIt = absEtaLowEdges_.size() - 1; bIt > -1; bIt--) {
0203         if (cEta > absEtaLowEdges_.at(bIt)) {
0204           iEA = bIt;
0205           break;
0206         }
0207       }
0208 
0209       sum = sum - rho * effectiveAreas_.at(iEA);
0210     }
0211 
0212     recoCandMap.insert(candRef, sum);
0213   }
0214 
0215   iEvent.put(std::make_unique<T1IsolationMap>(recoCandMap));
0216 }
0217 
0218 typedef HLTHcalPFClusterIsolationProducer<reco::RecoEcalCandidate> EgammaHLTHcalPFClusterIsolationProducer;
0219 typedef HLTHcalPFClusterIsolationProducer<reco::RecoChargedCandidate> MuonHLTHcalPFClusterIsolationProducer;
0220 
0221 DEFINE_FWK_MODULE(EgammaHLTHcalPFClusterIsolationProducer);
0222 DEFINE_FWK_MODULE(MuonHLTHcalPFClusterIsolationProducer);