EgammaIsolationAlgos/plugins/EgammaEcalRecHitIsolationProducer.cc

0001 //*****************************************************************************
0002 // File:      EgammaEcalRecHitIsolationProducer.cc
0003 // ----------------------------------------------------------------------------
0004 // OrigAuth:  Matthias Mozer, adapted from EgammaHcalIsolationProducer by S. Harper
0005 // Institute: IIHE-VUB, RAL
0006 //=============================================================================
0007 //*****************************************************************************
0008
0009 #include "DataFormats/Candidate/interface/Candidate.h"
0010 #include "DataFormats/Candidate/interface/CandAssociation.h"
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/DetId/interface/DetId.h"
0013 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0014 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0015 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0016 #include "FWCore/Framework/interface/ConsumesCollector.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/global/EDProducer.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/Utilities/interface/ESGetToken.h"
0023 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
0026 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0027 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0028
0029 class EgammaEcalRecHitIsolationProducer : public edm::global::EDProducer<> {
0030 public:
0031   explicit EgammaEcalRecHitIsolationProducer(const edm::ParameterSet&);
0032
0033   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0034
0035 private:
0036   const edm::EDGetTokenT<edm::View<reco::Candidate>> emObjectProducer_;
0037   const edm::EDGetTokenT<EcalRecHitCollection> ecalBarrelRecHitCollection_;
0038   const edm::EDGetTokenT<EcalRecHitCollection> ecalEndcapRecHitCollection_;
0039
0040   double egIsoPtMinBarrel_;       //minimum Et noise cut
0041   double egIsoEMinBarrel_;        //minimum E noise cut
0042   double egIsoPtMinEndcap_;       //minimum Et noise cut
0043   double egIsoEMinEndcap_;        //minimum E noise cut
0044   double egIsoConeSizeOut_;       //outer cone size
0045   double egIsoConeSizeInBarrel_;  //inner cone size
0046   double egIsoConeSizeInEndcap_;  //inner cone size
0047   double egIsoJurassicWidth_;     // exclusion strip width for jurassic veto
0048
0049   bool useIsolEt_;  //switch for isolEt rather than isolE
0050   bool tryBoth_;    // use rechits from barrel + endcap
0051   bool subtract_;   // subtract SC energy (allows veto cone of zero size)
0052
0053   bool useNumCrystals_;  // veto on number of crystals
0054   bool vetoClustered_;   // veto all clusterd rechits
0055
0056   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> sevLvToken_;
0057   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometrytoken_;
0058 };
0059
0060 #include "FWCore/Framework/interface/MakerMacros.h"
0061 DEFINE_FWK_MODULE(EgammaEcalRecHitIsolationProducer);
0062
0063 EgammaEcalRecHitIsolationProducer::EgammaEcalRecHitIsolationProducer(const edm::ParameterSet& config)
0064     //inputs
0065     : emObjectProducer_{consumes(config.getParameter<edm::InputTag>("emObjectProducer"))},
0066       ecalBarrelRecHitCollection_{consumes(config.getParameter<edm::InputTag>("ecalBarrelRecHitCollection"))},
0067       ecalEndcapRecHitCollection_{consumes(config.getParameter<edm::InputTag>("ecalEndcapRecHitCollection"))} {
0068   //vetos
0069   egIsoPtMinBarrel_ = config.getParameter<double>("etMinBarrel");
0070   egIsoEMinBarrel_ = config.getParameter<double>("eMinBarrel");
0071   egIsoPtMinEndcap_ = config.getParameter<double>("etMinEndcap");
0072   egIsoEMinEndcap_ = config.getParameter<double>("eMinEndcap");
0073   egIsoConeSizeInBarrel_ = config.getParameter<double>("intRadiusBarrel");
0074   egIsoConeSizeInEndcap_ = config.getParameter<double>("intRadiusEndcap");
0075   egIsoConeSizeOut_ = config.getParameter<double>("extRadius");
0076   egIsoJurassicWidth_ = config.getParameter<double>("jurassicWidth");
0077
0078   // options
0079   useIsolEt_ = config.getParameter<bool>("useIsolEt");
0080   tryBoth_ = config.getParameter<bool>("tryBoth");
0081   subtract_ = config.getParameter<bool>("subtract");
0082   useNumCrystals_ = config.getParameter<bool>("useNumCrystals");
0083   vetoClustered_ = config.getParameter<bool>("vetoClustered");
0084
0085   //EventSetup Tokens
0086   sevLvToken_ = esConsumes();
0087   caloGeometrytoken_ = esConsumes();
0088
0089   //register your products
0090   produces<edm::ValueMap<double>>();
0091 }
0092
0093 // ------------ method called to produce the data  ------------
0094 void EgammaEcalRecHitIsolationProducer::produce(edm::StreamID,
0095                                                 edm::Event& iEvent,
0096                                                 const edm::EventSetup& iSetup) const {
0097   // Get the  filtered objects
0098   auto emObjectHandle = iEvent.getHandle(emObjectProducer_);
0099
0100   // Next get Ecal hits barrel
0101   auto ecalBarrelRecHitHandle = iEvent.getHandle(ecalBarrelRecHitCollection_);
0102
0103   // Next get Ecal hits endcap
0104   auto ecalEndcapRecHitHandle = iEvent.getHandle(ecalEndcapRecHitCollection_);
0105
0106   edm::ESHandle<EcalSeverityLevelAlgo> sevlv = iSetup.getHandle(sevLvToken_);
0107   const EcalSeverityLevelAlgo* sevLevel = sevlv.product();
0108
0109   //Get Calo Geometry
0110   edm::ESHandle<CaloGeometry> pG = iSetup.getHandle(caloGeometrytoken_);
0111   const CaloGeometry* caloGeom = pG.product();
0112
0113   //reco::CandViewDoubleAssociations* isoMap = new reco::CandViewDoubleAssociations( reco::CandidateBaseRefProd( emObjectHandle ) );
0114   auto isoMap = std::make_unique<edm::ValueMap<double>>();
0115   edm::ValueMap<double>::Filler filler(*isoMap);
0116   std::vector<double> retV(emObjectHandle->size(), 0);
0117
0118   EgammaRecHitIsolation ecalBarrelIsol(egIsoConeSizeOut_,
0119                                        egIsoConeSizeInBarrel_,
0120                                        egIsoJurassicWidth_,
0121                                        egIsoPtMinBarrel_,
0122                                        egIsoEMinBarrel_,
0123                                        caloGeom,
0124                                        *ecalBarrelRecHitHandle,
0125                                        sevLevel,
0126                                        DetId::Ecal);
0127   ecalBarrelIsol.setUseNumCrystals(useNumCrystals_);
0128   ecalBarrelIsol.setVetoClustered(vetoClustered_);
0129
0130   EgammaRecHitIsolation ecalEndcapIsol(egIsoConeSizeOut_,
0131                                        egIsoConeSizeInEndcap_,
0132                                        egIsoJurassicWidth_,
0133                                        egIsoPtMinEndcap_,
0134                                        egIsoEMinEndcap_,
0135                                        caloGeom,
0136                                        *ecalEndcapRecHitHandle,
0137                                        sevLevel,
0138                                        DetId::Ecal);
0139   ecalEndcapIsol.setUseNumCrystals(useNumCrystals_);
0140   ecalEndcapIsol.setVetoClustered(vetoClustered_);
0141
0142   for (size_t i = 0; i < emObjectHandle->size(); ++i) {
0143     //i need to know if its in the barrel/endcap so I get the supercluster handle to find out the detector eta
0144     //this might not be the best way, are we guaranteed that eta<1.5 is barrel
0145     //this can be safely replaced by another method which determines where the emobject is
0146     //then we either get the isolation Et or isolation Energy depending on user selection
0147     double isoValue = 0.;
0148
0149     reco::SuperClusterRef superClus = emObjectHandle->at(i).get<reco::SuperClusterRef>();
0150
0151     if (tryBoth_) {  //barrel + endcap
0152       if (useIsolEt_)
0153         isoValue =
0154             ecalBarrelIsol.getEtSum(&(emObjectHandle->at(i))) + ecalEndcapIsol.getEtSum(&(emObjectHandle->at(i)));
0155       else
0156         isoValue = ecalBarrelIsol.getEnergySum(&(emObjectHandle->at(i))) +
0157                    ecalEndcapIsol.getEnergySum(&(emObjectHandle->at(i)));
0158     } else if (fabs(superClus->eta()) < 1.479) {  //barrel
0159       if (useIsolEt_)
0160         isoValue = ecalBarrelIsol.getEtSum(&(emObjectHandle->at(i)));
0161       else
0162         isoValue = ecalBarrelIsol.getEnergySum(&(emObjectHandle->at(i)));
0163     } else {  //endcap
0164       if (useIsolEt_)
0165         isoValue = ecalEndcapIsol.getEtSum(&(emObjectHandle->at(i)));
0166       else
0167         isoValue = ecalEndcapIsol.getEnergySum(&(emObjectHandle->at(i)));
0168     }
0169
0170     //we subtract off the electron energy here as well
0171     double subtractVal = 0;
0172
0173     if (useIsolEt_)
0174       subtractVal = superClus.get()->rawEnergy() * sin(2 * atan(exp(-superClus.get()->eta())));
0175     else
0176       subtractVal = superClus.get()->rawEnergy();
0177
0178     if (subtract_)
0179       isoValue -= subtractVal;
0180
0181     retV[i] = isoValue;
0182     //all done, isolation is now in the map
0183
0184   }  //end of loop over em objects
0185
0186   filler.insert(emObjectHandle, retV.begin(), retV.end());
0187   filler.fill();
0188
0189   iEvent.put(std::move(isoMap));
0190 }
0191
0192 //define this as a plug-in
0193 //DEFINE_FWK_MODULE(EgammaRecHitIsolation,Producer);