CMSSW Reference Manual

00001 //
00002 // $Id: TtSemiEvtSolutionMaker.cc,v 1.40.2.1.2.1 2009/04/21 19:50:27 rwolf Exp $
00003 //
00004 
00005 #include "TopQuarkAnalysis/TopEventProducers/interface/TtSemiEvtSolutionMaker.h"
00006 
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008 
00009 #include "PhysicsTools/Utilities/interface/deltaR.h"
00010 
00011 #include "TopQuarkAnalysis/TopTools/interface/JetPartonMatching.h"
00012 #include "AnalysisDataFormats/TopObjects/interface/TtSemiEvtSolution.h"
00013 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiSimpleBestJetComb.h"
00014 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiLRJetCombObservables.h"
00015 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiLRJetCombCalc.h"
00016 #include "TopQuarkAnalysis/TopEventSelection/interface/TtSemiLRSignalSelObservables.h"
00017 #include "TopQuarkAnalysis/TopEventSelection/interface/TtSemiLRSignalSelCalc.h"
00018 
00019 #include <memory>
00020 
00021 
00023 TtSemiEvtSolutionMaker::TtSemiEvtSolutionMaker(const edm::ParameterSet & iConfig) {
00024   // configurables
00025   electronSrc_     = iConfig.getParameter<edm::InputTag>    ("electronSource");
00026   muonSrc_         = iConfig.getParameter<edm::InputTag>    ("muonSource");
00027   metSrc_          = iConfig.getParameter<edm::InputTag>    ("metSource");
00028   jetSrc_          = iConfig.getParameter<edm::InputTag>    ("jetSource");
00029   leptonFlavour_   = iConfig.getParameter<std::string>      ("leptonFlavour");
00030   jetCorrScheme_   = iConfig.getParameter<int>              ("jetCorrectionScheme");
00031   nrCombJets_      = iConfig.getParameter<unsigned int>     ("nrCombJets");
00032   doKinFit_        = iConfig.getParameter<bool>             ("doKinFit");
00033   addLRSignalSel_  = iConfig.getParameter<bool>             ("addLRSignalSel");
00034   lrSignalSelObs_  = iConfig.getParameter<std::vector<int> >("lrSignalSelObs");
00035   lrSignalSelFile_ = iConfig.getParameter<std::string>      ("lrSignalSelFile");
00036   addLRJetComb_    = iConfig.getParameter<bool>             ("addLRJetComb");
00037   lrJetCombObs_    = iConfig.getParameter<std::vector<int> >("lrJetCombObs");
00038   lrJetCombFile_   = iConfig.getParameter<std::string>      ("lrJetCombFile");
00039   maxNrIter_       = iConfig.getParameter<int>              ("maxNrIter");
00040   maxDeltaS_       = iConfig.getParameter<double>           ("maxDeltaS");
00041   maxF_            = iConfig.getParameter<double>           ("maxF");
00042   jetParam_        = iConfig.getParameter<int>              ("jetParametrisation");
00043   lepParam_        = iConfig.getParameter<int>              ("lepParametrisation");
00044   metParam_        = iConfig.getParameter<int>              ("metParametrisation");
00045   constraints_     = iConfig.getParameter<std::vector<unsigned> >("constraints");
00046   matchToGenEvt_   = iConfig.getParameter<bool>             ("matchToGenEvt");
00047   matchingAlgo_    = iConfig.getParameter<int>              ("matchingAlgorithm");
00048   useMaxDist_      = iConfig.getParameter<bool>             ("useMaximalDistance");
00049   useDeltaR_       = iConfig.getParameter<bool>             ("useDeltaR");
00050   maxDist_         = iConfig.getParameter<double>           ("maximalDistance");
00051 
00052   // define kinfitter
00053   if(doKinFit_){
00054     myKinFitter = new TtSemiLepKinFitter(param(jetParam_), param(lepParam_), param(metParam_), maxNrIter_, maxDeltaS_, maxF_, constraints(constraints_));
00055   }
00056 
00057   // define jet combinations related calculators
00058   mySimpleBestJetComb                    = new TtSemiSimpleBestJetComb();
00059   myLRSignalSelObservables               = new TtSemiLRSignalSelObservables();
00060   myLRJetCombObservables                 = new TtSemiLRJetCombObservables();
00061   myLRJetCombObservables -> jetSource(jetSrc_);
00062   if (addLRJetComb_)   myLRJetCombCalc   = new TtSemiLRJetCombCalc(edm::FileInPath(lrJetCombFile_).fullPath(), lrJetCombObs_);
00063 
00064   // instantiate signal selection calculator
00065   if (addLRSignalSel_) myLRSignalSelCalc = new TtSemiLRSignalSelCalc(edm::FileInPath(lrSignalSelFile_).fullPath(), lrSignalSelObs_);
00066 
00067   // define what will be produced
00068   produces<std::vector<TtSemiEvtSolution> >();
00069 }
00070 
00071 
00073 TtSemiEvtSolutionMaker::~TtSemiEvtSolutionMaker() {
00074   if (doKinFit_)      delete myKinFitter;
00075   delete mySimpleBestJetComb;
00076   delete myLRSignalSelObservables;
00077   delete myLRJetCombObservables;
00078   if(addLRSignalSel_) delete myLRSignalSelCalc;
00079   if(addLRJetComb_)   delete myLRJetCombCalc;
00080 }
00081 
00082 
00083 void TtSemiEvtSolutionMaker::produce(edm::Event & iEvent, const edm::EventSetup & iSetup) {
00084 
00085   //
00086   //  TopObject Selection
00087   //
00088 
00089   // select lepton (the TtLepton vectors are, for the moment, sorted on pT)
00090   bool leptonFound = false;
00091   edm::Handle<std::vector<pat::Muon> > muons;
00092   if(leptonFlavour_ == "muon"){
00093     iEvent.getByLabel(muonSrc_, muons);
00094     if (muons->size() > 0) leptonFound = true;
00095   }
00096   edm::Handle<std::vector<pat::Electron> > electrons;
00097   if(leptonFlavour_ == "electron"){
00098     iEvent.getByLabel(electronSrc_, electrons);
00099     if (electrons->size() > 0) leptonFound = true;
00100   }  
00101 
00102   // select MET (TopMET vector is sorted on ET)
00103   bool metFound = false;
00104   edm::Handle<std::vector<pat::MET> > mets;
00105   iEvent.getByLabel(metSrc_, mets);
00106   if (mets->size() > 0) metFound = true;
00107 
00108   // select Jets
00109   bool jetsFound = false;
00110   edm::Handle<std::vector<pat::Jet> > jets;
00111   iEvent.getByLabel(jetSrc_, jets);
00112   if (jets->size() >=4) jetsFound = true;
00113 
00114   //
00115   // Build Event solutions according to the ambiguity in the jet combination
00116   //
00117   std::vector<TtSemiEvtSolution> * evtsols = new std::vector<TtSemiEvtSolution>();
00118   if(leptonFound && metFound && jetsFound){
00119     // protect against reading beyond array boundaries
00120     unsigned int nrCombJets = nrCombJets_; // do not overwrite nrCombJets_
00121     if (jets->size() < nrCombJets) nrCombJets = jets->size();
00122     // loop over all jets
00123     for (unsigned int p=0; p<nrCombJets; p++) {
00124       for (unsigned int q=0; q<nrCombJets; q++) {
00125         for (unsigned int bh=0; bh<nrCombJets; bh++) {
00126           if (q>p && !(bh==p || bh==q)) {
00127             for (unsigned int bl=0; bl<nrCombJets; bl++) {
00128               if (!(bl==p || bl==q || bl==bh)) {
00129                 TtSemiEvtSolution asol;
00130                 asol.setJetCorrectionScheme(jetCorrScheme_);
00131                 if(leptonFlavour_ == "muon")     asol.setMuon(muons, 0);
00132                 if(leptonFlavour_ == "electron") asol.setElectron(electrons, 0);
00133                 asol.setNeutrino(mets, 0);
00134                 asol.setHadp(jets, p);
00135                 asol.setHadq(jets, q);
00136                 asol.setHadb(jets, bh);
00137                 asol.setLepb(jets, bl);
00138                 if (doKinFit_) {
00139                   asol = myKinFitter->addKinFitInfo(&asol);
00140                   // just to keep a record in the event (drop? -> present in provenance anyway...)
00141                   asol.setJetParametrisation(jetParam_);
00142                   asol.setLeptonParametrisation(lepParam_);
00143                   asol.setNeutrinoParametrisation(metParam_);
00144                 }
00145                 if(matchToGenEvt_){
00146                   edm::Handle<TtGenEvent> genEvt;
00147                   iEvent.getByLabel ("genEvt",genEvt); 
00148                   if( genEvt->isSemiLeptonic(false) ){
00149                     asol.setGenEvt(genEvt);   
00150                     
00151                   }
00152                 }
00153                 // these lines calculate the observables to be used in the TtSemiSignalSelection LR
00154                 (*myLRSignalSelObservables)(asol, *jets);
00155 
00156                 // if asked for, calculate with these observable values the LRvalue and 
00157                 // (depending on the configuration) probability this event is signal
00158                 // FIXME: DO WE NEED TO DO THIS FOR EACH SOLUTION??? (S.L.20/8/07)
00159                 if(addLRSignalSel_) (*myLRSignalSelCalc)(asol);
00160 
00161                 // these lines calculate the observables to be used in the TtSemiJetCombination LR
00162                 //(*myLRJetCombObservables)(asol);
00163                 (*myLRJetCombObservables)(asol, iEvent);
00164 
00165                 // if asked for, calculate with these observable values the LRvalue and 
00166                 // (depending on the configuration) probability a jet combination is correct
00167                 if(addLRJetComb_) (*myLRJetCombCalc)(asol);
00168 
00169                 //std::cout<<"SignalSelLRval = "<<asol.getLRSignalEvtLRval()<<"  JetCombProb = "<<asol.getLRSignalEvtProb()<<std::endl;
00170                 //std::cout<<"JetCombLRval = "<<asol.getLRJetCombLRval()<<"  JetCombProb = "<<asol.getLRJetCombProb()<<std::endl;
00171 
00172                 // fill solution to vector
00173                 asol.setupHyp();
00174                 evtsols->push_back(asol);
00175               } 
00176             }
00177           }
00178         } 
00179       }
00180     }
00181     // if asked for, match the event solutions to the gen Event
00182     if(matchToGenEvt_){
00183       int bestSolution = -999; 
00184       int bestSolutionChangeWQ = -999;
00185       edm::Handle<TtGenEvent> genEvt;
00186       iEvent.getByLabel ("genEvt",genEvt); 
00187       if( genEvt->isSemiLeptonic(false) ){
00188         vector<const reco::Candidate*> quarks;
00189         const reco::Candidate & genp  = *(genEvt->hadronicDecayQuark());
00190         const reco::Candidate & genq  = *(genEvt->hadronicDecayQuarkBar());
00191         const reco::Candidate & genbh = *(genEvt->hadronicDecayB(false));
00192         const reco::Candidate & genbl = *(genEvt->leptonicDecayB(false));
00193         quarks.push_back( &genp );
00194         quarks.push_back( &genq );
00195         quarks.push_back( &genbh );
00196         quarks.push_back( &genbl );
00197         vector<const reco::Candidate*> recjets;  
00198         for(size_t s=0; s<evtsols->size(); s++) {
00199           recjets.clear();
00200           const reco::Candidate & jetp  = (*evtsols)[s].getRecHadp();
00201           const reco::Candidate & jetq  = (*evtsols)[s].getRecHadq();
00202           const reco::Candidate & jetbh = (*evtsols)[s].getRecHadb();
00203           const reco::Candidate & jetbl = (*evtsols)[s].getRecLepb();
00204           recjets.push_back( &jetp );
00205           recjets.push_back( &jetq );
00206           recjets.push_back( &jetbh );
00207           recjets.push_back( &jetbl );
00208 
00209           JetPartonMatching aMatch(quarks, recjets, matchingAlgo_, useMaxDist_, useDeltaR_, maxDist_);   
00210           (*evtsols)[s].setGenEvt(genEvt);   
00211           (*evtsols)[s].setMCBestSumAngles(aMatch.getSumDistances());
00212           (*evtsols)[s].setMCBestAngleHadp(aMatch.getDistanceForParton(0));
00213           (*evtsols)[s].setMCBestAngleHadq(aMatch.getDistanceForParton(1));
00214           (*evtsols)[s].setMCBestAngleHadb(aMatch.getDistanceForParton(2));
00215           (*evtsols)[s].setMCBestAngleLepb(aMatch.getDistanceForParton(3));
00216 
00217           if(aMatch.getMatchForParton(2) == 2 && aMatch.getMatchForParton(3) == 3){
00218             if(aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) {
00219               bestSolution = s;
00220               bestSolutionChangeWQ = 0;
00221             } else if(aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0) {
00222               bestSolution = s;
00223               bestSolutionChangeWQ = 1;
00224             }
00225           }
00226         }
00227       }
00228       for(size_t s=0; s<evtsols->size(); s++) {
00229         (*evtsols)[s].setMCBestJetComb(bestSolution);
00230         (*evtsols)[s].setMCChangeWQ(bestSolutionChangeWQ);     
00231       }
00232     }
00233 
00234     // add TtSemiSimpleBestJetComb to solutions
00235     int simpleBestJetComb = (*mySimpleBestJetComb)(*evtsols);
00236     for(size_t s=0; s<evtsols->size(); s++) (*evtsols)[s].setSimpleBestJetComb(simpleBestJetComb);
00237     // choose the best jet combination according to LR value
00238     if (addLRJetComb_ && evtsols->size()>0) {
00239       float bestLRVal = -1000000;
00240       int bestSol = (*evtsols)[0].getLRBestJetComb(); // duplicate the default
00241       for(size_t s=0; s<evtsols->size(); s++) {
00242         if ((*evtsols)[s].getLRJetCombLRval() > bestLRVal) {
00243           bestLRVal = (*evtsols)[s].getLRJetCombLRval();
00244           bestSol = s;
00245         }
00246       }
00247       for(size_t s=0; s<evtsols->size(); s++) {
00248         (*evtsols)[s].setLRBestJetComb(bestSol);
00249       }
00250     }
00251     //store the vector of solutions to the event     
00252     std::auto_ptr<std::vector<TtSemiEvtSolution> > pOut(evtsols);
00253     iEvent.put(pOut);
00254 
00255   } else {
00256 
00257     /*
00258     std::cout<<"No calibrated solutions built, because:  ";
00259     if(jets->size()<4)                                            std::cout<<"nr sel jets < 4"<<std::endl;
00260     if(leptonFlavour_ == "muon" && muons->size() == 0)            std::cout<<"no good muon candidate"<<std::endl;
00261     if(leptonFlavour_ == "electron" && electrons->size() == 0)   std::cout<<"no good electron candidate"<<std::endl;
00262     if(mets->size() == 0)                                         std::cout<<"no MET reconstruction"<<std::endl;
00263     */
00264 //    TtSemiEvtSolution asol;
00265 //    evtsols->push_back(asol);
00266     std::auto_ptr<std::vector<TtSemiEvtSolution> > pOut(evtsols);
00267     iEvent.put(pOut);
00268   }
00269 }
00270 
00271 TtSemiLepKinFitter::Param TtSemiEvtSolutionMaker::param(unsigned val) 
00272 {
00273   TtSemiLepKinFitter::Param result;
00274   switch(val){
00275   case TtSemiLepKinFitter::kEMom       : result=TtSemiLepKinFitter::kEMom;       break;
00276   case TtSemiLepKinFitter::kEtEtaPhi   : result=TtSemiLepKinFitter::kEtEtaPhi;   break;
00277   case TtSemiLepKinFitter::kEtThetaPhi : result=TtSemiLepKinFitter::kEtThetaPhi; break;
00278   default: 
00279     throw cms::Exception("WrongConfig") 
00280       << "Chosen jet parametrization is not supported: " << val << "\n";
00281     break;
00282   }
00283   return result;
00284 } 
00285 
00286 TtSemiLepKinFitter::Constraint TtSemiEvtSolutionMaker::constraint(unsigned val) 
00287 {
00288   TtSemiLepKinFitter::Constraint result;
00289   switch(val){
00290   case TtSemiLepKinFitter::kWHadMass     : result=TtSemiLepKinFitter::kWHadMass;     break;
00291   case TtSemiLepKinFitter::kWLepMass     : result=TtSemiLepKinFitter::kWLepMass;     break;
00292   case TtSemiLepKinFitter::kTopHadMass   : result=TtSemiLepKinFitter::kTopHadMass;   break;
00293   case TtSemiLepKinFitter::kTopLepMass   : result=TtSemiLepKinFitter::kTopLepMass;   break;
00294   case TtSemiLepKinFitter::kNeutrinoMass : result=TtSemiLepKinFitter::kNeutrinoMass; break;
00295   default: 
00296     throw cms::Exception("WrongConfig") 
00297       << "Chosen fit constraint is not supported: " << val << "\n";
00298     break;
00299   }
00300   return result;
00301 } 
00302 
00303 std::vector<TtSemiLepKinFitter::Constraint> TtSemiEvtSolutionMaker::constraints(std::vector<unsigned>& val)
00304 {
00305   std::vector<TtSemiLepKinFitter::Constraint> result;
00306   for(unsigned i=0; i<val.size(); ++i){
00307     result.push_back(constraint(val[i]));
00308   } 
00309   return result;
00310 }
00311