TopEventSelection/src/TtDilepLRSignalSelObservables.cc

0001 #include "TopQuarkAnalysis/TopEventSelection/interface/TtDilepLRSignalSelObservables.h"
0002 #include "DataFormats/PatCandidates/interface/Jet.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005 #include "DataFormats/Common/interface/Handle.h"
0006 #include "AnalysisDataFormats/TopObjects/interface/TtGenEvent.h"
0007
0008 /************** Definition of the functions of the class ***************/
0009
0010 //Constructor
0011 TtDilepLRSignalSelObservables::TtDilepLRSignalSelObservables(
0012     edm::ConsumesCollector &&iC, const edm::EDGetTokenT<std::vector<pat::Jet> > &jetSourceToken)
0013     : jetSourceToken_(jetSourceToken), genEvtToken_(iC.consumes<TtGenEvent>(edm::InputTag("genEvt"))) {
0014   count1 = 0;
0015   count2 = 0;
0016   count3 = 0;
0017   count4 = 0;
0018   count5 = 0;
0019   count3 = 0;
0020 }
0021
0022 // Destructor
0023 TtDilepLRSignalSelObservables::~TtDilepLRSignalSelObservables() {}
0024
0025 std::vector<TtDilepLRSignalSelObservables::IntBoolPair> TtDilepLRSignalSelObservables::operator()(
0026     TtDilepEvtSolution &solution, const edm::Event &iEvent, bool matchOnly) {
0027   evtselectVarVal.clear();
0028   evtselectVarMatch.clear();
0029
0030   // Check whether the objects are matched:
0031   bool matchB1 = false;
0032   bool matchB2 = false;
0033   bool matchB = false;
0034   bool matchBbar = false;
0035   bool matchLeptPos = false;
0036   bool matchLeptNeg = false;
0037
0038   edm::Handle<TtGenEvent> genEvent;
0039   iEvent.getByToken(genEvtToken_, genEvent);
0040
0041   if (!genEvent.failedToGet() && genEvent.isValid()) {
0042     // std::cout <<std::endl;
0043     double dr, dr1, dr2;
0044
0045     if (genEvent->isFullLeptonic()) {
0046       // Match the leptons, by type and deltaR
0047       dr = DeltaR<reco::Particle, reco::GenParticle>()(solution.getLeptPos(), *(solution.getGenLepp()));
0048       matchLeptPos = ((((solution.getWpDecay() == "electron") && (std::abs(solution.getGenLepp()->pdgId()) == 11)) ||
0049                        ((solution.getWpDecay() == "muon") && (std::abs(solution.getGenLepp()->pdgId()) == 13))) &&
0050                       (dr < 0.1));
0051
0052       dr = DeltaR<reco::Particle, reco::GenParticle>()(solution.getLeptNeg(), *(solution.getGenLepm()));
0053       matchLeptNeg = ((((solution.getWmDecay() == "electron") && (std::abs(solution.getGenLepm()->pdgId()) == 11)) ||
0054                        ((solution.getWmDecay() == "muon") && (std::abs(solution.getGenLepm()->pdgId()) == 13))) &&
0055                       (dr < 0.1));
0056     }
0057
0058     if (genEvent->isSemiLeptonic()) {
0059       int id = genEvent->singleLepton()->pdgId();
0060
0061       if (id > 0) {
0062         dr = DeltaR<reco::Particle, reco::GenParticle>()(solution.getLeptNeg(), *(genEvent->singleLepton()));
0063         matchLeptNeg = ((((solution.getWmDecay() == "electron") && (id == 11)) ||
0064                          ((solution.getWmDecay() == "muon") && (id == 13))) &&
0065                         (dr < 0.1));
0066       } else {
0067         dr = DeltaR<reco::Particle, reco::GenParticle>()(solution.getLeptPos(), *(genEvent->singleLepton()));
0068         matchLeptPos = ((((solution.getWpDecay() == "electron") && (id == -11)) ||
0069                          ((solution.getWpDecay() == "muon") && (id == -13))) &&
0070                         (dr < 0.1));
0071       }
0072     }
0073
0074     if (genEvent->isTtBar() && genEvent->numberOfBQuarks() > 1) {
0075       if (solution.getJetB().partonFlavour() == 5)
0076         ++count1;
0077       if (solution.getJetBbar().partonFlavour() == 5)
0078         ++count1;
0079
0080       dr1 = DeltaR<pat::Jet, reco::GenParticle>()(solution.getCalJetB(), *(genEvent->b()));
0081       dr2 = DeltaR<pat::Jet, reco::GenParticle>()(solution.getCalJetB(), *(genEvent->bBar()));
0082
0083       matchB1 = ((dr1 < 0.4) || (dr2 < 0.4));
0084       matchB = ((solution.getJetB().partonFlavour() == 5) && (dr1 < 0.4));
0085       if (matchB)
0086         ++count3;
0087       matchB = ((dr1 < 0.4));
0088       if (dr1 < 0.5)
0089         ++count2;
0090       if (dr1 < 0.4)
0091         ++count4;
0092       if (dr1 < 0.3)
0093         ++count5;
0094
0095       dr1 = DeltaR<pat::Jet, reco::GenParticle>()(solution.getCalJetBbar(), *(genEvent->b()));
0096       dr2 = DeltaR<pat::Jet, reco::GenParticle>()(solution.getCalJetBbar(), *(genEvent->bBar()));
0097
0098       matchBbar = ((solution.getJetBbar().partonFlavour() == 5) && (dr2 < 0.4));
0099       if (matchBbar)
0100         ++count3;
0101       matchBbar = ((dr2 < 0.4));
0102       matchB2 = ((dr1 < 0.4) || (dr2 < 0.4));
0103       if (dr2 < 0.5)
0104         ++count2;
0105       if (dr2 < 0.4)
0106         ++count4;
0107       if (dr2 < 0.3)
0108         ++count5;
0109     }
0110   }
0111
0112   edm::Handle<std::vector<pat::Jet> > jets;
0113   iEvent.getByToken(jetSourceToken_, jets);
0114
0115   //  Lower / Higher of both jet angles
0116
0117   double v1 = std::abs(solution.getJetB().p4().theta() - M_PI / 2);
0118   double v2 = std::abs(solution.getJetBbar().p4().theta() - M_PI / 2);
0119   fillMinMax(v1, v2, 1, evtselectVarVal, matchB1, matchB2, evtselectVarMatch);
0120
0121   //  Lower / Higher of both jet pT
0122
0123   double pt1 = solution.getJetB().p4().pt();
0124   double pt2 = solution.getJetBbar().p4().pt();
0125   fillMinMax(pt1, pt2, 3, evtselectVarVal, matchB1, matchB2, evtselectVarMatch);
0126
0127   //  Lower / Higher of both lepton pT
0128
0129   pt1 = solution.getLeptPos().p4().pt();
0130   pt2 = solution.getLeptNeg().p4().pt();
0131   fillMinMax(pt1, pt2, 5, evtselectVarVal, matchLeptPos, matchLeptNeg, evtselectVarMatch);
0132
0133   // delta theta btw the b-jets
0134
0135   double deltaPhi = std::abs(delta(solution.getJetB().p4().phi(), solution.getJetBbar().p4().phi()));
0136
0137   evtselectVarVal.push_back(IntDblPair(7, deltaPhi));
0138   evtselectVarMatch.push_back(IntBoolPair(7, matchB1 && matchB2));
0139
0140   // delta phi btw the b-jets
0141
0142   double deltaTheta = std::abs(delta(solution.getJetBbar().p4().theta(), solution.getJetB().p4().theta()));
0143
0144   evtselectVarVal.push_back(IntDblPair(8, deltaTheta));
0145   evtselectVarMatch.push_back(IntBoolPair(8, matchB1 && matchB2));
0146
0147   //  Lower / Higher of phi difference between the b and associated lepton
0148
0149   double deltaPhi1 = std::abs(delta(solution.getJetB().p4().phi(), solution.getLeptPos().p4().phi()));
0150   double deltaPhi2 = std::abs(delta(solution.getJetBbar().p4().phi(), solution.getLeptNeg().p4().phi()));
0151
0152   fillMinMax(
0153       deltaPhi1, deltaPhi2, 9, evtselectVarVal, matchB && matchLeptPos, matchBbar && matchLeptNeg, evtselectVarMatch);
0154
0155   //  Lower / Higher of theta difference between the b and associated lepton
0156
0157   double deltaTheta1 = std::abs(solution.getJetB().p4().theta() - solution.getLeptPos().p4().theta());
0158   double deltaTheta2 = std::abs(solution.getJetBbar().p4().theta() - solution.getLeptNeg().p4().theta());
0159   fillMinMax(deltaTheta1,
0160              deltaTheta2,
0161              11,
0162              evtselectVarVal,
0163              matchB && matchLeptPos,
0164              matchBbar && matchLeptNeg,
0165              evtselectVarMatch);
0166
0167   // Invariant Mass of lepton pair
0168
0169   math::XYZTLorentzVector pp = solution.getLeptPos().p4() + solution.getLeptNeg().p4();
0170   double mass = pp.mass();
0171   evtselectVarVal.push_back(IntDblPair(13, mass));
0172   evtselectVarMatch.push_back(IntBoolPair(13, matchLeptNeg && matchLeptPos));
0173
0174   evtselectVarVal.push_back(IntDblPair(13, mass));
0175   evtselectVarMatch.push_back(IntBoolPair(13, matchLeptNeg && matchLeptPos));
0176
0177   std::vector<pat::Jet> jet3;
0178   for (int i = 0; i < (int)jets->size(); ++i) {
0179     if (((*jets)[i].et() < solution.getJetB().et()) && ((*jets)[i].et() < solution.getJetBbar().et())) {
0180       jet3.push_back((*jets)[i]);
0181     }
0182   }
0183   double jet1Ratio = 0., jet2Ratio = 0.;
0184   if (!jet3.empty()) {
0185     jet1Ratio = jet3[0].et() / solution.getJetB().et();
0186     jet2Ratio = jet3[0].et() / solution.getJetBbar().et();
0187   }
0188   fillMinMax(jet1Ratio, jet2Ratio, 14, evtselectVarVal, matchB1, matchB2, evtselectVarMatch);
0189
0190   evtselectVarVal.push_back(IntDblPair(16, jets->size()));
0191   evtselectVarMatch.push_back(IntBoolPair(16, matchB && matchBbar));
0192
0193   if (!matchOnly)
0194     solution.setLRSignalEvtObservables(evtselectVarVal);
0195   return evtselectVarMatch;
0196 }
0197
0198 void TtDilepLRSignalSelObservables::fillMinMax(double v1,
0199                                                double v2,
0200                                                int obsNbr,
0201                                                std::vector<IntDblPair> &varList,
0202                                                bool match1,
0203                                                bool match2,
0204                                                std::vector<IntBoolPair> &matchList) {
0205   if (v1 < v2) {
0206     varList.push_back(IntDblPair(obsNbr, v1));
0207     varList.push_back(IntDblPair(obsNbr + 1, v2));
0208     matchList.push_back(IntBoolPair(obsNbr, match1));
0209     matchList.push_back(IntBoolPair(obsNbr + 1, match2));
0210
0211   } else {
0212     varList.push_back(IntDblPair(obsNbr, v2));
0213     varList.push_back(IntDblPair(obsNbr + 1, v1));
0214     matchList.push_back(IntBoolPair(obsNbr, match2));
0215     matchList.push_back(IntBoolPair(obsNbr + 1, match1));
0216   }
0217 }
0218
0219 double TtDilepLRSignalSelObservables::delta(double phi1, double phi2) {
0220   double deltaPhi = phi1 - phi2;
0221   while (deltaPhi > M_PI)
0222     deltaPhi -= 2 * M_PI;
0223   while (deltaPhi <= -M_PI)
0224     deltaPhi += 2 * M_PI;
0225   return deltaPhi;
0226 }