Reference Manual

Go to the documentation of this file.
00001 // $Id: TtHadEvtSolutionMaker.cc,v 1.21 2010/03/30 14:04:44 snaumann Exp $
00002 
00003 #include "TopQuarkAnalysis/TopEventProducers/interface/TtHadEvtSolutionMaker.h"
00004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00005 #include "DataFormats/Math/interface/deltaR.h"
00006 #include "AnalysisDataFormats/TopObjects/interface/TtHadEvtSolution.h"
00007 #include "TopQuarkAnalysis/TopTools/interface/JetPartonMatching.h"
00008 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullHadKinFitter.h"
00009 #include "TopQuarkAnalysis/TopJetCombination/interface/TtHadSimpleBestJetComb.h"
00010 #include "TopQuarkAnalysis/TopJetCombination/interface/TtHadLRJetCombObservables.h"
00011 #include "TopQuarkAnalysis/TopJetCombination/interface/TtHadLRJetCombCalc.h"
00012 #include "TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelObservables.h"
00013 #include "TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelCalc.h"
00014 
00015 #include <memory>
00016 
00017 
00019 TtHadEvtSolutionMaker::TtHadEvtSolutionMaker(const edm::ParameterSet & iConfig) {
00020   // configurables
00021   jetSrc_          = iConfig.getParameter<edm::InputTag>    ("jetSource");
00022   jetCorrScheme_   = iConfig.getParameter<int>              ("jetCorrectionScheme");
00023   doKinFit_        = iConfig.getParameter<bool>             ("doKinFit");
00024   addLRSignalSel_  = iConfig.getParameter<bool>             ("addLRSignalSel");
00025   lrSignalSelObs_  = iConfig.getParameter<std::vector<int> >("lrSignalSelObs");
00026   lrSignalSelFile_ = iConfig.getParameter<std::string>      ("lrSignalSelFile");
00027   addLRJetComb_    = iConfig.getParameter<bool>             ("addLRJetComb");
00028   lrJetCombObs_    = iConfig.getParameter<std::vector<int> >("lrJetCombObs");
00029   lrJetCombFile_   = iConfig.getParameter<std::string>      ("lrJetCombFile");
00030   maxNrIter_       = iConfig.getParameter<int>              ("maxNrIter");
00031   maxDeltaS_       = iConfig.getParameter<double>           ("maxDeltaS");
00032   maxF_            = iConfig.getParameter<double>           ("maxF");
00033   jetParam_        = iConfig.getParameter<int>              ("jetParametrisation");
00034   constraints_     = iConfig.getParameter<std::vector<unsigned int> >("constraints");
00035   matchToGenEvt_   = iConfig.getParameter<bool>             ("matchToGenEvt");
00036   matchingAlgo_    = iConfig.getParameter<bool>             ("matchingAlgorithm");
00037   useMaxDist_      = iConfig.getParameter<bool>             ("useMaximalDistance");
00038   useDeltaR_       = iConfig.getParameter<bool>             ("useDeltaR");
00039   maxDist_         = iConfig.getParameter<double>           ("maximalDistance");
00040 
00041   // define kinfitter
00042   if(doKinFit_){
00043     myKinFitter = new TtFullHadKinFitter(jetParam_, maxNrIter_, maxDeltaS_, maxF_, constraints_);
00044   }
00045   
00046   
00047   // define jet combinations related calculators
00048   mySimpleBestJetComb                    = new TtHadSimpleBestJetComb();
00049   myLRSignalSelObservables               = new TtHadLRSignalSelObservables();
00050   myLRJetCombObservables                 = new TtHadLRJetCombObservables();
00051   if (addLRJetComb_)   myLRJetCombCalc   = new TtHadLRJetCombCalc(lrJetCombFile_, lrJetCombObs_);
00052   
00053   // instantiate signal selection calculator
00054   if (addLRSignalSel_) myLRSignalSelCalc = new TtHadLRSignalSelCalc(lrSignalSelFile_, lrSignalSelObs_);
00055  
00056   // define what will be produced
00057   produces<std::vector<TtHadEvtSolution> >();
00058   
00059 }
00060 
00061 
00063 TtHadEvtSolutionMaker::~TtHadEvtSolutionMaker()
00064 {
00065   if (doKinFit_) {
00066     delete myKinFitter;
00067   }
00068   delete mySimpleBestJetComb;
00069   delete myLRSignalSelObservables;
00070   delete myLRJetCombObservables;
00071   if(addLRSignalSel_) delete myLRSignalSelCalc;
00072   if(addLRJetComb_)   delete myLRJetCombCalc;
00073 }
00074 
00075 
00076 void TtHadEvtSolutionMaker::produce(edm::Event & iEvent, const edm::EventSetup & iSetup) {
00077   // TopObject Selection
00078   // Select Jets
00079 
00080   bool jetsFound = false;
00081   edm::Handle<std::vector<pat::Jet> > jets;
00082   iEvent.getByLabel(jetSrc_, jets);
00083 
00084   if (jets->size() >= 6) jetsFound = true;
00085   
00086   // Build Event solutions according to the ambiguity in the jet combination
00087   // Note, hardcoded to only run through the 6 most energetic jets - could be changed ....
00088   
00089   std::vector<TtHadEvtSolution> * evtsols = new std::vector<TtHadEvtSolution>();
00090   if(jetsFound){
00091     for (unsigned int p=0; p<3; p++) {  // loop over light jet p
00092       for (unsigned int q=p+1; q<4; q++) { // loop over light jet q
00093         for (unsigned int j=q+1; j<5; j++) {  // loop over light jet j
00094           for (unsigned int k=j+1; k<6; k++) { // loop over light jet k
00095             for (unsigned int bh=0; bh!=jets->size(); bh++) { //loop over hadronic b-jet1
00096               if(!(bh==p || bh==q || bh==j || bh==k)) { 
00097                 for (unsigned int bbarh=0; bbarh!=jets->size(); bbarh++) { //loop over hadronic b-jet2
00098                   if (!(bbarh==p || bbarh==q || bbarh==j || bbarh==k) && !(bbarh==bh)) {
00099                     // Make event solutions for all possible combinations of the 4 light
00100                     // jets and 2 possible b-jets, not including the option of the b's being swapped. 
00101                     // Hadp,Hadq is one pair, Hadj,Hadk the other
00102                     std::vector<TtHadEvtSolution> asol;
00103                     asol.resize(3);
00104                     //[p][q][b] and [j][k][bbar]
00105                     asol[0].setJetCorrectionScheme(jetCorrScheme_);
00106                     asol[0].setHadp(jets, p);
00107                     asol[0].setHadq(jets, q);
00108                     asol[0].setHadj(jets, j);
00109                     asol[0].setHadk(jets, k);
00110                     asol[0].setHadb(jets, bh);
00111                     asol[0].setHadbbar(jets, bbarh);
00112 
00113                     //[p][j][b] and [q][k][bbar]
00114                     asol[1].setJetCorrectionScheme(jetCorrScheme_);
00115                     asol[1].setHadp(jets, p);
00116                     asol[1].setHadq(jets, j);
00117                     asol[1].setHadj(jets, q);
00118                     asol[1].setHadk(jets, k);
00119                     asol[1].setHadb(jets, bh);
00120                     asol[1].setHadbbar(jets, bbarh);
00121 
00122                     //[p][k][b] and [j][q][bbar]
00123                     asol[2].setJetCorrectionScheme(jetCorrScheme_);
00124                     asol[2].setHadp(jets, p);
00125                     asol[2].setHadq(jets, k);
00126                     asol[2].setHadj(jets, j);
00127                     asol[2].setHadk(jets, q);
00128                     asol[2].setHadb(jets, bh);
00129                     asol[2].setHadbbar(jets, bbarh);
00130                    
00131                     if(doKinFit_){
00132                       for(unsigned int i=0;i!=asol.size();i++){
00133                         asol[i] = myKinFitter->addKinFitInfo(&(asol[i]));
00134                         asol[i].setJetParametrisation(jetParam_);
00135                       }
00136                      
00137                     }else{
00138                       std::cout<<"Fitting needed to decide on best solution, enable fitting!"<<std::endl;
00139                     }
00140                     // these lines calculate the observables to be used in the TtHadSignalSelection LR
00141 
00142                     for(unsigned int i=0;i!=asol.size();i++){ 
00143                     (*myLRSignalSelObservables)(asol[i]);
00144                     }
00145                     // if asked for, calculate with these observable values the LRvalue and 
00146                     // (depending on the configuration) probability this event is signal
00147                     if(addLRSignalSel_){
00148                       for(unsigned int i=0;i!=asol.size();i++){
00149                         (*myLRSignalSelCalc)(asol[i]);
00150                       }
00151                     }
00152                     
00153                     // these lines calculate the observables to be used in the TtHadJetCombination LR
00154                     for(unsigned int i=0;i!=asol.size();i++){
00155                     (*myLRJetCombObservables)(asol[i]);
00156                     } 
00157                     // if asked for, calculate with these observable values the LRvalue and 
00158                     // (depending on the configuration) probability a jet combination is correct
00159                     if(addLRJetComb_){
00160                       for(unsigned int i=0;i!=asol.size();i++){
00161                         (*myLRJetCombCalc)(asol[i]);
00162                       }
00163                     }
00164                     //std::cout<<"SignalSelLRval = "<<asol.getLRSignalEvtLRval()<<"  JetCombProb = "<<asol.getLRSignalEvtProb()<<std::endl;
00165                     //std::cout<<"JetCombLRval = "<<asol.getLRJetCombLRval()<<"  JetCombProb = "<<asol.getLRJetCombProb()<<std::endl;
00166                     
00167                     // fill solution to vector with all possible solutions 
00168                     for(unsigned int i=0;i!=asol.size();i++){
00169                       evtsols->push_back(asol[i]);
00170                     }
00171                   }
00172                 }
00173               }
00174             }
00175           }
00176         }
00177       }
00178     }
00179 
00180     
00181     // add TtHadSimpleBestJetComb to solutions
00182     int simpleBestJetComb = (*mySimpleBestJetComb)(*evtsols);
00183 
00184     for(size_t s=0; s<evtsols->size(); s++){
00185       (*evtsols)[s].setSimpleBestJetComb(simpleBestJetComb);    
00186       // if asked for, match the event solutions to the gen Event
00187       if(matchToGenEvt_){
00188         int bestSolution = -999; 
00189         int bestSolutionChangeW1Q = -999;
00190         int bestSolutionChangeW2Q = -999;
00191         edm::Handle<TtGenEvent> genEvt;
00192         iEvent.getByLabel ("genEvt",genEvt); 
00193         std::vector<const reco::Candidate*> quarks;
00194         const reco::Candidate & genp  = *(genEvt->daughterQuarkOfWPlus());
00195         const reco::Candidate & genq  = *(genEvt->daughterQuarkBarOfWPlus());
00196         const reco::Candidate & genb  = *(genEvt->b());
00197         const reco::Candidate & genj  = *(genEvt->daughterQuarkOfWMinus());
00198         const reco::Candidate & genk  = *(genEvt->daughterQuarkBarOfWMinus());
00199         const reco::Candidate & genbbar = *(genEvt->bBar());
00200         quarks.push_back( &genp );       
00201         quarks.push_back( &genq );   
00202         quarks.push_back( &genb );
00203         quarks.push_back( &genj );       
00204         quarks.push_back( &genk );   
00205         quarks.push_back( &genbbar );
00206         std::vector<const reco::Candidate*> jets;         
00207         for(size_t s=0; s<evtsols->size(); s++) {
00208           jets.clear();     
00209           const reco::Candidate & jetp  = (*evtsols)[s].getRecHadp();
00210           const reco::Candidate & jetq  = (*evtsols)[s].getRecHadq();
00211           const reco::Candidate & jetbh = (*evtsols)[s].getRecHadb();
00212           const reco::Candidate & jetj  = (*evtsols)[s].getRecHadj();
00213           const reco::Candidate & jetk  = (*evtsols)[s].getRecHadk();
00214           const reco::Candidate & jetbbar = (*evtsols)[s].getRecHadbbar();
00215           jets.push_back( &jetp );      
00216           jets.push_back( &jetq );        
00217           jets.push_back( &jetbh );
00218           jets.push_back( &jetj );
00219           jets.push_back( &jetk );
00220           jets.push_back( &jetbbar );
00221           JetPartonMatching aMatch(quarks, jets, matchingAlgo_, useMaxDist_, useDeltaR_, maxDist_);  
00222           (*evtsols)[s].setGenEvt(genEvt);   
00223           (*evtsols)[s].setMCBestSumAngles(aMatch.getSumDistances());
00224           (*evtsols)[s].setMCBestAngleHadp(aMatch.getDistanceForParton(0));
00225           (*evtsols)[s].setMCBestAngleHadq(aMatch.getDistanceForParton(1));
00226           (*evtsols)[s].setMCBestAngleHadb(aMatch.getDistanceForParton(2));
00227           (*evtsols)[s].setMCBestAngleHadb(aMatch.getDistanceForParton(2));
00228           (*evtsols)[s].setMCBestAngleHadj(aMatch.getDistanceForParton(3));
00229           (*evtsols)[s].setMCBestAngleHadk(aMatch.getDistanceForParton(4));
00230           (*evtsols)[s].setMCBestAngleHadbbar(aMatch.getDistanceForParton(5));
00231           
00232           // Check match - checking if two light quarks are swapped wrt matched gen particle
00233           if((aMatch.getMatchForParton(2) == 2 && aMatch.getMatchForParton(5) == 5)
00234              || (aMatch.getMatchForParton(2) == 5 && aMatch.getMatchForParton(5) == 2)){ // check b-jets
00235             
00236             if(aMatch.getMatchForParton(3) == 3 && aMatch.getMatchForParton(4) == 4){ //check light jets
00237               bestSolutionChangeW2Q = 0;
00238               if(aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) { 
00239                 bestSolution = s;
00240                 bestSolutionChangeW1Q = 0;
00241               }else{
00242                 if(aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0){
00243                   bestSolution = s;
00244                   bestSolutionChangeW1Q = 1;
00245                 }
00246               }
00247             }else{
00248               if(aMatch.getMatchForParton(2) == 3 && aMatch.getMatchForParton(3) == 2){ // or check if swapped 
00249                 bestSolutionChangeW2Q = 1;
00250                 if(aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0){
00251                   bestSolution = s;
00252                   bestSolutionChangeW1Q = 1;
00253                 }else{
00254                   if(aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) { 
00255                     bestSolution = s;
00256                     bestSolutionChangeW1Q = 0;
00257                   }
00258                 }
00259               }
00260               if(aMatch.getMatchForParton(2) == 2 && aMatch.getMatchForParton(3) == 3){
00261                 bestSolutionChangeW2Q = 0;
00262                 if(aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) {
00263                   bestSolution = s; 
00264                   bestSolutionChangeW1Q = 0;
00265                 } else if(aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0) {
00266                   bestSolution = s;             
00267                   bestSolutionChangeW1Q = 1;           
00268                 }
00269               }
00270             }
00271           }
00272           for(size_t s=0; s<evtsols->size(); s++) {
00273             (*evtsols)[s].setMCBestJetComb(bestSolution);
00274             (*evtsols)[s].setMCChangeW1Q(bestSolutionChangeW1Q);
00275             (*evtsols)[s].setMCChangeW2Q(bestSolutionChangeW2Q);
00276           }
00277         }
00278       } // end matchEvt
00279     }       
00280     //store the vector of solutions to the event  
00281  
00282     std::auto_ptr<std::vector<TtHadEvtSolution> > pOut(evtsols);
00283     iEvent.put(pOut);
00284   }else {     //end loop jet/MET found
00285     std::cout<<"No calibrated solutions built, because only "<<jets->size()<<" were present";
00286      
00287     std::auto_ptr<std::vector<TtHadEvtSolution> > pOut(evtsols);
00288     iEvent.put(pOut);
00289   }
00290 }  
00291
/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/TopQuarkAnalysis/TopEventProducers/src/TtHadEvtSolutionMaker.cc
