Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-25 05:52:09

0001 #include "AnalysisDataFormats/TopObjects/interface/TtDilepEvtSolution.h"
0002 #include "DataFormats/Math/interface/deltaR.h"
0003 #include "DataFormats/Candidate/interface/Candidate.h"
0004 
0005 #include "FWCore/Framework/interface/EDProducer.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010 
0011 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullLepKinSolver.h"
0012 #include "TopQuarkAnalysis/TopEventSelection/interface/TtDilepLRSignalSelObservables.h"
0013 
0014 #include <memory>
0015 #include <string>
0016 #include <vector>
0017 
0018 class TtDilepEvtSolutionMaker : public edm::EDProducer {
0019 public:
0020   explicit TtDilepEvtSolutionMaker(const edm::ParameterSet& iConfig);
0021   ~TtDilepEvtSolutionMaker() override;
0022 
0023   void beginJob() override;
0024   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0025 
0026 private:
0027   // next methods are avoidable but they make the code legible
0028   inline bool PTComp(const reco::Candidate*, const reco::Candidate*) const;
0029   inline bool LepDiffCharge(const reco::Candidate*, const reco::Candidate*) const;
0030   inline bool HasPositiveCharge(const reco::Candidate*) const;
0031 
0032 private:
0033   edm::EDGetTokenT<std::vector<pat::Electron> > electronSourceToken_;
0034   edm::EDGetTokenT<std::vector<pat::Muon> > muonSourceToken_;
0035   edm::EDGetTokenT<std::vector<pat::Tau> > tauSourceToken_;
0036   edm::EDGetTokenT<std::vector<pat::MET> > metSourceToken_;
0037   edm::EDGetTokenT<std::vector<pat::Jet> > jetSourceToken_;
0038   edm::EDGetTokenT<TtGenEvent> evtSourceToken_;
0039   int jetCorrScheme_;
0040   unsigned int nrCombJets_;
0041   bool matchToGenEvt_, calcTopMass_, useMCforBest_;
0042   bool eeChannel_, emuChannel_, mumuChannel_, etauChannel_, mutauChannel_, tautauChannel_;
0043   double tmassbegin_, tmassend_, tmassstep_;
0044   std::vector<double> nupars_;
0045 
0046   TtDilepLRSignalSelObservables* myLRSignalSelObservables;
0047   TtFullLepKinSolver* solver;
0048 };
0049 
0050 inline bool TtDilepEvtSolutionMaker::PTComp(const reco::Candidate* l1, const reco::Candidate* l2) const {
0051   return (l1->pt() > l2->pt());
0052 }
0053 
0054 inline bool TtDilepEvtSolutionMaker::LepDiffCharge(const reco::Candidate* l1, const reco::Candidate* l2) const {
0055   return (l1->charge() != l2->charge());
0056 }
0057 
0058 inline bool TtDilepEvtSolutionMaker::HasPositiveCharge(const reco::Candidate* l) const { return (l->charge() > 0); }
0059 
0060 TtDilepEvtSolutionMaker::TtDilepEvtSolutionMaker(const edm::ParameterSet& iConfig) {
0061   // configurables
0062   electronSourceToken_ = consumes<std::vector<pat::Electron> >(iConfig.getParameter<edm::InputTag>("electronSource"));
0063   muonSourceToken_ = consumes<std::vector<pat::Muon> >(iConfig.getParameter<edm::InputTag>("muonSource"));
0064   tauSourceToken_ = consumes<std::vector<pat::Tau> >(iConfig.getParameter<edm::InputTag>("tauSource"));
0065   metSourceToken_ = consumes<std::vector<pat::MET> >(iConfig.getParameter<edm::InputTag>("metSource"));
0066   jetSourceToken_ = consumes<std::vector<pat::Jet> >(iConfig.getParameter<edm::InputTag>("jetSource"));
0067   jetCorrScheme_ = iConfig.getParameter<int>("jetCorrectionScheme");
0068   evtSourceToken_ = mayConsume<TtGenEvent>(iConfig.getParameter<edm::InputTag>("evtSource"));
0069   nrCombJets_ = iConfig.getParameter<unsigned int>("nrCombJets");
0070   matchToGenEvt_ = iConfig.getParameter<bool>("matchToGenEvt");
0071   calcTopMass_ = iConfig.getParameter<bool>("calcTopMass");
0072   useMCforBest_ = iConfig.getParameter<bool>("bestSolFromMC");
0073   eeChannel_ = iConfig.getParameter<bool>("eeChannel");
0074   emuChannel_ = iConfig.getParameter<bool>("emuChannel");
0075   mumuChannel_ = iConfig.getParameter<bool>("mumuChannel");
0076   mutauChannel_ = iConfig.getParameter<bool>("mutauChannel");
0077   etauChannel_ = iConfig.getParameter<bool>("etauChannel");
0078   tautauChannel_ = iConfig.getParameter<bool>("tautauChannel");
0079   tmassbegin_ = iConfig.getParameter<double>("tmassbegin");
0080   tmassend_ = iConfig.getParameter<double>("tmassend");
0081   tmassstep_ = iConfig.getParameter<double>("tmassstep");
0082   nupars_ = iConfig.getParameter<std::vector<double> >("neutrino_parameters");
0083 
0084   // define what will be produced
0085   produces<std::vector<TtDilepEvtSolution> >();
0086 
0087   myLRSignalSelObservables = new TtDilepLRSignalSelObservables(consumesCollector(), jetSourceToken_);
0088 }
0089 
0090 TtDilepEvtSolutionMaker::~TtDilepEvtSolutionMaker() {}
0091 
0092 void TtDilepEvtSolutionMaker::beginJob() {
0093   solver = new TtFullLepKinSolver(tmassbegin_, tmassend_, tmassstep_, nupars_);
0094 }
0095 
0096 void TtDilepEvtSolutionMaker::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0097   edm::Handle<std::vector<pat::Tau> > taus;
0098   iEvent.getByToken(tauSourceToken_, taus);
0099   edm::Handle<std::vector<pat::Muon> > muons;
0100   iEvent.getByToken(muonSourceToken_, muons);
0101   edm::Handle<std::vector<pat::Electron> > electrons;
0102   iEvent.getByToken(electronSourceToken_, electrons);
0103   edm::Handle<std::vector<pat::MET> > mets;
0104   iEvent.getByToken(metSourceToken_, mets);
0105   edm::Handle<std::vector<pat::Jet> > jets;
0106   iEvent.getByToken(jetSourceToken_, jets);
0107 
0108   int selMuonp = -1, selMuonm = -1;
0109   int selElectronp = -1, selElectronm = -1;
0110   int selTaup = -1, selTaum = -1;
0111   bool leptonFound = false;
0112   bool mumu = false;
0113   bool emu = false;
0114   bool ee = false;
0115   bool etau = false;
0116   bool mutau = false;
0117   bool tautau = false;
0118   bool leptonFoundEE = false;
0119   bool leptonFoundMM = false;
0120   bool leptonFoundTT = false;
0121   bool leptonFoundEpMm = false;
0122   bool leptonFoundEmMp = false;
0123   bool leptonFoundEpTm = false;
0124   bool leptonFoundEmTp = false;
0125   bool leptonFoundMpTm = false;
0126   bool leptonFoundMmTp = false;
0127   bool jetsFound = false;
0128   bool METFound = false;
0129   std::vector<int> JetVetoByTaus;
0130 
0131   //select MET (TopMET vector is sorted on ET)
0132   if (!mets->empty()) {
0133     METFound = true;
0134   }
0135 
0136   // If we have electrons and muons available,
0137   // build a solutions with electrons and muons.
0138   if (muons->size() + electrons->size() >= 2) {
0139     // select leptons
0140     if (electrons->empty())
0141       mumu = true;
0142     else if (muons->empty())
0143       ee = true;
0144     else if (electrons->size() == 1) {
0145       if (muons->size() == 1)
0146         emu = true;
0147       else if (PTComp(&(*electrons)[0], &(*muons)[1]))
0148         emu = true;
0149       else
0150         mumu = true;
0151     } else if (electrons->size() > 1) {
0152       if (PTComp(&(*electrons)[1], &(*muons)[0]))
0153         ee = true;
0154       else if (muons->size() == 1)
0155         emu = true;
0156       else if (PTComp(&(*electrons)[0], &(*muons)[1]))
0157         emu = true;
0158       else
0159         mumu = true;
0160     }
0161     if (ee) {
0162       if (LepDiffCharge(&(*electrons)[0], &(*electrons)[1])) {
0163         leptonFound = true;
0164         leptonFoundEE = true;
0165         if (HasPositiveCharge(&(*electrons)[0])) {
0166           selElectronp = 0;
0167           selElectronm = 1;
0168         } else {
0169           selElectronp = 1;
0170           selElectronm = 0;
0171         }
0172       }
0173     } else if (emu) {
0174       if (LepDiffCharge(&(*electrons)[0], &(*muons)[0])) {
0175         leptonFound = true;
0176         if (HasPositiveCharge(&(*electrons)[0])) {
0177           leptonFoundEpMm = true;
0178           selElectronp = 0;
0179           selMuonm = 0;
0180         } else {
0181           leptonFoundEmMp = true;
0182           selMuonp = 0;
0183           selElectronm = 0;
0184         }
0185       }
0186     } else if (mumu) {
0187       if (LepDiffCharge(&(*muons)[0], &(*muons)[1])) {
0188         leptonFound = true;
0189         leptonFoundMM = true;
0190         if (HasPositiveCharge(&(*muons)[0])) {
0191           selMuonp = 0;
0192           selMuonm = 1;
0193         } else {
0194           selMuonp = 1;
0195           selMuonm = 0;
0196         }
0197       }
0198     }
0199     //select Jets (TopJet vector is sorted on ET)
0200     if (jets->size() >= 2) {
0201       jetsFound = true;
0202     }
0203   }
0204   // If a tau is needed to have two leptons, then only consider the taus.
0205   // This is the minimal modification of the dilept selection that includes taus,
0206   // since we are considering taus only when no other solution exist.
0207   else if (muons->size() + electrons->size() == 1 && !taus->empty()) {
0208     // select leptons
0209     if (muons->size() == 1) {
0210       mutau = true;
0211       // depending on the muon charge, set the right muon index and specify channel
0212       int expectedCharge = -muons->begin()->charge();
0213       int* tauIdx = nullptr;
0214       if (expectedCharge < 0) {
0215         selMuonp = 0;
0216         tauIdx = &selTaum;
0217         leptonFoundMpTm = true;
0218       } else {
0219         selMuonm = 0;
0220         tauIdx = &selTaup;
0221         leptonFoundMmTp = true;
0222       }
0223       // loop over the vector of taus to find the ones
0224       // that have the charge opposite to the muon one, and do not match in eta-phi
0225       std::vector<std::vector<pat::Tau>::const_iterator> subset1;
0226       for (std::vector<pat::Tau>::const_iterator tau = taus->begin(); tau < taus->end(); ++tau) {
0227         if (tau->charge() * expectedCharge >= 0 && DeltaR<pat::Particle>()(*tau, *(muons->begin())) > 0.1) {
0228           *tauIdx = tau - taus->begin();
0229           leptonFound = true;
0230           subset1.push_back(tau);
0231         }
0232       }
0233       // if there are more than one tau with ecalIsol==0, take the smallest E/P
0234       float iso = 999.;
0235       for (std::vector<std::vector<pat::Tau>::const_iterator>::const_iterator tau = subset1.begin();
0236            tau < subset1.end();
0237            ++tau) {
0238         if ((*tau)->isPFTau() && (*tau)->isolationPFChargedHadrCandsPtSum() < iso) {
0239           *tauIdx = *tau - taus->begin();
0240           iso = (*tau)->isolationPFChargedHadrCandsPtSum();
0241         }
0242       }
0243 
0244       // check that one combination has been found
0245       if (!leptonFound) {
0246         leptonFoundMpTm = false;
0247         leptonFoundMmTp = false;
0248       }
0249       // discard the jet that matches the tau (if one)
0250       if (leptonFound) {
0251         for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet < jets->end(); ++jet) {
0252           if (DeltaR<pat::Particle, pat::Jet>()(*(taus->begin() + *tauIdx), *jet) < 0.1) {
0253             JetVetoByTaus.push_back(jet - jets->begin());
0254           }
0255         }
0256       }
0257     } else {
0258       etau = true;
0259       // depending on the electron charge, set the right electron index and specify channel
0260       int expectedCharge = -electrons->begin()->charge();
0261       int* tauIdx = nullptr;
0262       if (expectedCharge < 0) {
0263         selElectronp = 0;
0264         tauIdx = &selTaum;
0265         leptonFoundEpTm = true;
0266       } else {
0267         selElectronm = 0;
0268         tauIdx = &selTaup;
0269         leptonFoundEmTp = true;
0270       }
0271       // loop over the vector of taus to find the ones
0272       // that have the charge opposite to the muon one, and do not match in eta-phi
0273       std::vector<std::vector<pat::Tau>::const_iterator> subset1;
0274       for (std::vector<pat::Tau>::const_iterator tau = taus->begin(); tau < taus->end(); ++tau) {
0275         if (tau->charge() * expectedCharge >= 0 && DeltaR<pat::Particle>()(*tau, *(electrons->begin())) > 0.1) {
0276           *tauIdx = tau - taus->begin();
0277           leptonFound = true;
0278           subset1.push_back(tau);
0279         }
0280       }
0281       // if there are more than one tau with ecalIsol==0, take the smallest E/P
0282       float iso = 999.;
0283       for (std::vector<std::vector<pat::Tau>::const_iterator>::const_iterator tau = subset1.begin();
0284            tau < subset1.end();
0285            ++tau) {
0286         if ((*tau)->isPFTau() && (*tau)->isolationPFChargedHadrCandsPtSum() < iso) {
0287           *tauIdx = *tau - taus->begin();
0288           iso = (*tau)->isolationPFChargedHadrCandsPtSum();
0289         }
0290       }
0291 
0292       // check that one combination has been found
0293       if (!leptonFound) {
0294         leptonFoundEpTm = false;
0295         leptonFoundEmTp = false;
0296       }
0297       // discard the jet that matches the tau (if one)
0298       if (leptonFound) {
0299         for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet < jets->end(); ++jet) {
0300           if (DeltaR<pat::Particle, pat::Jet>()(*(taus->begin() + *tauIdx), *jet) < 0.1) {
0301             JetVetoByTaus.push_back(jet - jets->begin());
0302           }
0303         }
0304       }
0305     }
0306     // select Jets (TopJet vector is sorted on ET)
0307     jetsFound = ((jets->size() - JetVetoByTaus.size()) >= 2);
0308   } else if (taus->size() > 1) {
0309     tautau = true;
0310     if (LepDiffCharge(&(*taus)[0], &(*taus)[1])) {
0311       leptonFound = true;
0312       leptonFoundTT = true;
0313       if (HasPositiveCharge(&(*taus)[0])) {
0314         selTaup = 0;
0315         selTaum = 1;
0316       } else {
0317         selTaup = 1;
0318         selTaum = 0;
0319       }
0320     }
0321     for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet < jets->end(); ++jet) {
0322       if (DeltaR<pat::Particle, pat::Jet>()((*taus)[0], *jet) < 0.1 ||
0323           DeltaR<pat::Particle, pat::Jet>()((*taus)[1], *jet) < 0.1) {
0324         JetVetoByTaus.push_back(jet - jets->begin());
0325       }
0326     }
0327     // select Jets (TopJet vector is sorted on ET)
0328     jetsFound = ((jets->size() - JetVetoByTaus.size()) >= 2);
0329   }
0330 
0331   // Check that the above work makes sense
0332   if (int(ee) + int(emu) + int(mumu) + int(etau) + int(mutau) + int(tautau) > 1)
0333     std::cout << "[TtDilepEvtSolutionMaker]: "
0334               << "Lepton selection criteria uncorrectly defined" << std::endl;
0335 
0336   bool correctLepton = (leptonFoundEE && eeChannel_) || ((leptonFoundEmMp || leptonFoundEpMm) && emuChannel_) ||
0337                        (leptonFoundMM && mumuChannel_) || ((leptonFoundMmTp || leptonFoundMpTm) && mutauChannel_) ||
0338                        ((leptonFoundEmTp || leptonFoundEpTm) && etauChannel_) || (leptonFoundTT && tautauChannel_);
0339 
0340   std::vector<TtDilepEvtSolution>* evtsols = new std::vector<TtDilepEvtSolution>();
0341   if (correctLepton && METFound && jetsFound) {
0342     // protect against reading beyond array boundaries while discounting vetoed jets
0343     unsigned int nrCombJets = 0;
0344     unsigned int numberOfJets = 0;
0345     for (; nrCombJets < jets->size() && numberOfJets < nrCombJets_; ++nrCombJets) {
0346       if (find(JetVetoByTaus.begin(), JetVetoByTaus.end(), int(nrCombJets)) == JetVetoByTaus.end())
0347         ++numberOfJets;
0348     }
0349     // consider all permutations
0350     for (unsigned int ib = 0; ib < nrCombJets; ib++) {
0351       // skipped jet vetoed during components-flagging.
0352       if (find(JetVetoByTaus.begin(), JetVetoByTaus.end(), int(ib)) != JetVetoByTaus.end())
0353         continue;
0354       // second loop of the permutations
0355       for (unsigned int ibbar = 0; ibbar < nrCombJets; ibbar++) {
0356         // avoid the diagonal: b and bbar must be distinct jets
0357         if (ib == ibbar)
0358           continue;
0359         // skipped jet vetoed during components-flagging.
0360         if (find(JetVetoByTaus.begin(), JetVetoByTaus.end(), int(ibbar)) != JetVetoByTaus.end())
0361           continue;
0362         // Build and save a solution
0363         TtDilepEvtSolution asol;
0364         asol.setJetCorrectionScheme(jetCorrScheme_);
0365         double xconstraint = 0, yconstraint = 0;
0366         // Set e+ in the event
0367         if (leptonFoundEE || leptonFoundEpMm || leptonFoundEpTm) {
0368           asol.setElectronp(electrons, selElectronp);
0369           xconstraint += (*electrons)[selElectronp].px();
0370           yconstraint += (*electrons)[selElectronp].py();
0371         }
0372         // Set e- in the event
0373         if (leptonFoundEE || leptonFoundEmMp || leptonFoundEmTp) {
0374           asol.setElectronm(electrons, selElectronm);
0375           xconstraint += (*electrons)[selElectronm].px();
0376           yconstraint += (*electrons)[selElectronm].py();
0377         }
0378         // Set mu+ in the event
0379         if (leptonFoundMM || leptonFoundEmMp || leptonFoundMpTm) {
0380           asol.setMuonp(muons, selMuonp);
0381           xconstraint += (*muons)[selMuonp].px();
0382           yconstraint += (*muons)[selMuonp].py();
0383         }
0384         // Set mu- in the event
0385         if (leptonFoundMM || leptonFoundEpMm || leptonFoundMmTp) {
0386           asol.setMuonm(muons, selMuonm);
0387           xconstraint += (*muons)[selMuonm].px();
0388           yconstraint += (*muons)[selMuonm].py();
0389         }
0390         // Set tau- in the event
0391         if (leptonFoundEpTm || leptonFoundMpTm || leptonFoundTT) {
0392           asol.setTaum(taus, selTaum);
0393           xconstraint += (*taus)[selTaum].px();
0394           yconstraint += (*taus)[selTaum].py();
0395         }
0396         // Set tau+ in the event
0397         if (leptonFoundEmTp || leptonFoundMmTp || leptonFoundTT) {
0398           asol.setTaup(taus, selTaup);
0399           xconstraint += (*taus)[selTaup].px();
0400           yconstraint += (*taus)[selTaup].py();
0401         }
0402         // Set Jets/MET in the event
0403         asol.setB(jets, ib);
0404         asol.setBbar(jets, ibbar);
0405         asol.setMET(mets, 0);
0406         xconstraint += (*jets)[ib].px() + (*jets)[ibbar].px() + (*mets)[0].px();
0407         yconstraint += (*jets)[ib].py() + (*jets)[ibbar].py() + (*mets)[0].py();
0408         // if asked for, match the event solutions to the gen Event
0409         if (matchToGenEvt_) {
0410           edm::Handle<TtGenEvent> genEvt;
0411           iEvent.getByToken(evtSourceToken_, genEvt);
0412           asol.setGenEvt(genEvt);
0413         }
0414         // If asked, use the kin fitter to compute the top mass
0415         if (calcTopMass_) {
0416           solver->SetConstraints(xconstraint, yconstraint);
0417           solver->useWeightFromMC(useMCforBest_);
0418           asol = solver->addKinSolInfo(&asol);
0419         }
0420 
0421         // these lines calculate the observables to be used in the TtDilepSignalSelection LR
0422         (*myLRSignalSelObservables)(asol, iEvent);
0423 
0424         evtsols->push_back(asol);
0425       }
0426     }
0427     // flag the best solution (MC matching)
0428     if (matchToGenEvt_) {
0429       double bestSolDR = 9999.;
0430       int bestSol = -1;
0431       double dR = 0.;
0432       for (size_t s = 0; s < evtsols->size(); s++) {
0433         dR = (*evtsols)[s].getJetResidual();
0434         if (dR < bestSolDR) {
0435           bestSolDR = dR;
0436           bestSol = s;
0437         }
0438       }
0439       if (bestSol != -1)
0440         (*evtsols)[bestSol].setBestSol(true);
0441     }
0442     // put the result in the event
0443     std::unique_ptr<std::vector<TtDilepEvtSolution> > pOut(evtsols);
0444     iEvent.put(std::move(pOut));
0445   } else {
0446     // no solution: put a dummy solution in the event
0447     TtDilepEvtSolution asol;
0448     evtsols->push_back(asol);
0449     std::unique_ptr<std::vector<TtDilepEvtSolution> > pOut(evtsols);
0450     iEvent.put(std::move(pOut));
0451   }
0452 }
0453 
0454 #include "FWCore/Framework/interface/MakerMacros.h"
0455 DEFINE_FWK_MODULE(TtDilepEvtSolutionMaker);