Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-04 22:55:07 

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         leptonFoundEE = true;
0161         if (HasPositiveCharge(&(*electrons)[0])) {
0162           selElectronp = 0;
0163           selElectronm = 1;
0164         } else {
0165           selElectronp = 1;
0166           selElectronm = 0;
0167         }
0168       }
0169     } else if (emu) {
0170       if (LepDiffCharge(&(*electrons)[0], &(*muons)[0])) {
0171         if (HasPositiveCharge(&(*electrons)[0])) {
0172           leptonFoundEpMm = true;
0173           selElectronp = 0;
0174           selMuonm = 0;
0175         } else {
0176           leptonFoundEmMp = true;
0177           selMuonp = 0;
0178           selElectronm = 0;
0179         }
0180       }
0181     } else if (mumu) {
0182       if (LepDiffCharge(&(*muons)[0], &(*muons)[1])) {
0183         leptonFoundMM = true;
0184         if (HasPositiveCharge(&(*muons)[0])) {
0185           selMuonp = 0;
0186           selMuonm = 1;
0187         } else {
0188           selMuonp = 1;
0189           selMuonm = 0;
0190         }
0191       }
0192     }
0193     //select Jets (TopJet vector is sorted on ET)
0194     if (jets->size() >= 2) {
0195       jetsFound = true;
0196     }
0197   }
0198   // If a tau is needed to have two leptons, then only consider the taus.
0199   // This is the minimal modification of the dilept selection that includes taus,
0200   // since we are considering taus only when no other solution exist.
0201   else if (muons->size() + electrons->size() == 1 && !taus->empty()) {
0202     // select leptons
0203     if (muons->size() == 1) {
0204       mutau = true;
0205       // depending on the muon charge, set the right muon index and specify channel
0206       int expectedCharge = -muons->begin()->charge();
0207       int* tauIdx = nullptr;
0208       if (expectedCharge < 0) {
0209         selMuonp = 0;
0210         tauIdx = &selTaum;
0211         leptonFoundMpTm = true;
0212       } else {
0213         selMuonm = 0;
0214         tauIdx = &selTaup;
0215         leptonFoundMmTp = true;
0216       }
0217       // loop over the vector of taus to find the ones
0218       // that have the charge opposite to the muon one, and do not match in eta-phi
0219       std::vector<std::vector<pat::Tau>::const_iterator> subset1;
0220       for (std::vector<pat::Tau>::const_iterator tau = taus->begin(); tau < taus->end(); ++tau) {
0221         if (tau->charge() * expectedCharge >= 0 && DeltaR<pat::Particle>()(*tau, *(muons->begin())) > 0.1) {
0222           *tauIdx = tau - taus->begin();
0223           leptonFound = true;
0224           subset1.push_back(tau);
0225         }
0226       }
0227       // if there are more than one tau with ecalIsol==0, take the smallest E/P
0228       float iso = 999.;
0229       for (std::vector<std::vector<pat::Tau>::const_iterator>::const_iterator tau = subset1.begin();
0230            tau < subset1.end();
0231            ++tau) {
0232         if ((*tau)->isPFTau() && (*tau)->isolationPFChargedHadrCandsPtSum() < iso) {
0233           *tauIdx = *tau - taus->begin();
0234           iso = (*tau)->isolationPFChargedHadrCandsPtSum();
0235         }
0236       }
0237 
0238       // check that one combination has been found
0239       if (!leptonFound) {
0240         leptonFoundMpTm = false;
0241         leptonFoundMmTp = false;
0242       }
0243       // discard the jet that matches the tau (if one)
0244       if (leptonFound) {
0245         for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet < jets->end(); ++jet) {
0246           if (DeltaR<pat::Particle, pat::Jet>()(*(taus->begin() + *tauIdx), *jet) < 0.1) {
0247             JetVetoByTaus.push_back(jet - jets->begin());
0248           }
0249         }
0250       }
0251     } else {
0252       etau = true;
0253       // depending on the electron charge, set the right electron index and specify channel
0254       int expectedCharge = -electrons->begin()->charge();
0255       int* tauIdx = nullptr;
0256       if (expectedCharge < 0) {
0257         selElectronp = 0;
0258         tauIdx = &selTaum;
0259         leptonFoundEpTm = true;
0260       } else {
0261         selElectronm = 0;
0262         tauIdx = &selTaup;
0263         leptonFoundEmTp = true;
0264       }
0265       // loop over the vector of taus to find the ones
0266       // that have the charge opposite to the muon one, and do not match in eta-phi
0267       std::vector<std::vector<pat::Tau>::const_iterator> subset1;
0268       for (std::vector<pat::Tau>::const_iterator tau = taus->begin(); tau < taus->end(); ++tau) {
0269         if (tau->charge() * expectedCharge >= 0 && DeltaR<pat::Particle>()(*tau, *(electrons->begin())) > 0.1) {
0270           *tauIdx = tau - taus->begin();
0271           leptonFound = true;
0272           subset1.push_back(tau);
0273         }
0274       }
0275       // if there are more than one tau with ecalIsol==0, take the smallest E/P
0276       float iso = 999.;
0277       for (std::vector<std::vector<pat::Tau>::const_iterator>::const_iterator tau = subset1.begin();
0278            tau < subset1.end();
0279            ++tau) {
0280         if ((*tau)->isPFTau() && (*tau)->isolationPFChargedHadrCandsPtSum() < iso) {
0281           *tauIdx = *tau - taus->begin();
0282           iso = (*tau)->isolationPFChargedHadrCandsPtSum();
0283         }
0284       }
0285 
0286       // check that one combination has been found
0287       if (!leptonFound) {
0288         leptonFoundEpTm = false;
0289         leptonFoundEmTp = false;
0290       }
0291       // discard the jet that matches the tau (if one)
0292       if (leptonFound) {
0293         for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet < jets->end(); ++jet) {
0294           if (DeltaR<pat::Particle, pat::Jet>()(*(taus->begin() + *tauIdx), *jet) < 0.1) {
0295             JetVetoByTaus.push_back(jet - jets->begin());
0296           }
0297         }
0298       }
0299     }
0300     // select Jets (TopJet vector is sorted on ET)
0301     jetsFound = ((jets->size() - JetVetoByTaus.size()) >= 2);
0302   } else if (taus->size() > 1) {
0303     tautau = true;
0304     if (LepDiffCharge(&(*taus)[0], &(*taus)[1])) {
0305       leptonFoundTT = true;
0306       if (HasPositiveCharge(&(*taus)[0])) {
0307         selTaup = 0;
0308         selTaum = 1;
0309       } else {
0310         selTaup = 1;
0311         selTaum = 0;
0312       }
0313     }
0314     for (std::vector<pat::Jet>::const_iterator jet = jets->begin(); jet < jets->end(); ++jet) {
0315       if (DeltaR<pat::Particle, pat::Jet>()((*taus)[0], *jet) < 0.1 ||
0316           DeltaR<pat::Particle, pat::Jet>()((*taus)[1], *jet) < 0.1) {
0317         JetVetoByTaus.push_back(jet - jets->begin());
0318       }
0319     }
0320     // select Jets (TopJet vector is sorted on ET)
0321     jetsFound = ((jets->size() - JetVetoByTaus.size()) >= 2);
0322   }
0323 
0324   // Check that the above work makes sense
0325   if (int(ee) + int(emu) + int(mumu) + int(etau) + int(mutau) + int(tautau) > 1)
0326     std::cout << "[TtDilepEvtSolutionMaker]: "
0327               << "Lepton selection criteria uncorrectly defined" << std::endl;
0328 
0329   bool correctLepton = (leptonFoundEE && eeChannel_) || ((leptonFoundEmMp || leptonFoundEpMm) && emuChannel_) ||
0330                        (leptonFoundMM && mumuChannel_) || ((leptonFoundMmTp || leptonFoundMpTm) && mutauChannel_) ||
0331                        ((leptonFoundEmTp || leptonFoundEpTm) && etauChannel_) || (leptonFoundTT && tautauChannel_);
0332 
0333   std::vector<TtDilepEvtSolution>* evtsols = new std::vector<TtDilepEvtSolution>();
0334   if (correctLepton && METFound && jetsFound) {
0335     // protect against reading beyond array boundaries while discounting vetoed jets
0336     unsigned int nrCombJets = 0;
0337     unsigned int numberOfJets = 0;
0338     for (; nrCombJets < jets->size() && numberOfJets < nrCombJets_; ++nrCombJets) {
0339       if (find(JetVetoByTaus.begin(), JetVetoByTaus.end(), int(nrCombJets)) == JetVetoByTaus.end())
0340         ++numberOfJets;
0341     }
0342     // consider all permutations
0343     for (unsigned int ib = 0; ib < nrCombJets; ib++) {
0344       // skipped jet vetoed during components-flagging.
0345       if (find(JetVetoByTaus.begin(), JetVetoByTaus.end(), int(ib)) != JetVetoByTaus.end())
0346         continue;
0347       // second loop of the permutations
0348       for (unsigned int ibbar = 0; ibbar < nrCombJets; ibbar++) {
0349         // avoid the diagonal: b and bbar must be distinct jets
0350         if (ib == ibbar)
0351           continue;
0352         // skipped jet vetoed during components-flagging.
0353         if (find(JetVetoByTaus.begin(), JetVetoByTaus.end(), int(ibbar)) != JetVetoByTaus.end())
0354           continue;
0355         // Build and save a solution
0356         TtDilepEvtSolution asol;
0357         asol.setJetCorrectionScheme(jetCorrScheme_);
0358         double xconstraint = 0, yconstraint = 0;
0359         // Set e+ in the event
0360         if (leptonFoundEE || leptonFoundEpMm || leptonFoundEpTm) {
0361           asol.setElectronp(electrons, selElectronp);
0362           xconstraint += (*electrons)[selElectronp].px();
0363           yconstraint += (*electrons)[selElectronp].py();
0364         }
0365         // Set e- in the event
0366         if (leptonFoundEE || leptonFoundEmMp || leptonFoundEmTp) {
0367           asol.setElectronm(electrons, selElectronm);
0368           xconstraint += (*electrons)[selElectronm].px();
0369           yconstraint += (*electrons)[selElectronm].py();
0370         }
0371         // Set mu+ in the event
0372         if (leptonFoundMM || leptonFoundEmMp || leptonFoundMpTm) {
0373           asol.setMuonp(muons, selMuonp);
0374           xconstraint += (*muons)[selMuonp].px();
0375           yconstraint += (*muons)[selMuonp].py();
0376         }
0377         // Set mu- in the event
0378         if (leptonFoundMM || leptonFoundEpMm || leptonFoundMmTp) {
0379           asol.setMuonm(muons, selMuonm);
0380           xconstraint += (*muons)[selMuonm].px();
0381           yconstraint += (*muons)[selMuonm].py();
0382         }
0383         // Set tau- in the event
0384         if (leptonFoundEpTm || leptonFoundMpTm || leptonFoundTT) {
0385           asol.setTaum(taus, selTaum);
0386           xconstraint += (*taus)[selTaum].px();
0387           yconstraint += (*taus)[selTaum].py();
0388         }
0389         // Set tau+ in the event
0390         if (leptonFoundEmTp || leptonFoundMmTp || leptonFoundTT) {
0391           asol.setTaup(taus, selTaup);
0392           xconstraint += (*taus)[selTaup].px();
0393           yconstraint += (*taus)[selTaup].py();
0394         }
0395         // Set Jets/MET in the event
0396         asol.setB(jets, ib);
0397         asol.setBbar(jets, ibbar);
0398         asol.setMET(mets, 0);
0399         xconstraint += (*jets)[ib].px() + (*jets)[ibbar].px() + (*mets)[0].px();
0400         yconstraint += (*jets)[ib].py() + (*jets)[ibbar].py() + (*mets)[0].py();
0401         // if asked for, match the event solutions to the gen Event
0402         if (matchToGenEvt_) {
0403           edm::Handle<TtGenEvent> genEvt;
0404           iEvent.getByToken(evtSourceToken_, genEvt);
0405           asol.setGenEvt(genEvt);
0406         }
0407         // If asked, use the kin fitter to compute the top mass
0408         if (calcTopMass_) {
0409           solver->SetConstraints(xconstraint, yconstraint);
0410           solver->useWeightFromMC(useMCforBest_);
0411           asol = solver->addKinSolInfo(&asol);
0412         }
0413 
0414         // these lines calculate the observables to be used in the TtDilepSignalSelection LR
0415         (*myLRSignalSelObservables)(asol, iEvent);
0416 
0417         evtsols->push_back(asol);
0418       }
0419     }
0420     // flag the best solution (MC matching)
0421     if (matchToGenEvt_) {
0422       double bestSolDR = 9999.;
0423       int bestSol = -1;
0424       double dR = 0.;
0425       for (size_t s = 0; s < evtsols->size(); s++) {
0426         dR = (*evtsols)[s].getJetResidual();
0427         if (dR < bestSolDR) {
0428           bestSolDR = dR;
0429           bestSol = s;
0430         }
0431       }
0432       if (bestSol != -1)
0433         (*evtsols)[bestSol].setBestSol(true);
0434     }
0435     // put the result in the event
0436     std::unique_ptr<std::vector<TtDilepEvtSolution> > pOut(evtsols);
0437     iEvent.put(std::move(pOut));
0438   } else {
0439     // no solution: put a dummy solution in the event
0440     TtDilepEvtSolution asol;
0441     evtsols->push_back(asol);
0442     std::unique_ptr<std::vector<TtDilepEvtSolution> > pOut(evtsols);
0443     iEvent.put(std::move(pOut));
0444   }
0445 }
0446 
0447 #include "FWCore/Framework/interface/MakerMacros.h"
0448 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!