Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:03

0001 #include "PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h"
0002 
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0005 #include "DataFormats/Candidate/interface/Candidate.h"
0006 #include "DataFormats/MuonReco/interface/Muon.h"
0007 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0008 #include "DataFormats/VertexReco/interface/Vertex.h"
0009 #include "DataFormats/JetReco/interface/PFJet.h"
0010 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0011 #include "DataFormats/PatCandidates/interface/Muon.h"
0012 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0013 
0014 using namespace pat;
0015 
0016 MuonMvaEstimator::MuonMvaEstimator(const edm::FileInPath& weightsfile, float dRmax) : dRmax_(dRmax) {
0017   gbrForest_ = createGBRForest(weightsfile);
0018 }
0019 
0020 MuonMvaEstimator::~MuonMvaEstimator() {}
0021 
0022 namespace {
0023 
0024   enum inputIndexes {
0025     kPt,
0026     kEta,
0027     kJetNDauCharged,
0028     kMiniRelIsoCharged,
0029     kMiniRelIsoNeutral,
0030     kJetPtRel,
0031     kJetBTagCSV,
0032     kJetPtRatio,
0033     kLog_abs_dxyBS,
0034     kSip,
0035     kLog_abs_dzPV,
0036     kSegmentCompatibility,
0037     kLast
0038   };
0039 
0040   float ptRel(const reco::Candidate::LorentzVector& muP4,
0041               const reco::Candidate::LorentzVector& jetP4,
0042               bool subtractMuon = true) {
0043     reco::Candidate::LorentzVector jp4 = jetP4;
0044     if (subtractMuon)
0045       jp4 -= muP4;
0046     float dot = muP4.Vect().Dot(jp4.Vect());
0047     float ptrel = muP4.P2() - dot * dot / jp4.P2();
0048     ptrel = ptrel > 0 ? sqrt(ptrel) : 0.0;
0049     return ptrel;
0050   }
0051 }  // namespace
0052 
0053 float MuonMvaEstimator::computeMva(const pat::Muon& muon,
0054                                    const reco::Vertex& vertex,
0055                                    const reco::JetTagCollection& bTags,
0056                                    float& jetPtRatio,
0057                                    float& jetPtRel,
0058                                    float& miniIsoValue,
0059                                    const reco::JetCorrector* correctorL1,
0060                                    const reco::JetCorrector* correctorL1L2L3Res) const {
0061   float var[kLast]{};
0062 
0063   var[kPt] = muon.pt();
0064   var[kEta] = muon.eta();
0065   var[kSegmentCompatibility] = muon.segmentCompatibility();
0066   var[kMiniRelIsoCharged] = muon.miniPFIsolation().chargedHadronIso() / muon.pt();
0067   var[kMiniRelIsoNeutral] = miniIsoValue - var[kMiniRelIsoCharged];
0068 
0069   double dB2D = fabs(muon.dB(pat::Muon::PV2D));
0070   double dB3D = muon.dB(pat::Muon::PV3D);
0071   double edB3D = muon.edB(pat::Muon::PV3D);
0072   double dz = fabs(muon.muonBestTrack()->dz(vertex.position()));
0073   var[kSip] = edB3D > 0 ? fabs(dB3D / edB3D) : 0.0;
0074   var[kLog_abs_dxyBS] = dB2D > 0 ? log(dB2D) : 0;
0075   var[kLog_abs_dzPV] = dz > 0 ? log(dz) : 0;
0076 
0077   //Initialise loop variables
0078   double minDr = 9999;
0079   double jecL1L2L3Res = 1.;
0080   double jecL1 = 1.;
0081 
0082   // Compute corrected isolation variables
0083   double chIso = muon.pfIsolationR04().sumChargedHadronPt;
0084   double nIso = muon.pfIsolationR04().sumNeutralHadronEt;
0085   double phoIso = muon.pfIsolationR04().sumPhotonEt;
0086   double puIso = muon.pfIsolationR04().sumPUPt;
0087   double dbCorrectedIsolation = chIso + std::max(nIso + phoIso - .5 * puIso, 0.);
0088   double dbCorrectedRelIso = dbCorrectedIsolation / muon.pt();
0089 
0090   var[kJetPtRatio] = 1. / (1 + dbCorrectedRelIso);
0091   var[kJetPtRel] = 0;
0092   var[kJetBTagCSV] = -999;
0093   var[kJetNDauCharged] = -1;
0094 
0095   for (const auto& tagI : bTags) {
0096     // for each muon with the lepton
0097     double dr = deltaR(*(tagI.first), muon);
0098     if (dr > minDr)
0099       continue;
0100     minDr = dr;
0101 
0102     const reco::Candidate::LorentzVector& muP4(muon.p4());
0103     reco::Candidate::LorentzVector jetP4(tagI.first->p4());
0104 
0105     if (correctorL1 && correctorL1L2L3Res) {
0106       jecL1L2L3Res = correctorL1L2L3Res->correction(*(tagI.first));
0107       jecL1 = correctorL1->correction(*(tagI.first));
0108     }
0109 
0110     // Get b-jet info
0111     var[kJetBTagCSV] = tagI.second;
0112     var[kJetNDauCharged] = 0;
0113     for (auto jet : tagI.first->getJetConstituentsQuick()) {
0114       const reco::PFCandidate* pfcand = dynamic_cast<const reco::PFCandidate*>(jet);
0115       if (pfcand == nullptr)
0116         throw cms::Exception("ConfigurationError") << "Cannot get jet constituents";
0117       if (pfcand->charge() == 0)
0118         continue;
0119       auto bestTrackPtr = pfcand->bestTrack();
0120       if (!bestTrackPtr)
0121         continue;
0122       if (!bestTrackPtr->quality(reco::Track::highPurity))
0123         continue;
0124       if (bestTrackPtr->pt() < 1.)
0125         continue;
0126       if (bestTrackPtr->hitPattern().numberOfValidHits() < 8)
0127         continue;
0128       if (bestTrackPtr->hitPattern().numberOfValidPixelHits() < 2)
0129         continue;
0130       if (bestTrackPtr->normalizedChi2() >= 5)
0131         continue;
0132 
0133       if (std::fabs(bestTrackPtr->dxy(vertex.position())) > 0.2)
0134         continue;
0135       if (std::fabs(bestTrackPtr->dz(vertex.position())) > 17)
0136         continue;
0137       var[kJetNDauCharged]++;
0138     }
0139 
0140     if (minDr < dRmax_) {
0141       if ((jetP4 - muP4).Rho() < 0.0001) {
0142         var[kJetPtRel] = 0;
0143         var[kJetPtRatio] = 1;
0144       } else {
0145         jetP4 -= muP4 / jecL1;
0146         jetP4 *= jecL1L2L3Res;
0147         jetP4 += muP4;
0148 
0149         var[kJetPtRatio] = muP4.pt() / jetP4.pt();
0150         var[kJetPtRel] = ptRel(muP4, jetP4);
0151       }
0152     }
0153   }
0154 
0155   if (var[kJetPtRatio] > 1.5)
0156     var[kJetPtRatio] = 1.5;
0157   if (var[kJetBTagCSV] < 0)
0158     var[kJetBTagCSV] = 0;
0159   jetPtRatio = var[kJetPtRatio];
0160   jetPtRel = var[kJetPtRel];
0161   return gbrForest_->GetClassifier(var);
0162 };