Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HitCountEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HitCountEval.cc
1 
17 #include "HitCountEval.h"
18 
19 #include <phool/getClass.h>
20 #include <fun4all/Fun4AllReturnCodes.h>
21 #include <g4main/PHG4HitContainer.h>
22 #include <g4main/PHG4TruthInfoContainer.h>
23 #include <g4main/PHG4Particle.h>
24 #include <g4main/PHG4Hit.h>
25 #include <fun4all/PHTFileServer.h>
26 
27 #include <TTree.h>
28 #include <TH1F.h>
29 #include <iostream>
30 
31 const char *const HitCountEval::hit_containers[10] {"G4HIT_EGEM_0", "G4HIT_EGEM_1", "G4HIT_EGEM_3",
32  "G4HIT_FGEM_0", "G4HIT_FGEM_1", "G4HIT_FGEM_2", "G4HIT_FGEM_3","G4HIT_FGEM_4", "G4HIT_MAPS",
33  "G4HIT_SVTX"};
34 
35 
36 #define NELEMS(a) (sizeof(a)/sizeof(a[0]))
37 
38 HitCountEval::HitCountEval(const std::string &name, const std::string &file_name):
39  SubsysReco {name},
40  output_file_name {file_name},
41  event {0}
42 {
43 }
44 
45 
46 //----------------------------------------------------------------------------//
47 //-- Init():
48 //-- Intialize all histograms, trees, and ntuples
49 //----------------------------------------------------------------------------//
50 int HitCountEval::Init(PHCompositeNode *topNode) {
51  std::cout << PHWHERE << " Openning file " << output_file_name << std::endl;
52  PHTFileServer::get().open(output_file_name, "RECREATE");
53 
54  // create TTree
55  for (size_t i {0}; i < NELEMS(hit_containers); ++i) {
56  hit_count[i] = new TH1F {(std::string(hit_containers[i]) + "-hits").c_str(),
57  "Hit Count", 100, -4.5, 4.5};
58  eta_count[i] = new TH1F {(std::string(hit_containers[i]) + "-etas").c_str(),
59  "Eta Count", 100, -4.5, 4.5};
60 
61  hits[i] = new TTree {hit_containers[i], ""};
62  hits[i]->Branch("hit_id", &hit_id, "hit_id/l");
63  hits[i]->Branch("x", &hit_x, "x/F");
64  hits[i]->Branch("y", &hit_y, "y/F");
65  hits[i]->Branch("z", &hit_z, "z/F");
66  hits[i]->Branch("px", &hit_px, "px/F");
67  hits[i]->Branch("py", &hit_py, "py/F");
68  hits[i]->Branch("pz", &hit_pz, "pz/F");
69  hits[i]->Branch("particle_initial_px", &particle_initial_px, "particle_initial_px/F");
70  hits[i]->Branch("particle_initial_py", &particle_initial_py, "particle_initial_py/F");
71  hits[i]->Branch("particle_initial_pz", &particle_initial_pz, "particle_initial_pz/F");
72  hits[i]->Branch("layer", &hit_layer, "hit_layer/i");
73  hits[i]->Branch("event", &hit_event, "event/i");
74 
75  normalized_hits[i] = new TTree {(std::string(hit_containers[i]) + "_normalized").c_str(), ""};
76  normalized_hits[i]->Branch("eta", &normalized_hit_eta, "eta/D");
77  normalized_hits[i]->Branch("hit_count", &normalized_hit_count, "eta/D");
78  }
79 
80  return Fun4AllReturnCodes::EVENT_OK;
81 }
82 
83 //----------------------------------------------------------------------------//
84 //-- process_event():
85 //-- Call user instructions for every event.
86 //-- This function contains the analysis structure.
87 //----------------------------------------------------------------------------//
88 int HitCountEval::process_event(PHCompositeNode *const topNode) {
89  event++;
90  if (verbosity >= 2 and event % 1000 == 0)
91  std::cout << PHWHERE << "Events processed: " << event << std::endl;
92 
93  const auto truth = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
94  const auto particle_range = truth->GetPrimaryParticleRange();
95  for (auto it = particle_range.first; it != particle_range.second; ++it) {
96  PHG4Particle const* g4particle{it->second};
97  if(!g4particle) {
98  std::cerr << "Missing PHG4Particle\n";
99  continue;
100  }
101 
102  particle_initial_px = g4particle->get_px();
103  particle_initial_py = g4particle->get_py();
104  particle_initial_pz = g4particle->get_pz();
105 
106  const double mom {std::sqrt(particle_initial_px * particle_initial_px
107  + particle_initial_py * particle_initial_py
108  + particle_initial_pz * particle_initial_pz)};
109  const double eta {0.5 *std::log((mom + particle_initial_pz) / (mom - particle_initial_pz))};
110  for (size_t i {0}; i < NELEMS(hit_containers); ++i) {
111  eta_count[i]->Fill(eta);
112  const auto h = findNode::getClass<PHG4HitContainer>(topNode, hit_containers[i]);
113  if (!h) {
114  std::cerr << "Couldn't find " << hit_containers[i] << '\n';
115  continue;
116  }
117 
118  const auto hit_range = h->getHits();
119  for (auto hp = hit_range.first; hp != hit_range.second; ++hp) {
120  PHG4Hit *const hit {hp->second};
121  for (int j {0}; j < 2; ++j) {
122  hit_count[i]->Fill(eta);
123  hit_id = hit->get_hit_id();
124  hit_x = hit->get_x(j);
125  hit_y = hit->get_y(j);
126  hit_z = hit->get_z(j);
127  hit_px = hit->get_px(j);
128  hit_py = hit->get_py(j);
129  hit_pz = hit->get_pz(j);
130  hits[i]->Fill();
131  }
132  }
133  }
134  }
135 
136  return Fun4AllReturnCodes::EVENT_OK;
137 }
138 
139 int HitCountEval::End(PHCompositeNode *topNode) {
140  PHTFileServer::get().cd(output_file_name);
141 
142  for (size_t i {0}; i < NELEMS(hit_containers); ++i) {
143  const Long64_t nbins {hit_count[i]->GetSize() - 2}; /* -2 for underflow and overflow */
144  for (Long64_t j {0}; j < nbins; ++j) {
145  if (eta_count[i]->GetBinContent(j + 1) != 0) {
146  const Double_t n {hit_count[i]->GetBinContent(j + 1)};
147  const Double_t N {eta_count[i]->GetBinContent(j + 1)};
148 
149  normalized_hit_eta = hit_count[i]->GetBinCenter(j + 1);
150  normalized_hit_count = n / N;
151  normalized_hits[i]->Fill();
152 
153  }
154  normalized_hits[i]->Write();
155  hits[i]->Write();
156  }
157  }
158 
159 
160  return Fun4AllReturnCodes::EVENT_OK;
161 }
162 
163 void HitCountEval::set_filename(const char *const file) {
164  if (file)
165  output_file_name = file;
166 }