Physics/src/TopDQMHelpers.cc

0001 #include "DQM/Physics/interface/TopDQMHelpers.h"
0002
0003 Calculate::Calculate(int maxNJets, double wMass)
0004     : failed_(false),
0005       maxNJets_(maxNJets),
0006       wMass_(wMass),
0007       massWBoson_(-1.),
0008       massTopQuark_(-1.),
0009       massBTopQuark_(-1.),
0010       tmassWBoson_(-1),
0011       tmassTopQuark_(-1) {}
0012
0013 double Calculate::massWBoson(const std::vector<reco::Jet>& jets) {
0014   if (!failed_ && massWBoson_ < 0)
0015     operator()(jets);
0016   return massWBoson_;
0017 }
0018
0019 double Calculate::massTopQuark(const std::vector<reco::Jet>& jets) {
0020   if (!failed_ && massTopQuark_ < 0)
0021     operator()(jets);
0022   return massTopQuark_;
0023 }
0024
0025 double Calculate::massBTopQuark(const std::vector<reco::Jet>& jets, std::vector<double> VbtagWP, double btagWP_) {
0026   if (!failed_ && massBTopQuark_ < 0)
0027     operator2(jets, VbtagWP, btagWP_);
0028   return massBTopQuark_;
0029 }
0030
0031 double Calculate::tmassWBoson(reco::RecoCandidate* mu, const reco::MET& met, const reco::Jet& b) {
0032   if (tmassWBoson_ < 0)
0033     operator()(b, mu, met);
0034   return tmassWBoson_;
0035 }
0036
0037 double Calculate::tmassTopQuark(reco::RecoCandidate* lepton, const reco::MET& met, const reco::Jet& b) {
0038   if (tmassTopQuark_ < 0)
0039     operator()(b, lepton, met);
0040   return tmassTopQuark_;
0041 }
0042
0043 void Calculate::operator()(const reco::Jet& bJet, reco::RecoCandidate* lepton, const reco::MET& met) {
0044   double metT = sqrt(pow(met.px(), 2) + pow(met.py(), 2));
0045   double lepT = sqrt(pow(lepton->px(), 2) + pow(lepton->py(), 2));
0046   double bT = sqrt(pow(bJet.px(), 2) + pow(bJet.py(), 2));
0047   reco::Particle::LorentzVector WT = lepton->p4() + met.p4();
0048   tmassWBoson_ = sqrt(pow(metT + lepT, 2) - (WT.px() * WT.px()) - (WT.py() * WT.py()));
0049   reco::Particle::LorentzVector topT = WT + bJet.p4();
0050   tmassTopQuark_ = sqrt(pow((metT + lepT + bT), 2) - (topT.px() * topT.px()) - (topT.py() * topT.py()));
0051 }
0052
0053 void Calculate::operator()(const std::vector<reco::Jet>& jets) {
0054   if (maxNJets_ < 0)
0055     maxNJets_ = jets.size();
0056   failed_ = jets.size() < (unsigned int)maxNJets_;
0057   if (failed_) {
0058     return;
0059   }
0060
0061   // associate those jets with maximum pt of the vectorial
0062   // sum to the hadronic decay chain
0063   double maxPt = -1.;
0064   std::vector<int> maxPtIndices;
0065   maxPtIndices.push_back(-1);
0066   maxPtIndices.push_back(-1);
0067   maxPtIndices.push_back(-1);
0068
0069   for (int idx = 0; idx < maxNJets_; ++idx) {
0070     for (int jdx = idx + 1; jdx < maxNJets_; ++jdx) {
0071       //if (jdx <= idx) continue;
0072       for (int kdx = 0; kdx < maxNJets_; ++kdx) {
0073         if (kdx == idx || kdx == jdx)
0074           continue;
0075         reco::Particle::LorentzVector sum = jets[idx].p4() + jets[jdx].p4() + jets[kdx].p4();
0076         if (maxPt < 0. || maxPt < sum.pt()) {
0077           maxPt = sum.pt();
0078           maxPtIndices.clear();
0079           maxPtIndices.push_back(idx);
0080           maxPtIndices.push_back(jdx);
0081           maxPtIndices.push_back(kdx);
0082         }
0083       }
0084     }
0085   }
0086   massTopQuark_ = (jets[maxPtIndices[0]].p4() + jets[maxPtIndices[1]].p4() + jets[maxPtIndices[2]].p4()).mass();
0087
0088   // associate those jets that get closest to the W mass
0089   // with their invariant mass to the W boson
0090   double wDist = -1.;
0091   std::vector<int> wMassIndices;
0092   wMassIndices.push_back(-1);
0093   wMassIndices.push_back(-1);
0094   for (unsigned idx = 0; idx < maxPtIndices.size(); ++idx) {
0095     for (unsigned jdx = 0; jdx < maxPtIndices.size(); ++jdx) {
0096       if (jdx == idx || maxPtIndices[idx] > maxPtIndices[jdx])
0097         continue;
0098       reco::Particle::LorentzVector sum = jets[maxPtIndices[idx]].p4() + jets[maxPtIndices[jdx]].p4();
0099       if (wDist < 0. || wDist > fabs(sum.mass() - wMass_)) {
0100         wDist = fabs(sum.mass() - wMass_);
0101         wMassIndices.clear();
0102         wMassIndices.push_back(maxPtIndices[idx]);
0103         wMassIndices.push_back(maxPtIndices[jdx]);
0104       }
0105     }
0106   }
0107   massWBoson_ = (jets[wMassIndices[0]].p4() + jets[wMassIndices[1]].p4()).mass();
0108 }
0109
0110 void Calculate::operator2(const std::vector<reco::Jet>& jets, std::vector<double> bjet, double btagWP) {
0111   if (maxNJets_ < 0)
0112     maxNJets_ = jets.size();
0113   failed_ = jets.size() < (unsigned int)maxNJets_;
0114   if (failed_) {
0115     return;
0116   }
0117   if (jets.size() != bjet.size()) {
0118     return;
0119   }
0120
0121   // associate those jets with maximum pt of the vectorial
0122   // sum to the hadronic decay chain. Require ONLY 1 btagged jet
0123   double maxBPt = -1.;
0124   std::vector<int> maxBPtIndices;
0125   maxBPtIndices.push_back(-1);
0126   maxBPtIndices.push_back(-1);
0127   maxBPtIndices.push_back(-1);
0128   for (int idx = 0; idx < maxNJets_; ++idx) {
0129     for (int jdx = idx + 1; jdx < maxNJets_; ++jdx) {
0130       //if (jdx <= idx) continue;
0131       for (int kdx = 0; kdx < maxNJets_; ++kdx) {
0132         if (kdx == idx || kdx == jdx)
0133           continue;
0134         // require only 1b-jet
0135         if ((bjet[idx] > btagWP && bjet[jdx] <= btagWP && bjet[kdx] <= btagWP) ||
0136             (bjet[idx] <= btagWP && bjet[jdx] > btagWP && bjet[kdx] <= btagWP) ||
0137             (bjet[idx] <= btagWP && bjet[jdx] <= btagWP && bjet[kdx] > btagWP)) {
0138           reco::Particle::LorentzVector sum = jets[idx].p4() + jets[jdx].p4() + jets[kdx].p4();
0139           if (maxBPt < 0. || maxBPt < sum.pt()) {
0140             maxBPt = sum.pt();
0141             maxBPtIndices.clear();
0142             maxBPtIndices.push_back(idx);
0143             maxBPtIndices.push_back(jdx);
0144             maxBPtIndices.push_back(kdx);
0145           }
0146         }
0147       }
0148     }
0149   }
0150   if (maxBPtIndices[0] < 0 || maxBPtIndices[1] < 0 || maxBPtIndices[2] < 0)
0151     return;
0152   massBTopQuark_ = (jets[maxBPtIndices[0]].p4() + jets[maxBPtIndices[1]].p4() + jets[maxBPtIndices[2]].p4()).mass();
0153 }
0154
0155 Calculate_miniAOD::Calculate_miniAOD(int maxNJets, double wMass)
0156     : failed_(false),
0157       maxNJets_(maxNJets),
0158       wMass_(wMass),
0159       massWBoson_(-1.),
0160       massTopQuark_(-1.),
0161       massBTopQuark_(-1.),
0162       tmassWBoson_(-1),
0163       tmassTopQuark_(-1) {}
0164
0165 double Calculate_miniAOD::massWBoson(const std::vector<pat::Jet>& jets) {
0166   if (!failed_ && massWBoson_ < 0)
0167     operator()(jets);
0168   return massWBoson_;
0169 }
0170
0171 double Calculate_miniAOD::massTopQuark(const std::vector<pat::Jet>& jets) {
0172   if (!failed_ && massTopQuark_ < 0)
0173     operator()(jets);
0174   return massTopQuark_;
0175 }
0176
0177 double Calculate_miniAOD::massBTopQuark(const std::vector<pat::Jet>& jets,
0178                                         std::vector<double> VbtagWP,
0179                                         double btagWP_) {
0180   if (!failed_ && massBTopQuark_ < 0)
0181     operator2(jets, VbtagWP, btagWP_);
0182   return massBTopQuark_;
0183 }
0184
0185 double Calculate_miniAOD::tmassWBoson(pat::Muon* mu, const pat::MET& met, const pat::Jet& b) {
0186   if (tmassWBoson_ < 0)
0187     operator()(b, mu, met);
0188   return tmassWBoson_;
0189 }
0190
0191 double Calculate_miniAOD::tmassWBoson(pat::Electron* mu, const pat::MET& met, const pat::Jet& b) {
0192   if (tmassWBoson_ < 0)
0193     operator()(b, mu, met);
0194   return tmassWBoson_;
0195 }
0196
0197 double Calculate_miniAOD::tmassTopQuark(pat::Electron* lepton, const pat::MET& met, const pat::Jet& b) {
0198   if (tmassTopQuark_ < 0)
0199     operator()(b, lepton, met);
0200   return tmassTopQuark_;
0201 }
0202
0203 double Calculate_miniAOD::tmassTopQuark(pat::Muon* lepton, const pat::MET& met, const pat::Jet& b) {
0204   if (tmassTopQuark_ < 0)
0205     operator()(b, lepton, met);
0206   return tmassTopQuark_;
0207 }
0208
0209 void Calculate_miniAOD::operator()(const pat::Jet& bJet, pat::Muon* lepton, const pat::MET& met) {
0210   double metT = sqrt(pow(met.px(), 2) + pow(met.py(), 2));
0211   double lepT = sqrt(pow(lepton->px(), 2) + pow(lepton->py(), 2));
0212   double bT = sqrt(pow(bJet.px(), 2) + pow(bJet.py(), 2));
0213   reco::Particle::LorentzVector WT = lepton->p4() + met.p4();
0214   tmassWBoson_ = sqrt(pow(metT + lepT, 2) - (WT.px() * WT.px()) - (WT.py() * WT.py()));
0215   reco::Particle::LorentzVector topT = WT + bJet.p4();
0216   tmassTopQuark_ = sqrt(pow((metT + lepT + bT), 2) - (topT.px() * topT.px()) - (topT.py() * topT.py()));
0217 }
0218
0219 void Calculate_miniAOD::operator()(const pat::Jet& bJet, pat::Electron* lepton, const pat::MET& met) {
0220   double metT = sqrt(pow(met.px(), 2) + pow(met.py(), 2));
0221   double lepT = sqrt(pow(lepton->px(), 2) + pow(lepton->py(), 2));
0222   double bT = sqrt(pow(bJet.px(), 2) + pow(bJet.py(), 2));
0223   reco::Particle::LorentzVector WT = lepton->p4() + met.p4();
0224   tmassWBoson_ = sqrt(pow(metT + lepT, 2) - (WT.px() * WT.px()) - (WT.py() * WT.py()));
0225   reco::Particle::LorentzVector topT = WT + bJet.p4();
0226   tmassTopQuark_ = sqrt(pow((metT + lepT + bT), 2) - (topT.px() * topT.px()) - (topT.py() * topT.py()));
0227 }
0228
0229 void Calculate_miniAOD::operator()(const std::vector<pat::Jet>& jets) {
0230   if (maxNJets_ < 0)
0231     maxNJets_ = jets.size();
0232   failed_ = jets.size() < (unsigned int)maxNJets_;
0233   if (failed_) {
0234     return;
0235   }
0236
0237   // associate those jets with maximum pt of the vectorial
0238   // sum to the hadronic decay chain
0239   double maxPt = -1.;
0240   std::vector<int> maxPtIndices;
0241   maxPtIndices.push_back(-1);
0242   maxPtIndices.push_back(-1);
0243   maxPtIndices.push_back(-1);
0244
0245   for (int idx = 0; idx < maxNJets_; ++idx) {
0246     for (int jdx = idx + 1; jdx < maxNJets_; ++jdx) {
0247       //if (jdx <= idx) continue;
0248       for (int kdx = 0; kdx < maxNJets_; ++kdx) {
0249         if (kdx == idx || kdx == jdx)
0250           continue;
0251         reco::Particle::LorentzVector sum = jets[idx].p4() + jets[jdx].p4() + jets[kdx].p4();
0252         if (maxPt < 0. || maxPt < sum.pt()) {
0253           maxPt = sum.pt();
0254           maxPtIndices.clear();
0255           maxPtIndices.push_back(idx);
0256           maxPtIndices.push_back(jdx);
0257           maxPtIndices.push_back(kdx);
0258         }
0259       }
0260     }
0261   }
0262   massTopQuark_ = (jets[maxPtIndices[0]].p4() + jets[maxPtIndices[1]].p4() + jets[maxPtIndices[2]].p4()).mass();
0263
0264   // associate those jets that get closest to the W mass
0265   // with their invariant mass to the W boson
0266   double wDist = -1.;
0267   std::vector<int> wMassIndices;
0268   wMassIndices.push_back(-1);
0269   wMassIndices.push_back(-1);
0270   for (unsigned idx = 0; idx < maxPtIndices.size(); ++idx) {
0271     for (unsigned jdx = 0; jdx < maxPtIndices.size(); ++jdx) {
0272       if (jdx == idx || maxPtIndices[idx] > maxPtIndices[jdx])
0273         continue;
0274       reco::Particle::LorentzVector sum = jets[maxPtIndices[idx]].p4() + jets[maxPtIndices[jdx]].p4();
0275       if (wDist < 0. || wDist > fabs(sum.mass() - wMass_)) {
0276         wDist = fabs(sum.mass() - wMass_);
0277         wMassIndices.clear();
0278         wMassIndices.push_back(maxPtIndices[idx]);
0279         wMassIndices.push_back(maxPtIndices[jdx]);
0280       }
0281     }
0282   }
0283   massWBoson_ = (jets[wMassIndices[0]].p4() + jets[wMassIndices[1]].p4()).mass();
0284 }
0285
0286 void Calculate_miniAOD::operator2(const std::vector<pat::Jet>& jets, std::vector<double> bjet, double btagWP) {
0287   if (maxNJets_ < 0)
0288     maxNJets_ = jets.size();
0289   failed_ = jets.size() < (unsigned int)maxNJets_;
0290
0291   if (failed_) {
0292     return;
0293   }
0294
0295   if (jets.size() != bjet.size()) {
0296     return;
0297   }
0298
0299   // associate those jets with maximum pt of the vectorial
0300   // sum to the hadronic decay chain. Require ONLY 1 btagged jet
0301   double maxBPt = -1.;
0302   std::vector<int> maxBPtIndices;
0303   maxBPtIndices.push_back(-1);
0304   maxBPtIndices.push_back(-1);
0305   maxBPtIndices.push_back(-1);
0306   for (int idx = 0; idx < maxNJets_; ++idx) {
0307     for (int jdx = idx + 1; jdx < maxNJets_; ++jdx) {
0308       //if (jdx <= idx) continue;
0309       for (int kdx = 0; kdx < maxNJets_; ++kdx) {
0310         if (kdx == idx || kdx == jdx)
0311           continue;
0312         // require only 1b-jet
0313         if ((bjet[idx] > btagWP && bjet[jdx] <= btagWP && bjet[kdx] <= btagWP) ||
0314             (bjet[idx] <= btagWP && bjet[jdx] > btagWP && bjet[kdx] <= btagWP) ||
0315             (bjet[idx] <= btagWP && bjet[jdx] <= btagWP && bjet[kdx] > btagWP)) {
0316           reco::Particle::LorentzVector sum = jets[idx].p4() + jets[jdx].p4() + jets[kdx].p4();
0317           if (maxBPt < 0. || maxBPt < sum.pt()) {
0318             maxBPt = sum.pt();
0319             maxBPtIndices.clear();
0320             maxBPtIndices.push_back(idx);
0321             maxBPtIndices.push_back(jdx);
0322             maxBPtIndices.push_back(kdx);
0323           }
0324         }
0325       }
0326     }
0327   }
0328   if (maxBPtIndices[0] < 0 || maxBPtIndices[1] < 0 || maxBPtIndices[2] < 0)
0329     return;
0330   massBTopQuark_ = (jets[maxBPtIndices[0]].p4() + jets[maxBPtIndices[1]].p4() + jets[maxBPtIndices[2]].p4()).mass();
0331 }