plugins/cuts/PhoAnyPFIsoWithEACut.cc

0001 #include "PhysicsTools/SelectorUtils/interface/CutApplicatorWithEventContentBase.h"
0002 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0003 #include "CommonTools/Egamma/interface/EffectiveAreas.h"
0004
0005 class PhoAnyPFIsoWithEACut : public CutApplicatorWithEventContentBase {
0006 public:
0007   PhoAnyPFIsoWithEACut(const edm::ParameterSet& c);
0008
0009   result_type operator()(const reco::PhotonPtr&) const final;
0010
0011   void setConsumes(edm::ConsumesCollector&) final;
0012   void getEventContent(const edm::EventBase&) final;
0013
0014   double value(const reco::CandidatePtr& cand) const final;
0015
0016   CandidateType candidateType() const final { return PHOTON; }
0017
0018 private:
0019   // Cut values
0020   float _C1_EB;
0021   float _C2_EB;
0022   float _C1_EE;
0023   float _C2_EE;
0024   // Configuration
0025   float _barrelCutOff;
0026   bool _useRelativeIso;
0027   // Effective area constants
0028   EffectiveAreas _effectiveAreas;
0029   // The isolations computed upstream
0030   edm::Handle<edm::ValueMap<float> > _anyPFIsoMap;
0031   // The rho
0032   edm::Handle<double> _rhoHandle;
0033
0034   constexpr static char anyPFIsoWithEA_[] = "anyPFIsoWithEA";
0035   constexpr static char rhoString_[] = "rho";
0036 };
0037
0038 constexpr char PhoAnyPFIsoWithEACut::anyPFIsoWithEA_[];
0039 constexpr char PhoAnyPFIsoWithEACut::rhoString_[];
0040
0041 DEFINE_EDM_PLUGIN(CutApplicatorFactory, PhoAnyPFIsoWithEACut, "PhoAnyPFIsoWithEACut");
0042
0043 PhoAnyPFIsoWithEACut::PhoAnyPFIsoWithEACut(const edm::ParameterSet& c)
0044     : CutApplicatorWithEventContentBase(c),
0045       _C1_EB(c.getParameter<double>("C1_EB")),
0046       _C2_EB(c.getParameter<double>("C2_EB")),
0047       _C1_EE(c.getParameter<double>("C1_EE")),
0048       _C2_EE(c.getParameter<double>("C2_EE")),
0049       _barrelCutOff(c.getParameter<double>("barrelCutOff")),
0050       _useRelativeIso(c.getParameter<bool>("useRelativeIso")),
0051       _effectiveAreas((c.getParameter<edm::FileInPath>("effAreasConfigFile")).fullPath()) {
0052   edm::InputTag maptag = c.getParameter<edm::InputTag>("anyPFIsoMap");
0053   contentTags_.emplace(anyPFIsoWithEA_, maptag);
0054
0055   edm::InputTag rhoTag = c.getParameter<edm::InputTag>("rho");
0056   contentTags_.emplace(rhoString_, rhoTag);
0057 }
0058
0059 void PhoAnyPFIsoWithEACut::setConsumes(edm::ConsumesCollector& cc) {
0060   auto anyPFIsoWithEA = cc.consumes<edm::ValueMap<float> >(contentTags_[anyPFIsoWithEA_]);
0061   contentTokens_.emplace(anyPFIsoWithEA_, anyPFIsoWithEA);
0062
0063   auto rho = cc.consumes<double>(contentTags_[rhoString_]);
0064   contentTokens_.emplace(rhoString_, rho);
0065 }
0066
0067 void PhoAnyPFIsoWithEACut::getEventContent(const edm::EventBase& ev) {
0068   ev.getByLabel(contentTags_[anyPFIsoWithEA_], _anyPFIsoMap);
0069   ev.getByLabel(contentTags_[rhoString_], _rhoHandle);
0070 }
0071
0072 CutApplicatorBase::result_type PhoAnyPFIsoWithEACut::operator()(const reco::PhotonPtr& cand) const {
0073   // in case we are by-value
0074   const std::string& inst_name = contentTags_.find(anyPFIsoWithEA_)->second.instance();
0075   edm::Ptr<pat::Photon> pat(cand);
0076   float anyisoval = -1.0;
0077   if (_anyPFIsoMap.isValid() && _anyPFIsoMap->contains(cand.id())) {
0078     anyisoval = (*_anyPFIsoMap)[cand];
0079   } else if (_anyPFIsoMap.isValid() && _anyPFIsoMap->idSize() == 1 && cand.id() == edm::ProductID()) {
0080     // in case we have spoofed a ptr
0081     //note this must be a 1:1 valuemap (only one product input)
0082     anyisoval = _anyPFIsoMap->begin()[cand.key()];
0083   } else if (_anyPFIsoMap.isValid()) {  // throw an exception
0084     anyisoval = (*_anyPFIsoMap)[cand];
0085   }
0086
0087   // Figure out the cut value
0088   // The value is generally pt-dependent: C1 + pt * C2
0089   const double absEta = std::abs(cand->superCluster()->eta());
0090   const float anyPFIsoWithEACutValue =
0091       (absEta < _barrelCutOff ? _C1_EB + cand->pt() * _C2_EB : _C1_EE + cand->pt() * _C2_EE);
0092
0093   // Retrieve the variable value for this particle
0094   float anyPFIso = _anyPFIsoMap.isValid() ? anyisoval : pat->userFloat(inst_name);
0095
0096   // Apply pile-up correction
0097   const double eA = _effectiveAreas.getEffectiveArea(absEta);
0098   const double rho = _rhoHandle.isValid() ? *_rhoHandle : 0;
0099   const float anyPFIsoWithEA = std::max(0.0, anyPFIso - rho * eA);
0100
0101   // Apply the cut and return the result
0102   // Scale by pT if the relative isolation is requested but avoid division by 0
0103   return anyPFIsoWithEA < anyPFIsoWithEACutValue * (_useRelativeIso ? cand->pt() : 1.);
0104 }
0105
0106 double PhoAnyPFIsoWithEACut::value(const reco::CandidatePtr& cand) const {
0107   reco::PhotonPtr pho(cand);
0108
0109   // in case we are by-value
0110   const std::string& inst_name = contentTags_.find(anyPFIsoWithEA_)->second.instance();
0111   edm::Ptr<pat::Photon> pat(cand);
0112   float anyisoval = -1.0;
0113   if (_anyPFIsoMap.isValid() && _anyPFIsoMap->contains(cand.id())) {
0114     anyisoval = (*_anyPFIsoMap)[cand];
0115   } else if (_anyPFIsoMap.isValid() && _anyPFIsoMap->idSize() == 1 && cand.id() == edm::ProductID()) {
0116     // in case we have spoofed a ptr
0117     //note this must be a 1:1 valuemap (only one product input)
0118     anyisoval = _anyPFIsoMap->begin()[cand.key()];
0119   } else if (_anyPFIsoMap.isValid()) {  // throw an exception
0120     anyisoval = (*_anyPFIsoMap)[cand];
0121   }
0122
0123   // Figure out the cut value
0124   // The value is generally pt-dependent: C1 + pt * C2
0125   double absEta = std::abs(pho->superCluster()->eta());
0126
0127   // Retrieve the variable value for this particle
0128   float anyPFIso = _anyPFIsoMap.isValid() ? anyisoval : pat->userFloat(inst_name);
0129
0130   // Apply pile-up correction
0131   double eA = _effectiveAreas.getEffectiveArea(absEta);
0132   double rho = *_rhoHandle;
0133   float anyPFIsoWithEA = std::max(0.0, anyPFIso - rho * eA);
0134
0135   // Divide by pT if the relative isolation is requested
0136   if (_useRelativeIso)
0137     anyPFIsoWithEA /= pho->pt();
0138
0139   // Apply the cut and return the result
0140   return anyPFIsoWithEA;
0141 }