ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MyJetAnalysis.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MyJetAnalysis.cc
1 #include "MyJetAnalysis.h"
2 
6 
8 #include <phool/getClass.h>
9 
10 #include <g4eval/JetEvalStack.h>
11 
13 
14 #include <g4jets/JetMap.h>
15 
16 #include <TFile.h>
17 #include <TH1F.h>
18 #include <TH2F.h>
19 #include <TString.h>
20 #include <TTree.h>
21 #include <TVector3.h>
22 
23 #include <algorithm>
24 #include <cassert>
25 #include <cmath>
26 #include <iostream>
27 #include <limits>
28 #include <stdexcept>
29 
30 using namespace std;
31 
32 MyJetAnalysis::MyJetAnalysis(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename)
33  : SubsysReco("MyJetAnalysis_" + recojetname + "_" + truthjetname)
34  , m_recoJetName(recojetname)
35  , m_truthJetName(truthjetname)
36  , m_outputFileName(outputfilename)
37  , m_etaRange(-1, 1)
38  , m_ptRange(5, 100)
39  , m_trackJetMatchingRadius(.7)
40  , m_hInclusiveE(nullptr)
41  , m_hInclusiveEta(nullptr)
42  , m_hInclusivePhi(nullptr)
43  , m_T(nullptr)
44  , m_event(-1)
45  , m_id(-1)
46  , m_nComponent(-1)
47  , m_eta(numeric_limits<float>::signaling_NaN())
48  , m_phi(numeric_limits<float>::signaling_NaN())
49  , m_e(numeric_limits<float>::signaling_NaN())
50  , m_pt(numeric_limits<float>::signaling_NaN())
51  , m_truthID(-1)
52  , m_truthNComponent(-1)
53  , m_truthEta(numeric_limits<float>::signaling_NaN())
54  , m_truthPhi(numeric_limits<float>::signaling_NaN())
55  , m_truthE(numeric_limits<float>::signaling_NaN())
56  , m_truthPt(numeric_limits<float>::signaling_NaN())
57  , m_nMatchedTrack(-1)
58 {
59  m_trackdR.fill(numeric_limits<float>::signaling_NaN());
60  m_trackpT.fill(numeric_limits<float>::signaling_NaN());
61 }
62 
64 {
65 }
66 
68 {
70  cout << "MyJetAnalysis::Init - Outoput to " << m_outputFileName << endl;
71 
73 
74  // Histograms
75  m_hInclusiveE = new TH1F(
76  "hInclusive_E", //
77  TString(m_recoJetName) + " inclusive jet E;Total jet energy (GeV)", 100, 0, 100);
78 
80  new TH1F(
81  "hInclusive_eta", //
82  TString(m_recoJetName) + " inclusive jet #eta;#eta;Jet energy density", 50, -1, 1);
84  new TH1F(
85  "hInclusive_phi", //
86  TString(m_recoJetName) + " inclusive jet #phi;#phi;Jet energy density", 50, -M_PI, M_PI);
87 
88  //Trees
89  m_T = new TTree("T", "MyJetAnalysis Tree");
90 
91  // int m_event;
92  m_T->Branch("m_event", &m_event, "event/I");
93  // int m_id;
94  m_T->Branch("id", &m_id, "id/I");
95  // int m_nComponent;
96  m_T->Branch("nComponent", &m_nComponent, "nComponent/I");
97  // float m_eta;
98  m_T->Branch("eta", &m_eta, "eta/F");
99  // float m_phi;
100  m_T->Branch("phi", &m_phi, "phi/F");
101  // float m_e;
102  m_T->Branch("e", &m_e, "e/F");
103  // float m_pt;
104  m_T->Branch("pt", &m_pt, "pt/F");
105  //
106  // int m_truthID;
107  m_T->Branch("truthID", &m_truthID, "truthID/I");
108  // int m_truthNComponent;
109  m_T->Branch("truthNComponent", &m_truthNComponent, "truthNComponent/I");
110  // float m_truthEta;
111  m_T->Branch("truthEta", &m_truthEta, "truthEta/F");
112  // float m_truthPhi;
113  m_T->Branch("truthPhi", &m_truthPhi, "truthPhi/F");
114  // float m_truthE;
115  m_T->Branch("truthE", &m_truthE, "truthE/F");
116  // float m_truthPt;
117  m_T->Branch("truthPt", &m_truthPt, "truthPt/F");
118  //
119  // //! number of matched tracks
120  // int m_nMatchedTrack;
121  m_T->Branch("nMatchedTrack", &m_nMatchedTrack, "nMatchedTrack/I");
122  // std::array<float, kMaxMatchedTrack> m_trackdR;
123  m_T->Branch("id", m_trackdR.data(), "trackdR[nMatchedTrack]/F");
124  // std::array<float, kMaxMatchedTrack> m_trackpT;
125  m_T->Branch("id", m_trackpT.data(), "trackpT[nMatchedTrack]/F");
126 
128 }
129 
131 {
132  cout << "MyJetAnalysis::End - Output to " << m_outputFileName << endl;
134 
135  m_hInclusiveE->Write();
136  m_hInclusiveEta->Write();
137  m_hInclusivePhi->Write();
138  m_T->Write();
139 
141 }
142 
144 {
146  m_jetEvalStack->get_stvx_eval_stack()->set_use_initial_vertex(initial_vertex);
148 }
149 
151 {
153  cout << "MyJetAnalysis::process_event() entered" << endl;
154 
155  m_jetEvalStack->next_event(topNode);
156  JetRecoEval* recoeval = m_jetEvalStack->get_reco_eval();
157  ++m_event;
158 
159  // interface to jets
160  JetMap* jets = findNode::getClass<JetMap>(topNode, m_recoJetName);
161  if (!jets)
162  {
163  cout
164  << "MyJetAnalysis::process_event - Error can not find DST JetMap node "
165  << m_recoJetName << endl;
166  exit(-1);
167  }
168 
169  // interface to tracks
170  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
171  if (!trackmap)
172  {
173  trackmap = findNode::getClass<SvtxTrackMap>(topNode, "TrackMap");
174  if (!trackmap)
175  {
176  cout
177  << "MyJetAnalysis::process_event - Error can not find DST trackmap node SvtxTrackMap" << endl;
178  exit(-1);
179  }
180  }
181  for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
182  {
183  Jet* jet = iter->second;
184  assert(jet);
185 
186  bool eta_cut = (jet->get_eta() >= m_etaRange.first) and (jet->get_eta() <= m_etaRange.second);
187  bool pt_cut = (jet->get_pt() >= m_ptRange.first) and (jet->get_pt() <= m_ptRange.second);
188  if ((not eta_cut) or (not pt_cut))
189  {
191  {
192  cout << "MyJetAnalysis::process_event() - jet failed acceptance cut: ";
193  cout << "eta cut: " << eta_cut << ", ptcut: " << pt_cut << endl;
194  cout << "jet eta: " << jet->get_eta() << ", jet pt: " << jet->get_pt() << endl;
195  jet->identify();
196  }
197  continue;
198  }
199 
200  // fill histograms
201  assert(m_hInclusiveE);
202  m_hInclusiveE->Fill(jet->get_e());
203  assert(m_hInclusiveEta);
204  m_hInclusiveEta->Fill(jet->get_eta());
205  assert(m_hInclusivePhi);
206  m_hInclusivePhi->Fill(jet->get_phi());
207 
208  // fill trees - jet spectrum
209  Jet* truthjet = recoeval->max_truth_jet_by_energy(jet);
210 
211  m_id = jet->get_id();
212  m_nComponent = jet->size_comp();
213  m_eta = jet->get_eta();
214  m_phi = jet->get_phi();
215  m_e = jet->get_e();
216  m_pt = jet->get_pt();
217 
218  m_truthID = -1;
219  m_truthNComponent = -1;
220  m_truthEta = NAN;
221  m_truthPhi = NAN;
222  m_truthE = NAN;
223  m_truthPt = NAN;
224 
225  if (truthjet)
226  {
227  m_truthID = truthjet->get_id();
228  m_truthNComponent = truthjet->size_comp();
229  m_truthEta = truthjet->get_eta();
230  m_truthPhi = truthjet->get_phi();
231  m_truthE = truthjet->get_e();
232  m_truthPt = truthjet->get_pt();
233  }
234 
235  // fill trees - jet track matching
236  m_nMatchedTrack = 0;
237 
238  for (SvtxTrackMap::Iter iter = trackmap->begin();
239  iter != trackmap->end();
240  ++iter)
241  {
242  SvtxTrack* track = iter->second;
243 
244  TVector3 v(track->get_px(), track->get_py(), track->get_pz());
245  const double dEta = v.Eta() - m_eta;
246  const double dPhi = v.Phi() - m_phi;
247  const double dR = sqrt(dEta * dEta + dPhi * dPhi);
248 
249  if (dR < m_trackJetMatchingRadius)
250  {
251  //matched track to jet
252 
254 
256  m_trackpT[m_nMatchedTrack] = v.Perp();
257 
258  ++m_nMatchedTrack;
259  }
260 
262  {
263  cout << "MyJetAnalysis::process_event() - reached max track that matching a jet. Quit iterating tracks" << endl;
264  break;
265  }
266 
267  } // for (SvtxTrackMap::Iter iter = trackmap->begin();
268 
269  m_T->Fill();
270  } // for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
271 
273 }