plugins/cuts/PhoAnyPFIsoWithEAAndExpoScalingCut.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 PhoAnyPFIsoWithEAAndExpoScalingCut : public CutApplicatorWithEventContentBase {
0006 public:
0007   PhoAnyPFIsoWithEAAndExpoScalingCut(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 _C3_EB;
0023   float _C1_EE;
0024   float _C2_EE;
0025   float _C3_EE;
0026   // Configuration
0027   float _barrelCutOff;
0028   bool _useRelativeIso;
0029   // Effective area constants
0030   EffectiveAreas _effectiveAreas;
0031   // The isolations computed upstream
0032   edm::Handle<edm::ValueMap<float> > _anyPFIsoMap;
0033   // The rho
0034   edm::Handle<double> _rhoHandle;
0035
0036   constexpr static char anyPFIsoWithEA_[] = "anyPFIsoWithEA";
0037   constexpr static char rhoString_[] = "rho";
0038 };
0039
0040 constexpr char PhoAnyPFIsoWithEAAndExpoScalingCut::anyPFIsoWithEA_[];
0041 constexpr char PhoAnyPFIsoWithEAAndExpoScalingCut::rhoString_[];
0042
0043 DEFINE_EDM_PLUGIN(CutApplicatorFactory, PhoAnyPFIsoWithEAAndExpoScalingCut, "PhoAnyPFIsoWithEAAndExpoScalingCut");
0044
0045 PhoAnyPFIsoWithEAAndExpoScalingCut::PhoAnyPFIsoWithEAAndExpoScalingCut(const edm::ParameterSet& c)
0046     : CutApplicatorWithEventContentBase(c),
0047       _C1_EB(c.getParameter<double>("C1_EB")),
0048       _C2_EB(c.getParameter<double>("C2_EB")),
0049       _C3_EB(c.getParameter<double>("C3_EB")),
0050       _C1_EE(c.getParameter<double>("C1_EE")),
0051       _C2_EE(c.getParameter<double>("C2_EE")),
0052       _C3_EE(c.getParameter<double>("C3_EE")),
0053       _barrelCutOff(c.getParameter<double>("barrelCutOff")),
0054       _useRelativeIso(c.getParameter<bool>("useRelativeIso")),
0055       _effectiveAreas((c.getParameter<edm::FileInPath>("effAreasConfigFile")).fullPath()) {
0056   edm::InputTag maptag = c.getParameter<edm::InputTag>("anyPFIsoMap");
0057   contentTags_.emplace(anyPFIsoWithEA_, maptag);
0058
0059   edm::InputTag rhoTag = c.getParameter<edm::InputTag>("rho");
0060   contentTags_.emplace(rhoString_, rhoTag);
0061 }
0062
0063 void PhoAnyPFIsoWithEAAndExpoScalingCut::setConsumes(edm::ConsumesCollector& cc) {
0064   auto anyPFIsoWithEA = cc.consumes<edm::ValueMap<float> >(contentTags_[anyPFIsoWithEA_]);
0065   contentTokens_.emplace(anyPFIsoWithEA_, anyPFIsoWithEA);
0066
0067   auto rho = cc.consumes<double>(contentTags_[rhoString_]);
0068   contentTokens_.emplace(rhoString_, rho);
0069 }
0070
0071 void PhoAnyPFIsoWithEAAndExpoScalingCut::getEventContent(const edm::EventBase& ev) {
0072   ev.getByLabel(contentTags_[anyPFIsoWithEA_], _anyPFIsoMap);
0073   ev.getByLabel(contentTags_[rhoString_], _rhoHandle);
0074 }
0075
0076 CutApplicatorBase::result_type PhoAnyPFIsoWithEAAndExpoScalingCut::operator()(const reco::PhotonPtr& cand) const {
0077   // in case we are by-value
0078   const std::string& inst_name = contentTags_.find(anyPFIsoWithEA_)->second.instance();
0079   edm::Ptr<pat::Photon> pat(cand);
0080   float anyisoval = -1.0;
0081   if (_anyPFIsoMap.isValid() && _anyPFIsoMap->contains(cand.id())) {
0082     anyisoval = (*_anyPFIsoMap)[cand];
0083   } else if (_anyPFIsoMap.isValid() && _anyPFIsoMap->idSize() == 1 && cand.id() == edm::ProductID()) {
0084     // in case we have spoofed a ptr
0085     //note this must be a 1:1 valuemap (only one product input)
0086     anyisoval = _anyPFIsoMap->begin()[cand.key()];
0087   } else if (_anyPFIsoMap.isValid()) {  // throw an exception
0088     anyisoval = (*_anyPFIsoMap)[cand];
0089   }
0090
0091   // Figure out the cut value
0092   // The value is generally pt-dependent: C1 + pt * C2
0093   const float pt = cand->pt();
0094
0095   // In this version of the isolation cut we apply
0096   // exponential pt scaling to the barrel isolation cut,
0097   // and linear pt scaling to the endcap isolation cut.
0098   double absEta = std::abs(cand->superCluster()->eta());
0099   const float isolationCutValue =
0100       (absEta < _barrelCutOff ? _C1_EB + exp(pt * _C2_EB + _C3_EB) : _C1_EE + exp(pt * _C2_EE + _C3_EE));
0101
0102   // Retrieve the variable value for this particle
0103   float anyPFIso = _anyPFIsoMap.isValid() ? anyisoval : pat->userFloat(inst_name);
0104
0105   // Apply pile-up correction
0106   double eA = _effectiveAreas.getEffectiveArea(absEta);
0107   double rho = *_rhoHandle;
0108   float anyPFIsoWithEA = std::max(0.0, anyPFIso - rho * eA);
0109
0110   // Divide by pT if the relative isolation is requested
0111   if (_useRelativeIso)
0112     anyPFIsoWithEA /= pt;
0113
0114   // Apply the cut and return the result
0115   return anyPFIsoWithEA < isolationCutValue;
0116 }
0117
0118 double PhoAnyPFIsoWithEAAndExpoScalingCut::value(const reco::CandidatePtr& cand) const {
0119   reco::PhotonPtr pho(cand);
0120
0121   // in case we are by-value
0122   const std::string& inst_name = contentTags_.find(anyPFIsoWithEA_)->second.instance();
0123   edm::Ptr<pat::Photon> pat(cand);
0124   float anyisoval = -1.0;
0125   if (_anyPFIsoMap.isValid() && _anyPFIsoMap->contains(cand.id())) {
0126     anyisoval = (*_anyPFIsoMap)[cand];
0127   } else if (_anyPFIsoMap.isValid() && _anyPFIsoMap->idSize() == 1 && cand.id() == edm::ProductID()) {
0128     // in case we have spoofed a ptr
0129     //note this must be a 1:1 valuemap (only one product input)
0130     anyisoval = _anyPFIsoMap->begin()[cand.key()];
0131   } else if (_anyPFIsoMap.isValid()) {  // throw an exception
0132     anyisoval = (*_anyPFIsoMap)[cand];
0133   }
0134
0135   // Figure out the cut value
0136   // The value is generally pt-dependent: C1 + pt * C2
0137   const float pt = pho->pt();
0138
0139   // In this version of the isolation cut we apply
0140   // exponential pt scaling to the barrel isolation cut,
0141   // and linear pt scaling to the endcap isolation cut.
0142   double absEta = std::abs(pho->superCluster()->eta());
0143
0144   // Retrieve the variable value for this particle
0145   float anyPFIso = _anyPFIsoMap.isValid() ? anyisoval : pat->userFloat(inst_name);
0146
0147   // Apply pile-up correction
0148   double eA = _effectiveAreas.getEffectiveArea(absEta);
0149   double rho = *_rhoHandle;
0150   float anyPFIsoWithEA = std::max(0.0, anyPFIso - rho * eA);
0151
0152   // Divide by pT if the relative isolation is requested
0153   if (_useRelativeIso)
0154     anyPFIsoWithEA /= pt;
0155
0156   // Apply the cut and return the result
0157   return anyPFIsoWithEA;
0158 }