Reference Manual

Go to the documentation of this file.
00001 //
00002 // $Id: TtSemiEvtSolutionMaker.cc,v 1.43 2010/03/30 14:04:44 snaumann Exp $
00003 //
00004 
00005 #include "TopQuarkAnalysis/TopEventProducers/interface/TtSemiEvtSolutionMaker.h"
00006 
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008 
00009 #include "DataFormats/Math/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->numberOfBQuarks() == 2 &&  // FIXME: in rare cases W->bc decay, resulting in a wrong filled genEvt leading to a segmentation fault 
00149                       genEvt->numberOfLeptons() == 1) {  // FIXME: temporary solution to avoid crash in JetPartonMatching for non semi-leptonic events
00150                     asol.setGenEvt(genEvt);   
00151                   
00152                   }
00153                 }
00154                 // these lines calculate the observables to be used in the TtSemiSignalSelection LR
00155                 (*myLRSignalSelObservables)(asol, *jets);
00156 
00157                 // if asked for, calculate with these observable values the LRvalue and 
00158                 // (depending on the configuration) probability this event is signal
00159                 // FIXME: DO WE NEED TO DO THIS FOR EACH SOLUTION??? (S.L.20/8/07)
00160                 if(addLRSignalSel_) (*myLRSignalSelCalc)(asol);
00161 
00162                 // these lines calculate the observables to be used in the TtSemiJetCombination LR
00163                 //(*myLRJetCombObservables)(asol);
00164                 
00165                 (*myLRJetCombObservables)(asol, iEvent);
00166                 
00167                 // if asked for, calculate with these observable values the LRvalue and 
00168                 // (depending on the configuration) probability a jet combination is correct
00169                 if(addLRJetComb_) (*myLRJetCombCalc)(asol);
00170 
00171                 //std::cout<<"SignalSelLRval = "<<asol.getLRSignalEvtLRval()<<"  JetCombProb = "<<asol.getLRSignalEvtProb()<<std::endl;
00172                 //std::cout<<"JetCombLRval = "<<asol.getLRJetCombLRval()<<"  JetCombProb = "<<asol.getLRJetCombProb()<<std::endl;
00173 
00174                 // fill solution to vector
00175                 asol.setupHyp();
00176                 evtsols->push_back(asol);
00177               } 
00178             }
00179           }
00180         } 
00181       }
00182     }
00183 
00184     // if asked for, match the event solutions to the gen Event
00185     if(matchToGenEvt_){
00186       int bestSolution = -999; 
00187       int bestSolutionChangeWQ = -999;
00188       edm::Handle<TtGenEvent> genEvt;
00189       iEvent.getByLabel ("genEvt",genEvt); 
00190       if (genEvt->numberOfBQuarks() == 2 &&   // FIXME: in rare cases W->bc decay, resulting in a wrong filled genEvt leading to a segmentation fault
00191           genEvt->numberOfLeptons() == 1) {   // FIXME: temporary solution to avoid crash in JetPartonMatching for non semi-leptonic events
00192         std::vector<const reco::Candidate*> quarks;
00193         const reco::Candidate & genp  = *(genEvt->hadronicDecayQuark());
00194         const reco::Candidate & genq  = *(genEvt->hadronicDecayQuarkBar());
00195         const reco::Candidate & genbh = *(genEvt->hadronicDecayB());
00196         const reco::Candidate & genbl = *(genEvt->leptonicDecayB());
00197         quarks.push_back( &genp );
00198         quarks.push_back( &genq );
00199         quarks.push_back( &genbh );
00200         quarks.push_back( &genbl );
00201         std::vector<const reco::Candidate*> recjets;  
00202         for(size_t s=0; s<evtsols->size(); s++) {
00203           recjets.clear();
00204           const reco::Candidate & jetp  = (*evtsols)[s].getRecHadp();
00205           const reco::Candidate & jetq  = (*evtsols)[s].getRecHadq();
00206           const reco::Candidate & jetbh = (*evtsols)[s].getRecHadb();
00207           const reco::Candidate & jetbl = (*evtsols)[s].getRecLepb();
00208           recjets.push_back( &jetp );
00209           recjets.push_back( &jetq );
00210           recjets.push_back( &jetbh );
00211           recjets.push_back( &jetbl );
00212           JetPartonMatching aMatch(quarks, recjets, matchingAlgo_, useMaxDist_, useDeltaR_, maxDist_);   
00213           (*evtsols)[s].setGenEvt(genEvt);   
00214           (*evtsols)[s].setMCBestSumAngles(aMatch.getSumDistances());
00215           (*evtsols)[s].setMCBestAngleHadp(aMatch.getDistanceForParton(0));
00216           (*evtsols)[s].setMCBestAngleHadq(aMatch.getDistanceForParton(1));
00217           (*evtsols)[s].setMCBestAngleHadb(aMatch.getDistanceForParton(2));
00218           (*evtsols)[s].setMCBestAngleLepb(aMatch.getDistanceForParton(3));
00219           if(aMatch.getMatchForParton(2) == 2 && aMatch.getMatchForParton(3) == 3){
00220             if(aMatch.getMatchForParton(0) == 0 && aMatch.getMatchForParton(1) == 1) {
00221               bestSolution = s;
00222               bestSolutionChangeWQ = 0;
00223             } else if(aMatch.getMatchForParton(0) == 1 && aMatch.getMatchForParton(1) == 0) {
00224               bestSolution = s;
00225               bestSolutionChangeWQ = 1;
00226             }
00227           }
00228         }
00229       }
00230       for(size_t s=0; s<evtsols->size(); s++) {
00231         (*evtsols)[s].setMCBestJetComb(bestSolution);
00232         (*evtsols)[s].setMCChangeWQ(bestSolutionChangeWQ);     
00233       }
00234     }
00235 
00236     // add TtSemiSimpleBestJetComb to solutions
00237     int simpleBestJetComb = (*mySimpleBestJetComb)(*evtsols);
00238     for(size_t s=0; s<evtsols->size(); s++) (*evtsols)[s].setSimpleBestJetComb(simpleBestJetComb);
00239 
00240     // choose the best jet combination according to LR value
00241     if (addLRJetComb_ && evtsols->size()>0) {
00242       float bestLRVal = -1000000;
00243       int bestSol = (*evtsols)[0].getLRBestJetComb(); // duplicate the default
00244       for(size_t s=0; s<evtsols->size(); s++) {
00245         if ((*evtsols)[s].getLRJetCombLRval() > bestLRVal) {
00246           bestLRVal = (*evtsols)[s].getLRJetCombLRval();
00247           bestSol = s;
00248         }
00249       }
00250       for(size_t s=0; s<evtsols->size(); s++) {
00251         (*evtsols)[s].setLRBestJetComb(bestSol);
00252       }
00253     }
00254 
00255     //store the vector of solutions to the event     
00256     std::auto_ptr<std::vector<TtSemiEvtSolution> > pOut(evtsols);
00257     iEvent.put(pOut);
00258 
00259   } else {
00260 
00261     /*
00262     std::cout<<"No calibrated solutions built, because:  ";
00263     if(jets->size()<4)                                            std::cout<<"nr sel jets < 4"<<std::endl;
00264     if(leptonFlavour_ == "muon" && muons->size() == 0)            std::cout<<"no good muon candidate"<<std::endl;
00265     if(leptonFlavour_ == "electron" && electrons->size() == 0)   std::cout<<"no good electron candidate"<<std::endl;
00266     if(mets->size() == 0)                                         std::cout<<"no MET reconstruction"<<std::endl;
00267     */
00268 //    TtSemiEvtSolution asol;
00269 //    evtsols->push_back(asol);
00270     std::auto_ptr<std::vector<TtSemiEvtSolution> > pOut(evtsols);
00271     iEvent.put(pOut);
00272   }
00273 }
00274 
00275 TtSemiLepKinFitter::Param TtSemiEvtSolutionMaker::param(unsigned val) 
00276 {
00277   TtSemiLepKinFitter::Param result;
00278   switch(val){
00279   case TtSemiLepKinFitter::kEMom       : result=TtSemiLepKinFitter::kEMom;       break;
00280   case TtSemiLepKinFitter::kEtEtaPhi   : result=TtSemiLepKinFitter::kEtEtaPhi;   break;
00281   case TtSemiLepKinFitter::kEtThetaPhi : result=TtSemiLepKinFitter::kEtThetaPhi; break;
00282   default: 
00283     throw cms::Exception("WrongConfig") 
00284       << "Chosen jet parametrization is not supported: " << val << "\n";
00285     break;
00286   }
00287   return result;
00288 } 
00289 
00290 TtSemiLepKinFitter::Constraint TtSemiEvtSolutionMaker::constraint(unsigned val) 
00291 {
00292   TtSemiLepKinFitter::Constraint result;
00293   switch(val){
00294   case TtSemiLepKinFitter::kWHadMass     : result=TtSemiLepKinFitter::kWHadMass;     break;
00295   case TtSemiLepKinFitter::kWLepMass     : result=TtSemiLepKinFitter::kWLepMass;     break;
00296   case TtSemiLepKinFitter::kTopHadMass   : result=TtSemiLepKinFitter::kTopHadMass;   break;
00297   case TtSemiLepKinFitter::kTopLepMass   : result=TtSemiLepKinFitter::kTopLepMass;   break;
00298   case TtSemiLepKinFitter::kNeutrinoMass : result=TtSemiLepKinFitter::kNeutrinoMass; break;
00299   default: 
00300     throw cms::Exception("WrongConfig") 
00301       << "Chosen fit constraint is not supported: " << val << "\n";
00302     break;
00303   }
00304   return result;
00305 } 
00306 
00307 std::vector<TtSemiLepKinFitter::Constraint> TtSemiEvtSolutionMaker::constraints(std::vector<unsigned>& val)
00308 {
00309   std::vector<TtSemiLepKinFitter::Constraint> result;
00310   for(unsigned i=0; i<val.size(); ++i){
00311     result.push_back(constraint(val[i]));
00312   } 
00313   return result;
00314 }
00315
/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/TopQuarkAnalysis/TopEventProducers/src/TtSemiEvtSolutionMaker.cc
