Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include <memory>
0002 #include <string>
0003 #include <vector>
0004 #include "TLorentzVector.h"
0005 #include "DataFormats/Math/interface/deltaR.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Framework/interface/stream/EDProducer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0012 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullLepKinSolver.h"
0013 
0014 class TtFullLepKinSolutionProducer : public edm::stream::EDProducer<> {
0015 public:
0016   explicit TtFullLepKinSolutionProducer(const edm::ParameterSet& iConfig);
0017   ~TtFullLepKinSolutionProducer() override = default;
0018 
0019   void produce(edm::Event& evt, const edm::EventSetup& iSetup) override;
0020 
0021 private:
0022   // next methods are avoidable but they make the code legible
0023   inline bool PTComp(const reco::Candidate*, const reco::Candidate*) const;
0024   inline bool LepDiffCharge(const reco::Candidate*, const reco::Candidate*) const;
0025   inline bool HasPositiveCharge(const reco::Candidate*) const;
0026 
0027   struct Compare {
0028     bool operator()(std::pair<double, int> a, std::pair<double, int> b) { return a.first > b.first; }
0029   };
0030 
0031   edm::EDGetTokenT<std::vector<pat::Jet> > jetsToken_;
0032   edm::EDGetTokenT<std::vector<pat::Electron> > electronsToken_;
0033   edm::EDGetTokenT<std::vector<pat::Muon> > muonsToken_;
0034   edm::EDGetTokenT<std::vector<pat::MET> > metsToken_;
0035 
0036   std::string jetCorrLevel_;
0037   int maxNJets_, maxNComb_;
0038   bool eeChannel_, emuChannel_, mumuChannel_, searchWrongCharge_;
0039   double tmassbegin_, tmassend_, tmassstep_;
0040   std::vector<double> nupars_;
0041 
0042   std::unique_ptr<TtFullLepKinSolver> solver;
0043 };
0044 
0045 inline bool TtFullLepKinSolutionProducer::PTComp(const reco::Candidate* l1, const reco::Candidate* l2) const {
0046   return (l1->pt() > l2->pt());
0047 }
0048 
0049 inline bool TtFullLepKinSolutionProducer::LepDiffCharge(const reco::Candidate* l1, const reco::Candidate* l2) const {
0050   return (l1->charge() != l2->charge());
0051 }
0052 
0053 inline bool TtFullLepKinSolutionProducer::HasPositiveCharge(const reco::Candidate* l) const {
0054   return (l->charge() > 0);
0055 }
0056 
0057 inline TtFullLepKinSolutionProducer::TtFullLepKinSolutionProducer(const edm::ParameterSet& iConfig) {
0058   // configurables
0059   jetsToken_ = consumes<std::vector<pat::Jet> >(iConfig.getParameter<edm::InputTag>("jets"));
0060   electronsToken_ = consumes<std::vector<pat::Electron> >(iConfig.getParameter<edm::InputTag>("electrons"));
0061   muonsToken_ = consumes<std::vector<pat::Muon> >(iConfig.getParameter<edm::InputTag>("muons"));
0062   metsToken_ = consumes<std::vector<pat::MET> >(iConfig.getParameter<edm::InputTag>("mets"));
0063   jetCorrLevel_ = iConfig.getParameter<std::string>("jetCorrectionLevel");
0064   maxNJets_ = iConfig.getParameter<int>("maxNJets");
0065   maxNComb_ = iConfig.getParameter<int>("maxNComb");
0066   eeChannel_ = iConfig.getParameter<bool>("eeChannel");
0067   emuChannel_ = iConfig.getParameter<bool>("emuChannel");
0068   mumuChannel_ = iConfig.getParameter<bool>("mumuChannel");
0069   searchWrongCharge_ = iConfig.getParameter<bool>("searchWrongCharge");
0070   tmassbegin_ = iConfig.getParameter<double>("tmassbegin");
0071   tmassend_ = iConfig.getParameter<double>("tmassend");
0072   tmassstep_ = iConfig.getParameter<double>("tmassstep");
0073   nupars_ = iConfig.getParameter<std::vector<double> >("neutrino_parameters");
0074 
0075   // define what will be produced
0076   produces<std::vector<std::vector<int> > >();  // vector of the particle inices (b, bbar, e1, e2, mu1, mu2)
0077   produces<std::vector<reco::LeafCandidate> >("fullLepNeutrinos");
0078   produces<std::vector<reco::LeafCandidate> >("fullLepNeutrinoBars");
0079   produces<std::vector<double> >("solWeight");  //weight for a specific kin solution
0080   produces<bool>("isWrongCharge");              //true if leptons have the same charge
0081 
0082   solver = std::make_unique<TtFullLepKinSolver>(tmassbegin_, tmassend_, tmassstep_, nupars_);
0083 }
0084 
0085 inline void TtFullLepKinSolutionProducer::produce(edm::Event& evt, const edm::EventSetup& iSetup) {
0086   //create vectors fo runsorted output
0087   std::vector<std::vector<int> > idcsV;
0088   std::vector<reco::LeafCandidate> nusV;
0089   std::vector<reco::LeafCandidate> nuBarsV;
0090   std::vector<std::pair<double, int> > weightsV;
0091 
0092   //create pointer for products
0093   std::unique_ptr<std::vector<std::vector<int> > > pIdcs(new std::vector<std::vector<int> >);
0094   std::unique_ptr<std::vector<reco::LeafCandidate> > pNus(new std::vector<reco::LeafCandidate>);
0095   std::unique_ptr<std::vector<reco::LeafCandidate> > pNuBars(new std::vector<reco::LeafCandidate>);
0096   std::unique_ptr<std::vector<double> > pWeight(new std::vector<double>);
0097   std::unique_ptr<bool> pWrongCharge(new bool);
0098 
0099   const edm::Handle<std::vector<pat::Jet> >& jets = evt.getHandle(jetsToken_);
0100   const edm::Handle<std::vector<pat::Electron> >& electrons = evt.getHandle(electronsToken_);
0101   const edm::Handle<std::vector<pat::Muon> >& muons = evt.getHandle(muonsToken_);
0102   const edm::Handle<std::vector<pat::MET> >& mets = evt.getHandle(metsToken_);
0103 
0104   int selMuon1 = -1, selMuon2 = -1;
0105   int selElectron1 = -1, selElectron2 = -1;
0106   bool ee = false;
0107   bool emu = false;
0108   bool mumu = false;
0109   bool isWrongCharge = false;
0110   bool jetsFound = false;
0111   bool METFound = false;
0112   bool electronsFound = false;
0113   bool electronMuonFound = false;
0114   bool muonsFound = false;
0115 
0116   //select Jets (TopJet vector is sorted on ET)
0117   if (jets->size() >= 2) {
0118     jetsFound = true;
0119   }
0120 
0121   //select MET (TopMET vector is sorted on ET)
0122   if (!mets->empty()) {
0123     METFound = true;
0124   }
0125 
0126   // If we have electrons and muons available,
0127   // build a solutions with electrons and muons.
0128   if (muons->size() + electrons->size() >= 2) {
0129     // select leptons
0130     if (electrons->empty())
0131       mumu = true;
0132     else if (muons->empty())
0133       ee = true;
0134     else if (electrons->size() == 1) {
0135       if (muons->size() == 1)
0136         emu = true;
0137       else if (PTComp(&(*electrons)[0], &(*muons)[1]))
0138         emu = true;
0139       else
0140         mumu = true;
0141     } else if (electrons->size() > 1) {
0142       if (PTComp(&(*electrons)[1], &(*muons)[0]))
0143         ee = true;
0144       else if (muons->size() == 1)
0145         emu = true;
0146       else if (PTComp(&(*electrons)[0], &(*muons)[1]))
0147         emu = true;
0148       else
0149         mumu = true;
0150     }
0151     if (ee && eeChannel_) {
0152       if (LepDiffCharge(&(*electrons)[0], &(*electrons)[1]) || searchWrongCharge_) {
0153         if (HasPositiveCharge(&(*electrons)[0]) || !LepDiffCharge(&(*electrons)[0], &(*electrons)[1])) {
0154           selElectron1 = 0;
0155           selElectron2 = 1;
0156         } else {
0157           selElectron1 = 1;
0158           selElectron2 = 0;
0159         }
0160         electronsFound = true;
0161         if (!LepDiffCharge(&(*electrons)[0], &(*electrons)[1]))
0162           isWrongCharge = true;
0163       }
0164     } else if (emu && emuChannel_) {
0165       if (LepDiffCharge(&(*electrons)[0], &(*muons)[0]) || searchWrongCharge_) {
0166         selElectron1 = 0;
0167         selMuon1 = 0;
0168         electronMuonFound = true;
0169         if (!LepDiffCharge(&(*electrons)[0], &(*muons)[0]))
0170           isWrongCharge = true;
0171       }
0172     } else if (mumu && mumuChannel_) {
0173       if (LepDiffCharge(&(*muons)[0], &(*muons)[1]) || searchWrongCharge_) {
0174         if (HasPositiveCharge(&(*muons)[0]) || !LepDiffCharge(&(*muons)[0], &(*muons)[1])) {
0175           selMuon1 = 0;
0176           selMuon2 = 1;
0177         } else {
0178           selMuon1 = 1;
0179           selMuon2 = 0;
0180         }
0181         muonsFound = true;
0182         if (!LepDiffCharge(&(*muons)[0], &(*muons)[1]))
0183           isWrongCharge = true;
0184       }
0185     }
0186   }
0187 
0188   *pWrongCharge = isWrongCharge;
0189 
0190   // Check that the above work makes sense
0191   if (int(ee) + int(emu) + int(mumu) > 1) {
0192     edm::LogWarning("TtFullLepKinSolutionProducer") << "Lepton selection criteria uncorrectly defined";
0193   }
0194 
0195   // Check if the leptons for the required Channel are available
0196   bool correctLeptons =
0197       ((electronsFound && eeChannel_) || (muonsFound && mumuChannel_) || (electronMuonFound && emuChannel_));
0198   // Check for equally charged leptons if for wrong charge combinations is searched
0199   if (isWrongCharge) {
0200     correctLeptons = correctLeptons && searchWrongCharge_;
0201   }
0202 
0203   if (correctLeptons && METFound && jetsFound) {
0204     // run over all jets if input parameter maxNJets is -1 or
0205     // adapt maxNJets if number of present jets is smaller than selected
0206     // number of jets
0207     int stop = maxNJets_;
0208     if (jets->size() < static_cast<unsigned int>(stop) || stop < 0)
0209       stop = jets->size();
0210 
0211     // counter for number of found kinematic solutions
0212     int nSol = 0;
0213 
0214     // consider all permutations
0215     for (int ib = 0; ib < stop; ib++) {
0216       // second loop of the permutations
0217       for (int ibbar = 0; ibbar < stop; ibbar++) {
0218         // avoid the diagonal: b and bbar must be distinct jets
0219         if (ib == ibbar)
0220           continue;
0221 
0222         std::vector<int> idcs;
0223 
0224         // push back the indices of the jets
0225         idcs.push_back(ib);
0226         idcs.push_back(ibbar);
0227 
0228         TLorentzVector LV_l1;
0229         TLorentzVector LV_l2;
0230         TLorentzVector LV_b = TLorentzVector((*jets)[ib].correctedJet(jetCorrLevel_, "bottom").px(),
0231                                              (*jets)[ib].correctedJet(jetCorrLevel_, "bottom").py(),
0232                                              (*jets)[ib].correctedJet(jetCorrLevel_, "bottom").pz(),
0233                                              (*jets)[ib].correctedJet(jetCorrLevel_, "bottom").energy());
0234         TLorentzVector LV_bbar = TLorentzVector((*jets)[ibbar].correctedJet(jetCorrLevel_, "bottom").px(),
0235                                                 (*jets)[ibbar].correctedJet(jetCorrLevel_, "bottom").py(),
0236                                                 (*jets)[ibbar].correctedJet(jetCorrLevel_, "bottom").pz(),
0237                                                 (*jets)[ibbar].correctedJet(jetCorrLevel_, "bottom").energy());
0238 
0239         double xconstraint = 0, yconstraint = 0;
0240 
0241         if (ee) {
0242           idcs.push_back(selElectron1);
0243           LV_l1.SetXYZT((*electrons)[selElectron1].px(),
0244                         (*electrons)[selElectron1].py(),
0245                         (*electrons)[selElectron1].pz(),
0246                         (*electrons)[selElectron1].energy());
0247           xconstraint += (*electrons)[selElectron1].px();
0248           yconstraint += (*electrons)[selElectron1].py();
0249 
0250           idcs.push_back(selElectron2);
0251           LV_l2.SetXYZT((*electrons)[selElectron2].px(),
0252                         (*electrons)[selElectron2].py(),
0253                         (*electrons)[selElectron2].pz(),
0254                         (*electrons)[selElectron2].energy());
0255           xconstraint += (*electrons)[selElectron2].px();
0256           yconstraint += (*electrons)[selElectron2].py();
0257 
0258           idcs.push_back(-1);
0259           idcs.push_back(-1);
0260         }
0261 
0262         else if (emu) {
0263           if (!isWrongCharge) {
0264             if (HasPositiveCharge(&(*electrons)[selElectron1])) {
0265               idcs.push_back(selElectron1);
0266               LV_l1.SetXYZT((*electrons)[selElectron1].px(),
0267                             (*electrons)[selElectron1].py(),
0268                             (*electrons)[selElectron1].pz(),
0269                             (*electrons)[selElectron1].energy());
0270               xconstraint += (*electrons)[selElectron1].px();
0271               yconstraint += (*electrons)[selElectron1].py();
0272 
0273               idcs.push_back(-1);
0274               idcs.push_back(-1);
0275 
0276               idcs.push_back(selMuon1);
0277               LV_l2.SetXYZT((*muons)[selMuon1].px(),
0278                             (*muons)[selMuon1].py(),
0279                             (*muons)[selMuon1].pz(),
0280                             (*muons)[selMuon1].energy());
0281               xconstraint += (*muons)[selMuon1].px();
0282               yconstraint += (*muons)[selMuon1].py();
0283             } else {
0284               idcs.push_back(-1);
0285 
0286               idcs.push_back(selMuon1);
0287               LV_l1.SetXYZT((*muons)[selMuon1].px(),
0288                             (*muons)[selMuon1].py(),
0289                             (*muons)[selMuon1].pz(),
0290                             (*muons)[selMuon1].energy());
0291               xconstraint += (*muons)[selMuon1].px();
0292               yconstraint += (*muons)[selMuon1].py();
0293 
0294               idcs.push_back(selElectron1);
0295               LV_l2.SetXYZT((*electrons)[selElectron1].px(),
0296                             (*electrons)[selElectron1].py(),
0297                             (*electrons)[selElectron1].pz(),
0298                             (*electrons)[selElectron1].energy());
0299               xconstraint += (*electrons)[selElectron1].px();
0300               yconstraint += (*electrons)[selElectron1].py();
0301 
0302               idcs.push_back(-1);
0303             }
0304           } else {                                                 // means "if wrong charge"
0305             if (HasPositiveCharge(&(*electrons)[selElectron1])) {  // both leps positive
0306               idcs.push_back(selElectron1);
0307               LV_l1.SetXYZT((*electrons)[selElectron1].px(),
0308                             (*electrons)[selElectron1].py(),
0309                             (*electrons)[selElectron1].pz(),
0310                             (*electrons)[selElectron1].energy());
0311               xconstraint += (*electrons)[selElectron1].px();
0312               yconstraint += (*electrons)[selElectron1].py();
0313 
0314               idcs.push_back(-1);
0315 
0316               idcs.push_back(selMuon1);
0317               LV_l2.SetXYZT((*muons)[selMuon1].px(),
0318                             (*muons)[selMuon1].py(),
0319                             (*muons)[selMuon1].pz(),
0320                             (*muons)[selMuon1].energy());
0321               xconstraint += (*muons)[selMuon1].px();
0322               yconstraint += (*muons)[selMuon1].py();
0323 
0324               idcs.push_back(-1);
0325             } else {  // both leps negative
0326               idcs.push_back(-1);
0327 
0328               idcs.push_back(selElectron1);
0329               LV_l2.SetXYZT((*electrons)[selElectron1].px(),
0330                             (*electrons)[selElectron1].py(),
0331                             (*electrons)[selElectron1].pz(),
0332                             (*electrons)[selElectron1].energy());
0333               xconstraint += (*electrons)[selElectron1].px();
0334               yconstraint += (*electrons)[selElectron1].py();
0335 
0336               idcs.push_back(-1);
0337 
0338               idcs.push_back(selMuon1);
0339               LV_l1.SetXYZT((*muons)[selMuon1].px(),
0340                             (*muons)[selMuon1].py(),
0341                             (*muons)[selMuon1].pz(),
0342                             (*muons)[selMuon1].energy());
0343               xconstraint += (*muons)[selMuon1].px();
0344               yconstraint += (*muons)[selMuon1].py();
0345             }
0346           }
0347         }
0348 
0349         else if (mumu) {
0350           idcs.push_back(-1);
0351           idcs.push_back(-1);
0352 
0353           idcs.push_back(selMuon1);
0354           LV_l1.SetXYZT(
0355               (*muons)[selMuon1].px(), (*muons)[selMuon1].py(), (*muons)[selMuon1].pz(), (*muons)[selMuon1].energy());
0356           xconstraint += (*muons)[selMuon1].px();
0357           yconstraint += (*muons)[selMuon1].py();
0358 
0359           idcs.push_back(selMuon2);
0360           LV_l2.SetXYZT(
0361               (*muons)[selMuon2].px(), (*muons)[selMuon2].py(), (*muons)[selMuon2].pz(), (*muons)[selMuon2].energy());
0362           xconstraint += (*muons)[selMuon2].px();
0363           yconstraint += (*muons)[selMuon2].py();
0364         }
0365 
0366         xconstraint += (*jets)[ib].px() + (*jets)[ibbar].px() + (*mets)[0].px();
0367         yconstraint += (*jets)[ib].py() + (*jets)[ibbar].py() + (*mets)[0].py();
0368 
0369         // calculate neutrino momenta and eventweight
0370         solver->SetConstraints(xconstraint, yconstraint);
0371         TtFullLepKinSolver::NeutrinoSolution nuSol = solver->getNuSolution(LV_l1, LV_l2, LV_b, LV_bbar);
0372 
0373         // add solution to the vectors of solutions if solution exists
0374         if (nuSol.weight > 0) {
0375           // add the leptons and jets indices to the vector of combinations
0376           idcsV.push_back(idcs);
0377 
0378           // add the neutrinos
0379           nusV.push_back(nuSol.neutrino);
0380           nuBarsV.push_back(nuSol.neutrinoBar);
0381 
0382           // add the solution weight
0383           weightsV.push_back(std::make_pair(nuSol.weight, nSol));
0384 
0385           nSol++;
0386         }
0387       }
0388     }
0389   }
0390 
0391   if (weightsV.empty()) {
0392     //create dmummy vector
0393     std::vector<int> idcs;
0394     idcs.reserve(6);
0395     for (int i = 0; i < 6; ++i)
0396       idcs.push_back(-1);
0397 
0398     idcsV.push_back(idcs);
0399     weightsV.push_back(std::make_pair(-1, 0));
0400     reco::LeafCandidate nu;
0401     nusV.push_back(nu);
0402     reco::LeafCandidate nuBar;
0403     nuBarsV.push_back(nuBar);
0404   }
0405 
0406   // check if all vectors have correct length
0407   int weightL = weightsV.size();
0408   int nuL = nusV.size();
0409   int nuBarL = nuBarsV.size();
0410   int idxL = idcsV.size();
0411 
0412   if (weightL != nuL || weightL != nuBarL || weightL != idxL) {
0413     edm::LogWarning("TtFullLepKinSolutionProducer")
0414         << "Output vectors are of different length:"
0415         << "\n weight: " << weightL << "\n     nu: " << nuL << "\n  nubar: " << nuBarL << "\n   idcs: " << idxL;
0416   }
0417 
0418   // sort vectors by weight in decreasing order
0419   if (weightsV.size() > 1) {
0420     sort(weightsV.begin(), weightsV.end(), Compare());
0421   }
0422 
0423   // determine the number of solutions which is written in the event
0424   int stop = weightL;
0425   if (maxNComb_ > 0 && maxNComb_ < stop)
0426     stop = maxNComb_;
0427 
0428   for (int i = 0; i < stop; ++i) {
0429     pWeight->push_back(weightsV[i].first);
0430     pNus->push_back(nusV[weightsV[i].second]);
0431     pNuBars->push_back(nuBarsV[weightsV[i].second]);
0432     pIdcs->push_back(idcsV[weightsV[i].second]);
0433   }
0434 
0435   // put the results in the event
0436   evt.put(std::move(pIdcs));
0437   evt.put(std::move(pNus), "fullLepNeutrinos");
0438   evt.put(std::move(pNuBars), "fullLepNeutrinoBars");
0439   evt.put(std::move(pWeight), "solWeight");
0440   evt.put(std::move(pWrongCharge), "isWrongCharge");
0441 }
0442 
0443 #include "FWCore/Framework/interface/MakerMacros.h"
0444 DEFINE_FWK_MODULE(TtFullLepKinSolutionProducer);
[ 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!