Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SnglNtuple.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4SnglNtuple.cc
1 #include "G4SnglNtuple.h"
2 
3 #include <g4main/PHG4Hit.h>
5 #include <g4main/PHG4Particle.h>
7 
9 #include <fun4all/SubsysReco.h> // for SubsysReco
10 
11 #include <phool/getClass.h>
12 
13 #include <TFile.h>
14 #include <TH1.h>
15 #include <TNtuple.h>
16 
17 #include <cmath> // for atan2, sqrt
18 #include <sstream>
19 #include <utility> // for pair
20 
21 using namespace std;
22 
24  : SubsysReco(name)
25  , nblocks(0)
26  , hm(nullptr)
27  , _filename(filename)
28  , ntup(nullptr)
29  , ntup_e(nullptr)
30  , outfile(nullptr)
31 {
32 }
33 
35 {
36  // delete ntup;
37  delete hm;
38 }
39 
41 {
42  hm = new Fun4AllHistoManager(Name());
43  outfile = new TFile(_filename.c_str(), "RECREATE");
44  ntup = new TNtuple("snglntup", "G4Sngls", "phi:theta:e:detid:layer:x0:y0:z0:x1:y1:z1:edep");
45  ntup_e = new TNtuple("sngl_e", "G4SnglEdeps", "phi:theta:e:detid:layer:edep");
46  // ntup->SetDirectory(0);
47  TH1 *h1 = new TH1F("edep1GeV", "edep 0-1GeV", 1000, 0, 1);
48  eloss.push_back(h1);
49  h1 = new TH1F("edep100GeV", "edep 0-100GeV", 1000, 0, 100);
50  eloss.push_back(h1);
51  return 0;
52 }
53 
55 {
56  // get the primary particle which did this to us....
57  PHG4TruthInfoContainer *truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
58  double px = NAN;
59  double py = NAN;
60  double pz = NAN;
61  double e = NAN;
62  double pt = NAN;
63  double phi = NAN;
64  double theta = NAN;
65  if (truthInfoList)
66  {
67  const PHG4TruthInfoContainer::Range primRange = truthInfoList->GetPrimaryParticleRange();
68  px = primRange.first->second->get_px();
69  py = primRange.first->second->get_py();
70  pz = primRange.first->second->get_pz();
71  e = primRange.first->second->get_e();
72  pt = sqrt(px * px + py * py);
73  phi = atan2(py, px);
74  theta = atan2(pt, pz);
75  }
76  ostringstream nodename;
77  set<string>::const_iterator iter;
78  vector<TH1 *>::const_iterator eiter;
79 
80  map<int, double> layer_edep_map;
81  map<int, double>::const_iterator edepiter;
82 
83  for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
84  {
85  layer_edep_map.clear();
86  int detid = (_detid.find(*iter))->second;
87  nodename.str("");
88  nodename << "G4HIT_" << *iter;
89  PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str());
90  if (hits)
91  {
92  double esum = 0;
93  PHG4HitContainer::ConstRange hit_range = hits->getHits();
94  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
95  {
96  layer_edep_map[hit_iter->second->get_layer()] += hit_iter->second->get_edep();
97  esum += hit_iter->second->get_edep();
98  ntup->Fill(phi, theta, e, detid,
99  hit_iter->second->get_detid(),
100  hit_iter->second->get_x(0),
101  hit_iter->second->get_y(0),
102  hit_iter->second->get_z(0),
103  hit_iter->second->get_x(1),
104  hit_iter->second->get_y(1),
105  hit_iter->second->get_z(1),
106  hit_iter->second->get_edep());
107  }
108  for (edepiter = layer_edep_map.begin(); edepiter != layer_edep_map.end(); ++edepiter)
109  {
110  ntup_e->Fill(phi, theta, e, detid, edepiter->first, edepiter->second);
111  }
112  ntup_e->Fill(phi, theta, e, detid, -1, esum); // fill sum over all layers for each detector
113 
114  for (eiter = eloss.begin(); eiter != eloss.end(); ++eiter)
115  {
116  (*eiter)->Fill(esum);
117  }
118  }
119  }
120  return 0;
121 }
122 
124 {
125  outfile->cd();
126  ntup->Write();
127  outfile->Write();
128  outfile->Close();
129  delete outfile;
130  hm->dumpHistos(_filename, "UPDATE");
131  return 0;
132 }
133 
134 void G4SnglNtuple::AddNode(const std::string &name, const int detid)
135 {
136  _node_postfix.insert(name);
137  _detid[name] = detid;
138  return;
139 }