TopEventProducers/src/TtHadEvtSolutionMaker.cc

0001 //
0002 // adapted TtSemiEvtSolutionMaker.h, v1.13 2007/07/06 02:49:42 lowette Exp $
0003 // for fully hadronic channel.
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/TopTools/interface/JetPartonMatching.h"
0012 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullHadKinFitter.h"
0013 #include "TopQuarkAnalysis/TopJetCombination/interface/TtHadSimpleBestJetComb.h"
0014 #include "TopQuarkAnalysis/TopJetCombination/interface/TtHadLRJetCombObservables.h"
0015 #include "TopQuarkAnalysis/TopJetCombination/interface/TtHadLRJetCombCalc.h"
0016 #include "TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelObservables.h"
0017 #include "TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelCalc.h"
0018
0019 #include "DataFormats/Math/interface/deltaR.h"
0020 #include "AnalysisDataFormats/TopObjects/interface/TtHadEvtSolution.h"
0021
0022 #include <memory>
0023 #include <string>
0024 #include <vector>
0025
0026 class TtHadEvtSolutionMaker : public edm::stream::EDProducer<> {
0027 public:
0028   explicit TtHadEvtSolutionMaker(const edm::ParameterSet& iConfig);
0029   ~TtHadEvtSolutionMaker() override;
0030
0031   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0032
0033 private:
0034   // configurables
0035
0036   edm::EDGetTokenT<std::vector<pat::Jet> > jetSrcToken_;
0037   int jetCorrScheme_;
0038   std::string lrSignalSelFile_, lrJetCombFile_;
0039   bool addLRSignalSel_, addLRJetComb_, doKinFit_, matchToGenEvt_;
0040   int matchingAlgo_;
0041   bool useMaxDist_, useDeltaR_;
0042   double maxDist_;
0043   int maxNrIter_;
0044   double maxDeltaS_, maxF_;
0045   int jetParam_;
0046   std::vector<int> lrSignalSelObs_, lrJetCombObs_;
0047   std::vector<unsigned int> constraints_;
0048   edm::EDGetTokenT<TtGenEvent> genEvtToken_;
0049   // tools
0050   TtFullHadKinFitter* myKinFitter;
0051   TtHadSimpleBestJetComb* mySimpleBestJetComb;
0052   TtHadLRJetCombObservables* myLRJetCombObservables;
0053   TtHadLRJetCombCalc* myLRJetCombCalc;
0054   TtHadLRSignalSelObservables* myLRSignalSelObservables;
0055   TtHadLRSignalSelCalc* myLRSignalSelCalc;
0056 };
0057
0058 /// constructor
0059 TtHadEvtSolutionMaker::TtHadEvtSolutionMaker(const edm::ParameterSet& iConfig) {
0060   // configurables
0061   jetSrcToken_ = consumes<std::vector<pat::Jet> >(iConfig.getParameter<edm::InputTag>("jetSource"));
0062   jetCorrScheme_ = iConfig.getParameter<int>("jetCorrectionScheme");
0063   doKinFit_ = iConfig.getParameter<bool>("doKinFit");
0064   addLRSignalSel_ = iConfig.getParameter<bool>("addLRSignalSel");
0065   lrSignalSelObs_ = iConfig.getParameter<std::vector<int> >("lrSignalSelObs");
0066   lrSignalSelFile_ = iConfig.getParameter<std::string>("lrSignalSelFile");
0067   addLRJetComb_ = iConfig.getParameter<bool>("addLRJetComb");
0068   lrJetCombObs_ = iConfig.getParameter<std::vector<int> >("lrJetCombObs");
0069   lrJetCombFile_ = iConfig.getParameter<std::string>("lrJetCombFile");
0070   maxNrIter_ = iConfig.getParameter<int>("maxNrIter");
0071   maxDeltaS_ = iConfig.getParameter<double>("maxDeltaS");
0072   maxF_ = iConfig.getParameter<double>("maxF");
0073   jetParam_ = iConfig.getParameter<int>("jetParametrisation");
0074   constraints_ = iConfig.getParameter<std::vector<unsigned int> >("constraints");
0075   matchToGenEvt_ = iConfig.getParameter<bool>("matchToGenEvt");
0076   matchingAlgo_ = iConfig.getParameter<bool>("matchingAlgorithm");
0077   useMaxDist_ = iConfig.getParameter<bool>("useMaximalDistance");
0078   useDeltaR_ = iConfig.getParameter<bool>("useDeltaR");
0079   maxDist_ = iConfig.getParameter<double>("maximalDistance");
0080   genEvtToken_ = mayConsume<TtGenEvent>(edm::InputTag("genEvt"));
0081
0082   // define kinfitter
0083   if (doKinFit_) {
0084     myKinFitter = new TtFullHadKinFitter(jetParam_, maxNrIter_, maxDeltaS_, maxF_, constraints_);
0085   }
0086
0087   // define jet combinations related calculators
0088   mySimpleBestJetComb = new TtHadSimpleBestJetComb();
0089   myLRSignalSelObservables = new TtHadLRSignalSelObservables();
0090   myLRJetCombObservables = new TtHadLRJetCombObservables();
0091   if (addLRJetComb_)
0092     myLRJetCombCalc = new TtHadLRJetCombCalc(lrJetCombFile_, lrJetCombObs_);
0093
0094   // instantiate signal selection calculator
0095   if (addLRSignalSel_)
0096     myLRSignalSelCalc = new TtHadLRSignalSelCalc(lrSignalSelFile_, lrSignalSelObs_);
0097
0098   // define what will be produced
0099   produces<std::vector<TtHadEvtSolution> >();
0100 }
0101
0102 /// destructor
0103 TtHadEvtSolutionMaker::~TtHadEvtSolutionMaker() {
0104   if (doKinFit_) {
0105     delete myKinFitter;
0106   }
0107   delete mySimpleBestJetComb;
0108   delete myLRSignalSelObservables;
0109   delete myLRJetCombObservables;
0110   if (addLRSignalSel_)
0111     delete myLRSignalSelCalc;
0112   if (addLRJetComb_)
0113     delete myLRJetCombCalc;
0114 }
0115
0116 void TtHadEvtSolutionMaker::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0117   // TopObject Selection
0118   // Select Jets
0119
0120   bool jetsFound = false;
0121   edm::Handle<std::vector<pat::Jet> > jets;
0122   iEvent.getByToken(jetSrcToken_, jets);
0123
0124   if (jets->size() >= 6)
0125     jetsFound = true;
0126
0127   // Build Event solutions according to the ambiguity in the jet combination
0128   // Note, hardcoded to only run through the 6 most energetic jets - could be changed ....
0129
0130   std::vector<TtHadEvtSolution>* evtsols = new std::vector<TtHadEvtSolution>();
0131   if (jetsFound) {
0132     for (unsigned int p = 0; p < 3; p++) {                         // loop over light jet p
0133       for (unsigned int q = p + 1; q < 4; q++) {                   // loop over light jet q
0134         for (unsigned int j = q + 1; j < 5; j++) {                 // loop over light jet j
0135           for (unsigned int k = j + 1; k < 6; k++) {               // loop over light jet k
0136             for (unsigned int bh = 0; bh != jets->size(); bh++) {  //loop over hadronic b-jet1
0137               if (!(bh == p || bh == q || bh == j || bh == k)) {
0138                 for (unsigned int bbarh = 0; bbarh != jets->size(); bbarh++) {  //loop over hadronic b-jet2
0139                   if (!(bbarh == p || bbarh == q || bbarh == j || bbarh == k) && !(bbarh == bh)) {
0140                     // Make event solutions for all possible combinations of the 4 light
0141                     // jets and 2 possible b-jets, not including the option of the b's being swapped.
0142                     // Hadp,Hadq is one pair, Hadj,Hadk the other
0143                     std::vector<TtHadEvtSolution> asol;
0144                     asol.resize(3);
0145                     //[p][q][b] and [j][k][bbar]
0146                     asol[0].setJetCorrectionScheme(jetCorrScheme_);
0147                     asol[0].setHadp(jets, p);
0148                     asol[0].setHadq(jets, q);
0149                     asol[0].setHadj(jets, j);
0150                     asol[0].setHadk(jets, k);
0151                     asol[0].setHadb(jets, bh);
0152                     asol[0].setHadbbar(jets, bbarh);
0153
0154                     //[p][j][b] and [q][k][bbar]
0155                     asol[1].setJetCorrectionScheme(jetCorrScheme_);
0156                     asol[1].setHadp(jets, p);
0157                     asol[1].setHadq(jets, j);
0158                     asol[1].setHadj(jets, q);
0159                     asol[1].setHadk(jets, k);
0160                     asol[1].setHadb(jets, bh);
0161                     asol[1].setHadbbar(jets, bbarh);
0162
0163                     //[p][k][b] and [j][q][bbar]
0164                     asol[2].setJetCorrectionScheme(jetCorrScheme_);
0165                     asol[2].setHadp(jets, p);
0166                     asol[2].setHadq(jets, k);
0167                     asol[2].setHadj(jets, j);
0168                     asol[2].setHadk(jets, q);
0169                     asol[2].setHadb(jets, bh);
0170                     asol[2].setHadbbar(jets, bbarh);
0171
0172                     if (doKinFit_) {
0173                       for (unsigned int i = 0; i != asol.size(); i++) {
0174                         asol[i] = myKinFitter->addKinFitInfo(&(asol[i]));
0175                         asol[i].setJetParametrisation(jetParam_);
0176                       }
0177
0178                     } else {
0179                       std::cout << "Fitting needed to decide on best solution, enable fitting!" << std::endl;
0180                     }
0181                     // these lines calculate the observables to be used in the TtHadSignalSelection LR
0182
0183                     for (unsigned int i = 0; i != asol.size(); i++) {
0184                       (*myLRSignalSelObservables)(asol[i]);
0185                     }
0186                     // if asked for, calculate with these observable values the LRvalue and
0187                     // (depending on the configuration) probability this event is signal
0188                     if (addLRSignalSel_) {
0189                       for (unsigned int i = 0; i != asol.size(); i++) {
0190                         (*myLRSignalSelCalc)(asol[i]);
0191                       }
0192                     }
0193
0194                     // these lines calculate the observables to be used in the TtHadJetCombination LR
0195                     for (unsigned int i = 0; i != asol.size(); i++) {
0196                       (*myLRJetCombObservables)(asol[i]);
0197                     }
0198                     // if asked for, calculate with these observable values the LRvalue and
0199                     // (depending on the configuration) probability a jet combination is correct
0200                     if (addLRJetComb_) {
0201                       for (unsigned int i = 0; i != asol.size(); i++) {
0202                         (*myLRJetCombCalc)(asol[i]);
0203                       }
0204                     }
0205                     //std::cout<<"SignalSelLRval = "<<asol.getLRSignalEvtLRval()<<"  JetCombProb = "<<asol.getLRSignalEvtProb()<<std::endl;
0206                     //std::cout<<"JetCombLRval = "<<asol.getLRJetCombLRval()<<"  JetCombProb = "<<asol.getLRJetCombProb()<<std::endl;
0207
0208                     // fill solution to vector with all possible solutions
0209                     for (unsigned int i = 0; i != asol.size(); i++) {
0210                       evtsols->push_back(asol[i]);
0211                     }
0212                   }
0213                 }
0214               }
0215             }
0216           }
0217         }
0218       }
0219     }
0220
0221     // add TtHadSimpleBestJetComb to solutions
0222     int simpleBestJetComb = (*mySimpleBestJetComb)(*evtsols);
0223
0224     for (size_t s = 0; s < evtsols->size(); s++) {
0225       (*evtsols)[s].setSimpleBestJetComb(simpleBestJetComb);
0226       // if asked for, match the event solutions to the gen Event
0227       if (matchToGenEvt_) {
0228         int bestSolution = -999;
0229         int bestSolutionChangeW1Q = -999;
0230         int bestSolutionChangeW2Q = -999;
0231         edm::Handle<TtGenEvent> genEvt;
0232         iEvent.getByToken(genEvtToken_, genEvt);
0233         std::vector<const reco::Candidate*> quarks;
0234         const reco::Candidate& genp = *(genEvt->daughterQuarkOfWPlus());
0235         const reco::Candidate& genq = *(genEvt->daughterQuarkBarOfWPlus());
0236         const reco::Candidate& genb = *(genEvt->b());
0237         const reco::Candidate& genj = *(genEvt->daughterQuarkOfWMinus());
0238         const reco::Candidate& genk = *(genEvt->daughterQuarkBarOfWMinus());
0239         const reco::Candidate& genbbar = *(genEvt->bBar());
0240         quarks.push_back(&genp);
0241         quarks.push_back(&genq);
0242         quarks.push_back(&genb);
0243         quarks.push_back(&genj);
0244         quarks.push_back(&genk);
0245         quarks.push_back(&genbbar);
0246         std::vector<const reco::Candidate*> jets;
0247         for (size_t s = 0; s < evtsols->size(); s++) {
0248           jets.clear();
0249           const reco::Candidate& jetp = (*evtsols)[s].getRecHadp();
0250           const reco::Candidate& jetq = (*evtsols)[s].getRecHadq();
0251           const reco::Candidate& jetbh = (*evtsols)[s].getRecHadb();
0252           const reco::Candidate& jetj = (*evtsols)[s].getRecHadj();
0253           const reco::Candidate& jetk = (*evtsols)[s].getRecHadk();
0254           const reco::Candidate& jetbbar = (*evtsols)[s].getRecHadbbar();
0255           jets.push_back(&jetp);
0256           jets.push_back(&jetq);
0257           jets.push_back(&jetbh);
0258           jets.push_back(&jetj);
0259           jets.push_back(&jetk);
0260           jets.push_back(&jetbbar);
0261           JetPartonMatching aMatch(quarks, jets, matchingAlgo_, useMaxDist_, useDeltaR_, maxDist_);
0262           (*evtsols)[s].setGenEvt(genEvt);
0263           (*evtsols)[s].setMCBestSumAngles(aMatch.getSumDistances());
0264           (*evtsols)[s].setMCBestAngleHadp(aMatch.getDistanceForParton(0));
0265           (*evtsols)[s].setMCBestAngleHadq(aMatch.getDistanceForParton(1));
0266           (*evtsols)[s].setMCBestAngleHadb(aMatch.getDistanceForParton(2));
0267           (*evtsols)[s].setMCBestAngleHadb(aMatch.getDistanceForParton(2));
0268           (*evtsols)[s].setMCBestAngleHadj(aMatch.getDistanceForParton(3));
0269           (*evtsols)[s].setMCBestAngleHadk(aMatch.getDistanceForParton(4));
0270           (*evtsols)[s].setMCBestAngleHadbbar(aMatch.getDistanceForParton(5));
0271
0272           // Check match - checking if two light quarks are swapped wrt matched gen particle
0273           if ((aMatch.getMatchForParton(2) == 2 && aMatch.getMatchForParton(5) == 5) ||
0274               (aMatch.getMatchForParton(2) == 5 && aMatch.getMatchForParton(5) == 2)) {  // check b-jets
0275
0276             if (aMatch.getMatchForParton(3) == 3 && aMatch.getMatchForParton(4) == 4) {  //check light jets
0277               bestSolutionChangeW2Q = 0;
0278               if (aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) {
0279                 bestSolution = s;
0280                 bestSolutionChangeW1Q = 0;
0281               } else {
0282                 if (aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0) {
0283                   bestSolution = s;
0284                   bestSolutionChangeW1Q = 1;
0285                 }
0286               }
0287             } else {
0288               if (aMatch.getMatchForParton(2) == 3 && aMatch.getMatchForParton(3) == 2) {  // or check if swapped
0289                 bestSolutionChangeW2Q = 1;
0290                 if (aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0) {
0291                   bestSolution = s;
0292                   bestSolutionChangeW1Q = 1;
0293                 } else {
0294                   if (aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) {
0295                     bestSolution = s;
0296                     bestSolutionChangeW1Q = 0;
0297                   }
0298                 }
0299               }
0300               if (aMatch.getMatchForParton(2) == 2 && aMatch.getMatchForParton(3) == 3) {
0301                 bestSolutionChangeW2Q = 0;
0302                 if (aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) {
0303                   bestSolution = s;
0304                   bestSolutionChangeW1Q = 0;
0305                 } else if (aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0) {
0306                   bestSolution = s;
0307                   bestSolutionChangeW1Q = 1;
0308                 }
0309               }
0310             }
0311           }
0312           for (size_t s = 0; s < evtsols->size(); s++) {
0313             (*evtsols)[s].setMCBestJetComb(bestSolution);
0314             (*evtsols)[s].setMCChangeW1Q(bestSolutionChangeW1Q);
0315             (*evtsols)[s].setMCChangeW2Q(bestSolutionChangeW2Q);
0316           }
0317         }
0318       }  // end matchEvt
0319     }
0320     //store the vector of solutions to the event
0321
0322     std::unique_ptr<std::vector<TtHadEvtSolution> > pOut(evtsols);
0323     iEvent.put(std::move(pOut));
0324   } else {  //end loop jet/MET found
0325     std::cout << "No calibrated solutions built, because only " << jets->size() << " were present";
0326
0327     std::unique_ptr<std::vector<TtHadEvtSolution> > pOut(evtsols);
0328     iEvent.put(std::move(pOut));
0329   }
0330 }
0331
0332 #include "FWCore/Framework/interface/MakerMacros.h"
0333 DEFINE_FWK_MODULE(TtHadEvtSolutionMaker);