File indexing completed on 2024-04-06 12:24:06
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
0009
0010
0011
0012
0013
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
0042 chiso += pt;
0043 } else if (!fromPV && pt > ptthresh && dr2 > deadcone_pu * deadcone_pu) {
0044
0045 puiso += pt;
0046 }
0047 }
0048
0049 if (std::abs(id) == 130 && pt > ptthresh && dr2 > deadcone_nh * deadcone_nh)
0050 nhiso += pt;
0051
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
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 }