Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:14 

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