METPUSubtraction/plugins/PFMETProducerMVA.cc

0001 #include "RecoMET/METPUSubtraction/plugins/PFMETProducerMVA.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 using namespace reco;
0005
0006 const double dR2Min = 0.01 * 0.01;
0007 const double dR2Max = 0.5 * 0.5;
0008 const double dPtMatch = 0.1;
0009
0010 PFMETProducerMVA::PFMETProducerMVA(const edm::ParameterSet& cfg)
0011     : mvaMEtAlgo_(cfg, consumesCollector()), mvaMEtAlgo_isInitialized_(false) {
0012   srcCorrJets_ = consumes<reco::PFJetCollection>(cfg.getParameter<edm::InputTag>("srcCorrJets"));
0013   srcUncorrJets_ = consumes<reco::PFJetCollection>(cfg.getParameter<edm::InputTag>("srcUncorrJets"));
0014   srcJetIds_ = consumes<edm::ValueMap<float> >(cfg.getParameter<edm::InputTag>("srcMVAPileupJetId"));
0015   srcPFCandidatesView_ = consumes<reco::CandidateView>(cfg.getParameter<edm::InputTag>("srcPFCandidates"));
0016   srcVertices_ = consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("srcVertices"));
0017   vInputTag srcLeptonsTags = cfg.getParameter<vInputTag>("srcLeptons");
0018   for (vInputTag::const_iterator it = srcLeptonsTags.begin(); it != srcLeptonsTags.end(); it++) {
0019     srcLeptons_.push_back(consumes<reco::CandidateView>(*it));
0020   }
0021   mJetCorrector_ = consumes<reco::JetCorrector>(cfg.getParameter<edm::InputTag>("corrector"));
0022   minNumLeptons_ = cfg.getParameter<int>("minNumLeptons");
0023
0024   globalThreshold_ = cfg.getParameter<double>("globalThreshold");
0025
0026   minCorrJetPt_ = cfg.getParameter<double>("minCorrJetPt");
0027   useType1_ = cfg.getParameter<bool>("useType1");
0028
0029   verbosity_ = (cfg.exists("verbosity")) ? cfg.getParameter<int>("verbosity") : 0;
0030
0031   produces<reco::PFMETCollection>();
0032 }
0033
0034 PFMETProducerMVA::~PFMETProducerMVA() {}
0035
0036 void PFMETProducerMVA::produce(edm::Event& evt, const edm::EventSetup& es) {
0037   // CV: check if the event is to be skipped
0038   if (minNumLeptons_ > 0) {
0039     int numLeptons = 0;
0040     for (std::vector<edm::EDGetTokenT<reco::CandidateView> >::const_iterator srcLeptons_i = srcLeptons_.begin();
0041          srcLeptons_i != srcLeptons_.end();
0042          ++srcLeptons_i) {
0043       edm::Handle<reco::CandidateView> leptons;
0044       evt.getByToken(*srcLeptons_i, leptons);
0045       numLeptons += leptons->size();
0046     }
0047     if (!(numLeptons >= minNumLeptons_)) {
0048       LogDebug("produce") << "<PFMETProducerMVA::produce>:" << std::endl
0049                           << "Run: " << evt.id().run() << ", LS: " << evt.luminosityBlock()
0050                           << ", Event: " << evt.id().event() << std::endl
0051                           << " numLeptons = " << numLeptons << ", minNumLeptons = " << minNumLeptons_
0052                           << " --> skipping !!" << std::endl;
0053
0054       reco::PFMET pfMEt;
0055       auto pfMEtCollection = std::make_unique<reco::PFMETCollection>();
0056       pfMEtCollection->push_back(pfMEt);
0057       evt.put(std::move(pfMEtCollection));
0058       return;
0059     }
0060   }
0061
0062   //get jet IDs
0063   edm::Handle<edm::ValueMap<float> > jetIds;
0064   evt.getByToken(srcJetIds_, jetIds);
0065
0066   // get jets (corrected and uncorrected)
0067   edm::Handle<reco::PFJetCollection> corrJets;
0068   evt.getByToken(srcCorrJets_, corrJets);
0069
0070   edm::Handle<reco::PFJetCollection> uncorrJets;
0071   evt.getByToken(srcUncorrJets_, uncorrJets);
0072
0073   edm::Handle<reco::JetCorrector> corrector;
0074   if (useType1_) {
0075     evt.getByToken(mJetCorrector_, corrector);
0076   }
0077
0078   edm::Handle<reco::CandidateView> pfCandidates_view;
0079   evt.getByToken(srcPFCandidatesView_, pfCandidates_view);
0080
0081   // get vertices
0082   edm::Handle<reco::VertexCollection> vertices;
0083   evt.getByToken(srcVertices_, vertices);
0084   // take vertex with highest sum(trackPt) as the vertex of the "hard scatter" interaction
0085   // (= first entry in vertex collection)
0086   const reco::Vertex* hardScatterVertex = (!vertices->empty()) ? &(vertices->front()) : nullptr;
0087
0088   // get leptons
0089   // (excluded from sum over PFCandidates when computing hadronic recoil)
0090   int lId = 0;
0091   bool lHasPhotons = false;
0092   std::vector<reco::PUSubMETCandInfo> leptonInfo =
0093       computeLeptonInfo(srcLeptons_, *pfCandidates_view, hardScatterVertex, lId, lHasPhotons, evt);
0094
0095   // initialize MVA MET algorithm
0096   // (this will load the BDTs, stored as GBRForrest objects;
0097   //  either in input ROOT files or in SQL-lite files/the Conditions Database)
0098   if (!mvaMEtAlgo_isInitialized_) {
0099     mvaMEtAlgo_.initialize(es);
0100     mvaMEtAlgo_isInitialized_ = true;
0101   }
0102
0103   // reconstruct "standard" particle-flow missing Et
0104   CommonMETData pfMEt_data = metAlgo_.run((*pfCandidates_view), globalThreshold_);
0105   SpecificPFMETData specificPfMET = pfMEtSpecificAlgo_.run((*pfCandidates_view));
0106   const reco::Candidate::LorentzVector p4(pfMEt_data.mex, pfMEt_data.mey, 0.0, pfMEt_data.met);
0107   const reco::Candidate::Point vtx(0.0, 0.0, 0.0);
0108   reco::PFMET pfMEt(specificPfMET, pfMEt_data.sumet, p4, vtx);
0109   reco::Candidate::LorentzVector pfMEtP4_original = pfMEt.p4();
0110
0111   // compute objects specific to MVA based MET reconstruction
0112   std::vector<reco::PUSubMETCandInfo> pfCandidateInfo = computePFCandidateInfo(*pfCandidates_view, hardScatterVertex);
0113   std::vector<reco::PUSubMETCandInfo> jetInfo = computeJetInfo(
0114       *uncorrJets, corrJets, *jetIds, *vertices, hardScatterVertex, *corrector, evt, es, leptonInfo, pfCandidateInfo);
0115   std::vector<reco::Vertex::Point> vertexInfo = computeVertexInfo(*vertices);
0116   // compute MVA based MET and estimate of its uncertainty
0117   mvaMEtAlgo_.setInput(leptonInfo, jetInfo, pfCandidateInfo, vertexInfo);
0118   mvaMEtAlgo_.setHasPhotons(lHasPhotons);
0119   mvaMEtAlgo_.evaluateMVA();
0120   pfMEt.setP4(mvaMEtAlgo_.getMEt());
0121   pfMEt.setSignificanceMatrix(mvaMEtAlgo_.getMEtCov());
0122
0123   LogDebug("produce") << "Run: " << evt.id().run() << ", LS: " << evt.luminosityBlock()
0124                       << ", Event: " << evt.id().event() << std::endl
0125                       << " PFMET: Pt = " << pfMEtP4_original.pt() << ", phi = " << pfMEtP4_original.phi() << " "
0126                       << "(Px = " << pfMEtP4_original.px() << ", Py = " << pfMEtP4_original.py() << ")" << std::endl
0127                       << " MVA MET: Pt = " << pfMEt.pt() << " phi = " << pfMEt.phi() << " (Px = " << pfMEt.px()
0128                       << ", Py = " << pfMEt.py() << ")" << std::endl
0129                       << " Cov:" << std::endl
0130                       << (mvaMEtAlgo_.getMEtCov())(0, 0) << "  " << (mvaMEtAlgo_.getMEtCov())(0, 1) << std::endl
0131                       << (mvaMEtAlgo_.getMEtCov())(1, 0) << "  " << (mvaMEtAlgo_.getMEtCov())(1, 1) << std::endl
0132                       << std::endl;
0133
0134   // add PFMET object to the event
0135   auto pfMEtCollection = std::make_unique<reco::PFMETCollection>();
0136   pfMEtCollection->push_back(pfMEt);
0137   evt.put(std::move(pfMEtCollection));
0138 }
0139
0140 std::vector<reco::PUSubMETCandInfo> PFMETProducerMVA::computeLeptonInfo(
0141     const std::vector<edm::EDGetTokenT<reco::CandidateView> >& srcLeptons_,
0142     const reco::CandidateView& pfCandidates_view,
0143     const reco::Vertex* hardScatterVertex,
0144     int& lId,
0145     bool& lHasPhotons,
0146     edm::Event& evt) {
0147   std::vector<reco::PUSubMETCandInfo> leptonInfo;
0148
0149   for (std::vector<edm::EDGetTokenT<reco::CandidateView> >::const_iterator srcLeptons_i = srcLeptons_.begin();
0150        srcLeptons_i != srcLeptons_.end();
0151        ++srcLeptons_i) {
0152     edm::Handle<reco::CandidateView> leptons;
0153     evt.getByToken(*srcLeptons_i, leptons);
0154     for (reco::CandidateView::const_iterator lepton1 = leptons->begin(); lepton1 != leptons->end(); ++lepton1) {
0155       bool pMatch = false;
0156       for (std::vector<edm::EDGetTokenT<reco::CandidateView> >::const_iterator srcLeptons_j = srcLeptons_.begin();
0157            srcLeptons_j != srcLeptons_.end();
0158            ++srcLeptons_j) {
0159         edm::Handle<reco::CandidateView> leptons2;
0160         evt.getByToken(*srcLeptons_j, leptons2);
0161         for (reco::CandidateView::const_iterator lepton2 = leptons2->begin(); lepton2 != leptons2->end(); ++lepton2) {
0162           if (&(*lepton1) == &(*lepton2)) {
0163             continue;
0164           }
0165           if (deltaR2(lepton1->p4(), lepton2->p4()) < dR2Max) {
0166             pMatch = true;
0167           }
0168           if (pMatch && !istau(&(*lepton1)) && istau(&(*lepton2))) {
0169             pMatch = false;
0170           }
0171           if (pMatch && ((istau(&(*lepton1)) && istau(&(*lepton2))) || (!istau(&(*lepton1)) && !istau(&(*lepton2)))) &&
0172               lepton1->pt() > lepton2->pt()) {
0173             pMatch = false;
0174           }
0175           if (pMatch && lepton1->pt() == lepton2->pt()) {
0176             pMatch = false;
0177             for (unsigned int i0 = 0; i0 < leptonInfo.size(); i0++) {
0178               if (std::abs(lepton1->pt() - leptonInfo[i0].p4().pt()) < dPtMatch) {
0179                 pMatch = true;
0180                 break;
0181               }
0182             }
0183           }
0184           if (pMatch)
0185             break;
0186         }
0187         if (pMatch)
0188           break;
0189       }
0190       if (pMatch)
0191         continue;
0192       reco::PUSubMETCandInfo pLeptonInfo;
0193       pLeptonInfo.setP4(lepton1->p4());
0194       pLeptonInfo.setChargedEnFrac(chargedEnFrac(&(*lepton1), pfCandidates_view, hardScatterVertex));
0195       leptonInfo.push_back(pLeptonInfo);
0196       if (lepton1->isPhoton()) {
0197         lHasPhotons = true;
0198       }
0199     }
0200     lId++;
0201   }
0202
0203   return leptonInfo;
0204 }
0205
0206 std::vector<reco::PUSubMETCandInfo> PFMETProducerMVA::computeJetInfo(const reco::PFJetCollection& uncorrJets,
0207                                                                      const edm::Handle<reco::PFJetCollection>& corrJets,
0208                                                                      const edm::ValueMap<float>& jetIds,
0209                                                                      const reco::VertexCollection& vertices,
0210                                                                      const reco::Vertex* hardScatterVertex,
0211                                                                      const reco::JetCorrector& iCorrector,
0212                                                                      edm::Event& iEvent,
0213                                                                      const edm::EventSetup& iSetup,
0214                                                                      std::vector<reco::PUSubMETCandInfo>& iLeptons,
0215                                                                      std::vector<reco::PUSubMETCandInfo>& iCands) {
0216   std::vector<reco::PUSubMETCandInfo> retVal;
0217   for (reco::PFJetCollection::const_iterator uncorrJet = uncorrJets.begin(); uncorrJet != uncorrJets.end();
0218        ++uncorrJet) {
0219     // for ( reco::PFJetCollection::const_iterator corrJet = corrJets.begin();
0220     //    corrJet != corrJets.end(); ++corrJet ) {
0221     auto corrJet = corrJets->begin();
0222     for (size_t cjIdx = 0; cjIdx < corrJets->size(); ++cjIdx, ++corrJet) {
0223       reco::PFJetRef corrJetRef(corrJets, cjIdx);
0224
0225       // match corrected and uncorrected jets
0226       if (uncorrJet->jetArea() != corrJet->jetArea())
0227         continue;
0228       if (deltaR2(corrJet->p4(), uncorrJet->p4()) > dR2Min)
0229         continue;
0230
0231       // check that jet passes loose PFJet id.
0232       if (!passPFLooseId(&(*uncorrJet)))
0233         continue;
0234
0235       // compute jet energy correction factor
0236       // (= ratio of corrected/uncorrected jet Pt)
0237       //double jetEnCorrFactor = corrJet->pt()/uncorrJet->pt();
0238       reco::PUSubMETCandInfo jetInfo;
0239
0240       // PH: apply jet energy corrections for all Jets ignoring recommendations
0241       jetInfo.setP4(corrJet->p4());
0242       double lType1Corr = 0;
0243       if (useType1_) {  //Compute the type 1 correction ===> This code is crap
0244         double pCorr = iCorrector.correction(*uncorrJet);
0245         lType1Corr = std::abs(corrJet->pt() - pCorr * uncorrJet->pt());
0246         TLorentzVector pVec;
0247         pVec.SetPtEtaPhiM(lType1Corr, 0, corrJet->phi(), 0);
0248         reco::Candidate::LorentzVector pType1Corr;
0249         pType1Corr.SetCoordinates(pVec.Px(), pVec.Py(), pVec.Pz(), pVec.E());
0250         //Filter to leptons
0251         bool pOnLepton = false;
0252         for (unsigned int i0 = 0; i0 < iLeptons.size(); i0++) {
0253           if (deltaR2(iLeptons[i0].p4(), corrJet->p4()) < dR2Max) {
0254             pOnLepton = true;
0255             break;
0256           }
0257         }
0258         //Add it to PF Collection
0259         if (corrJet->pt() > 10 && !pOnLepton) {
0260           reco::PUSubMETCandInfo pfCandidateInfo;
0261           pfCandidateInfo.setP4(pType1Corr);
0262           pfCandidateInfo.setDZ(-999);
0263           iCands.push_back(pfCandidateInfo);
0264         }
0265         //Scale
0266         lType1Corr = (pCorr * uncorrJet->pt() - uncorrJet->pt());
0267         lType1Corr /= corrJet->pt();
0268       }
0269
0270       // check that jet Pt used to compute MVA based jet id. is above threshold
0271       if (!(jetInfo.p4().pt() > minCorrJetPt_))
0272         continue;
0273
0274       jetInfo.setMvaVal(jetIds[corrJetRef]);
0275       float chEnF = (uncorrJet->chargedEmEnergy() + uncorrJet->chargedHadronEnergy() + uncorrJet->chargedMuEnergy()) /
0276                     uncorrJet->energy();
0277       if (useType1_)
0278         chEnF += lType1Corr * (1 - jetInfo.chargedEnFrac());
0279       jetInfo.setChargedEnFrac(chEnF);
0280       retVal.push_back(jetInfo);
0281       break;
0282     }
0283   }
0284
0285   //order jets per pt
0286   std::sort(retVal.begin(), retVal.end());
0287
0288   return retVal;
0289 }
0290
0291 std::vector<reco::PUSubMETCandInfo> PFMETProducerMVA::computePFCandidateInfo(const reco::CandidateView& pfCandidates,
0292                                                                              const reco::Vertex* hardScatterVertex) {
0293   std::vector<reco::PUSubMETCandInfo> retVal;
0294   for (reco::CandidateView::const_iterator pfCandidate = pfCandidates.begin(); pfCandidate != pfCandidates.end();
0295        ++pfCandidate) {
0296     double dZ = -999.;  // PH: If no vertex is reconstructed in the event
0297                         //     or PFCandidate has no track, set dZ to -999
0298     if (hardScatterVertex) {
0299       const reco::PFCandidate* pfc = dynamic_cast<const reco::PFCandidate*>(&(*pfCandidate));
0300       if (pfc != nullptr) {  //PF candidate for RECO and PAT levels
0301         if (pfc->trackRef().isNonnull())
0302           dZ = std::abs(pfc->trackRef()->dz(hardScatterVertex->position()));
0303         else if (pfc->gsfTrackRef().isNonnull())
0304           dZ = std::abs(pfc->gsfTrackRef()->dz(hardScatterVertex->position()));
0305       } else {  //if not, then packedCandidate for miniAOD level
0306         const pat::PackedCandidate* pfc = dynamic_cast<const pat::PackedCandidate*>(&(*pfCandidate));
0307         dZ = std::abs(pfc->dz(hardScatterVertex->position()));
0308         //exact dz=zero corresponds to the -999 case for pfcandidate
0309         if (dZ == 0) {
0310           dZ = -999;
0311         }
0312       }
0313     }
0314     reco::PUSubMETCandInfo pfCandidateInfo;
0315     pfCandidateInfo.setP4(pfCandidate->p4());
0316     pfCandidateInfo.setDZ(dZ);
0317     retVal.push_back(pfCandidateInfo);
0318   }
0319   return retVal;
0320 }
0321
0322 std::vector<reco::Vertex::Point> PFMETProducerMVA::computeVertexInfo(const reco::VertexCollection& vertices) {
0323   std::vector<reco::Vertex::Point> retVal;
0324   for (reco::VertexCollection::const_iterator vertex = vertices.begin(); vertex != vertices.end(); ++vertex) {
0325     if (std::abs(vertex->z()) > 24.)
0326       continue;
0327     if (vertex->ndof() < 4.)
0328       continue;
0329     if (vertex->position().Rho() > 2.)
0330       continue;
0331     retVal.push_back(vertex->position());
0332   }
0333   return retVal;
0334 }
0335 double PFMETProducerMVA::chargedEnFrac(const reco::Candidate* iCand,
0336                                        const reco::CandidateView& pfCandidates,
0337                                        const reco::Vertex* hardScatterVertex) {
0338   if (iCand->isMuon()) {
0339     return 1;
0340   }
0341   if (iCand->isElectron()) {
0342     return 1.;
0343   }
0344   if (iCand->isPhoton()) {
0345     return chargedFracInCone(iCand, pfCandidates, hardScatterVertex);
0346   }
0347   double lPtTot = 0;
0348   double lPtCharged = 0;
0349   const reco::PFTau* lPFTau = nullptr;
0350   lPFTau = dynamic_cast<const reco::PFTau*>(iCand);
0351   if (lPFTau != nullptr) {
0352     for (UInt_t i0 = 0; i0 < lPFTau->signalCands().size(); i0++) {
0353       lPtTot += (lPFTau->signalCands())[i0]->pt();
0354       if ((lPFTau->signalCands())[i0]->charge() == 0)
0355         continue;
0356       lPtCharged += (lPFTau->signalCands())[i0]->pt();
0357     }
0358   } else {
0359     const pat::Tau* lPatPFTau = nullptr;
0360     lPatPFTau = dynamic_cast<const pat::Tau*>(iCand);
0361     if (lPatPFTau != nullptr) {
0362       for (UInt_t i0 = 0; i0 < lPatPFTau->signalCands().size(); i0++) {
0363         lPtTot += (lPatPFTau->signalCands())[i0]->pt();
0364         if ((lPatPFTau->signalCands())[i0]->charge() == 0)
0365           continue;
0366         lPtCharged += (lPatPFTau->signalCands())[i0]->pt();
0367       }
0368     }
0369   }
0370   if (lPtTot == 0)
0371     lPtTot = 1.;
0372   return lPtCharged / lPtTot;
0373 }
0374 //Return tau id by process of elimination
0375 bool PFMETProducerMVA::istau(const reco::Candidate* iCand) {
0376   if (iCand->isMuon())
0377     return false;
0378   if (iCand->isElectron())
0379     return false;
0380   if (iCand->isPhoton())
0381     return false;
0382   return true;
0383 }
0384 bool PFMETProducerMVA::passPFLooseId(const PFJet* iJet) {
0385   if (iJet->energy() == 0)
0386     return false;
0387   if (iJet->neutralHadronEnergy() / iJet->energy() > 0.99)
0388     return false;
0389   if (iJet->neutralEmEnergy() / iJet->energy() > 0.99)
0390     return false;
0391   if (iJet->nConstituents() < 2)
0392     return false;
0393   if (iJet->chargedHadronEnergy() / iJet->energy() <= 0 && std::abs(iJet->eta()) < 2.4)
0394     return false;
0395   if (iJet->chargedEmEnergy() / iJet->energy() > 0.99 && std::abs(iJet->eta()) < 2.4)
0396     return false;
0397   if (iJet->chargedMultiplicity() < 1 && std::abs(iJet->eta()) < 2.4)
0398     return false;
0399   return true;
0400 }
0401
0402 double PFMETProducerMVA::chargedFracInCone(const reco::Candidate* iCand,
0403                                            const reco::CandidateView& pfCandidates,
0404                                            const reco::Vertex* hardScatterVertex,
0405                                            double iDRMax) {
0406   double iDR2Max = iDRMax * iDRMax;
0407   reco::Candidate::LorentzVector lVis(0, 0, 0, 0);
0408   for (reco::CandidateView::const_iterator pfCandidate = pfCandidates.begin(); pfCandidate != pfCandidates.end();
0409        ++pfCandidate) {
0410     if (deltaR2(iCand->p4(), pfCandidate->p4()) > iDR2Max)
0411       continue;
0412     double dZ = -999.;  // PH: If no vertex is reconstructed in the event
0413                         //     or PFCandidate has no track, set dZ to -999
0414     if (hardScatterVertex) {
0415       const reco::PFCandidate* pfc = dynamic_cast<const reco::PFCandidate*>((&(*pfCandidate)));
0416       if (pfc != nullptr) {  //PF candidate for RECO and PAT levels
0417         if (pfc->trackRef().isNonnull())
0418           dZ = std::abs(pfc->trackRef()->dz(hardScatterVertex->position()));
0419         else if (pfc->gsfTrackRef().isNonnull())
0420           dZ = std::abs(pfc->gsfTrackRef()->dz(hardScatterVertex->position()));
0421       } else {  //if not, then packedCandidate for miniAOD level
0422         const pat::PackedCandidate* pfc = dynamic_cast<const pat::PackedCandidate*>(&(*pfCandidate));
0423         dZ = std::abs(pfc->dz(hardScatterVertex->position()));
0424       }
0425     }
0426     if (std::abs(dZ) > 0.1)
0427       continue;
0428     lVis += pfCandidate->p4();
0429   }
0430   return lVis.pt() / iCand->pt();
0431 }
0432
0433 #include "FWCore/Framework/interface/MakerMacros.h"
0434
0435 DEFINE_FWK_MODULE(PFMETProducerMVA);