CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CalcTopMassExample.cc
Go to the documentation of this file.
1 //
2 // Example to calculated top mass from genjets corrected on the fly with parton calibration
3 // Top mass is calculated from:
4 // 1) Uncorrected jets
5 // 2) Corrected jets with flavour ttbar correction assuming jet mass = quark mass
6 // 3) The same as 3 assuming massless jets
7 // 4) Corrected jets with mixed ttbar correction assuming massless jets
8 //
9 // Author: Attilio Santocchia
10 // Date: 15.06.2009
11 // Tested on CMSSW_3_1_0
12 //
13 // user include files
14 #include "FWCore/Framework/interface/Frameworkfwd.h"
15 #include "FWCore/Framework/interface/EDAnalyzer.h"
16 #include "FWCore/Framework/interface/ESHandle.h"
17 #include "FWCore/Framework/interface/Event.h"
18 #include "FWCore/Framework/interface/EventSetup.h"
19 #include "FWCore/Framework/interface/MakerMacros.h"
20 #include "FWCore/ParameterSet/interface/ParameterSet.h"
21 #include "FWCore/ServiceRegistry/interface/Service.h"
22 
23 #include "DataFormats/Candidate/interface/Candidate.h"
24 #include "DataFormats/Candidate/interface/CandidateFwd.h"
25 #include "DataFormats/JetReco/interface/Jet.h"
26 #include "DataFormats/JetReco/interface/JetCollection.h"
27 #include "DataFormats/JetReco/interface/JetFloatAssociation.h"
28 
29 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
30 
31 #include "DataFormats/Common/interface/View.h"
32 #include "DataFormats/Common/interface/Ref.h"
33 #include "DataFormats/Common/interface/getRef.h"
34 
35 #include "SimDataFormats/JetMatching/interface/JetFlavour.h"
36 #include "SimDataFormats/JetMatching/interface/JetFlavourMatching.h"
37 #include "SimDataFormats/JetMatching/interface/MatchedPartons.h"
38 #include "SimDataFormats/JetMatching/interface/JetMatchedPartons.h"
39 
40 #include "CommonTools/UtilAlgos/interface/TFileService.h"
41 
42 #include "TFile.h"
43 #include "TH1.h"
44 #include "TF1.h"
45 #include "TLorentzVector.h"
46 
47 
48 // system include files
49 #include <memory>
50 #include <string>
51 #include <iostream>
52 #include <vector>
53 #include <Math/VectorUtil.h>
54 #include <TMath.h>
55 
56 class calcTopMass : public edm::EDAnalyzer {
57  public:
58  explicit calcTopMass(const edm::ParameterSet & );
59  ~calcTopMass() {};
60  virtual void analyze( const edm::Event& iEvent, const edm::EventSetup& iSetup);
61  private:
62  edm::InputTag sourcePartons_;
63  edm::InputTag sourceByRefer_;
64  edm::Handle<reco::JetMatchedPartonsCollection> theJetPartonMatch;
65 
66  std::string m_qT_CorrectorName;
67  std::string m_cT_CorrectorName;
68  std::string m_bT_CorrectorName;
69  std::string m_tT_CorrectorName;
70 
71  float bMass;
72  float cMass;
73  float qMass;
74 
75  TH1F *hMassNoCorr;
76  TH1F *hMassCorFl0;
77  TH1F *hMassCorFlM;
78  TH1F *hMassCorMix;
79 };
80 
81 using namespace std;
82 using namespace reco;
83 using namespace edm;
84 using namespace ROOT::Math::VectorUtil;
85 
86 calcTopMass::calcTopMass(const edm::ParameterSet& iConfig)
87 {
88 
89  sourceByRefer_ = iConfig.getParameter<InputTag> ("srcByReference");
90  m_qT_CorrectorName = iConfig.getParameter <std::string> ("qTopCorrector");
91  m_cT_CorrectorName = iConfig.getParameter <std::string> ("cTopCorrector");
92  m_bT_CorrectorName = iConfig.getParameter <std::string> ("bTopCorrector");
93  m_tT_CorrectorName = iConfig.getParameter <std::string> ("tTopCorrector");
94 
95  bMass = 4.5;
96  cMass = 1.5;
97  qMass = 0.3;
98 
99  Service<TFileService> fs;
100  hMassNoCorr = fs->make<TH1F>("hMassNoCorr","",100,100,300);
101  hMassCorFl0 = fs->make<TH1F>("hMassCorFl0","",100,100,300);
102  hMassCorFlM = fs->make<TH1F>("hMassCorFlM","",100,100,300);
103  hMassCorMix = fs->make<TH1F>("hMassCorMix","",100,100,300);
104 
105 }
106 
107 void calcTopMass::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
108 {
109  cout << "[calcTopMass] analysing event " << iEvent.id() << endl;
110  try {
111  iEvent.getByLabel (sourceByRefer_ , theJetPartonMatch );
112  } catch(std::exception& ce) {
113  cerr << "[calcTopMass] caught std::exception " << ce.what() << endl;
114  return;
115  }
116 
117  // get all correctors from top events
118  const JetCorrector* qTopCorrector = JetCorrector::getJetCorrector (m_qT_CorrectorName, iSetup);
119  const JetCorrector* cTopCorrector = JetCorrector::getJetCorrector (m_cT_CorrectorName, iSetup);
120  const JetCorrector* bTopCorrector = JetCorrector::getJetCorrector (m_bT_CorrectorName, iSetup);
121  const JetCorrector* tTopCorrector = JetCorrector::getJetCorrector (m_tT_CorrectorName, iSetup);
122 
123  TLorentzVector jetsNoCorr[6];
124  TLorentzVector jetsCorFl0[6];
125  TLorentzVector jetsCorFlM[6];
126  TLorentzVector jetsCorMix[6];
127 
128  bool isQuarkFound[6] = {false};
129 
130  for ( JetMatchedPartonsCollection::const_iterator j = theJetPartonMatch->begin();
131  j != theJetPartonMatch->end();
132  j ++ ) {
133 
134  const math::XYZTLorentzVector theJet = (*j).first.get()->p4();
135  const MatchedPartons aMatch = (*j).second;
136  const GenParticleRef thePhyDef = aMatch.physicsDefinitionParton() ;
137 
138  if(thePhyDef.isNonnull()) {
139 
140  int particPDG = thePhyDef.get()->pdgId();
141 
142  if( particPDG == 5 ) { //b from t
143  double bTcorr = bTopCorrector->correction (theJet);
144  double tTcorr = tTopCorrector->correction (theJet);
145  jetsNoCorr[0].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
146  jetsCorFl0[0].SetPtEtaPhiM( theJet.pt()*bTcorr , theJet.eta() , theJet.phi() , 0 );
147  jetsCorFlM[0].SetPtEtaPhiM( theJet.pt()*bTcorr , theJet.eta() , theJet.phi() , bMass );
148  jetsCorMix[0].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
149  isQuarkFound[0]=true;
150  } else if ( particPDG == -5 ) { //bbar from tbar
151  double bTcorr = bTopCorrector->correction (theJet);
152  double tTcorr = tTopCorrector->correction (theJet);
153  jetsNoCorr[3].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
154  jetsCorFl0[3].SetPtEtaPhiM( theJet.pt()*bTcorr , theJet.eta() , theJet.phi() , 0 );
155  jetsCorFlM[3].SetPtEtaPhiM( theJet.pt()*bTcorr , theJet.eta() , theJet.phi() , bMass );
156  jetsCorMix[3].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
157  isQuarkFound[3]=true;
158  } else if ( particPDG == 2 ) { //W+ from t
159  double qTcorr = qTopCorrector->correction (theJet);
160  double tTcorr = tTopCorrector->correction (theJet);
161  jetsNoCorr[1].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
162  jetsCorFl0[1].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0 );
163  jetsCorFlM[1].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass );
164  jetsCorMix[1].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
165  isQuarkFound[1]=true;
166  } else if ( particPDG == 4 ) { //W+ from t
167  double cTcorr = cTopCorrector->correction (theJet);
168  double tTcorr = tTopCorrector->correction (theJet);
169  jetsNoCorr[1].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
170  jetsCorFl0[1].SetPtEtaPhiM( theJet.pt()*cTcorr , theJet.eta() , theJet.phi() , 0 );
171  jetsCorFlM[1].SetPtEtaPhiM( theJet.pt()*cTcorr , theJet.eta() , theJet.phi() , cMass );
172  jetsCorMix[1].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
173  isQuarkFound[1]=true;
174  } else if ( particPDG == -1 ) { //W+ from t
175  double qTcorr = qTopCorrector->correction (theJet);
176  double tTcorr = tTopCorrector->correction (theJet);
177  jetsNoCorr[2].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
178  jetsCorFl0[2].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0 );
179  jetsCorFlM[2].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass );
180  jetsCorMix[2].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
181  isQuarkFound[2]=true;
182  } else if ( particPDG == -3 ) { //W+ from t
183  double qTcorr = qTopCorrector->correction (theJet);
184  double tTcorr = tTopCorrector->correction (theJet);
185  jetsNoCorr[2].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
186  jetsCorFl0[2].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0 );
187  jetsCorFlM[2].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass );
188  jetsCorMix[2].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
189  isQuarkFound[2]=true;
190  } else if ( particPDG == -2 ) { //W- from tbar
191  double qTcorr = qTopCorrector->correction (theJet);
192  double tTcorr = tTopCorrector->correction (theJet);
193  jetsNoCorr[4].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
194  jetsCorFl0[4].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0 );
195  jetsCorFlM[4].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass );
196  jetsCorMix[4].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
197  isQuarkFound[4]=true;
198  } else if ( particPDG == -4 ) { //W- from tbar
199  double qTcorr = qTopCorrector->correction (theJet);
200  double tTcorr = tTopCorrector->correction (theJet);
201  jetsNoCorr[4].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
202  jetsCorFl0[4].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0 );
203  jetsCorFlM[4].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass );
204  jetsCorMix[4].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
205  isQuarkFound[4]=true;
206  } else if ( particPDG == 1 ) { //W- from tbar
207  double qTcorr = qTopCorrector->correction (theJet);
208  double tTcorr = tTopCorrector->correction (theJet);
209  jetsNoCorr[5].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
210  jetsCorFl0[5].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , 0 );
211  jetsCorFlM[5].SetPtEtaPhiM( theJet.pt()*qTcorr , theJet.eta() , theJet.phi() , qMass );
212  jetsCorMix[5].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
213  isQuarkFound[5]=true;
214  } else if ( particPDG == 3 ) { //W- from tbar
215  double cTcorr = cTopCorrector->correction (theJet);
216  double tTcorr = tTopCorrector->correction (theJet);
217  jetsNoCorr[5].SetPtEtaPhiM( theJet.pt() , theJet.eta() , theJet.phi() , theJet.mass() );
218  jetsCorFl0[5].SetPtEtaPhiM( theJet.pt()*cTcorr , theJet.eta() , theJet.phi() , 0 );
219  jetsCorFlM[5].SetPtEtaPhiM( theJet.pt()*cTcorr , theJet.eta() , theJet.phi() , cMass );
220  jetsCorMix[5].SetPtEtaPhiM( theJet.pt()*tTcorr , theJet.eta() , theJet.phi() , 0 );
221  isQuarkFound[5]=true;
222  }
223  }
224  }
225 
226  if( isQuarkFound[0] && isQuarkFound[1] && isQuarkFound[2] ) {
227  TLorentzVector *theTopPNoCorr = new TLorentzVector( jetsNoCorr[0] + jetsNoCorr[1] + jetsNoCorr[2] );
228  TLorentzVector *theTopPCorFl0 = new TLorentzVector( jetsCorFl0[0] + jetsCorFl0[1] + jetsCorFl0[2] );
229  TLorentzVector *theTopPCorFlM = new TLorentzVector( jetsCorFlM[0] + jetsCorFlM[1] + jetsCorFlM[2] );
230  TLorentzVector *theTopPCorMix = new TLorentzVector( jetsCorMix[0] + jetsCorMix[1] + jetsCorMix[2] );
231  hMassNoCorr->Fill( theTopPNoCorr->M() );
232  hMassCorFl0->Fill( theTopPCorFl0->M() );
233  hMassCorFlM->Fill( theTopPCorFlM->M() );
234  hMassCorMix->Fill( theTopPCorMix->M() );
235  }
236 
237  if( isQuarkFound[3] && isQuarkFound[4] && isQuarkFound[5] ) {
238  TLorentzVector *theTopMNoCorr = new TLorentzVector( jetsNoCorr[3] + jetsNoCorr[4] + jetsNoCorr[5] );
239  TLorentzVector *theTopMCorFl0 = new TLorentzVector( jetsCorFl0[3] + jetsCorFl0[4] + jetsCorFl0[5] );
240  TLorentzVector *theTopMCorFlM = new TLorentzVector( jetsCorFlM[3] + jetsCorFlM[4] + jetsCorFlM[5] );
241  TLorentzVector *theTopMCorMix = new TLorentzVector( jetsCorMix[3] + jetsCorMix[4] + jetsCorMix[5] );
242  hMassNoCorr->Fill( theTopMNoCorr->M() );
243  hMassCorFl0->Fill( theTopMCorFl0->M() );
244  hMassCorFlM->Fill( theTopMCorFlM->M() );
245  hMassCorMix->Fill( theTopMCorMix->M() );
246  }
247 
248 }
249 
250 DEFINE_FWK_MODULE( calcTopMass );
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
calcTopMass::theJetPartonMatch
edm::Handle< reco::JetMatchedPartonsCollection > theJetPartonMatch
Definition: CalcTopMassExample.cc:64
calcTopMass::m_cT_CorrectorName
std::string m_cT_CorrectorName
Definition: CalcTopMassExample.cc:67
JetCorrector::correction
virtual double correction(const LorentzVector &fJet) const =0
get correction using Jet information only
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
calcTopMass::sourceByRefer_
edm::InputTag sourceByRefer_
Definition: CalcTopMassExample.cc:63
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:30
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
iEvent
int iEvent
Definition: GenABIO.cc:243
L1Trigger_dataformats.reco
dictionary reco
Definition: L1Trigger_dataformats.py:9
j
int j
Definition: DBlmapReader.cc:9
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
calcTopMass::calcTopMass
calcTopMass(const edm::ParameterSet &)
Definition: CalcTopMassExample.cc:86
calcTopMass::m_qT_CorrectorName
std::string m_qT_CorrectorName
Definition: CalcTopMassExample.cc:66
reco::MatchedPartons::physicsDefinitionParton
const GenParticleRef & physicsDefinitionParton() const
Definition: MatchedPartons.h:41
JetCorrector::getJetCorrector
static const JetCorrector * getJetCorrector(const std::string &fName, const edm::EventSetup &fSetup)
retrieve corrector from the event setup. troughs exception if something is missing ...
Definition: JetCorrector.cc:51
calcTopMass::m_tT_CorrectorName
std::string m_tT_CorrectorName
Definition: CalcTopMassExample.cc:69
dtNoiseDBValidation_cfg.cerr
tuple cerr
Definition: dtNoiseDBValidation_cfg.py:22
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
calcTopMass::sourcePartons_
edm::InputTag sourcePartons_
Definition: CalcTopMassExample.cc:62
TFileDirectory::make
T * make() const
make new ROOT object
Definition: TFileDirectory.h:55
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:242
calcTopMass::analyze
virtual void analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup)
Definition: CalcTopMassExample.cc:107
calcTopMass::m_bT_CorrectorName
std::string m_bT_CorrectorName
Definition: CalcTopMassExample.cc:68
Generated for CMSSW Reference Manual by doxygen 1.8.5