PatUtils/src/MiniIsolation.cc

0001 #include "DataFormats/Candidate/interface/Candidate.h"
0002 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004 #include "PhysicsTools/PatUtils/interface/MiniIsolation.h"
0005
0006 namespace pat {
0007
0008   // Computes MiniIsolation given a 4-vector of the object in question
0009   // and a collection of packed PF cands. Computed as a sum of PFCand pT
0010   // inside a cone of radius dR = max(mindir, min(maxdr, kt_scale/pT))
0011   // Excludes PFCands inside of "deadcone" radius.
0012   // For nh, ph, pu, only include particles with pT > ptthresh
0013   // Some documentation can be found here: https://hypernews.cern.ch/HyperNews/CMS/get/susy/1991.html
0014
0015   float miniIsoDr(const reco::Candidate::PolarLorentzVector &p4, float mindr, float maxdr, float kt_scale) {
0016     return std::max(mindr, std::min(maxdr, float(kt_scale / p4.pt())));
0017   }
0018
0019   PFIsolation getMiniPFIsolation(const pat::PackedCandidateCollection *pfcands,
0020                                  const reco::Candidate::PolarLorentzVector &p4,
0021                                  float mindr,
0022                                  float maxdr,
0023                                  float kt_scale,
0024                                  float ptthresh,
0025                                  float deadcone_ch,
0026                                  float deadcone_pu,
0027                                  float deadcone_ph,
0028                                  float deadcone_nh,
0029                                  float dZ_cut) {
0030     float chiso = 0, nhiso = 0, phiso = 0, puiso = 0;
0031     float drcut = miniIsoDr(p4, mindr, maxdr, kt_scale);
0032     for (auto const &pc : *pfcands) {
0033       float dr2 = deltaR2(p4, pc);
0034       if (dr2 > drcut * drcut)
0035         continue;
0036       float pt = pc.p4().pt();
0037       int id = pc.pdgId();
0038       if (std::abs(id) == 211) {
0039         bool fromPV = (pc.fromPV() > 1 || fabs(pc.dz()) < dZ_cut);
0040         if (fromPV && dr2 > deadcone_ch * deadcone_ch) {
0041           // if charged hadron and from primary vertex, add to charged hadron isolation
0042           chiso += pt;
0043         } else if (!fromPV && pt > ptthresh && dr2 > deadcone_pu * deadcone_pu) {
0044           // if charged hadron and NOT from primary vertex, add to pileup isolation
0045           puiso += pt;
0046         }
0047       }
0048       // if neutral hadron, add to neutral hadron isolation
0049       if (std::abs(id) == 130 && pt > ptthresh && dr2 > deadcone_nh * deadcone_nh)
0050         nhiso += pt;
0051       // if photon, add to photon isolation
0052       if (std::abs(id) == 22 && pt > ptthresh && dr2 > deadcone_ph * deadcone_ph)
0053         phiso += pt;
0054     }
0055
0056     return pat::PFIsolation(chiso, nhiso, phiso, puiso);
0057   }
0058
0059   double muonRelMiniIsoPUCorrected(const PFIsolation &iso,
0060                                    const reco::Candidate::PolarLorentzVector &p4,
0061                                    double dr,
0062                                    double rho,
0063                                    const std::vector<double> &area) {
0064     double absEta = std::abs(p4.eta());
0065     double ea = 0;
0066     //Eta dependent effective area
0067     if (absEta < 0.800)
0068       ea = area.at(0);
0069     else if (absEta < 1.300)
0070       ea = area.at(1);
0071     else if (absEta < 2.000)
0072       ea = area.at(2);
0073     else if (absEta < 2.200)
0074       ea = area.at(3);
0075     else if (absEta < 2.500)
0076       ea = area.at(4);
0077
0078     double correction = rho * ea * (dr / 0.3) * (dr / 0.3);
0079     double correctedIso = iso.chargedHadronIso() + std::max(0.0, iso.neutralHadronIso() + iso.photonIso() - correction);
0080     return correctedIso / p4.pt();
0081   }
0082
0083 }  // namespace pat