Heppy/src/IsolationComputer.cc

0001 #include "PhysicsTools/Heppy/interface/IsolationComputer.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003
0004 #include <algorithm>
0005 namespace {
0006   struct ByEta {
0007     bool operator()(const pat::PackedCandidate *c1, const pat::PackedCandidate *c2) const {
0008       return c1->eta() < c2->eta();
0009     }
0010     bool operator()(float c1eta, const pat::PackedCandidate *c2) const { return c1eta < c2->eta(); }
0011     bool operator()(const pat::PackedCandidate *c1, float c2eta) const { return c1->eta() < c2eta; }
0012   };
0013 }  // namespace
0014 void heppy::IsolationComputer::setPackedCandidates(const std::vector<pat::PackedCandidate> &all,
0015                                                    int fromPV_thresh,
0016                                                    float dz_thresh,
0017                                                    float dxy_thresh,
0018                                                    bool also_leptons) {
0019   allcands_ = &all;
0020   charged_.clear();
0021   neutral_.clear();
0022   pileup_.clear();
0023
0024   for (const pat::PackedCandidate &p : all) {
0025     if (p.charge() == 0) {
0026       neutral_.push_back(&p);
0027     } else {
0028       if ((abs(p.pdgId()) == 211) || (also_leptons && ((abs(p.pdgId()) == 11) || (abs(p.pdgId()) == 13)))) {
0029         if (p.fromPV() > fromPV_thresh && fabs(p.dz()) < dz_thresh && fabs(p.dxy()) < dxy_thresh) {
0030           charged_.push_back(&p);
0031         } else {
0032           pileup_.push_back(&p);
0033         }
0034       }
0035     }
0036   }
0037   if (weightCone_ > 0)
0038     weights_.resize(neutral_.size());
0039   std::fill(weights_.begin(), weights_.end(), -1.f);
0040   std::sort(charged_.begin(), charged_.end(), ByEta());
0041   std::sort(neutral_.begin(), neutral_.end(), ByEta());
0042   std::sort(pileup_.begin(), pileup_.end(), ByEta());
0043   clearVetos();
0044 }
0045
0046 /// veto footprint from this candidate
0047 void heppy::IsolationComputer::addVetos(const reco::Candidate &cand) {
0048   for (unsigned int i = 0, n = cand.numberOfSourceCandidatePtrs(); i < n; ++i) {
0049     const reco::CandidatePtr &cp = cand.sourceCandidatePtr(i);
0050     if (cp.isNonnull() && cp.isAvailable())
0051       vetos_.push_back(&*cp);
0052   }
0053 }
0054
0055 /// clear all vetos
0056 void heppy::IsolationComputer::clearVetos() { vetos_.clear(); }
0057
0058 /// Isolation from charged from the PV
0059 float heppy::IsolationComputer::chargedAbsIso(
0060     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto) const {
0061   return isoSumRaw(charged_, cand, dR, innerR, threshold, selfVeto);
0062 }
0063
0064 /// Isolation from charged from PU
0065 float heppy::IsolationComputer::puAbsIso(
0066     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto) const {
0067   return isoSumRaw(pileup_, cand, dR, innerR, threshold, selfVeto);
0068 }
0069
0070 float heppy::IsolationComputer::neutralAbsIsoRaw(
0071     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto) const {
0072   return isoSumRaw(neutral_, cand, dR, innerR, threshold, selfVeto);
0073 }
0074
0075 float heppy::IsolationComputer::neutralAbsIsoWeighted(
0076     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto) const {
0077   return isoSumNeutralsWeighted(cand, dR, innerR, threshold, selfVeto);
0078 }
0079
0080 float heppy::IsolationComputer::neutralHadAbsIsoRaw(
0081     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto) const {
0082   return isoSumRaw(neutral_, cand, dR, innerR, threshold, selfVeto, 130);
0083 }
0084
0085 float heppy::IsolationComputer::neutralHadAbsIsoWeighted(
0086     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto) const {
0087   return isoSumNeutralsWeighted(cand, dR, innerR, threshold, selfVeto, 130);
0088 }
0089
0090 float heppy::IsolationComputer::photonAbsIsoRaw(
0091     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto) const {
0092   return isoSumRaw(neutral_, cand, dR, innerR, threshold, selfVeto, 22);
0093 }
0094
0095 float heppy::IsolationComputer::photonAbsIsoWeighted(
0096     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto) const {
0097   return isoSumNeutralsWeighted(cand, dR, innerR, threshold, selfVeto, 22);
0098 }
0099
0100 float heppy::IsolationComputer::isoSumRaw(const std::vector<const pat::PackedCandidate *> &cands,
0101                                           const reco::Candidate &cand,
0102                                           float dR,
0103                                           float innerR,
0104                                           float threshold,
0105                                           SelfVetoPolicy selfVeto,
0106                                           int pdgId) const {
0107   float dR2 = dR * dR, innerR2 = innerR * innerR;
0108
0109   std::vector<const reco::Candidate *> vetos(vetos_);
0110   for (unsigned int i = 0, n = cand.numberOfSourceCandidatePtrs(); i < n; ++i) {
0111     if (selfVeto == selfVetoNone)
0112       break;
0113     const reco::CandidatePtr &cp = cand.sourceCandidatePtr(i);
0114     if (cp.isNonnull() && cp.isAvailable()) {
0115       vetos.push_back(&*cp);
0116       if (selfVeto == selfVetoFirst)
0117         break;
0118     }
0119   }
0120
0121   typedef std::vector<const pat::PackedCandidate *>::const_iterator IT;
0122   IT candsbegin = std::lower_bound(cands.begin(), cands.end(), cand.eta() - dR, ByEta());
0123   IT candsend = std::upper_bound(candsbegin, cands.end(), cand.eta() + dR, ByEta());
0124
0125   double isosum = 0;
0126   for (IT icharged = candsbegin; icharged < candsend; ++icharged) {
0127     // pdgId
0128     if (pdgId > 0 && abs((*icharged)->pdgId()) != pdgId)
0129       continue;
0130     // threshold
0131     if (threshold > 0 && (*icharged)->pt() < threshold)
0132       continue;
0133     // cone
0134     float mydr2 = reco::deltaR2(**icharged, cand);
0135     if (mydr2 >= dR2 || mydr2 < innerR2)
0136       continue;
0137     // veto
0138     if (std::find(vetos.begin(), vetos.end(), *icharged) != vetos.end()) {
0139       continue;
0140     }
0141     // add to sum
0142     isosum += (*icharged)->pt();
0143   }
0144   return isosum;
0145 }
0146
0147 float heppy::IsolationComputer::isoSumNeutralsWeighted(
0148     const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto, int pdgId) const {
0149   if (weightCone_ <= 0)
0150     throw cms::Exception("LogicError", "you must set a valid weight cone to use this method");
0151   float dR2 = dR * dR, innerR2 = innerR * innerR, weightCone2 = weightCone_ * weightCone_;
0152
0153   std::vector<const reco::Candidate *> vetos(vetos_);
0154   for (unsigned int i = 0, n = cand.numberOfSourceCandidatePtrs(); i < n; ++i) {
0155     if (selfVeto == selfVetoNone)
0156       break;
0157     const reco::CandidatePtr &cp = cand.sourceCandidatePtr(i);
0158     if (cp.isNonnull() && cp.isAvailable()) {
0159       vetos.push_back(&*cp);
0160       if (selfVeto == selfVetoFirst)
0161         break;
0162     }
0163   }
0164
0165   typedef std::vector<const pat::PackedCandidate *>::const_iterator IT;
0166   IT charged_begin = std::lower_bound(charged_.begin(), charged_.end(), cand.eta() - dR - weightCone_, ByEta());
0167   IT charged_end = std::upper_bound(charged_begin, charged_.end(), cand.eta() + dR + weightCone_, ByEta());
0168   IT pileup_begin = std::lower_bound(pileup_.begin(), pileup_.end(), cand.eta() - dR - weightCone_, ByEta());
0169   IT pileup_end = std::upper_bound(pileup_begin, pileup_.end(), cand.eta() + dR + weightCone_, ByEta());
0170   IT neutral_begin = std::lower_bound(neutral_.begin(), neutral_.end(), cand.eta() - dR, ByEta());
0171   IT neutral_end = std::upper_bound(neutral_begin, neutral_.end(), cand.eta() + dR, ByEta());
0172
0173   double isosum = 0.0;
0174   for (IT ineutral = neutral_begin; ineutral < neutral_end; ++ineutral) {
0175     // pdgId
0176     if (pdgId > 0 && abs((*ineutral)->pdgId()) != pdgId)
0177       continue;
0178     // threshold
0179     if (threshold > 0 && (*ineutral)->pt() < threshold)
0180       continue;
0181     // cone
0182     float mydr2 = reco::deltaR2(**ineutral, cand);
0183     if (mydr2 >= dR2 || mydr2 < innerR2)
0184       continue;
0185     // veto
0186     if (std::find(vetos.begin(), vetos.end(), *ineutral) != vetos.end()) {
0187       continue;
0188     }
0189     // weight
0190     float &w = weights_[ineutral - neutral_.begin()];
0191     if (w == -1.f) {
0192       double sumc = 0, sump = 0.0;
0193       for (IT icharged = charged_begin; icharged < charged_end; ++icharged) {
0194         float hisdr2 = std::max<float>(reco::deltaR2(**icharged, **ineutral), 0.01f);
0195         if (hisdr2 > weightCone2)
0196           continue;
0197         if (std::find(vetos_.begin(), vetos_.end(), *icharged) != vetos_.end()) {
0198           continue;
0199         }
0200         sumc += std::log((*icharged)->pt() / std::sqrt(hisdr2));
0201       }
0202       for (IT ipileup = pileup_begin; ipileup < pileup_end; ++ipileup) {
0203         float hisdr2 = std::max<float>(reco::deltaR2(**ipileup, **ineutral), 0.01f);
0204         if (hisdr2 > weightCone2)
0205           continue;
0206         if (std::find(vetos_.begin(), vetos_.end(), *ipileup) != vetos_.end()) {
0207           continue;
0208         }
0209         sumc += std::log((*ipileup)->pt() / std::sqrt(hisdr2));
0210       }
0211       w = (sump == 0 ? 1 : sumc / (sump + sumc));
0212     }
0213     // add to sum
0214     isosum += w * (*ineutral)->pt();
0215   }
0216   return isosum;
0217 }