EgammaHLTProducers/plugins/HLTHGCalLayerClusterIsolationProducer.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/Math/interface/deltaR.h"
0015
0016 #include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
0017
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/stream/EDProducer.h"
0020
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025
0026 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0027 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h"
0028
0029 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0030 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateIsolation.h"
0031
0032 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0033 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0034
0035 #include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h"
0036
0037 template <typename T1>
0038 class HLTHGCalLayerClusterIsolationProducer : public edm::stream::EDProducer<> {
0039   typedef std::vector<T1> T1Collection;
0040   typedef edm::Ref<T1Collection> T1Ref;
0041   typedef edm::AssociationMap<edm::OneToValue<std::vector<T1>, float>> T1IsolationMap;
0042
0043 public:
0044   explicit HLTHGCalLayerClusterIsolationProducer(const edm::ParameterSet&);
0045   ~HLTHGCalLayerClusterIsolationProducer() override = default;
0046
0047   void produce(edm::Event&, const edm::EventSetup&) override;
0048   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0049
0050 private:
0051   edm::EDGetTokenT<T1Collection> recoCandidateProducer_;
0052   const edm::EDGetTokenT<reco::CaloClusterCollection> layerClusterProducer_;
0053   const edm::EDGetTokenT<double> rhoProducer_;
0054
0055   const double drMax_;
0056   const double drVetoEM_;
0057   const double drVetoHad_;
0058   const double minEnergyEM_;
0059   const double minEnergyHad_;
0060   const double minEtEM_;
0061   const double minEtHad_;
0062   const bool useEt_;
0063   const bool doRhoCorrection_;
0064   const double rhoMax_;
0065   const double rhoScale_;
0066   const std::vector<double> effectiveAreas_;
0067 };
0068
0069 template <typename T1>
0070 HLTHGCalLayerClusterIsolationProducer<T1>::HLTHGCalLayerClusterIsolationProducer(const edm::ParameterSet& config)
0071     : layerClusterProducer_(
0072           consumes<reco::CaloClusterCollection>(config.getParameter<edm::InputTag>("layerClusterProducer"))),
0073       rhoProducer_(consumes<double>(config.getParameter<edm::InputTag>("rhoProducer"))),
0074       drMax_(config.getParameter<double>("drMax")),
0075       drVetoEM_(config.getParameter<double>("drVetoEM")),
0076       drVetoHad_(config.getParameter<double>("drVetoHad")),
0077       minEnergyEM_(config.getParameter<double>("minEnergyEM")),
0078       minEnergyHad_(config.getParameter<double>("minEnergyHad")),
0079       minEtEM_(config.getParameter<double>("minEtEM")),
0080       minEtHad_(config.getParameter<double>("minEtHad")),
0081       useEt_(config.getParameter<bool>("useEt")),
0082       doRhoCorrection_(config.getParameter<bool>("doRhoCorrection")),
0083       rhoMax_(config.getParameter<double>("rhoMax")),
0084       rhoScale_(config.getParameter<double>("rhoScale")),
0085       effectiveAreas_(config.getParameter<std::vector<double>>("effectiveAreas")) {
0086   if (doRhoCorrection_) {
0087     if (effectiveAreas_.size() != 2)
0088       throw cms::Exception("IncompatibleVects")
0089           << "effectiveAreas should have two elements for em and had components. \n";
0090   }
0091
0092   std::string recoCandidateProducerName = "recoCandidateProducer";
0093   if ((typeid(HLTHGCalLayerClusterIsolationProducer<T1>) ==
0094        typeid(HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate>)))
0095     recoCandidateProducerName = "recoEcalCandidateProducer";
0096
0097   recoCandidateProducer_ = consumes<T1Collection>(config.getParameter<edm::InputTag>(recoCandidateProducerName));
0098   produces<T1IsolationMap>();
0099   produces<T1IsolationMap>("em");
0100   produces<T1IsolationMap>("had");
0101 }
0102
0103 template <typename T1>
0104 void HLTHGCalLayerClusterIsolationProducer<T1>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0105   std::string recoCandidateProducerName = "recoCandidateProducer";
0106   if ((typeid(HLTHGCalLayerClusterIsolationProducer<T1>) ==
0107        typeid(HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate>)))
0108     recoCandidateProducerName = "recoEcalCandidateProducer";
0109
0110   edm::ParameterSetDescription desc;
0111   desc.add<edm::InputTag>(recoCandidateProducerName, edm::InputTag("hltL1SeededRecoEcalCandidatePF"));
0112   desc.add<edm::InputTag>("layerClusterProducer", edm::InputTag("hltParticleFlowClusterECAL"));
0113   desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAllCalo"));
0114   desc.add<bool>("doRhoCorrection", false);
0115   desc.add<bool>("useEt", false);
0116   desc.add<double>("rhoMax", 9.9999999E7);
0117   desc.add<double>("rhoScale", 1.0);
0118   desc.add<double>("drMax", 0.3);
0119   desc.add<double>("drVetoEM", 0.0);
0120   desc.add<double>("drVetoHad", 0.0);
0121   desc.add<double>("minEnergyEM", 0.0);
0122   desc.add<double>("minEnergyHad", 0.0);
0123   desc.add<double>("minEtEM", 0.0);
0124   desc.add<double>("minEtHad", 0.0);
0125   desc.add<std::vector<double>>("effectiveAreas", {0.0, 0.0});  // for em and had components
0126   descriptions.add(defaultModuleLabel<HLTHGCalLayerClusterIsolationProducer<T1>>(), desc);
0127 }
0128
0129 template <typename T1>
0130 void HLTHGCalLayerClusterIsolationProducer<T1>::produce(edm::Event& iEvent, const edm::EventSetup&) {
0131   edm::Handle<double> rhoHandle;
0132   double rho = 0.0;
0133   if (doRhoCorrection_) {
0134     iEvent.getByToken(rhoProducer_, rhoHandle);
0135     rho = *(rhoHandle.product());
0136   }
0137
0138   rho = std::min(rho, rhoMax_);
0139   rho = rho * rhoScale_;
0140
0141   edm::Handle<T1Collection> recoCandHandle;
0142   edm::Handle<reco::CaloClusterCollection> clusterHandle;
0143
0144   iEvent.getByToken(recoCandidateProducer_, recoCandHandle);
0145   iEvent.getByToken(layerClusterProducer_, clusterHandle);
0146
0147   const std::vector<reco::CaloCluster> layerClusters = *(clusterHandle.product());
0148
0149   T1IsolationMap recoCandMap(recoCandHandle);
0150   T1IsolationMap recoCandMapEm(recoCandHandle);
0151   T1IsolationMap recoCandMapHad(recoCandHandle);
0152
0153   for (unsigned int iReco = 0; iReco < recoCandHandle->size(); iReco++) {
0154     T1Ref candRef(recoCandHandle, iReco);
0155
0156     float sumEm =
0157         HGCalClusterTools::emEnergyInCone(candRef->eta(),
0158                                           candRef->phi(),
0159                                           layerClusters,
0160                                           drVetoEM_,
0161                                           drMax_,
0162                                           minEtEM_,
0163                                           minEnergyEM_,
0164                                           useEt_ ? HGCalClusterTools::EType::ET : HGCalClusterTools::EType::ENERGY);
0165
0166     float sumHad =
0167         HGCalClusterTools::hadEnergyInCone(candRef->eta(),
0168                                            candRef->phi(),
0169                                            layerClusters,
0170                                            drVetoHad_,
0171                                            drMax_,
0172                                            minEtHad_,
0173                                            minEnergyHad_,
0174                                            useEt_ ? HGCalClusterTools::EType::ET : HGCalClusterTools::EType::ENERGY);
0175
0176     if (doRhoCorrection_) {
0177       sumEm = sumEm - rho * effectiveAreas_.at(0);
0178       sumHad = sumHad - rho * effectiveAreas_.at(1);
0179     }
0180
0181     float sum = sumEm + sumHad;
0182
0183     recoCandMap.insert(candRef, sum);
0184     recoCandMapEm.insert(candRef, sumEm);
0185     recoCandMapHad.insert(candRef, sumHad);
0186   }
0187
0188   iEvent.put(std::make_unique<T1IsolationMap>(recoCandMap));
0189   iEvent.put(std::make_unique<T1IsolationMap>(recoCandMapEm), "em");
0190   iEvent.put(std::make_unique<T1IsolationMap>(recoCandMapHad), "had");
0191 }
0192
0193 typedef HLTHGCalLayerClusterIsolationProducer<reco::RecoEcalCandidate> EgammaHLTHGCalLayerClusterIsolationProducer;
0194 typedef HLTHGCalLayerClusterIsolationProducer<reco::RecoChargedCandidate> MuonHLTHGCalLayerClusterIsolationProducer;
0195
0196 DEFINE_FWK_MODULE(EgammaHLTHGCalLayerClusterIsolationProducer);
0197 DEFINE_FWK_MODULE(MuonHLTHGCalLayerClusterIsolationProducer);