Reference Manual

Go to the documentation of this file.
00001 //
00002 // Example to calculated top mass from genjets corrected on the fly with parton calibration 
00003 // Top mass is calculated from: 
00004 // 1) Uncorrected jets
00005 // 2) Corrected jets with flavour ttbar correction assuming jet mass = quark mass
00006 // 3) The same as 3 assuming massless jets
00007 // 4) Corrected jets with mixed ttbar correction assuming massless jets
00008 // 
00009 // Author: Attilio Santocchia
00010 // Date: 15.06.2009
00011 // Tested on CMSSW_3_1_0
00012 //
00013 // user include files
00014 #include "FWCore/Framework/interface/Frameworkfwd.h"
00015 #include "FWCore/Framework/interface/EDAnalyzer.h"
00016 #include "FWCore/Framework/interface/ESHandle.h"
00017 #include "FWCore/Framework/interface/Event.h"
00018 #include "FWCore/Framework/interface/EventSetup.h"
00019 #include "FWCore/Framework/interface/MakerMacros.h"
00020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00021 #include "FWCore/ServiceRegistry/interface/Service.h"
00022 
00023 #include "DataFormats/Candidate/interface/Candidate.h"
00024 #include "DataFormats/Candidate/interface/CandidateFwd.h"
00025 #include "DataFormats/JetReco/interface/Jet.h"
00026 #include "DataFormats/JetReco/interface/JetCollection.h"
00027 #include "DataFormats/JetReco/interface/JetFloatAssociation.h"
00028 
00029 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
00030 
00031 #include "DataFormats/Common/interface/View.h"
00032 #include "DataFormats/Common/interface/Ref.h"
00033 #include "DataFormats/Common/interface/getRef.h"
00034 
00035 #include "SimDataFormats/JetMatching/interface/JetFlavour.h"
00036 #include "SimDataFormats/JetMatching/interface/JetFlavourMatching.h"
00037 #include "SimDataFormats/JetMatching/interface/MatchedPartons.h"
00038 #include "SimDataFormats/JetMatching/interface/JetMatchedPartons.h"
00039 
00040 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00041 
00042 #include "TFile.h"
00043 #include "TH1.h"
00044 #include "TF1.h"
00045 #include "TLorentzVector.h"
00046 
00047 
00048 // system include files
00049 #include <memory>
00050 #include <string>
00051 #include <iostream>
00052 #include <vector>
00053 #include <Math/VectorUtil.h>
00054 #include <TMath.h>
00055 
00056 class calcTopMass : public edm::EDAnalyzer {
00057   public:
00058     explicit calcTopMass(const edm::ParameterSet & );
00059     ~calcTopMass() {};
00060     virtual void analyze( const edm::Event& iEvent, const edm::EventSetup& iSetup);
00061   private:
00062     edm::InputTag sourcePartons_;
00063     edm::InputTag sourceByRefer_;
00064     edm::Handle<reco::JetMatchedPartonsCollection> theJetPartonMatch;
00065 
00066     std::string m_qT_CorrectorName;
00067     std::string m_cT_CorrectorName;
00068     std::string m_bT_CorrectorName;
00069     std::string m_tT_CorrectorName;
00070 
00071     float bMass;
00072     float cMass;
00073     float qMass;
00074     
00075     TH1F *hMassNoCorr;
00076     TH1F *hMassCorFl0;
00077     TH1F *hMassCorFlM;
00078     TH1F *hMassCorMix;
00079 };
00080 
00081 using namespace std;
00082 using namespace reco;
00083 using namespace edm;
00084 using namespace ROOT::Math::VectorUtil;
00085 
00086 calcTopMass::calcTopMass(const edm::ParameterSet& iConfig)
00087 {
00088 
00089   sourceByRefer_ = iConfig.getParameter<InputTag> ("srcByReference");
00090   m_qT_CorrectorName = iConfig.getParameter <std::string> ("qTopCorrector");
00091   m_cT_CorrectorName = iConfig.getParameter <std::string> ("cTopCorrector");
00092   m_bT_CorrectorName = iConfig.getParameter <std::string> ("bTopCorrector");
00093   m_tT_CorrectorName = iConfig.getParameter <std::string> ("tTopCorrector");
00094 
00095   bMass = 4.5;
00096   cMass = 1.5;
00097   qMass = 0.3;
00098 
00099   Service<TFileService> fs;
00100   hMassNoCorr = fs->make<TH1F>("hMassNoCorr","",100,100,300);
00101   hMassCorFl0 = fs->make<TH1F>("hMassCorFl0","",100,100,300);
00102   hMassCorFlM = fs->make<TH1F>("hMassCorFlM","",100,100,300);
00103   hMassCorMix = fs->make<TH1F>("hMassCorMix","",100,100,300);
00104 
00105 }
00106 
00107 void calcTopMass::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00108 {
00109   cout << "[calcTopMass] analysing event " << iEvent.id() << endl;
00110   try {
00111     iEvent.getByLabel (sourceByRefer_ , theJetPartonMatch   );
00112   } catch(std::exception& ce) {
00113     cerr << "[calcTopMass] caught std::exception " << ce.what() << endl;
00114     return;
00115   }
00116 
00117   // get all correctors from top events
00118   const JetCorrector* qTopCorrector = JetCorrector::getJetCorrector (m_qT_CorrectorName, iSetup);
00119   const JetCorrector* cTopCorrector = JetCorrector::getJetCorrector (m_cT_CorrectorName, iSetup);
00120   const JetCorrector* bTopCorrector = JetCorrector::getJetCorrector (m_bT_CorrectorName, iSetup);
00121   const JetCorrector* tTopCorrector = JetCorrector::getJetCorrector (m_tT_CorrectorName, iSetup);
00122 
00123   TLorentzVector jetsNoCorr[6];
00124   TLorentzVector jetsCorFl0[6];
00125   TLorentzVector jetsCorFlM[6];
00126   TLorentzVector jetsCorMix[6];
00127 
00128   bool isQuarkFound[6] = {false};
00129 
00130   for ( JetMatchedPartonsCollection::const_iterator j  = theJetPartonMatch->begin();
00131                                      j != theJetPartonMatch->end();
00132                                      j ++ ) {
00133 
00134     const math::XYZTLorentzVector theJet = (*j).first.get()->p4();
00135     const MatchedPartons aMatch = (*j).second;
00136     const GenParticleRef thePhyDef = aMatch.physicsDefinitionParton() ;
00137     
00138     if(thePhyDef.isNonnull()) {
00139 
00140       int particPDG = thePhyDef.get()->pdgId();
00141 
00142       if(         particPDG ==  5 ) { //b from t
00143         double bTcorr = bTopCorrector->correction (theJet);
00144         double tTcorr = tTopCorrector->correction (theJet);
00145         jetsNoCorr[0].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00146         jetsCorFl0[0].SetPtEtaPhiM( theJet.pt()*bTcorr , theJet.eta() , theJet.phi() , 0             );
00147         jetsCorFlM[0].SetPtEtaPhiM( theJet.pt()*bTcorr , theJet.eta() , theJet.phi() , bMass         );
00148         jetsCorMix[0].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00149         isQuarkFound[0]=true;
00150       } else if ( particPDG == -5 ) { //bbar from tbar
00151         double bTcorr = bTopCorrector->correction (theJet);
00152         double tTcorr = tTopCorrector->correction (theJet);
00153         jetsNoCorr[3].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00154         jetsCorFl0[3].SetPtEtaPhiM( theJet.pt()*bTcorr , theJet.eta() , theJet.phi() , 0             );
00155         jetsCorFlM[3].SetPtEtaPhiM( theJet.pt()*bTcorr , theJet.eta() , theJet.phi() , bMass         );
00156         jetsCorMix[3].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00157         isQuarkFound[3]=true;
00158       } else if ( particPDG == 2 ) { //W+ from t
00159         double qTcorr = qTopCorrector->correction (theJet);
00160         double tTcorr = tTopCorrector->correction (theJet);
00161         jetsNoCorr[1].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00162         jetsCorFl0[1].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0             );
00163         jetsCorFlM[1].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass         );
00164         jetsCorMix[1].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00165         isQuarkFound[1]=true;
00166       } else if ( particPDG == 4 ) { //W+ from t
00167         double cTcorr = cTopCorrector->correction (theJet);
00168         double tTcorr = tTopCorrector->correction (theJet);
00169         jetsNoCorr[1].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00170         jetsCorFl0[1].SetPtEtaPhiM( theJet.pt()*cTcorr , theJet.eta() , theJet.phi() , 0             );
00171         jetsCorFlM[1].SetPtEtaPhiM( theJet.pt()*cTcorr , theJet.eta() , theJet.phi() , cMass         );
00172         jetsCorMix[1].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00173         isQuarkFound[1]=true;
00174       } else if ( particPDG == -1 ) { //W+ from t
00175         double qTcorr = qTopCorrector->correction (theJet);
00176         double tTcorr = tTopCorrector->correction (theJet);
00177         jetsNoCorr[2].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00178         jetsCorFl0[2].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0             );
00179         jetsCorFlM[2].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass         );
00180         jetsCorMix[2].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00181         isQuarkFound[2]=true;
00182       } else if ( particPDG == -3 ) { //W+ from t
00183         double qTcorr = qTopCorrector->correction (theJet);
00184         double tTcorr = tTopCorrector->correction (theJet);
00185         jetsNoCorr[2].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00186         jetsCorFl0[2].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0             );
00187         jetsCorFlM[2].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass         );
00188         jetsCorMix[2].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00189         isQuarkFound[2]=true;
00190       } else if ( particPDG == -2 ) { //W- from tbar
00191         double qTcorr = qTopCorrector->correction (theJet);
00192         double tTcorr = tTopCorrector->correction (theJet);
00193         jetsNoCorr[4].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00194         jetsCorFl0[4].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0             );
00195         jetsCorFlM[4].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass         );
00196         jetsCorMix[4].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00197         isQuarkFound[4]=true;
00198       } else if ( particPDG == -4 ) { //W- from tbar
00199         double qTcorr = qTopCorrector->correction (theJet);
00200         double tTcorr = tTopCorrector->correction (theJet);
00201         jetsNoCorr[4].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00202         jetsCorFl0[4].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0             );
00203         jetsCorFlM[4].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass         );
00204         jetsCorMix[4].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00205         isQuarkFound[4]=true;
00206       } else if ( particPDG == 1 ) { //W- from tbar
00207         double qTcorr = qTopCorrector->correction (theJet);
00208         double tTcorr = tTopCorrector->correction (theJet);
00209         jetsNoCorr[5].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00210         jetsCorFl0[5].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0             );
00211         jetsCorFlM[5].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass         );
00212         jetsCorMix[5].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00213         isQuarkFound[5]=true;
00214       } else if ( particPDG == 3 ) { //W- from tbar
00215         double cTcorr = cTopCorrector->correction (theJet);
00216         double tTcorr = tTopCorrector->correction (theJet);
00217         jetsNoCorr[5].SetPtEtaPhiM( theJet.pt()        , theJet.eta() , theJet.phi() , theJet.mass() );
00218         jetsCorFl0[5].SetPtEtaPhiM( theJet.pt()*cTcorr , theJet.eta() , theJet.phi() , 0             );
00219         jetsCorFlM[5].SetPtEtaPhiM( theJet.pt()*cTcorr , theJet.eta() , theJet.phi() , cMass         );
00220         jetsCorMix[5].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0             );
00221         isQuarkFound[5]=true;
00222       }
00223     }
00224   }
00225 
00226   if( isQuarkFound[0] && isQuarkFound[1] && isQuarkFound[2] ) {
00227     TLorentzVector *theTopPNoCorr = new TLorentzVector( jetsNoCorr[0] + jetsNoCorr[1] + jetsNoCorr[2] );
00228     TLorentzVector *theTopPCorFl0 = new TLorentzVector( jetsCorFl0[0] + jetsCorFl0[1] + jetsCorFl0[2] );
00229     TLorentzVector *theTopPCorFlM = new TLorentzVector( jetsCorFlM[0] + jetsCorFlM[1] + jetsCorFlM[2] );
00230     TLorentzVector *theTopPCorMix = new TLorentzVector( jetsCorMix[0] + jetsCorMix[1] + jetsCorMix[2] );
00231     hMassNoCorr->Fill( theTopPNoCorr->M() );
00232     hMassCorFl0->Fill( theTopPCorFl0->M() );
00233     hMassCorFlM->Fill( theTopPCorFlM->M() );
00234     hMassCorMix->Fill( theTopPCorMix->M() );
00235   }
00236 
00237   if( isQuarkFound[3] && isQuarkFound[4] && isQuarkFound[5] ) {
00238     TLorentzVector *theTopMNoCorr = new TLorentzVector( jetsNoCorr[3] + jetsNoCorr[4] + jetsNoCorr[5] );
00239     TLorentzVector *theTopMCorFl0 = new TLorentzVector( jetsCorFl0[3] + jetsCorFl0[4] + jetsCorFl0[5] );
00240     TLorentzVector *theTopMCorFlM = new TLorentzVector( jetsCorFlM[3] + jetsCorFlM[4] + jetsCorFlM[5] );
00241     TLorentzVector *theTopMCorMix = new TLorentzVector( jetsCorMix[3] + jetsCorMix[4] + jetsCorMix[5] );
00242     hMassNoCorr->Fill( theTopMNoCorr->M() );
00243     hMassCorFl0->Fill( theTopMCorFl0->M() );
00244     hMassCorFlM->Fill( theTopMCorFlM->M() );
00245     hMassCorMix->Fill( theTopMCorMix->M() );
00246   }
00247 
00248 }
00249 
00250 DEFINE_FWK_MODULE( calcTopMass );
/data/refman/pasoursint/CMSSW_5_3_3/src/RecoJets/JetAnalyzers/src/CalcTopMassExample.cc
