CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/TopQuarkAnalysis/TopJetCombination/plugins/TtSemiLepJetCombMVAComputer.cc

Go to the documentation of this file.
00001 #include "PhysicsTools/JetMCUtils/interface/combination.h"
00002 
00003 #include "DataFormats/PatCandidates/interface/Jet.h"
00004 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00005 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
00006 #include "TopQuarkAnalysis/TopJetCombination/interface/TtSemiLepJetCombEval.h"
00007 #include "TopQuarkAnalysis/TopJetCombination/plugins/TtSemiLepJetCombMVAComputer.h"
00008 
00009 TtSemiLepJetCombMVAComputer::TtSemiLepJetCombMVAComputer(const edm::ParameterSet& cfg):
00010   leps_    (cfg.getParameter<edm::InputTag>("leps")),
00011   jets_    (cfg.getParameter<edm::InputTag>("jets")),
00012   mets_    (cfg.getParameter<edm::InputTag>("mets")),
00013   maxNJets_(cfg.getParameter<int>("maxNJets")),
00014   maxNComb_(cfg.getParameter<int>("maxNComb"))
00015 {
00016   produces<std::vector<std::vector<int> > >();
00017   produces<std::vector<double>            >("Discriminators");
00018   produces<std::string                    >("Method");
00019   produces<int                            >("NumberOfConsideredJets");
00020 }
00021 
00022 TtSemiLepJetCombMVAComputer::~TtSemiLepJetCombMVAComputer()
00023 {
00024 }
00025 
00026 void
00027 TtSemiLepJetCombMVAComputer::produce(edm::Event& evt, const edm::EventSetup& setup)
00028 {
00029   std::auto_ptr<std::vector<std::vector<int> > > pOut    (new std::vector<std::vector<int> >);
00030   std::auto_ptr<std::vector<double>            > pOutDisc(new std::vector<double>);
00031   std::auto_ptr<std::string                    > pOutMeth(new std::string);
00032   std::auto_ptr<int                            > pJetsConsidered(new int);
00033 
00034   mvaComputer.update<TtSemiLepJetCombMVARcd>(setup, "ttSemiLepJetCombMVA");
00035 
00036   // read name of the processor that provides the MVA discriminator
00037   // (to be used as meta information)
00038   edm::ESHandle<PhysicsTools::Calibration::MVAComputerContainer> calibContainer;
00039   setup.get<TtSemiLepJetCombMVARcd>().get( calibContainer );
00040   std::vector<PhysicsTools::Calibration::VarProcessor*> processors
00041     = (calibContainer->find("ttSemiLepJetCombMVA")).getProcessors();
00042   *pOutMeth = ( processors[ processors.size()-3 ] )->getInstanceName();
00043   evt.put(pOutMeth, "Method");
00044 
00045   // get lepton, jets and mets
00046   edm::Handle< edm::View<reco::RecoCandidate> > leptons; 
00047   evt.getByLabel(leps_, leptons);
00048 
00049   edm::Handle< std::vector<pat::Jet> > jets;
00050   evt.getByLabel(jets_, jets);
00051 
00052   edm::Handle< std::vector<pat::MET> > mets;
00053   evt.getByLabel(mets_, mets);
00054 
00055   const unsigned int nPartons = 4;
00056 
00057   // skip events with no appropriate lepton candidate,
00058   // empty METs vector or less jets than partons
00059   if( leptons->empty() || mets->empty() || jets->size() < nPartons ) {
00060     std::vector<int> invalidCombi;
00061     for(unsigned int i = 0; i < nPartons; ++i) 
00062       invalidCombi.push_back( -1 );
00063     pOut->push_back( invalidCombi );
00064     evt.put(pOut);
00065     pOutDisc->push_back( 0. );
00066     evt.put(pOutDisc, "Discriminators");
00067     *pJetsConsidered = jets->size();
00068     evt.put(pJetsConsidered, "NumberOfConsideredJets");
00069     return;
00070   }
00071 
00072   const math::XYZTLorentzVector lepton = leptons->begin()->p4();
00073 
00074   const pat::MET *met = &(*mets)[0];
00075 
00076   // analyze jet combinations
00077   std::vector<int> jetIndices;
00078   for(unsigned int i=0; i<jets->size(); ++i){
00079     if(maxNJets_ >= (int) nPartons && maxNJets_ == (int) i) {
00080       *pJetsConsidered = i;
00081       break;
00082     }
00083     jetIndices.push_back(i);
00084   }
00085   
00086   std::vector<int> combi;
00087   for(unsigned int i=0; i<nPartons; ++i) 
00088     combi.push_back(i);
00089 
00090   typedef std::pair<double, std::vector<int> > discCombPair;
00091   std::list<discCombPair> discCombList;
00092 
00093   do{
00094     for(int cnt = 0; cnt < TMath::Factorial( combi.size() ); ++cnt){
00095       // take into account indistinguishability of the two jets from the hadr. W decay,
00096       // reduces combinatorics by a factor of 2
00097       if(combi[TtSemiLepEvtPartons::LightQ] < combi[TtSemiLepEvtPartons::LightQBar]) {
00098 
00099         TtSemiLepJetComb jetComb(*jets, combi, lepton, *met);
00100 
00101         // feed MVA input variables into a ValueList
00102         PhysicsTools::Variable::ValueList values;
00103         evaluateTtSemiLepJetComb(values, jetComb);
00104 
00105         // get discriminator from the MVAComputer
00106         double discrim = mvaComputer->eval( values );
00107 
00108         discCombList.push_back( std::make_pair(discrim, combi) );
00109 
00110       }
00111       next_permutation( combi.begin() , combi.end() );
00112     }
00113   }
00114   while(stdcomb::next_combination( jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end() ));
00115 
00116   // sort results w.r.t. discriminator values
00117   discCombList.sort();
00118 
00119   // write result into the event
00120   // (starting with the JetComb having the highest discriminator value -> reverse iterator)
00121   unsigned int iDiscComb = 0;
00122   typedef std::list<discCombPair>::reverse_iterator discCombIterator;
00123   for(discCombIterator discCombPair = discCombList.rbegin(); discCombPair != discCombList.rend(); ++discCombPair) {
00124     if(maxNComb_ >= 1 && iDiscComb == (unsigned int) maxNComb_) break;
00125     pOut    ->push_back( discCombPair->second );
00126     pOutDisc->push_back( discCombPair->first  );
00127     iDiscComb++;
00128   }
00129   evt.put(pOut);
00130   evt.put(pOutDisc, "Discriminators");
00131   evt.put(pJetsConsidered, "NumberOfConsideredJets");
00132 }
00133 
00134 void 
00135 TtSemiLepJetCombMVAComputer::beginJob()
00136 {
00137 }
00138 
00139 void 
00140 TtSemiLepJetCombMVAComputer::endJob()
00141 {
00142 }
00143 
00144 // implement the plugins for the computer container
00145 // -> register TtSemiLepJetCombMVARcd
00146 // -> define TtSemiLepJetCombMVAFileSource
00147 MVA_COMPUTER_CONTAINER_IMPLEMENT(TtSemiLepJetCombMVA);