Reference Manual

Go to the documentation of this file.
00001 #ifndef TtSemiLepKinFitProducer_h
00002 #define TtSemiLepKinFitProducer_h
00003 
00004 #include "FWCore/Framework/interface/Event.h"
00005 #include "FWCore/Framework/interface/EDProducer.h"
00006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00007 
00008 #include "PhysicsTools/JetMCUtils/interface/combination.h"
00009 #include "AnalysisDataFormats/TopObjects/interface/TtSemiLepEvtPartons.h"
00010 #include "TopQuarkAnalysis/TopKinFitter/interface/TtSemiLepKinFitter.h"
00011 
00012 template <typename LeptonCollection>
00013 class TtSemiLepKinFitProducer : public edm::EDProducer {
00014   
00015  public:
00016   
00017   explicit TtSemiLepKinFitProducer(const edm::ParameterSet&);
00018   ~TtSemiLepKinFitProducer();
00019   
00020  private:
00021   // produce
00022   virtual void produce(edm::Event&, const edm::EventSetup&);
00023 
00024   // convert unsigned to Param
00025   TtSemiLepKinFitter::Param param(unsigned);
00026   // convert unsigned to Param
00027   TtSemiLepKinFitter::Constraint constraint(unsigned);
00028   // convert unsigned to Param
00029   std::vector<TtSemiLepKinFitter::Constraint> constraints(std::vector<unsigned>&);
00030   // helper function for b-tagging
00031   bool doBTagging(bool& useBTag_, edm::Handle<std::vector<pat::Jet> >& jets, std::vector<int>& combi,
00032                   std::string& bTagAlgo_, double& minBTagValueBJets_, double& maxBTagValueNonBJets_);
00033 
00034   edm::InputTag jets_;
00035   edm::InputTag leps_;
00036   edm::InputTag mets_;
00037   
00038   edm::InputTag match_;
00040   bool useOnlyMatch_;
00042   std::string bTagAlgo_;
00044   double minBTagValueBJet_;
00046   double maxBTagValueNonBJet_;
00048   bool useBTag_;
00050   int maxNJets_;
00052   int maxNComb_;
00053 
00055   unsigned int maxNrIter_;
00057   double maxDeltaS_;
00059   double maxF_;
00060   unsigned int jetParam_;
00061   unsigned int lepParam_;
00062   unsigned int metParam_;
00064   std::vector<unsigned> constraints_;
00065   double mW_;
00066   double mTop_;
00068   double jetEnergyResolutionSmearFactor_;
00069 
00070   TtSemiLepKinFitter* fitter;
00071 
00072   struct KinFitResult {
00073     int Status;
00074     double Chi2;
00075     double Prob;
00076     pat::Particle HadB;
00077     pat::Particle HadP;
00078     pat::Particle HadQ;
00079     pat::Particle LepB;
00080     pat::Particle LepL;
00081     pat::Particle LepN;
00082     std::vector<int> JetCombi;
00083     bool operator< (const KinFitResult& rhs) { return Chi2 < rhs.Chi2; };
00084   };
00085 };
00086 
00087 template<typename LeptonCollection>
00088 TtSemiLepKinFitProducer<LeptonCollection>::TtSemiLepKinFitProducer(const edm::ParameterSet& cfg):
00089   jets_                    (cfg.getParameter<edm::InputTag>("jets")),
00090   leps_                    (cfg.getParameter<edm::InputTag>("leps")),
00091   mets_                    (cfg.getParameter<edm::InputTag>("mets")),
00092   match_                   (cfg.getParameter<edm::InputTag>("match")),
00093   useOnlyMatch_            (cfg.getParameter<bool>         ("useOnlyMatch"        )),
00094   bTagAlgo_                (cfg.getParameter<std::string>  ("bTagAlgo"            )),
00095   minBTagValueBJet_        (cfg.getParameter<double>       ("minBDiscBJets"       )),
00096   maxBTagValueNonBJet_     (cfg.getParameter<double>       ("maxBDiscLightJets"   )),
00097   useBTag_                 (cfg.getParameter<bool>         ("useBTagging"         )),
00098   maxNJets_                (cfg.getParameter<int>          ("maxNJets"            )),
00099   maxNComb_                (cfg.getParameter<int>          ("maxNComb"            )),
00100   maxNrIter_               (cfg.getParameter<unsigned>     ("maxNrIter"           )),
00101   maxDeltaS_               (cfg.getParameter<double>       ("maxDeltaS"           )),
00102   maxF_                    (cfg.getParameter<double>       ("maxF"                )),
00103   jetParam_                (cfg.getParameter<unsigned>     ("jetParametrisation"  )),
00104   lepParam_                (cfg.getParameter<unsigned>     ("lepParametrisation"  )),
00105   metParam_                (cfg.getParameter<unsigned>     ("metParametrisation"  )),
00106   constraints_             (cfg.getParameter<std::vector<unsigned> >("constraints")),
00107   mW_                      (cfg.getParameter<double>       ("mW"                  )),
00108   mTop_                    (cfg.getParameter<double>       ("mTop"                )),
00109   jetEnergyResolutionSmearFactor_(cfg.getParameter<double> ("jetEnergyResolutionSmearFactor"))
00110 
00111 {
00112   fitter = new TtSemiLepKinFitter(param(jetParam_), param(lepParam_), param(metParam_), maxNrIter_, maxDeltaS_, maxF_,
00113                                   constraints(constraints_), mW_, mTop_);
00114 
00115   produces< std::vector<pat::Particle> >("PartonsHadP");
00116   produces< std::vector<pat::Particle> >("PartonsHadQ");
00117   produces< std::vector<pat::Particle> >("PartonsHadB");
00118   produces< std::vector<pat::Particle> >("PartonsLepB");
00119   produces< std::vector<pat::Particle> >("Leptons");
00120   produces< std::vector<pat::Particle> >("Neutrinos");
00121 
00122   produces< std::vector<std::vector<int> > >();
00123   produces< std::vector<double> >("Chi2");
00124   produces< std::vector<double> >("Prob");
00125   produces< std::vector<int> >("Status");
00126 }
00127 
00128 template<typename LeptonCollection>
00129 TtSemiLepKinFitProducer<LeptonCollection>::~TtSemiLepKinFitProducer()
00130 {
00131   delete fitter;
00132 }
00133 
00134 template<typename LeptonCollection>
00135 bool TtSemiLepKinFitProducer<LeptonCollection>::doBTagging(bool& useBTag_, edm::Handle<std::vector<pat::Jet> >& jets, std::vector<int>& combi,
00136                                                            std::string& bTagAlgo_, double& minBTagValueBJet_, double& maxBTagValueNonBJet_){
00137   
00138   if( !useBTag_ ) {
00139     return true;
00140   }
00141   if( useBTag_ &&
00142       (*jets)[combi[TtSemiLepEvtPartons::HadB     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00143       (*jets)[combi[TtSemiLepEvtPartons::LepB     ]].bDiscriminator(bTagAlgo_) >= minBTagValueBJet_ &&
00144       (*jets)[combi[TtSemiLepEvtPartons::LightQ   ]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ &&
00145       (*jets)[combi[TtSemiLepEvtPartons::LightQBar]].bDiscriminator(bTagAlgo_) <  maxBTagValueNonBJet_ ) {
00146     return true;
00147   }
00148   else{
00149     return false;
00150   }
00151 }
00152 
00153 template<typename LeptonCollection>
00154 void TtSemiLepKinFitProducer<LeptonCollection>::produce(edm::Event& evt, const edm::EventSetup& setup)
00155 {
00156   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadP( new std::vector<pat::Particle> );
00157   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadQ( new std::vector<pat::Particle> );
00158   std::auto_ptr< std::vector<pat::Particle> > pPartonsHadB( new std::vector<pat::Particle> );
00159   std::auto_ptr< std::vector<pat::Particle> > pPartonsLepB( new std::vector<pat::Particle> );
00160   std::auto_ptr< std::vector<pat::Particle> > pLeptons    ( new std::vector<pat::Particle> );
00161   std::auto_ptr< std::vector<pat::Particle> > pNeutrinos  ( new std::vector<pat::Particle> );
00162 
00163   std::auto_ptr< std::vector<std::vector<int> > > pCombi ( new std::vector<std::vector<int> > );
00164   std::auto_ptr< std::vector<double>            > pChi2  ( new std::vector<double> );
00165   std::auto_ptr< std::vector<double>            > pProb  ( new std::vector<double> );
00166   std::auto_ptr< std::vector<int>               > pStatus( new std::vector<int> );
00167 
00168   edm::Handle<std::vector<pat::Jet> > jets;
00169   evt.getByLabel(jets_, jets);
00170 
00171   edm::Handle<std::vector<pat::MET> > mets;
00172   evt.getByLabel(mets_, mets);
00173 
00174   edm::Handle<LeptonCollection> leps;
00175   evt.getByLabel(leps_, leps);
00176 
00177   unsigned int nPartons = 4;
00178 
00179   std::vector<int> match;
00180   bool invalidMatch = false;
00181   if(useOnlyMatch_) {
00182     edm::Handle<std::vector<std::vector<int> > > matchHandle;
00183     evt.getByLabel(match_, matchHandle);
00184     match = *(matchHandle->begin());
00185     // check if match is valid
00186     if(match.size()!=nPartons) invalidMatch=true;
00187     else {
00188       for(unsigned int idx=0; idx<match.size(); ++idx) {
00189         if(match[idx]<0 || match[idx]>=(int)jets->size()) {
00190           invalidMatch=true;
00191           break;
00192         }
00193       }
00194     }
00195   }
00196 
00197   // -----------------------------------------------------
00198   // skip events with no appropriate lepton candidate in
00199   // or empty MET or less jets than partons or invalid match
00200   // -----------------------------------------------------
00201 
00202   if( leps->empty() || mets->empty() || jets->size()<nPartons || invalidMatch ) {
00203     // the kinFit getters return empty objects here
00204     pPartonsHadP->push_back( fitter->fittedHadP()     );
00205     pPartonsHadQ->push_back( fitter->fittedHadQ()     );
00206     pPartonsHadB->push_back( fitter->fittedHadB()     );
00207     pPartonsLepB->push_back( fitter->fittedLepB()     );
00208     pLeptons    ->push_back( fitter->fittedLepton()   );
00209     pNeutrinos  ->push_back( fitter->fittedNeutrino() );
00210     // indices referring to the jet combination
00211     std::vector<int> invalidCombi;
00212     for(unsigned int i = 0; i < nPartons; ++i) 
00213       invalidCombi.push_back( -1 );
00214     pCombi->push_back( invalidCombi );
00215     // chi2
00216     pChi2->push_back( -1. );
00217     // chi2 probability
00218     pProb->push_back( -1. );
00219     // status of the fitter
00220     pStatus->push_back( -1 );
00221     // feed out all products
00222     evt.put(pCombi);
00223     evt.put(pPartonsHadP, "PartonsHadP");
00224     evt.put(pPartonsHadQ, "PartonsHadQ");
00225     evt.put(pPartonsHadB, "PartonsHadB");
00226     evt.put(pPartonsLepB, "PartonsLepB");
00227     evt.put(pLeptons    , "Leptons"    );
00228     evt.put(pNeutrinos  , "Neutrinos"  );
00229     evt.put(pChi2       , "Chi2"       );
00230     evt.put(pProb       , "Prob"       );
00231     evt.put(pStatus     , "Status"     );
00232     return;
00233   }
00234 
00235   // -----------------------------------------------------
00236   // analyze different jet combinations using the KinFitter
00237   // (or only a given jet combination if useOnlyMatch=true)
00238   // -----------------------------------------------------
00239   
00240   std::vector<int> jetIndices;
00241   if(!useOnlyMatch_) {
00242     for(unsigned int i=0; i<jets->size(); ++i){
00243       if(maxNJets_ >= (int) nPartons && maxNJets_ == (int) i) break;
00244       jetIndices.push_back(i);
00245     }
00246   }
00247   
00248   std::vector<int> combi;
00249   for(unsigned int i=0; i<nPartons; ++i) {
00250     if(useOnlyMatch_) combi.push_back( match[i] );
00251     else combi.push_back(i);
00252   }
00253 
00254   std::list<KinFitResult> FitResultList;
00255 
00256   do{
00257     for(int cnt = 0; cnt < TMath::Factorial( combi.size() ); ++cnt){
00258       // take into account indistinguishability of the two jets from the hadr. W decay,
00259       // reduces combinatorics by a factor of 2
00260       if( (combi[TtSemiLepEvtPartons::LightQ] < combi[TtSemiLepEvtPartons::LightQBar]
00261          || useOnlyMatch_ ) && doBTagging(useBTag_, jets, combi, bTagAlgo_, minBTagValueBJet_, maxBTagValueNonBJet_) ){
00262         
00263         std::vector<pat::Jet> jetCombi;
00264         jetCombi.resize(nPartons);
00265         jetCombi[TtSemiLepEvtPartons::LightQ   ] = (*jets)[combi[TtSemiLepEvtPartons::LightQ   ]];
00266         jetCombi[TtSemiLepEvtPartons::LightQBar] = (*jets)[combi[TtSemiLepEvtPartons::LightQBar]];
00267         jetCombi[TtSemiLepEvtPartons::HadB     ] = (*jets)[combi[TtSemiLepEvtPartons::HadB     ]];
00268         jetCombi[TtSemiLepEvtPartons::LepB     ] = (*jets)[combi[TtSemiLepEvtPartons::LepB     ]];
00269 
00270         // do the kinematic fit
00271         int status = fitter->fit(jetCombi, (*leps)[0], (*mets)[0], jetEnergyResolutionSmearFactor_);
00272 
00273         if( status == 0 ) { // only take into account converged fits
00274           KinFitResult result;
00275           result.Status = status;
00276           result.Chi2 = fitter->fitS();
00277           result.Prob = fitter->fitProb();
00278           result.HadB = fitter->fittedHadB();
00279           result.HadP = fitter->fittedHadP();
00280           result.HadQ = fitter->fittedHadQ();
00281           result.LepB = fitter->fittedLepB();
00282           result.LepL = fitter->fittedLepton();
00283           result.LepN = fitter->fittedNeutrino();
00284           result.JetCombi = combi;
00285 
00286           FitResultList.push_back(result);
00287         }
00288 
00289       }
00290       if(useOnlyMatch_) break; // don't go through combinatorics if useOnlyMatch was chosen
00291       next_permutation( combi.begin(), combi.end() );
00292     }
00293     if(useOnlyMatch_) break; // don't go through combinatorics if useOnlyMatch was chosen
00294   }
00295   while(stdcomb::next_combination( jetIndices.begin(), jetIndices.end(), combi.begin(), combi.end() ));
00296 
00297   // sort results w.r.t. chi2 values
00298   FitResultList.sort();
00299   
00300   // -----------------------------------------------------
00301   // feed out result
00302   // starting with the JetComb having the smallest chi2
00303   // -----------------------------------------------------
00304 
00305   if( (unsigned)FitResultList.size() < 1 ) { // in case no fit results were stored in the list (all fits aborted)
00306     pPartonsHadP->push_back( fitter->fittedHadP()     );
00307     pPartonsHadQ->push_back( fitter->fittedHadQ()     );
00308     pPartonsHadB->push_back( fitter->fittedHadB()     );
00309     pPartonsLepB->push_back( fitter->fittedLepB()     );
00310     pLeptons    ->push_back( fitter->fittedLepton()   );
00311     pNeutrinos  ->push_back( fitter->fittedNeutrino() );
00312     // indices referring to the jet combination
00313     std::vector<int> invalidCombi;
00314     for(unsigned int i = 0; i < nPartons; ++i) 
00315       invalidCombi.push_back( -1 );
00316     pCombi->push_back( invalidCombi );
00317     // chi2
00318     pChi2->push_back( -1. );
00319     // chi2 probability
00320     pProb->push_back( -1. );
00321     // status of the fitter
00322     pStatus->push_back( -1 );
00323   }
00324   else {
00325     unsigned int iComb = 0;
00326     for(typename std::list<KinFitResult>::const_iterator result = FitResultList.begin(); result != FitResultList.end(); ++result) {
00327       if(maxNComb_ >= 1 && iComb == (unsigned int) maxNComb_) break;
00328       iComb++;
00329       // partons
00330       pPartonsHadP->push_back( result->HadP );
00331       pPartonsHadQ->push_back( result->HadQ );
00332       pPartonsHadB->push_back( result->HadB );
00333       pPartonsLepB->push_back( result->LepB );
00334       // lepton
00335       pLeptons->push_back( result->LepL );
00336       // neutrino
00337       pNeutrinos->push_back( result->LepN );
00338       // indices referring to the jet combination
00339       pCombi->push_back( result->JetCombi );
00340       // chi2
00341       pChi2->push_back( result->Chi2 );
00342       // chi2 probability
00343       pProb->push_back( result->Prob );
00344       // status of the fitter
00345       pStatus->push_back( result->Status );
00346     }
00347   }
00348   evt.put(pCombi);
00349   evt.put(pPartonsHadP, "PartonsHadP");
00350   evt.put(pPartonsHadQ, "PartonsHadQ");
00351   evt.put(pPartonsHadB, "PartonsHadB");
00352   evt.put(pPartonsLepB, "PartonsLepB");
00353   evt.put(pLeptons    , "Leptons"    );
00354   evt.put(pNeutrinos  , "Neutrinos"  );
00355   evt.put(pChi2       , "Chi2"       );
00356   evt.put(pProb       , "Prob"       );
00357   evt.put(pStatus     , "Status"     );
00358 }
00359  
00360 template<typename LeptonCollection>
00361 TtSemiLepKinFitter::Param TtSemiLepKinFitProducer<LeptonCollection>::param(unsigned val) 
00362 {
00363   TtSemiLepKinFitter::Param result;
00364   switch(val){
00365   case TtSemiLepKinFitter::kEMom       : result=TtSemiLepKinFitter::kEMom;       break;
00366   case TtSemiLepKinFitter::kEtEtaPhi   : result=TtSemiLepKinFitter::kEtEtaPhi;   break;
00367   case TtSemiLepKinFitter::kEtThetaPhi : result=TtSemiLepKinFitter::kEtThetaPhi; break;
00368   default: 
00369     throw cms::Exception("WrongConfig") 
00370       << "Chosen jet parametrization is not supported: " << val << "\n";
00371     break;
00372   }
00373   return result;
00374 } 
00375 
00376 template<typename LeptonCollection>
00377 TtSemiLepKinFitter::Constraint TtSemiLepKinFitProducer<LeptonCollection>::constraint(unsigned val) 
00378 {
00379   TtSemiLepKinFitter::Constraint result;
00380   switch(val){
00381   case TtSemiLepKinFitter::kWHadMass       : result=TtSemiLepKinFitter::kWHadMass;       break;
00382   case TtSemiLepKinFitter::kWLepMass       : result=TtSemiLepKinFitter::kWLepMass;       break;
00383   case TtSemiLepKinFitter::kTopHadMass     : result=TtSemiLepKinFitter::kTopHadMass;     break;
00384   case TtSemiLepKinFitter::kTopLepMass     : result=TtSemiLepKinFitter::kTopLepMass;     break;
00385   case TtSemiLepKinFitter::kNeutrinoMass   : result=TtSemiLepKinFitter::kNeutrinoMass;   break;
00386   case TtSemiLepKinFitter::kEqualTopMasses : result=TtSemiLepKinFitter::kEqualTopMasses; break;
00387   default: 
00388     throw cms::Exception("WrongConfig") 
00389       << "Chosen fit constraint is not supported: " << val << "\n";
00390     break;
00391   }
00392   return result;
00393 } 
00394 
00395 template<typename LeptonCollection>
00396 std::vector<TtSemiLepKinFitter::Constraint> TtSemiLepKinFitProducer<LeptonCollection>::constraints(std::vector<unsigned>& val)
00397 {
00398   std::vector<TtSemiLepKinFitter::Constraint> result;
00399   for(unsigned i=0; i<val.size(); ++i){
00400     result.push_back(constraint(val[i]));
00401   }
00402   return result; 
00403 }
00404 
00405 #endif
/data/refman/pasoursint/CMSSW_5_3_0/src/TopQuarkAnalysis/TopKinFitter/plugins/TtSemiLepKinFitProducer.h
