ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SnglTree.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4SnglTree.cc
1 #include "G4SnglTree.h"
2 
3 #include <g4main/PHG4Hit.h>
4 #include <g4main/PHG4HitContainer.h> // for PHG4HitContainer, PHG4Hit...
5 #include <g4main/PHG4Particle.h>
7 
8 #include <fun4all/SubsysReco.h> // for SubsysReco
9 
10 #include <phool/getClass.h>
11 
12 #include <TFile.h>
13 #include <TTree.h>
14 
15 #include <cmath> // for atan2, sqrt
16 #include <cstring> // for strcmp
17 #include <iostream> // for ostringstream, operator<<
18 #include <sstream>
19 #include <utility> // for pair
20 
21 using namespace std;
22 
23 G4SnglTree::G4SnglTree(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 
33 {
34  outfile = new TFile(_filename.c_str(), "RECREATE");
35  g4tree = new TTree("mG4EvtTree", "g4tree");
36  g4tree->SetAutoSave(1000000);
37 
38  cout << "Initialize Geant 4 Tree ... << " << endl;
39 
41  g4tree->Branch("energy", &mG4EvtTree.energy, "energy/F");
42  g4tree->Branch("theta", &mG4EvtTree.theta, "theta/F");
43  g4tree->Branch("phi", &mG4EvtTree.phi, "phi/F");
44  g4tree->Branch("px", &mG4EvtTree.px, "px/F");
45  g4tree->Branch("py", &mG4EvtTree.py, "py/F");
46  g4tree->Branch("pz", &mG4EvtTree.pz, "pz/F");
47  g4tree->Branch("cemcactLayers", &mG4EvtTree.cemcactLayers, "cemcactLayers/I");
48  g4tree->Branch("cemcabsLayers", &mG4EvtTree.cemcabsLayers, "cemcabsLayers/I");
49  g4tree->Branch("hcalactLayers", &mG4EvtTree.hcalactLayers, "hcalactLayers/I");
50  g4tree->Branch("hcalabsLayers", &mG4EvtTree.hcalabsLayers, "hcalabsLayers/I");
51  g4tree->Branch("cemcactESum", mG4EvtTree.cemcactESum, "cemcactESum[cemcactLayers]/F");
52  g4tree->Branch("cemcabsESum", mG4EvtTree.cemcabsESum, "cemcabsESum[cemcabsLayers]/F");
53  g4tree->Branch("hcalactESum", mG4EvtTree.hcalactESum, "hcalactESum[hcalactLayers]/F");
54  g4tree->Branch("hcalabsESum", mG4EvtTree.hcalabsESum, "hcalabsESum[hcalabsLayers]/F");
55 
57  g4tree->Branch("nhits", &mG4EvtTree.nhits, "nhits/I");
58  g4tree->Branch("detid", mG4EvtTree.detid, "detid[nhits]/I");
59  g4tree->Branch("layer", mG4EvtTree.layer, "layer[nhits]/I");
60  g4tree->Branch("hitid", mG4EvtTree.hitid, "hitid[nhits]/I");
61  g4tree->Branch("trkid", mG4EvtTree.trkid, "trkid[nhits]/I");
62  g4tree->Branch("scintid", mG4EvtTree.scintid, "scintid[nhits]/I");
63  g4tree->Branch("x0", mG4EvtTree.x0, "x0[nhits]/F");
64  g4tree->Branch("y0", mG4EvtTree.y0, "y0[nhits]/F");
65  g4tree->Branch("z0", mG4EvtTree.z0, "z0[nhits]/F");
66  g4tree->Branch("x1", mG4EvtTree.x1, "x1[nhits]/F");
67  g4tree->Branch("y1", mG4EvtTree.y1, "y1[nhits]/F");
68  g4tree->Branch("z1", mG4EvtTree.z1, "z1[nhits]/F");
69  g4tree->Branch("edep", mG4EvtTree.edep, "edep[nhits]/F");
70 
71  return 0;
72 }
73 
75 {
76  // get the primary particle which did this to us....
77  PHG4TruthInfoContainer *truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
78 
79  const PHG4TruthInfoContainer::Range primRange = truthInfoList->GetPrimaryParticleRange();
80  double px = primRange.first->second->get_px();
81  double py = primRange.first->second->get_py();
82  double pz = primRange.first->second->get_pz();
83  double e = primRange.first->second->get_e();
84  double pt = sqrt(px * px + py * py);
85  double phi = atan2(py, px);
86  double theta = atan2(pt, pz);
87 
90  mG4EvtTree.phi = phi;
91  mG4EvtTree.px = px;
92  mG4EvtTree.py = py;
93  mG4EvtTree.pz = pz;
94 
95  int nhits = 0;
96 
97  ostringstream nodename;
98  set<string>::const_iterator iter;
99 
100  for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
101  {
102  int detid = (_detid.find(*iter))->second;
103  nodename.str("");
104  nodename << "G4HIT_" << *iter;
105  PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
106 
107  if (!strcmp("G4HIT_CEMC", nodename.str().c_str())) //CEMC scintillator
108  {
109  mG4EvtTree.cemcactLayers = process_hit(hits, "G4HIT_CEMC", detid, nhits);
110  }
111  else if (!strcmp("G4HIT_ABSORBER_CEMC", nodename.str().c_str())) //CEMC Aabsorber G4_W
112  {
113  mG4EvtTree.cemcabsLayers = process_hit(hits, "G4HIT_ABSORBER_CEMC", detid, nhits);
114  }
115  else if (!strcmp("G4HIT_HCAL", nodename.str().c_str())) //HCAL Active scintilltor
116  {
117  mG4EvtTree.hcalactLayers = process_hit(hits, "G4HIT_HCAL", detid, nhits);
118  }
119  else if (!strcmp("G4HIT_ABSORBER_HCAL", nodename.str().c_str())) //HCAL Aabsorber steel
120  {
121  mG4EvtTree.hcalabsLayers = process_hit(hits, "G4HIT_ABSORBER_HCAL", detid, nhits);
122  }
123  }
124 
125  mG4EvtTree.nhits = nhits;
126 
127  if (g4tree) g4tree->Fill();
128 
129  return 0;
130 }
131 
133 {
134  outfile->cd();
135  g4tree->Write();
136  outfile->Write();
137  outfile->Close();
138  delete outfile;
139  return 0;
140 }
141 
142 void G4SnglTree::AddNode(const std::string &name, const int detid)
143 {
144  _node_postfix.insert(name);
145  _detid[name] = detid;
146  return;
147 }
148 
149 int G4SnglTree::process_hit(PHG4HitContainer *hits, const string &dName, int detid, int &nhits)
150 {
151  map<int, double> layer_edep_map;
152  map<int, double>::const_iterator edepiter;
153 
154  int nLayers = 0;
155  if (hits)
156  {
157  // double esum = 0;
158  PHG4HitContainer::ConstRange hit_range = hits->getHits();
159  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
160  {
161  layer_edep_map[hit_iter->second->get_layer()] += hit_iter->second->get_edep();
162  mG4EvtTree.detid[nhits] = detid;
163  mG4EvtTree.layer[nhits] = hit_iter->second->get_layer();
164  mG4EvtTree.hitid[nhits] = hit_iter->second->get_hit_id();
165  mG4EvtTree.trkid[nhits] = hit_iter->second->get_trkid();
166  mG4EvtTree.scintid[nhits] = hit_iter->second->get_scint_id();
167  mG4EvtTree.x0[nhits] = hit_iter->second->get_x(0);
168  mG4EvtTree.y0[nhits] = hit_iter->second->get_y(0);
169  mG4EvtTree.z0[nhits] = hit_iter->second->get_z(0);
170  mG4EvtTree.x1[nhits] = hit_iter->second->get_x(1);
171  mG4EvtTree.y1[nhits] = hit_iter->second->get_y(1);
172  mG4EvtTree.z1[nhits] = hit_iter->second->get_z(1);
173  mG4EvtTree.edep[nhits] = hit_iter->second->get_edep();
174  nhits++;
175  // esum += hit_iter->second->get_edep();
176  }
177  for (edepiter = layer_edep_map.begin(); edepiter != layer_edep_map.end(); ++edepiter)
178  {
179  nLayers = edepiter->first - 1;
180  if (!strcmp("G4HIT_CEMC", dName.c_str()))
181  mG4EvtTree.cemcactESum[nLayers] = edepiter->second;
182  else if (!strcmp("G4HIT_ABSORBER_CEMC", dName.c_str()))
183  mG4EvtTree.cemcabsESum[nLayers] = edepiter->second;
184  else if (!strcmp("G4HIT_HCAL", dName.c_str()))
185  mG4EvtTree.hcalactESum[nLayers] = edepiter->second;
186  else if (!strcmp("G4HIT_ABSORBER_HCAL", dName.c_str()))
187  mG4EvtTree.hcalabsESum[nLayers] = edepiter->second;
188  }
189  }
190 
191  return nLayers + 1;
192 }