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