CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
JetAnaPythia.cc
Go to the documentation of this file.
1 // Name: JetAnaPythia
2 // Description: Example of analysis of Pythia produced partons & jets
3 // Based on Kostas Kousouris' templated JetPlotsExample.
4 // Plots are tailored to needs of dijet mass and ratio analysis.
5 // Author: R. Harris
6 // Date: 28 - Oct - 2008
7 #include "RecoJets/JetAnalyzers/interface/JetAnaPythia.h"
8 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
9 #include "DataFormats/JetReco/interface/PFJetCollection.h"
10 #include "DataFormats/JetReco/interface/GenJetCollection.h"
11 #include "DataFormats/JetReco/interface/CaloJet.h"
12 #include "DataFormats/JetReco/interface/PFJet.h"
13 #include "DataFormats/JetReco/interface/GenJet.h"
14 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
15 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
16 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
17 #include "FWCore/Framework/interface/Event.h"
18 #include "FWCore/Framework/interface/Run.h"
19 #include "FWCore/ParameterSet/interface/ParameterSet.h"
20 #include <TFile.h>
21 #include <cmath>
22 using namespace edm;
23 using namespace reco;
24 using namespace std;
26 template <class Jet>
27 JetAnaPythia<Jet>::JetAnaPythia(edm::ParameterSet const& cfg) {
28  JetAlgorithm = cfg.getParameter<std::string>("JetAlgorithm");
29  HistoFileName = cfg.getParameter<std::string>("HistoFileName");
30  NJets = cfg.getParameter<int>("NJets");
31  debug = cfg.getParameter<bool>("debug");
32  eventsGen = cfg.getParameter<int>("eventsGen");
33  anaLevel = cfg.getParameter<std::string>("anaLevel");
34  xsecGen = cfg.getParameter<vector<double> >("xsecGen");
35  ptHatEdges = cfg.getParameter<vector<double> >("ptHatEdges");
36 }
38 template <class Jet>
39 void JetAnaPythia<Jet>::beginJob() {
40  TString hname;
41  m_file = new TFile(HistoFileName.c_str(), "RECREATE");
43  const int nMassBins = 103;
44  double massBoundaries[nMassBins + 1] = {
45  1, 3, 6, 10, 16, 23, 31, 40, 50, 61, 74, 88, 103, 119, 137,
46  156, 176, 197, 220, 244, 270, 296, 325, 354, 386, 419, 453, 489, 526, 565,
47  606, 649, 693, 740, 788, 838, 890, 944, 1000, 1058, 1118, 1181, 1246, 1313, 1383,
48  1455, 1530, 1607, 1687, 1770, 1856, 1945, 2037, 2132, 2231, 2332, 2438, 2546, 2659, 2775,
49  2895, 3019, 3147, 3279, 3416, 3558, 3704, 3854, 4010, 4171, 4337, 4509, 4686, 4869, 5058,
50  5253, 5455, 5663, 5877, 6099, 6328, 6564, 6808, 7060, 7320, 7589, 7866, 8152, 8447, 8752,
51  9067, 9391, 9726, 10072, 10430, 10798, 11179, 11571, 11977, 12395, 12827, 13272, 13732, 14000};
52 
53  hname = "JetPt";
54  m_HistNames1D[hname] = new TH1F(hname, hname, 500, 0, 5000);
55 
56  hname = "JetEta";
57  m_HistNames1D[hname] = new TH1F(hname, hname, 120, -6, 6);
58 
59  hname = "JetPhi";
60  m_HistNames1D[hname] = new TH1F(hname, hname, 100, -M_PI, M_PI);
61 
62  hname = "NumberOfJets";
63  m_HistNames1D[hname] = new TH1F(hname, hname, 100, 0, 100);
64 
65  hname = "DijetMass";
66  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
67 
68  hname = "DijetMassWt";
69  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
70  m_HistNames1D.find(hname)->second->Sumw2();
71 
72  hname = "DijetMassIn";
73  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
74 
75  hname = "DijetMassInWt";
76  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
77  m_HistNames1D.find(hname)->second->Sumw2();
78 
79  hname = "DijetMassOut";
80  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
81 
82  hname = "DijetMassOutWt";
83  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
84  m_HistNames1D.find(hname)->second->Sumw2();
85 
86  hname = "ResonanceMass";
87  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
88 
89  hname = "DipartonMass";
90  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
91 
92  hname = "DipartonMassWt";
93  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
94  m_HistNames1D.find(hname)->second->Sumw2();
95 
96  hname = "DipartonMassIn";
97  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
98 
99  hname = "DipartonMassInWt";
100  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
101  m_HistNames1D.find(hname)->second->Sumw2();
102 
103  hname = "DipartonMassOut";
104  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
105 
106  hname = "DipartonMassOutWt";
107  m_HistNames1D[hname] = new TH1F(hname, hname, nMassBins, massBoundaries);
108  m_HistNames1D.find(hname)->second->Sumw2();
109 
110  hname = "PtHat";
111  m_HistNames1D[hname] = new TH1F(hname, hname, 1000, 0, 5000);
112 
113  hname = "PtHatWt";
114  m_HistNames1D[hname] = new TH1F(hname, hname, 1000, 0, 5000);
115  m_HistNames1D.find(hname)->second->Sumw2();
116 
117  hname = "PtHatFine";
118  m_HistNames1D[hname] = new TH1F(hname, hname, 5000, 0, 5000);
119 
120  hname = "PtHatFineWt";
121  m_HistNames1D[hname] = new TH1F(hname, hname, 5000, 0, 5000);
122  m_HistNames1D.find(hname)->second->Sumw2();
123 
124  mcTruthTree_ = new TTree("mcTruthTree", "mcTruthTree");
125  mcTruthTree_->Branch("xsec", &xsec, "xsec/F");
126  mcTruthTree_->Branch("weight", &weight, "weight/F");
127  mcTruthTree_->Branch("pt_hat", &pt_hat, "pt_hat/F");
128  mcTruthTree_->Branch("nJets", &nJets, "nJets/I");
129  mcTruthTree_->Branch("etaJet1", &etaJet1, "etaJet1/F");
130  mcTruthTree_->Branch("etaJet2", &etaJet2, "etaJet2/F");
131  mcTruthTree_->Branch("ptJet1", &ptJet1, "ptJet1/F");
132  mcTruthTree_->Branch("ptJet2", &ptJet2, "ptJet2/F");
133  mcTruthTree_->Branch("diJetMass", &diJetMass, "diJetMass/F");
134  mcTruthTree_->Branch("etaPart1", &etaPart1, "etaPart1/F");
135  mcTruthTree_->Branch("etaPart2", &etaPart2, "etaPart2/F");
136  mcTruthTree_->Branch("ptPart1", &ptPart1, "ptPart1/F");
137  mcTruthTree_->Branch("ptPart2", &ptPart2, "ptPart2/F");
138  mcTruthTree_->Branch("diPartMass", &diPartMass, "diPartMass/F");
139 }
141 template <class Jet>
142 void JetAnaPythia<Jet>::analyze(edm::Event const& evt, edm::EventSetup const& iSetup) {
143  int notDone = 1;
144  while (notDone) { //while loop to allow us to tailor the analysis level for faster running.
145  TString hname;
146 
147  // Process Info
148 
149  //edm::Handle< double > genEventScale;
150  //evt.getByLabel("genEventScale", genEventScale );
151  //pt_hat = *genEventScale;
152 
153  edm::Handle<edm::HepMCProduct> MCevt;
154  evt.getByLabel("generatorSmeared", MCevt);
155  HepMC::GenEvent* myGenEvent = new HepMC::GenEvent(*(MCevt->GetEvent()));
156 
157  double pthat = myGenEvent->event_scale();
158  pt_hat = float(pthat);
159 
160  delete myGenEvent;
161 
162  if (anaLevel != "generating") { //We are not generating events, so xsec is there
163  //edm::Handle< GenRunInfoProduct > genInfoProduct;
165  //xsec = (double)genInfoProduct->externalXSecLO();
166  xsec = 0.0;
167  if (ptHatEdges.size() > xsecGen.size()) {
168  for (unsigned int i_pthat = 0; i_pthat < xsecGen.size(); ++i_pthat) {
169  if (pthat >= ptHatEdges[i_pthat] && pthat < ptHatEdges[i_pthat + 1])
170  xsec = float(xsecGen[i_pthat]);
171  }
172  } else {
173  std::cout << "Number of PtHat bin edges too small. Xsec set to zero" << std::endl;
174  }
175  } else {
176  xsec = xsecGen[0]; //Generating events, no xsec in event, get xsec from user input
177  }
178  if (debug)
179  std::cout << "cross section=" << xsec << " pb" << std::endl;
180  weight = xsec / eventsGen;
181 
182  if (debug)
183  std::cout << "pt_hat=" << pt_hat << std::endl;
184  hname = "PtHat";
185  FillHist1D(hname, pt_hat, 1.0);
186  hname = "PtHatFine";
187  FillHist1D(hname, pt_hat, 1.0);
188  hname = "PtHatWt";
189  FillHist1D(hname, pt_hat, weight);
190  hname = "PtHatFineWt";
191  FillHist1D(hname, pt_hat, weight);
192  if (anaLevel == "PtHatOnly")
193  break; //ptHatOnly should be very fast
194 
195  // Jet Info
196  math::XYZTLorentzVector p4jet[2];
197  float etajet[2];
199  Handle<JetCollection> jets;
200  evt.getByLabel(JetAlgorithm, jets);
201  typename JetCollection::const_iterator i_jet;
202  int index = 0;
203 
205  hname = "NumberOfJets";
206  nJets = jets->size();
207  FillHist1D(hname, nJets, 1.0);
208 
209  // Two Leading Jet Info
210  for (i_jet = jets->begin(); i_jet != jets->end() && index < 2; ++i_jet) {
211  hname = "JetPt";
212  FillHist1D(hname, i_jet->pt(), 1.0);
213  hname = "JetEta";
214  FillHist1D(hname, i_jet->eta(), 1.0);
215  hname = "JetPhi";
216  FillHist1D(hname, i_jet->phi(), 1.0);
217  p4jet[index] = i_jet->p4();
218  etajet[index] = i_jet->eta();
219  if (debug)
220  std::cout << "jet " << index + 1 << ": pt=" << i_jet->pt() << ", eta=" << etajet[index] << std::endl;
221  index++;
222  }
223 
224  // TTree variables //
225  etaJet1 = etajet[0];
226  etaJet2 = etajet[1];
227  ptJet1 = p4jet[0].pt();
228  ptJet2 = p4jet[1].pt();
229  diJetMass = (p4jet[0] + p4jet[1]).mass();
230 
232  if (index == 2 && abs(etaJet1) < 1.3 && abs(etaJet2) < 1.3) {
233  hname = "DijetMass";
234  FillHist1D(hname, diJetMass, 1.0);
235  hname = "DijetMassWt";
236  FillHist1D(hname, diJetMass, weight);
237  }
238 
240  if (index == 2 && abs(etaJet1) < 0.7 && abs(etaJet2) < 0.7) {
241  hname = "DijetMassIn";
242  FillHist1D(hname, diJetMass, 1.0);
243  hname = "DijetMassInWt";
244  FillHist1D(hname, diJetMass, weight);
245  }
247  if (index == 2 && (abs(etaJet1) > 0.7 && abs(etaJet1) < 1.3) && (abs(etaJet2) > 0.7 && abs(etaJet2) < 1.3)) {
248  hname = "DijetMassOut";
249  FillHist1D(hname, diJetMass, 1.0);
250  hname = "DijetMassOutWt";
251  FillHist1D(hname, diJetMass, weight);
252  }
253  if (anaLevel == "Jets")
254  break; //Jets level for samples without genParticles
255 
256  // Parton Info
257  edm::Handle<std::vector<reco::GenParticle> > genParticlesHandle_;
258  evt.getByLabel("genParticles", genParticlesHandle_);
259  if (debug)
260  for (size_t i = 0; i < genParticlesHandle_->size(); ++i) {
261  const reco::GenParticle& p = (*genParticlesHandle_)[i];
262  int id = p.pdgId();
263  int st = p.status();
264  const math::XYZTLorentzVector& genP4 = p.p4();
265  if (i >= 2 && i <= 8)
266  std::cout << "particle " << i << ": id=" << id << ", status=" << st << ", mass=" << genP4.mass()
267  << ", pt=" << genP4.pt() << ", eta=" << genP4.eta() << std::endl;
268  }
269  // Examine the 7th particle in pythia.
270  // It should be either a resonance (abs(id)>=32) or the first outgoing parton
271  // for the processes we will consider: dijet resonances, QCD, or QCD +contact interactions.
272  const reco::GenParticle& p = (*genParticlesHandle_)[6];
273  int id = p.pdgId();
274  math::XYZTLorentzVector resonance_p, parton1_p, parton2_p;
275  if (abs(id) >= 32) {
277  resonance_p = p.p4();
278  hname = "ResonanceMass";
279  FillHist1D(hname, resonance_p.mass(), 1.0);
280  const reco::GenParticle& q = (*genParticlesHandle_)[7];
281  parton1_p = q.p4();
282  const reco::GenParticle& r = (*genParticlesHandle_)[8];
283  parton2_p = r.p4();
284  if (debug)
285  std::cout << "Resonance mass=" << resonance_p.mass() << ", parton 1 pt=" << parton1_p.pt()
286  << ", parton 2 pt=" << parton2_p.pt() << ", diparton mass=" << (parton1_p + parton2_p).mass()
287  << std::endl;
288  } else {
290  parton1_p = p.p4();
291  const reco::GenParticle& q = (*genParticlesHandle_)[7];
292  parton2_p = q.p4();
293  if (debug)
294  std::cout << "parton 1 pt=" << parton1_p.pt() << ", parton 2 pt=" << parton2_p.pt()
295  << ", diparton mass=" << (parton1_p + parton2_p).mass() << std::endl;
296  }
297 
298  etaPart1 = parton1_p.eta();
299  etaPart2 = parton2_p.eta();
300  ptPart1 = parton1_p.pt();
301  ptPart2 = parton2_p.pt();
302  diPartMass = (parton1_p + parton2_p).mass();
304  if (abs(etaPart1) < 1.3 && abs(etaPart2) < 1.3) {
305  hname = "DipartonMass";
306  FillHist1D(hname, diPartMass, 1.0);
307  hname = "DipartonMassWt";
308  FillHist1D(hname, diPartMass, weight);
309  }
311  if (abs(etaPart1) < 0.7 && abs(etaPart2) < 0.7) {
312  hname = "DipartonMassIn";
313  FillHist1D(hname, diPartMass, 1.0);
314  hname = "DipartonMassInWt";
315  FillHist1D(hname, diPartMass, weight);
316  }
318  if ((abs(etaPart1) > 0.7 && abs(etaPart1) < 1.3) && (abs(etaPart2) > 0.7 && abs(etaPart2) < 1.3)) {
319  hname = "DipartonMassOut";
320  FillHist1D(hname, diPartMass, 1.0);
321  hname = "DipartonMassOutWt";
322  FillHist1D(hname, diPartMass, weight);
323  }
324 
325  // Fill the TTree //
326  mcTruthTree_->Fill();
327 
328  notDone = 0; //We are done, exit the while loop
329  } //end of while
330 }
332 template <class Jet>
333 void JetAnaPythia<Jet>::endJob() {
335  if (m_file != nullptr) {
336  m_file->cd();
337  mcTruthTree_->Write();
338  for (std::map<TString, TH1*>::iterator hid = m_HistNames1D.begin(); hid != m_HistNames1D.end(); hid++)
339  hid->second->Write();
340  delete m_file;
341  m_file = nullptr;
342  }
343 }
345 template <class Jet>
346 void JetAnaPythia<Jet>::FillHist1D(const TString& histName, const Double_t& value, const Double_t& wt) {
347  std::map<TString, TH1*>::iterator hid = m_HistNames1D.find(histName);
348  if (hid == m_HistNames1D.end())
349  std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
350  else
351  hid->second->Fill(value, wt);
352 }
354 #include "FWCore/Framework/interface/MakerMacros.h"
356 typedef JetAnaPythia<CaloJet> CaloJetAnaPythia;
357 DEFINE_FWK_MODULE(CaloJetAnaPythia);
359 typedef JetAnaPythia<GenJet> GenJetAnaPythia;
360 DEFINE_FWK_MODULE(GenJetAnaPythia);
362 typedef JetAnaPythia<PFJet> PFJetAnaPythia;
363 DEFINE_FWK_MODULE(PFJetAnaPythia);
JetAnaPythia::beginJob
void beginJob() override
Definition: JetAnaPythia.cc:39
looper.cfg
tuple cfg
Definition: looper.py:296
mps_fire.i
i
Definition: mps_fire.py:428
JetAnaPythia::JetAnaPythia
JetAnaPythia(edm::ParameterSet const &cfg)
Definition: JetAnaPythia.cc:27
JetAnaPythia::endJob
void endJob() override
Definition: JetAnaPythia.cc:333
CaloJetAnaPythia
JetAnaPythia< CaloJet > CaloJetAnaPythia
Definition: JetAnaPythia.cc:356
JetAnaPythia::FillHist1D
void FillHist1D(const TString &histName, const Double_t &x, const Double_t &wt)
Definition: JetAnaPythia.cc:346
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
reco::LeafCandidate::status
int status() const final
status word
Definition: LeafCandidate.h:180
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:111
reco::LeafCandidate::p4
const LorentzVector & p4() const final
four-momentum Lorentz vector
Definition: LeafCandidate.h:114
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
reco::LeafCandidate::pdgId
int pdgId() const final
PDG identifier.
Definition: LeafCandidate.h:176
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
JetAnaPythia::analyze
void analyze(edm::Event const &e, edm::EventSetup const &iSetup) override
Definition: JetAnaPythia.cc:142
fwrapper::jets
vector< PseudoJet > jets
Definition: fastjetfortran_madfks.cc:44
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:500
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:50
debug
#define debug
Definition: HDRShower.cc:19
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
GenJetAnaPythia
JetAnaPythia< GenJet > GenJetAnaPythia
Definition: JetAnaPythia.cc:359
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
histoStyle.weight
int weight
Definition: histoStyle.py:51
PFJetAnaPythia
JetAnaPythia< PFJet > PFJetAnaPythia
Definition: JetAnaPythia.cc:362
Generated for CMSSW Reference Manual by doxygen 1.8.5