Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JetValidation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JetValidation.cc
1 //module for producing a TTree with jet information for doing jet validation studies
2 // for questions/bugs please contact Virginia Bailey vbailey13@gsu.edu
3 #include <fun4all/Fun4AllBase.h>
4 #include "JetValidation.h"
5 #include <fun4all/Fun4AllReturnCodes.h>
6 #include <fun4all/PHTFileServer.h>
7 #include <phool/PHCompositeNode.h>
8 #include <phool/getClass.h>
9 #include <jetbase/JetMap.h>
10 #include <jetbase/JetContainer.h>
11 #include <jetbase/Jetv2.h>
12 #include <jetbase/Jetv1.h>
13 #include <centrality/CentralityInfo.h>
14 #include <globalvertex/GlobalVertex.h>
15 #include <globalvertex/GlobalVertexMap.h>
16 #include <calobase/RawTower.h>
17 #include <calobase/RawTowerContainer.h>
18 #include <calobase/RawTowerGeom.h>
19 #include <calobase/RawTowerGeomContainer.h>
20 #include <calobase/TowerInfoContainer.h>
21 #include <calobase/TowerInfo.h>
22 
23 #include <Pythia8/Pythia.h> // Include the Pythia header
24 #include <jetbackground/TowerBackground.h>
25 
26 #include <TTree.h>
27 #include <iostream>
28 #include <fastjet/PseudoJet.hh>
29 #include <sstream>
30 #include <iomanip>
31 #include <cmath>
32 #include <vector>
33 #include "fastjet/ClusterSequence.hh"
34 #include "fastjet/contrib/SoftDrop.hh" // In external code, this should be fastjet/contrib/SoftDrop.hh
35 
36 using namespace fastjet;
37 
38 // ROOT, for histogramming.
39 #include "TH1.h"
40 #include "TH2.h"
41 
42 // ROOT, for interactive graphics.
43 #include "TVirtualPad.h"
44 #include "TApplication.h"
45 
46 // ROOT, for saving file.
47 #include "TFile.h"
48 #include <TTree.h>
49 
50 //____________________________________________________________________________..
51 JetValidation::JetValidation(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename):
52  SubsysReco("JetValidation_" + recojetname + "_" + truthjetname)
53  , m_recoJetName(recojetname)
54  , m_truthJetName(truthjetname)
55  , m_outputFileName(outputfilename)
56  , m_etaRange(-1, 1)
57  , m_ptRange(5, 100)
58  , m_doTruthJets(0)
59  , m_doSeeds(0)
60  , m_doUnsubJet(0)
61  , m_T(nullptr)
62  , m_event(-1)
63  , m_nTruthJet(-1)
64  , m_nJet(-1)
65  , m_id()
66  , m_nComponent()
67  , m_eta()
68  , m_phi()
69  , m_e()
70  , m_pt()
71  , m_sub_et()
72  , m_truthID()
73  , m_truthNComponent()
74  , m_truthEta()
75  , m_truthPhi()
76  , m_truthE()
77  , m_truthPt()
78  , m_eta_rawseed()
79  , m_phi_rawseed()
80  , m_pt_rawseed()
81  , m_e_rawseed()
82  , m_rawseed_cut()
83  , m_eta_subseed()
84  , m_phi_subseed()
85  , m_pt_subseed()
86  , m_e_subseed()
87  , m_subseed_cut()
88 {
89  std::cout << "JetValidation::JetValidation(const std::string &name) Calling ctor" << std::endl;
90 }
91 
92 //____________________________________________________________________________..
93 JetValidation::~JetValidation()
94 {
95  std::cout << "JetValidation::~JetValidation() Calling dtor" << std::endl;
96 }
97 
98 //____________________________________________________________________________..
99 int JetValidation::Init(PHCompositeNode *topNode)
100 {
101  std::cout << "JetValidation::Init(PHCompositeNode *topNode) Initializing" << std::endl;
102  PHTFileServer::get().open(m_outputFileName, "RECREATE");
103  std::cout << "JetValidation::Init - Output to " << m_outputFileName << std::endl;
104 
105  // configure Tree
106  m_T = new TTree("T", "MyJetAnalysis Tree");
107  m_T->Branch("m_event", &m_event, "event/I");
108  m_T->Branch("nJet", &m_nJet, "nJet/I");
109  m_T->Branch("cent", &m_centrality);
110  m_T->Branch("zvtx", &m_zvtx);
111  m_T->Branch("b", &m_impactparam);
112  m_T->Branch("id", &m_id);
113  m_T->Branch("nComponent", &m_nComponent);
114 
115  m_T->Branch("eta", &m_eta);
116  m_T->Branch("phi", &m_phi);
117  m_T->Branch("e", &m_e);
118  m_T->Branch("pt", &m_pt);
119  if(m_doUnsubJet)
120  {
121  m_T->Branch("pt_unsub", &m_unsub_pt);
122  m_T->Branch("subtracted_et", &m_sub_et);
123  }
124  if(m_doTruthJets){
125  m_T->Branch("nTruthJet", &m_nTruthJet);
126  m_T->Branch("truthID", &m_truthID);
127  m_T->Branch("truthNComponent", &m_truthNComponent);
128  m_T->Branch("truthEta", &m_truthEta);
129  m_T->Branch("truthPhi", &m_truthPhi);
130  m_T->Branch("truthE", &m_truthE);
131  m_T->Branch("truthPt", &m_truthPt);
132  }
133 
134  if(m_doSeeds){
135  m_T->Branch("rawseedEta", &m_eta_rawseed);
136  m_T->Branch("rawseedPhi", &m_phi_rawseed);
137  m_T->Branch("rawseedPt", &m_pt_rawseed);
138  m_T->Branch("rawseedE", &m_e_rawseed);
139  m_T->Branch("rawseedCut", &m_rawseed_cut);
140  m_T->Branch("subseedEta", &m_eta_subseed);
141  m_T->Branch("subseedPhi", &m_phi_subseed);
142  m_T->Branch("subseedPt", &m_pt_subseed);
143  m_T->Branch("subseedE", &m_e_subseed);
144  m_T->Branch("subseedCut", &m_subseed_cut);
145  }
146 
147  return Fun4AllReturnCodes::EVENT_OK;
148 }
149 
150 //____________________________________________________________________________..
151 int JetValidation::InitRun(PHCompositeNode *topNode)
152 {
153  std::cout << "JetValidation::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
154  return Fun4AllReturnCodes::EVENT_OK;
155 }
156 
157 //____________________________________________________________________________..
158 int JetValidation::process_event(PHCompositeNode *topNode)
159 {
160  std::cout << "JetValidation::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
161  ++m_event;
162 
163  // interface to reco jets
164  JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
165  if (!jets)
166  {
167  std::cout
168  << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node "
169  << m_recoJetName << std::endl;
170  exit(-1);
171  }
172 
173  //interface to truth jets
174  JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
175  if (!jetsMC && m_doTruthJets)
176  {
177  std::cout
178  << "MyJetAnalysis::process_event - Error can not find DST Truth JetMap node "
179  << m_truthJetName << std::endl;
180  exit(-1);
181  }
182 
183  // interface to jet seeds
184  JetContainer* seedjetsraw = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
185  if (!seedjetsraw && m_doSeeds)
186  {
187  std::cout
188  << "MyJetAnalysis::process_event - Error can not find DST raw seed jets "
189  << std::endl;
190  exit(-1);
191  }
192 
193  JetContainer* seedjetssub = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_HIRecoSeedsSub_r02");
194  if (!seedjetssub && m_doSeeds)
195  {
196  std::cout
197  << "MyJetAnalysis::process_event - Error can not find DST subtracted seed jets "
198  << std::endl;
199  exit(-1);
200  }
201 
202  //centrality
203  CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
204  if (!cent_node)
205  {
206  std::cout
207  << "MyJetAnalysis::process_event - Error can not find centrality node "
208  << std::endl;
209  exit(-1);
210  }
211 
212  //zvertex
213  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
214  GlobalVertex *vtx = vertexmap->begin()->second;
215  m_zvtx = vtx->get_z();
216 
217 
218  //calorimeter towers
219  TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");//_SUB1
220  TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");//_SUB1
221  TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");//_SUB1 added at the ends but I am removing them for now
222  RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
223  RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
224  if(!towersEM3 || !towersIH3 || !towersOH3){
225  std::cout
226  <<"MyJetAnalysis::process_event - Error can not find raw tower node "
227  << std::endl;
228  exit(-1);
229  }
230 
231  if(!tower_geom || !tower_geomOH){
232  std::cout
233  <<"MyJetAnalysis::process_event - Error can not find raw tower geometry "
234  << std::endl;
235  exit(-1);
236  }
237 
238  //underlying event
239  TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
240  if(!background){
241  std::cout<<"Can't get background. Exiting"<<std::endl;
242  return Fun4AllReturnCodes::EVENT_OK;
243  }
244 
245  //get the event centrality/impact parameter from HIJING
246  m_centrality = cent_node->get_centile(CentralityInfo::PROP::bimp);
247  m_impactparam = cent_node->get_quantity(CentralityInfo::PROP::bimp);
248 
249  //get reco jets
250  m_nJet = 0;
251  float background_v2 = 0;
252  float background_Psi2 = 0;
253  if(m_doUnsubJet)
254  {
255  background_v2 = background->get_v2();
256  background_Psi2 = background->get_Psi2();
257  }
258  /* for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter) */
259  for (auto jet : *jets)
260  {
261 
262  if(jet->get_pt() < 1) continue; // to remove noise jets
263 
264  m_id.push_back(jet->get_id());
265  m_nComponent.push_back(jet->size_comp());
266  m_eta.push_back(jet->get_eta());
267  m_phi.push_back(jet->get_phi());
268  m_e.push_back(jet->get_e());
269  m_pt.push_back(jet->get_pt());
270  /*
271  if(m_doUnsubJet)
272  {
273  Jet* unsubjet = new Jetv1();
274  */
275  float totalPx = 0;
276  float totalPy = 0;
277  float totalPz = 0;
278  float totalE = 0;
279  int nconst = 0;
280 
281  // for (Jet::ConstIter comp = jet->begin_comp(); comp != jet->end_comp(); ++comp)
282  for (auto comp: jet->get_comp_vec())
283  {
284  TowerInfo *tower;
285  nconst++;
286  unsigned int channel = comp.second;
287 
288  if (comp.first == 15 || comp.first == 30)
289  {
290  tower = towersIH3->get_tower_at_channel(channel);
291  if(!tower || !tower_geom){
292  continue;
293  }
294  unsigned int calokey = towersIH3->encode_key(channel);
295  int ieta = towersIH3->getTowerEtaBin(calokey);
296  int iphi = towersIH3->getTowerPhiBin(calokey);
297  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
298  float UE = background->get_UE(1).at(ieta);
299  float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
300  float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
301 
302  UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
303  totalE += tower->get_energy() + UE;
304  double pt = (tower->get_energy() + UE) / cosh(tower_eta);
305  totalPx += pt * cos(tower_phi);
306  totalPy += pt * sin(tower_phi);
307  totalPz += pt * sinh(tower_eta);
308  }
309  else if (comp.first == 16 || comp.first == 31)
310  {
311  tower = towersOH3->get_tower_at_channel(channel);
312  if(!tower || !tower_geomOH)
313  {
314  continue;
315  }
316 
317  unsigned int calokey = towersOH3->encode_key(channel);
318  int ieta = towersOH3->getTowerEtaBin(calokey);
319  int iphi = towersOH3->getTowerPhiBin(calokey);
320  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
321  float UE = background->get_UE(2).at(ieta);
322  float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
323  float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
324 
325  UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
326  totalE +=tower->get_energy() + UE;
327  double pt = (tower->get_energy() + UE) / cosh(tower_eta);
328  totalPx += pt * cos(tower_phi);
329  totalPy += pt * sin(tower_phi);
330  totalPz += pt * sinh(tower_eta);
331  }
332  else if (comp.first == 14 || comp.first == 29)
333  {
334  tower = towersEM3->get_tower_at_channel(channel);
335  if(!tower || !tower_geom)
336  {
337  continue;
338  }
339 
340  unsigned int calokey = towersEM3->encode_key(channel);
341  int ieta = towersEM3->getTowerEtaBin(calokey);
342  int iphi = towersEM3->getTowerPhiBin(calokey);
343  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
344  float UE = background->get_UE(0).at(ieta);
345  float tower_phi = tower_geom->get_tower_geometry(key)->get_phi();
346  float tower_eta = tower_geom->get_tower_geometry(key)->get_eta();
347 
348  UE = UE * (1 + 2 * background_v2 * cos(2 * (tower_phi - background_Psi2)));
349  totalE +=tower->get_energy() + UE;
350  double pt = (tower->get_energy() + UE) / cosh(tower_eta);
351  totalPx += pt * cos(tower_phi);
352  totalPy += pt * sin(tower_phi);
353  totalPz += pt * sinh(tower_eta);
354 
355  }
356  }
357  /* //get unsubtracted jet
358  unsubjet->set_px(totalPx);
359  unsubjet->set_py(totalPy);
360  unsubjet->set_pz(totalPz);
361  unsubjet->set_e(totalE);
362  m_unsub_pt.push_back(unsubjet->get_pt());
363  m_sub_et.push_back(unsubjet->get_et() - jet->get_et());
364  */
365 
366 
367  m_nJet++;
368  }
369 
370  //get truth jets
371  if(m_doTruthJets)
372  {
373  m_nTruthJet = 0;
374  for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter)
375  {
376  Jet* truthjet = iter->second;
377 
378  bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
379  bool pt_cut = (truthjet->get_pt() >= m_ptRange.first) and (truthjet->get_pt() <= m_ptRange.second);
380  if ((not eta_cut) or (not pt_cut)) continue;
381  m_truthID.push_back(truthjet->get_id());
382  m_truthNComponent.push_back(truthjet->size_comp());
383  m_truthEta.push_back(truthjet->get_eta());
384  m_truthPhi.push_back(truthjet->get_phi());
385  m_truthE.push_back(truthjet->get_e());
386  m_truthPt.push_back(truthjet->get_pt());
387  m_nTruthJet++;
388  }
389  }
390 
391  //get seed jets
392  if(m_doSeeds)
393  {
394  for (auto jet : *seedjetsraw)
395  {
396  int passesCut = jet->get_property(Jet::PROPERTY::prop_SeedItr);
397  m_eta_rawseed.push_back(jet->get_eta());
398  m_phi_rawseed.push_back(jet->get_phi());
399  m_e_rawseed.push_back(jet->get_e());
400  m_pt_rawseed.push_back(jet->get_pt());
401  m_rawseed_cut.push_back(passesCut);
402  }
403 
404  for (auto jet : *seedjetssub) //JetMap::Iter iter = seedjetssub->begin(); iter != seedjetssub->end(); ++iter)
405  {
406  int passesCut = jet->get_property(Jet::PROPERTY::prop_SeedItr);
407  m_eta_subseed.push_back(jet->get_eta());
408  m_phi_subseed.push_back(jet->get_phi());
409  m_e_subseed.push_back(jet->get_e());
410  m_pt_subseed.push_back(jet->get_pt());
411  m_subseed_cut.push_back(passesCut);
412  }
413  }
414 
415  //fill the tree
416  m_T->Fill();
417 
418  return Fun4AllReturnCodes::EVENT_OK;
419 }
420 
421 //____________________________________________________________________________..
422 int JetValidation::ResetEvent(PHCompositeNode *topNode)
423 {
424  //std::cout << "JetValidation::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
425  m_id.clear();
426  m_nComponent.clear();
427  m_eta.clear();
428  m_phi.clear();
429  m_e.clear();
430  m_pt.clear();
431  m_unsub_pt.clear();
432  m_sub_et.clear();
433 
434  m_truthID.clear();
435  m_truthNComponent.clear();
436  m_truthEta.clear();
437  m_truthPhi.clear();
438  m_truthE.clear();
439  m_truthPt.clear();
440  m_truthdR.clear();
441 
442  m_eta_subseed.clear();
443  m_phi_subseed.clear();
444  m_e_subseed.clear();
445  m_pt_subseed.clear();
446  m_subseed_cut.clear();
447 
448  m_eta_rawseed.clear();
449  m_phi_rawseed.clear();
450  m_e_rawseed.clear();
451  m_pt_rawseed.clear();
452  m_rawseed_cut.clear();
453  return Fun4AllReturnCodes::EVENT_OK;
454 }
455 
456 //____________________________________________________________________________..
457 int JetValidation::EndRun(const int runnumber)
458 {
459  std::cout << "JetValidation::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
460  return Fun4AllReturnCodes::EVENT_OK;
461 }
462 
463 //____________________________________________________________________________..
464 int JetValidation::End(PHCompositeNode *topNode)
465 {
466  std::cout << "JetValidation::End - Output to " << m_outputFileName << std::endl;
467  PHTFileServer::get().cd(m_outputFileName);
468 
469  m_T->Write();
470  std::cout << "JetValidation::End(PHCompositeNode *topNode) This is the End..." << std::endl;
471  return Fun4AllReturnCodes::EVENT_OK;
472 }
473 
474 //____________________________________________________________________________..
475 int JetValidation::Reset(PHCompositeNode *topNode)
476 {
477  std::cout << "JetValidation::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
478  return Fun4AllReturnCodes::EVENT_OK;
479 }
480 
481 //____________________________________________________________________________..
482 void JetValidation::Print(const std::string &what) const
483 {
484  std::cout << "JetValidation::Print(const std::string &what) const Printing info for " << what << std::endl;
485 }