PatAlgos/src/CalculatePtRatioRel.cc

0001 #include "PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h"
0002
0003 #include "DataFormats/JetReco/interface/PFJet.h"
0004 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0005 #include "DataFormats/PatCandidates/interface/Muon.h"
0006 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0007
0008 using namespace pat;
0009
0010 CalculatePtRatioRel::CalculatePtRatioRel(float dR2max) : dR2max_(dR2max) {}
0011
0012 CalculatePtRatioRel::~CalculatePtRatioRel() {}
0013
0014 namespace {
0015
0016   float ptRel(const reco::Candidate::LorentzVector& muP4,
0017               const reco::Candidate::LorentzVector& jetP4,
0018               bool subtractMuon = true) {
0019     reco::Candidate::LorentzVector jp4 = jetP4;
0020     if (subtractMuon)
0021       jp4 -= muP4;
0022     float dot = muP4.Vect().Dot(jp4.Vect());
0023     float ptrel = muP4.P2() - dot * dot / jp4.P2();
0024     ptrel = ptrel > 0 ? sqrt(ptrel) : 0.0;
0025     return ptrel;
0026   }
0027 }  // namespace
0028
0029 std::array<float, 2> CalculatePtRatioRel::computePtRatioRel(const pat::Muon& muon,
0030                                                             const reco::JetTagCollection& bTags,
0031                                                             const reco::JetCorrector* correctorL1,
0032                                                             const reco::JetCorrector* correctorL1L2L3Res) const {
0033   //Initialise loop variables
0034   float minDr2 = 9999;
0035   double jecL1L2L3Res = 1.;
0036   double jecL1 = 1.;
0037
0038   // Compute corrected isolation variables
0039   const float chIso = muon.pfIsolationR04().sumChargedHadronPt;
0040   const float nIso = muon.pfIsolationR04().sumNeutralHadronEt;
0041   const float phoIso = muon.pfIsolationR04().sumPhotonEt;
0042   const float puIso = muon.pfIsolationR04().sumPUPt;
0043   const float dbCorrectedIsolation = chIso + std::max(nIso + phoIso - .5f * puIso, 0.f);
0044   const float dbCorrectedRelIso = dbCorrectedIsolation / muon.pt();
0045
0046   float JetPtRatio = 1. / (1 + dbCorrectedRelIso);
0047   float JetPtRel = 0.;
0048
0049   const reco::Candidate::LorentzVector& muP4(muon.p4());
0050   for (const auto& tagI : bTags) {
0051     // for each muon with the lepton
0052     float dr2 = deltaR2(*(tagI.first), muon);
0053     if (dr2 > minDr2)
0054       continue;
0055     minDr2 = dr2;
0056
0057     reco::Candidate::LorentzVector jetP4(tagI.first->p4());
0058
0059     if (correctorL1 && correctorL1L2L3Res) {
0060       jecL1L2L3Res = correctorL1L2L3Res->correction(*(tagI.first));
0061       jecL1 = correctorL1->correction(*(tagI.first));
0062     }
0063
0064     if (minDr2 < dR2max_) {
0065       if ((jetP4 - muP4).Rho() < 0.0001) {
0066         JetPtRel = 0.;
0067         JetPtRatio = 1.;
0068       } else {
0069         jetP4 -= muP4 / jecL1;
0070         jetP4 *= jecL1L2L3Res;
0071         jetP4 += muP4;
0072
0073         JetPtRatio = muP4.pt() / jetP4.pt();
0074         JetPtRel = ptRel(muP4, jetP4);
0075       }
0076     }
0077   }
0078
0079   if (JetPtRatio > 1.5)
0080     JetPtRatio = 1.5;
0081
0082   std::array<float, 2> outputs = {{JetPtRatio, JetPtRel}};
0083   return outputs;
0084 };