CMS 3D CMS Logo
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 // system include files
48 #include <memory>
49 #include <string>
50 #include <iostream>
51 #include <vector>
52 #include <Math/VectorUtil.h>
53 #include <TMath.h>
54 
55 class calcTopMass : public edm::EDAnalyzer {
56 public:
57  explicit calcTopMass(const edm::ParameterSet &);
58  ~calcTopMass() override{};
59  void analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) override;
60 
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  sourceByRefer_ = iConfig.getParameter<InputTag>("srcByReference");
88  m_qT_CorrectorName = iConfig.getParameter<std::string>("qTopCorrector");
89  m_cT_CorrectorName = iConfig.getParameter<std::string>("cTopCorrector");
90  m_bT_CorrectorName = iConfig.getParameter<std::string>("bTopCorrector");
91  m_tT_CorrectorName = iConfig.getParameter<std::string>("tTopCorrector");
92 
93  bMass = 4.5;
94  cMass = 1.5;
95  qMass = 0.3;
96 
97  Service<TFileService> fs;
98  hMassNoCorr = fs->make<TH1F>("hMassNoCorr", "", 100, 100, 300);
99  hMassCorFl0 = fs->make<TH1F>("hMassCorFl0", "", 100, 100, 300);
100  hMassCorFlM = fs->make<TH1F>("hMassCorFlM", "", 100, 100, 300);
101  hMassCorMix = fs->make<TH1F>("hMassCorMix", "", 100, 100, 300);
102 }
103 
104 void calcTopMass::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
105  cout << "[calcTopMass] analysing event " << iEvent.id() << endl;
106  try {
107  iEvent.getByLabel(sourceByRefer_, theJetPartonMatch);
108  } catch (std::exception &ce) {
109  cerr << "[calcTopMass] caught std::exception " << ce.what() << endl;
110  return;
111  }
112 
113  // get all correctors from top events
114  const JetCorrector *qTopCorrector = JetCorrector::getJetCorrector(m_qT_CorrectorName, iSetup);
115  const JetCorrector *cTopCorrector = JetCorrector::getJetCorrector(m_cT_CorrectorName, iSetup);
116  const JetCorrector *bTopCorrector = JetCorrector::getJetCorrector(m_bT_CorrectorName, iSetup);
117  const JetCorrector *tTopCorrector = JetCorrector::getJetCorrector(m_tT_CorrectorName, iSetup);
118 
119  TLorentzVector jetsNoCorr[6];
120  TLorentzVector jetsCorFl0[6];
121  TLorentzVector jetsCorFlM[6];
122  TLorentzVector jetsCorMix[6];
123 
124  bool isQuarkFound[6] = {false};
125 
126  for (JetMatchedPartonsCollection::const_iterator j = theJetPartonMatch->begin(); j != theJetPartonMatch->end(); j++) {
127  const math::XYZTLorentzVector theJet = (*j).first.get()->p4();
128  const MatchedPartons aMatch = (*j).second;
129  const GenParticleRef &thePhyDef = aMatch.physicsDefinitionParton();
130 
131  if (thePhyDef.isNonnull()) {
132  int particPDG = thePhyDef.get()->pdgId();
133 
134  if (particPDG == 5) { //b from t
135  double bTcorr = bTopCorrector->correction(theJet);
136  double tTcorr = tTopCorrector->correction(theJet);
137  jetsNoCorr[0].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
138  jetsCorFl0[0].SetPtEtaPhiM(theJet.pt() * bTcorr, theJet.eta(), theJet.phi(), 0);
139  jetsCorFlM[0].SetPtEtaPhiM(theJet.pt() * bTcorr, theJet.eta(), theJet.phi(), bMass);
140  jetsCorMix[0].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
141  isQuarkFound[0] = true;
142  } else if (particPDG == -5) { //bbar from tbar
143  double bTcorr = bTopCorrector->correction(theJet);
144  double tTcorr = tTopCorrector->correction(theJet);
145  jetsNoCorr[3].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
146  jetsCorFl0[3].SetPtEtaPhiM(theJet.pt() * bTcorr, theJet.eta(), theJet.phi(), 0);
147  jetsCorFlM[3].SetPtEtaPhiM(theJet.pt() * bTcorr, theJet.eta(), theJet.phi(), bMass);
148  jetsCorMix[3].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
149  isQuarkFound[3] = true;
150  } else if (particPDG == 2) { //W+ from t
151  double qTcorr = qTopCorrector->correction(theJet);
152  double tTcorr = tTopCorrector->correction(theJet);
153  jetsNoCorr[1].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
154  jetsCorFl0[1].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
155  jetsCorFlM[1].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
156  jetsCorMix[1].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
157  isQuarkFound[1] = true;
158  } else if (particPDG == 4) { //W+ from t
159  double cTcorr = cTopCorrector->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() * cTcorr, theJet.eta(), theJet.phi(), 0);
163  jetsCorFlM[1].SetPtEtaPhiM(theJet.pt() * cTcorr, theJet.eta(), theJet.phi(), cMass);
164  jetsCorMix[1].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
165  isQuarkFound[1] = true;
166  } else if (particPDG == -1) { //W+ from t
167  double qTcorr = qTopCorrector->correction(theJet);
168  double tTcorr = tTopCorrector->correction(theJet);
169  jetsNoCorr[2].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
170  jetsCorFl0[2].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
171  jetsCorFlM[2].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
172  jetsCorMix[2].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
173  isQuarkFound[2] = true;
174  } else if (particPDG == -3) { //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 == -2) { //W- from tbar
183  double qTcorr = qTopCorrector->correction(theJet);
184  double tTcorr = tTopCorrector->correction(theJet);
185  jetsNoCorr[4].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
186  jetsCorFl0[4].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
187  jetsCorFlM[4].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
188  jetsCorMix[4].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
189  isQuarkFound[4] = true;
190  } else if (particPDG == -4) { //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 == 1) { //W- from tbar
199  double qTcorr = qTopCorrector->correction(theJet);
200  double tTcorr = tTopCorrector->correction(theJet);
201  jetsNoCorr[5].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
202  jetsCorFl0[5].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
203  jetsCorFlM[5].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
204  jetsCorMix[5].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
205  isQuarkFound[5] = true;
206  } else if (particPDG == 3) { //W- from tbar
207  double cTcorr = cTopCorrector->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() * cTcorr, theJet.eta(), theJet.phi(), 0);
211  jetsCorFlM[5].SetPtEtaPhiM(theJet.pt() * cTcorr, theJet.eta(), theJet.phi(), cMass);
212  jetsCorMix[5].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
213  isQuarkFound[5] = true;
214  }
215  }
216  }
217 
218  if (isQuarkFound[0] && isQuarkFound[1] && isQuarkFound[2]) {
219  TLorentzVector *theTopPNoCorr = new TLorentzVector(jetsNoCorr[0] + jetsNoCorr[1] + jetsNoCorr[2]);
220  TLorentzVector *theTopPCorFl0 = new TLorentzVector(jetsCorFl0[0] + jetsCorFl0[1] + jetsCorFl0[2]);
221  TLorentzVector *theTopPCorFlM = new TLorentzVector(jetsCorFlM[0] + jetsCorFlM[1] + jetsCorFlM[2]);
222  TLorentzVector *theTopPCorMix = new TLorentzVector(jetsCorMix[0] + jetsCorMix[1] + jetsCorMix[2]);
223  hMassNoCorr->Fill(theTopPNoCorr->M());
224  hMassCorFl0->Fill(theTopPCorFl0->M());
225  hMassCorFlM->Fill(theTopPCorFlM->M());
226  hMassCorMix->Fill(theTopPCorMix->M());
227  }
228 
229  if (isQuarkFound[3] && isQuarkFound[4] && isQuarkFound[5]) {
230  TLorentzVector *theTopMNoCorr = new TLorentzVector(jetsNoCorr[3] + jetsNoCorr[4] + jetsNoCorr[5]);
231  TLorentzVector *theTopMCorFl0 = new TLorentzVector(jetsCorFl0[3] + jetsCorFl0[4] + jetsCorFl0[5]);
232  TLorentzVector *theTopMCorFlM = new TLorentzVector(jetsCorFlM[3] + jetsCorFlM[4] + jetsCorFlM[5]);
233  TLorentzVector *theTopMCorMix = new TLorentzVector(jetsCorMix[3] + jetsCorMix[4] + jetsCorMix[5]);
234  hMassNoCorr->Fill(theTopMNoCorr->M());
235  hMassCorFl0->Fill(theTopMCorFl0->M());
236  hMassCorFlM->Fill(theTopMCorFlM->M());
237  hMassCorMix->Fill(theTopMCorMix->M());
238  }
239 }
240 
241 DEFINE_FWK_MODULE(calcTopMass);
reco::MatchedPartons::physicsDefinitionParton
const GenParticleRef & physicsDefinitionParton() const
Definition: MatchedPartons.h:28
JetMatchedPartons.h
ESHandle.h
JetCorrector.h
calcTopMass::hMassCorFlM
TH1F * hMassCorFlM
Definition: CalcTopMassExample.cc:77
edm
HLT enums.
Definition: AlignableModifier.h:19
JetFloatAssociation.h
calcTopMass::m_cT_CorrectorName
std::string m_cT_CorrectorName
Definition: CalcTopMassExample.cc:67
gather_cfg.cout
cout
Definition: gather_cfg.py:144
reco::JetCorrector
Definition: JetCorrector.h:35
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:232
getRef.h
Jet.h
reco::JetCorrector::correction
double correction(const LorentzVector &fJet) const
get correction using Jet information only
Definition: JetCorrector.h:47
EDAnalyzer.h
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
edm::Handle< reco::JetMatchedPartonsCollection >
edm::Ref< GenParticleCollection >
edm::EDAnalyzer
Definition: EDAnalyzer.h:29
CandidateFwd.h
calcTopMass::hMassNoCorr
TH1F * hMassNoCorr
Definition: CalcTopMassExample.cc:75
MakerMacros.h
calcTopMass::sourceByRefer_
edm::InputTag sourceByRefer_
Definition: CalcTopMassExample.cc:63
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
MatchedPartons.h
Service.h
calcTopMass::calcTopMass
calcTopMass(const edm::ParameterSet &)
Definition: CalcTopMassExample.cc:86
calcTopMass::m_qT_CorrectorName
std::string m_qT_CorrectorName
Definition: CalcTopMassExample.cc:66
JetCollection.h
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
TFileService.h
cppFunctionSkipper.exception
exception
Definition: cppFunctionSkipper.py:10
calcTopMass::hMassCorFl0
TH1F * hMassCorFl0
Definition: CalcTopMassExample.cc:76
edm::ParameterSet
Definition: ParameterSet.h:36
Event.h
calcTopMass::analyze
void analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) override
Definition: CalcTopMassExample.cc:104
beam_dqm_sourceclient-live_cfg.cerr
cerr
Definition: beam_dqm_sourceclient-live_cfg.py:17
calcTopMass::m_tT_CorrectorName
std::string m_tT_CorrectorName
Definition: CalcTopMassExample.cc:69
JetFlavour.h
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
edm::Service< TFileService >
iEvent
int iEvent
Definition: GenABIO.cc:224
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:48
edm::EventSetup
Definition: EventSetup.h:57
calcTopMass::~calcTopMass
~calcTopMass() override
Definition: CalcTopMassExample.cc:58
calcTopMass::qMass
float qMass
Definition: CalcTopMassExample.cc:73
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
calcTopMass::hMassCorMix
TH1F * hMassCorMix
Definition: CalcTopMassExample.cc:78
std
Definition: JetResolutionObject.h:76
Ref.h
calcTopMass::bMass
float bMass
Definition: CalcTopMassExample.cc:71
calcTopMass
Definition: CalcTopMassExample.cc:55
Frameworkfwd.h
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
EventSetup.h
calcTopMass::m_bT_CorrectorName
std::string m_bT_CorrectorName
Definition: CalcTopMassExample.cc:68
Candidate.h
View.h
ParameterSet.h
calcTopMass::cMass
float cMass
Definition: CalcTopMassExample.cc:72
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
edm::Event
Definition: Event.h:73
JetFlavourMatching.h
edm::InputTag
Definition: InputTag.h:15
calcTopMass::theJetPartonMatch
edm::Handle< reco::JetMatchedPartonsCollection > theJetPartonMatch
Definition: CalcTopMassExample.cc:64
TFileService::make
T * make(const Args &... args) const
make new ROOT object
Definition: TFileService.h:64
reco::MatchedPartons
Definition: MatchedPartons.h:12
calcTopMass::sourcePartons_
edm::InputTag sourcePartons_
Definition: CalcTopMassExample.cc:62
Generated for CMSSW Reference Manual by doxygen 1.8.16