Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:48 

0001 #include "RecoMET/METPUSubtraction/interface/MvaMEtUtilities.h"
0002 
0003 #include "DataFormats/Math/interface/deltaR.h"
0004 
0005 #include <algorithm>
0006 #include <cmath>
0007 
0008 MvaMEtUtilities::MvaMEtUtilities(const edm::ParameterSet& cfg) {
0009   // jet id
0010   /* ===> this code parses an xml it uses parameter set 
0011   edm::ParameterSet jetConfig = iConfig.getParameter<edm::ParameterSet>("JetIdParams");
0012   for(int i0 = 0; i0 < 3; i0++) { 
0013     std::string lCutType                            = "Tight";
0014     if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium";
0015     if(i0 == PileupJetIdentifier::kLoose)  lCutType = "Loose";
0016     std::vector<double> pt010  = jetConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
0017     std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
0018     std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
0019     std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
0020     for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][0][i1] = pt010 [i1];
0021     for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][1][i1] = pt1020[i1];
0022     for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][2][i1] = pt2030[i1];
0023     for(int i1 = 0; i1 < 4; i1++) mvacut_[i0][3][i1] = pt3050[i1];
0024     }
0025   */
0026   //Tight Id => not used
0027   mvaCut_[0][0][0] = 0.5;
0028   mvaCut_[0][0][1] = 0.6;
0029   mvaCut_[0][0][2] = 0.6;
0030   mvaCut_[0][0][3] = 0.9;
0031   mvaCut_[0][1][0] = -0.2;
0032   mvaCut_[0][1][1] = 0.2;
0033   mvaCut_[0][1][2] = 0.2;
0034   mvaCut_[0][1][3] = 0.6;
0035   mvaCut_[0][2][0] = 0.3;
0036   mvaCut_[0][2][1] = 0.4;
0037   mvaCut_[0][2][2] = 0.7;
0038   mvaCut_[0][2][3] = 0.8;
0039   mvaCut_[0][3][0] = 0.5;
0040   mvaCut_[0][3][1] = 0.4;
0041   mvaCut_[0][3][2] = 0.8;
0042   mvaCut_[0][3][3] = 0.9;
0043   //Medium id => not used
0044   mvaCut_[1][0][0] = 0.2;
0045   mvaCut_[1][0][1] = 0.4;
0046   mvaCut_[1][0][2] = 0.2;
0047   mvaCut_[1][0][3] = 0.6;
0048   mvaCut_[1][1][0] = -0.3;
0049   mvaCut_[1][1][1] = 0.;
0050   mvaCut_[1][1][2] = 0.;
0051   mvaCut_[1][1][3] = 0.5;
0052   mvaCut_[1][2][0] = 0.2;
0053   mvaCut_[1][2][1] = 0.2;
0054   mvaCut_[1][2][2] = 0.5;
0055   mvaCut_[1][2][3] = 0.7;
0056   mvaCut_[1][3][0] = 0.3;
0057   mvaCut_[1][3][1] = 0.2;
0058   mvaCut_[1][3][2] = 0.7;
0059   mvaCut_[1][3][3] = 0.8;
0060   //Met Id => used
0061   mvaCut_[2][0][0] = -0.2;
0062   mvaCut_[2][0][1] = -0.3;
0063   mvaCut_[2][0][2] = -0.5;
0064   mvaCut_[2][0][3] = -0.5;
0065   mvaCut_[2][1][0] = -0.2;
0066   mvaCut_[2][1][1] = -0.2;
0067   mvaCut_[2][1][2] = -0.5;
0068   mvaCut_[2][1][3] = -0.3;
0069   mvaCut_[2][2][0] = -0.2;
0070   mvaCut_[2][2][1] = -0.2;
0071   mvaCut_[2][2][2] = -0.2;
0072   mvaCut_[2][2][3] = 0.1;
0073   mvaCut_[2][3][0] = -0.2;
0074   mvaCut_[2][3][1] = -0.2;
0075   mvaCut_[2][3][2] = 0.;
0076   mvaCut_[2][3][3] = 0.2;
0077 
0078   dzCut_ = cfg.getParameter<double>("dZcut");
0079   ptThreshold_ = (cfg.exists("ptThreshold")) ? cfg.getParameter<int>("ptThreshold") : -1000;
0080 }
0081 
0082 MvaMEtUtilities::~MvaMEtUtilities() {
0083   // nothing to be done yet...
0084 }
0085 
0086 bool MvaMEtUtilities::passesMVA(const reco::Candidate::LorentzVector& jetP4, double mvaJetId) {
0087   int ptBin = 0;
0088   if (jetP4.pt() >= 10. && jetP4.pt() < 20.)
0089     ptBin = 1;
0090   if (jetP4.pt() >= 20. && jetP4.pt() < 30.)
0091     ptBin = 2;
0092   if (jetP4.pt() >= 30.)
0093     ptBin = 3;
0094 
0095   int etaBin = 0;
0096   if (std::abs(jetP4.eta()) >= 2.5 && std::abs(jetP4.eta()) < 2.75)
0097     etaBin = 1;
0098   if (std::abs(jetP4.eta()) >= 2.75 && std::abs(jetP4.eta()) < 3.0)
0099     etaBin = 2;
0100   if (std::abs(jetP4.eta()) >= 3.0 && std::abs(jetP4.eta()) < 5.0)
0101     etaBin = 3;
0102 
0103   return (mvaJetId > mvaCut_[2][ptBin][etaBin]);
0104 }
0105 
0106 reco::Candidate::LorentzVector MvaMEtUtilities::leadJetP4(const std::vector<reco::PUSubMETCandInfo>& jets) {
0107   return jetP4(jets, 0);
0108 }
0109 
0110 reco::Candidate::LorentzVector MvaMEtUtilities::subleadJetP4(const std::vector<reco::PUSubMETCandInfo>& jets) {
0111   return jetP4(jets, 1);
0112 }
0113 
0114 reco::Candidate::LorentzVector MvaMEtUtilities::jetP4(const std::vector<reco::PUSubMETCandInfo>& jets, unsigned idx) {
0115   reco::Candidate::LorentzVector retVal(0., 0., 0., 0.);
0116   if (idx < jets.size()) {
0117     std::vector<reco::PUSubMETCandInfo> jets_sorted = jets;
0118     std::sort(jets_sorted.rbegin(), jets_sorted.rend());
0119     retVal = jets_sorted[idx].p4();
0120   }
0121   return retVal;
0122 }
0123 unsigned MvaMEtUtilities::numJetsAboveThreshold(const std::vector<reco::PUSubMETCandInfo>& jets, double ptThreshold) {
0124   unsigned retVal = 0;
0125   for (std::vector<reco::PUSubMETCandInfo>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
0126     if (jet->p4().pt() > ptThreshold)
0127       ++retVal;
0128   }
0129   return retVal;
0130 }
0131 std::vector<reco::PUSubMETCandInfo> MvaMEtUtilities::cleanJets(const std::vector<reco::PUSubMETCandInfo>& jets,
0132                                                                const std::vector<reco::PUSubMETCandInfo>& leptons,
0133                                                                double ptThreshold,
0134                                                                double dRmatch) {
0135   double dR2match = dRmatch * dRmatch;
0136   std::vector<reco::PUSubMETCandInfo> retVal;
0137   for (std::vector<reco::PUSubMETCandInfo>::const_iterator jet = jets.begin(); jet != jets.end(); ++jet) {
0138     bool isOverlap = false;
0139     for (std::vector<reco::PUSubMETCandInfo>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
0140          ++lepton) {
0141       if (deltaR2(jet->p4(), lepton->p4()) < dR2match)
0142         isOverlap = true;
0143     }
0144     if (jet->p4().pt() > ptThreshold && !isOverlap)
0145       retVal.push_back(*jet);
0146   }
0147   return retVal;
0148 }
0149 
0150 std::vector<reco::PUSubMETCandInfo> MvaMEtUtilities::cleanPFCands(
0151     const std::vector<reco::PUSubMETCandInfo>& pfCandidates,
0152     const std::vector<reco::PUSubMETCandInfo>& leptons,
0153     double dRmatch,
0154     bool invert) {
0155   double dR2match = dRmatch * dRmatch;
0156   std::vector<reco::PUSubMETCandInfo> retVal;
0157   for (std::vector<reco::PUSubMETCandInfo>::const_iterator pfCandidate = pfCandidates.begin();
0158        pfCandidate != pfCandidates.end();
0159        ++pfCandidate) {
0160     bool isOverlap = false;
0161     for (std::vector<reco::PUSubMETCandInfo>::const_iterator lepton = leptons.begin(); lepton != leptons.end();
0162          ++lepton) {
0163       if (deltaR2(pfCandidate->p4(), lepton->p4()) < dR2match)
0164         isOverlap = true;
0165     }
0166     if ((!isOverlap && !invert) || (isOverlap && invert))
0167       retVal.push_back(*pfCandidate);
0168   }
0169   return retVal;
0170 }
0171 void MvaMEtUtilities::finalize(CommonMETData& metData) {
0172   metData.met = sqrt(metData.mex * metData.mex + metData.mey * metData.mey);
0173   metData.mez = 0.;
0174   metData.phi = atan2(metData.mey, metData.mex);
0175 }
0176 
0177 CommonMETData MvaMEtUtilities::computeCandSum(int compKey,
0178                                               double dZmax,
0179                                               int dZflag,
0180                                               bool iCharged,
0181                                               bool mvaPassFlag,
0182                                               const std::vector<reco::PUSubMETCandInfo>& objects) {
0183   CommonMETData retVal;
0184   retVal.mex = 0.;
0185   retVal.mey = 0.;
0186   retVal.sumet = 0.;
0187 
0188   for (std::vector<reco::PUSubMETCandInfo>::const_iterator object = objects.begin(); object != objects.end();
0189        ++object) {
0190     double pFrac = 1;
0191 
0192     //pf candidates
0193     // dZcut
0194     //   maximum distance within which tracks are
0195     //considered to be associated to hard scatter vertex
0196     // dZflag
0197     //   0 : select charged PFCandidates originating from hard scatter vertex
0198     //   1 : select charged PFCandidates originating from pile-up vertices
0199     //   2 : select all PFCandidates
0200     if (compKey == MvaMEtUtilities::kPFCands) {
0201       if (object->dZ() < 0. && dZflag != 2)
0202         continue;
0203       if (object->dZ() > dZmax && dZflag == 0)
0204         continue;
0205       if (object->dZ() < dZmax && dZflag == 1)
0206         continue;
0207     }
0208 
0209     //leptons
0210     if (compKey == MvaMEtUtilities::kLeptons) {
0211       if (iCharged)
0212         pFrac = object->chargedEnFrac();
0213     }
0214 
0215     //jets
0216     if (compKey == MvaMEtUtilities::kJets) {
0217       bool passesMVAjetId = passesMVA(object->p4(), object->mva());
0218 
0219       if (passesMVAjetId && !mvaPassFlag)
0220         continue;
0221       if (!passesMVAjetId && mvaPassFlag)
0222         continue;
0223 
0224       pFrac = 1 - object->chargedEnFrac();  //neutral energy fraction
0225     }
0226 
0227     retVal.mex += object->p4().px() * pFrac;
0228     retVal.mey += object->p4().py() * pFrac;
0229     retVal.sumet += object->p4().pt() * pFrac;
0230   }
0231 
0232   finalize(retVal);
0233   return retVal;
0234 }
0235 
0236 CommonMETData MvaMEtUtilities::computeRecoil(int metType) {
0237   CommonMETData retVal;
0238 
0239   if (metType == MvaMEtUtilities::kPF) {
0240     //MET = pfMET = - all candidates
0241     // MET (1) in JME-13-003
0242     retVal.mex = leptonsSum_.mex - pfCandSum_.mex;
0243     retVal.mey = leptonsSum_.mey - pfCandSum_.mey;
0244     retVal.sumet = pfCandSum_.sumet - leptonsSum_.sumet;
0245   }
0246   if (metType == MvaMEtUtilities::kChHS) {
0247     //MET = - charged HS
0248     // MET (2) in JME-13-003
0249     retVal.mex = leptonsChSum_.mex - pfCandChHSSum_.mex;
0250     retVal.mey = leptonsChSum_.mey - pfCandChHSSum_.mey;
0251     retVal.sumet = pfCandChHSSum_.sumet - leptonsChSum_.sumet;
0252   }
0253   if (metType == MvaMEtUtilities::kHS) {
0254     //MET = - charged HS - neutral HS in jets
0255     // MET (3) in JME-13-003
0256     retVal.mex = leptonsChSum_.mex - (pfCandChHSSum_.mex + neutralJetHSSum_.mex);
0257     retVal.mey = leptonsChSum_.mey - (pfCandChHSSum_.mey + neutralJetHSSum_.mey);
0258     retVal.sumet = pfCandChHSSum_.sumet + neutralJetHSSum_.sumet - leptonsChSum_.sumet;
0259   }
0260   if (metType == MvaMEtUtilities::kPU) {
0261     //MET = - charged PU - neutral PU in jets
0262     //MET = -recoil in that particular case, - sign not useful for the MVA and then discarded
0263     //motivated as PU IS its own recoil
0264     // MET (4) in JME-13-003
0265     retVal.mex = -(pfCandChPUSum_.mex + neutralJetPUSum_.mex);
0266     retVal.mey = -(pfCandChPUSum_.mey + neutralJetPUSum_.mey);
0267     retVal.sumet = pfCandChPUSum_.sumet + neutralJetPUSum_.sumet;
0268   }
0269   if (metType == MvaMEtUtilities::kHSMinusNeutralPU) {
0270     //MET = all candidates - charged PU - neutral PU in jets
0271     // = all charged HS + all neutrals - neutral PU in jets
0272     // MET (5) in JME-13-003
0273     retVal.mex = leptonsSum_.mex - (pfCandSum_.mex - pfCandChPUSum_.mex - neutralJetPUSum_.mex);
0274     retVal.mey = leptonsSum_.mey - (pfCandSum_.mey - pfCandChPUSum_.mey - neutralJetPUSum_.mey);
0275     retVal.sumet = (pfCandSum_.sumet - pfCandChPUSum_.sumet - neutralJetPUSum_.sumet) - leptonsSum_.sumet;
0276   }
0277 
0278   finalize(retVal);
0279   return retVal;
0280 }
0281 
0282 void MvaMEtUtilities::computeAllSums(const std::vector<reco::PUSubMETCandInfo>& jets,
0283                                      const std::vector<reco::PUSubMETCandInfo>& leptons,
0284                                      const std::vector<reco::PUSubMETCandInfo>& pfCandidates) {
0285   cleanedJets_ = cleanJets(jets, leptons, ptThreshold_, 0.5);
0286 
0287   leptonsSum_ = computeCandSum(kLeptons, 0., 0, false, false, leptons);
0288   leptonsChSum_ = computeCandSum(kLeptons, 0., 0, true, false, leptons);
0289   pfCandSum_ = computeCandSum(kPFCands, dzCut_, 2, false, false, pfCandidates);
0290   pfCandChHSSum_ = computeCandSum(kPFCands, dzCut_, 0, false, false, pfCandidates);
0291   pfCandChPUSum_ = computeCandSum(kPFCands, dzCut_, 1, false, false, pfCandidates);
0292   neutralJetHSSum_ = computeCandSum(kJets, 0., 0, false, true, jets);
0293   neutralJetPUSum_ = computeCandSum(kJets, 0., 0, false, false, jets);
0294 }
0295 
0296 double MvaMEtUtilities::getLeptonsSumMEX() const { return leptonsSum_.mex; }
0297 
0298 double MvaMEtUtilities::getLeptonsSumMEY() const { return leptonsSum_.mey; }
0299 
0300 double MvaMEtUtilities::getLeptonsChSumMEX() const { return leptonsChSum_.mex; }
0301 
0302 double MvaMEtUtilities::getLeptonsChSumMEY() const { return leptonsChSum_.mey; }
0303 
0304 const std::vector<reco::PUSubMETCandInfo>& MvaMEtUtilities::getCleanedJets() const { return cleanedJets_; }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!