ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DIRCTree.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DIRCTree.cc
1 #include "G4DIRCTree.h"
2 
3 #include "PrtHit.h"
4 
5 #include <g4main/PHG4Hit.h>
6 #include <g4main/PHG4HitContainer.h> // for PHG4HitContainer, PHG4Hit...
7 #include <g4main/PHG4Particle.h>
8 #include <g4main/PHG4TruthInfoContainer.h>
9 
10 #include <fun4all/SubsysReco.h> // for SubsysReco
11 
12 #include <phool/getClass.h>
13 
14 #include <TFile.h>
15 #include <TTree.h>
16 
17 #include <cmath> // for atan2, sqrt
18 #include <cstring> // for strcmp
19 #include <iostream> // for ostringstream, operator<<
20 #include <sstream>
21 #include <utility> // for pair
22 
23 G4DIRCTree::G4DIRCTree(const std::string &name, const std::string &filename)
24  : SubsysReco(name)
25  , nblocks(0)
26  , _filename(filename)
27  , g4tree(nullptr)
28  , outfile(nullptr)
29 {
30 }
31 
32 int G4DIRCTree::Init(PHCompositeNode *)
33 {
34  outfile = new TFile(_filename.c_str(), "RECREATE");
35  g4tree = new TTree("mG4EvtTree", "g4tree");
36  g4tree->SetAutoSave(1000000);
37 
38  std::cout << "Initialize Geant 4 Tree ... << " << std::endl;
39 
41  g4tree->Branch("momentum", &mG4EvtTree.momentum, "p/D");
42  g4tree->Branch("theta", &mG4EvtTree.theta, "theta/D");
43  g4tree->Branch("phi", &mG4EvtTree.phi, "phi/D");
44  g4tree->Branch("px", &mG4EvtTree.px, "px/D");
45  g4tree->Branch("py", &mG4EvtTree.py, "py/D");
46  g4tree->Branch("pz", &mG4EvtTree.pz, "pz/D");
47  g4tree->Branch("pid", &mG4EvtTree.pid, "pid/I");
48 
49  //g4tree->Branch("track_mom_bar", mG4EvtTree.track_mom_bar, "track_mom_bar[3]/D");
50  //g4tree->Branch("track_hit_pos_bar", mG4EvtTree.track_hit_pos_bar, "track_hit_pos_bar[3]/D");
51 
53  g4tree->Branch("nhits", &mG4EvtTree.nhits, "nhits/I");
54  g4tree->Branch("detid", mG4EvtTree.detid, "detid[nhits]/I");
55  g4tree->Branch("hitid", mG4EvtTree.hitid, "hitid[nhits]/I");
56  g4tree->Branch("trkid", mG4EvtTree.trkid, "trkid[nhits]/I");
57  g4tree->Branch("x0", mG4EvtTree.x0, "x0[nhits]/D");
58  g4tree->Branch("y0", mG4EvtTree.y0, "y0[nhits]/D");
59  g4tree->Branch("z0", mG4EvtTree.z0, "z0[nhits]/D");
60  g4tree->Branch("x1", mG4EvtTree.x1, "x1[nhits]/D");
61  g4tree->Branch("y1", mG4EvtTree.y1, "y1[nhits]/D");
62  g4tree->Branch("z1", mG4EvtTree.z1, "z1[nhits]/D");
63  g4tree->Branch("edep", mG4EvtTree.edep, "edep[nhits]/D");
64 
65  g4tree->Branch("mcp_id", mG4EvtTree.mcp_id, "mcp_id[nhits]/I");
66  g4tree->Branch("pixel_id", mG4EvtTree.pixel_id, "pixel_id[nhits]/I");
67  g4tree->Branch("lead_time", mG4EvtTree.lead_time, "lead_time[nhits]/D");
68  g4tree->Branch("wavelength", mG4EvtTree.wavelength, "wavelength[nhits]/D");
69  g4tree->Branch("hit_pathId", mG4EvtTree.hit_pathId, "hit_pathId[nhits]/L");
70  g4tree->Branch("nrefl", mG4EvtTree.nrefl, "nrefl[nhits]/I");
71 
72  g4tree->Branch("hit_globalPos", mG4EvtTree.hit_globalPos, "hit_globalPos[nhits][3]/D");
73  g4tree->Branch("hit_localPos", mG4EvtTree.hit_localPos, "hit_localPos[nhits][3]/D");
74  g4tree->Branch("hit_digiPos", mG4EvtTree.hit_digiPos, "hit_digiPos[nhits][3]/D");
75  g4tree->Branch("hit_mom", mG4EvtTree.hit_mom, "hit_mom[nhits][3]/D");
76  g4tree->Branch("hit_pos", mG4EvtTree.hit_pos, "hit_pos[nhits][3]/D");
77  //g4tree->Branch("track_mom_bar", mG4EvtTree.track_mom_bar, "track_mom_bar[nhits][3]/D");
78  //g4tree->Branch("track_hit_pos_bar", mG4EvtTree.track_hit_pos_bar, "track_hit_pos_bar[nhits][3]/D");
79 
80  return 0;
81 }
82 
83 int G4DIRCTree::process_event(PHCompositeNode *topNode)
84 {
85  // get the primary particle which did this to us....
86  PHG4TruthInfoContainer *truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
87 
88  const PHG4TruthInfoContainer::Range primRange = truthInfoList->GetPrimaryParticleRange();
89  double px = primRange.first->second->get_px();
90  double py = primRange.first->second->get_py();
91  double pz = primRange.first->second->get_pz();
92  double p = sqrt(px * px + py * py + pz * pz);
93  double pt = sqrt(px * px + py * py);
94  double phi = atan2(py, px);
95  double theta = atan2(pt, pz);
96  double pid = primRange.first->second->get_pid();
97 
98  mG4EvtTree.momentum = p;
99  mG4EvtTree.theta = theta;
100  mG4EvtTree.phi = phi;
101  mG4EvtTree.px = px;
102  mG4EvtTree.py = py;
103  mG4EvtTree.pz = pz;
104  mG4EvtTree.pid = pid;
105 
106  /*for(int i=0; i < 3; i++)
107  {
108  mG4EvtTree.track_mom_bar[i] = bar_vectors::get_p_bar()[evt_num](i);
109  mG4EvtTree.track_hit_pos_bar[i] = bar_vectors::get_pos_bar()[evt_num](i);
110  }*/
111 
112  int nhits = 0;
113 
114  std::ostringstream nodename;
115  std::set<std::string>::const_iterator iter;
116 
117  for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
118  {
119  int detid = (_detid.find(*iter))->second;
120  nodename.str("");
121  nodename << "G4HIT_" << *iter;
122  PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
123 
124  if (!strcmp("G4HIT_DIRC", nodename.str().c_str())) // DIRC
125  {
126  process_hit(hits, "G4HIT_DIRC", detid, nhits);
127  }
128  }
129 
130  mG4EvtTree.nhits = nhits;
131 
132  if (g4tree) g4tree->Fill();
133 
134  evt_num++;
135 
136  return 0;
137 }
138 
139 int G4DIRCTree::End(PHCompositeNode *topNode)
140 {
141  outfile->cd();
142  g4tree->Write();
143  outfile->Write();
144  outfile->Close();
145  delete outfile;
146  return 0;
147 }
148 
149 void G4DIRCTree::AddNode(const std::string &name, const int detid)
150 {
151  _node_postfix.insert(name);
152  _detid[name] = detid;
153  return;
154 }
155 
156 int G4DIRCTree::process_hit(PHG4HitContainer *hits, const std::string &dName, int detid, int &nhits)
157 {
158  if (hits)
159  {
160  PHG4HitContainer::ConstRange hit_range_0 = hits->getHits();
161  for (PHG4HitContainer::ConstIterator hit_iter_0 = hit_range_0.first; hit_iter_0 != hit_range_0.second; hit_iter_0++)
162  {
163  PrtHit *dirc_hit = dynamic_cast<PrtHit *>(hit_iter_0->second);
164 
165  mG4EvtTree.detid[nhits] = detid;
166  mG4EvtTree.hitid[nhits] = dirc_hit->get_hit_id();
167  mG4EvtTree.trkid[nhits] = dirc_hit->get_trkid();
168 
169  mG4EvtTree.x0[nhits] = dirc_hit->get_x(0);
170  mG4EvtTree.y0[nhits] = dirc_hit->get_y(0);
171  mG4EvtTree.z0[nhits] = dirc_hit->get_z(0);
172  mG4EvtTree.x1[nhits] = dirc_hit->get_x(1);
173  mG4EvtTree.y1[nhits] = dirc_hit->get_y(1);
174  mG4EvtTree.z1[nhits] = dirc_hit->get_z(1);
175  mG4EvtTree.edep[nhits] = dirc_hit->get_edep();
176 
177  mG4EvtTree.mcp_id[nhits] = dirc_hit->GetMcpId();
178  mG4EvtTree.pixel_id[nhits] = dirc_hit->GetPixelId();
179  mG4EvtTree.lead_time[nhits] = dirc_hit->GetLeadTime();
180  mG4EvtTree.wavelength[nhits] = dirc_hit->GetTotTime();
181  mG4EvtTree.hit_pathId[nhits] = dirc_hit->GetPathInPrizm();
182  mG4EvtTree.nrefl[nhits] = dirc_hit->GetNreflectionsInPrizm();
183 
184  for (int i = 0; i < 3; i++)
185  {
186  mG4EvtTree.hit_globalPos[nhits][i] = dirc_hit->GetGlobalPos()(i);
187  mG4EvtTree.hit_localPos[nhits][i] = dirc_hit->GetLocalPos()(i);
188  mG4EvtTree.hit_digiPos[nhits][i] = dirc_hit->GetDigiPos()(i);
189  mG4EvtTree.hit_mom[nhits][i] = dirc_hit->GetMomentum()(i);
190  mG4EvtTree.hit_pos[nhits][i] = dirc_hit->GetPosition()(i);
191  }
192 
193  nhits++;
194  }
195 
196  /*PHG4HitContainer::ConstRange hit_range_1 = hits->getHits();
197  for (PHG4HitContainer::ConstIterator hit_iter_1 = hit_range_1.first; hit_iter_1 != hit_range_1.second; hit_iter_1++)
198  {
199  PrtHit *dirc_hit = dynamic_cast<PrtHit*>(hit_iter_1->second);
200 
201  mG4EvtTree.detid[nhits] = detid;
202 
203  for(int i=0; i < 3; i++)
204  {
205  mG4EvtTree.track_mom_bar[nhits][i] = dirc_hit->GetMomentumAtBar()(i);
206  mG4EvtTree.track_hit_pos_bar[nhits][i] = dirc_hit->GetPositionAtBar()(i);
207  }
208 
209  nhits++;
210  }*/
211  }
212 
213  return 0;
214 }