ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
getVectors.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file getVectors.cc
1 #include "getVectors.h"
2 
3 #include <fun4all/Fun4AllReturnCodes.h>
4 
5 #include <phool/PHCompositeNode.h>
6 #include <phool/getClass.h>
7 
8 #include <CLHEP/Vector/ThreeVector.h>
9 
10 #include <math.h>
11 #include <float.h>
12 #include <numeric>
13 
14 #include <g4main/PHG4HitContainer.h>
15 //____________________________________________________________________________..
16 getVectors::getVectors(const std::string &name)
17  : SubsysReco(name)
18  , m_truth_info(nullptr)
19  , m_g4particle(nullptr)
20  , m_outfile_name("outputData.root")
21  , m_outfile(nullptr)
22  , m_tree(nullptr)
23 {}
24 
25 //____________________________________________________________________________..
26 getVectors::~getVectors() {}
27 
28 //____________________________________________________________________________..
29 int getVectors::Init(PHCompositeNode *topNode)
30 {
31  assert(topNode);
32  initializeBranches();
33  return Fun4AllReturnCodes::EVENT_OK;
34 }
35 
36 //____________________________________________________________________________..
37 int getVectors::process_event(PHCompositeNode *topNode)
38 {
39  m_truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
40  if (!m_truth_info)
41  {
42  std::cout << "DecayFinder: Missing node G4TruthInfo" << std::endl;
43  return Fun4AllReturnCodes::ABORTEVENT;
44  }
45 
46  CLHEP::Hep3Vector *momentum_hit = NULL;
47  CLHEP::Hep3Vector *the_first_hit = NULL;
48 
49  PHG4HitContainer* hit_container = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MySuperDetector");
50  if (!hit_container)
51  {
52  std::cout << __FILE__ << ": Missing node G4TruthInfo" << std::endl;
53  return Fun4AllReturnCodes::ABORTEVENT;
54  }
55 
56  PHG4TruthInfoContainer::ConstRange range = m_truth_info->GetParticleRange();
57  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
58  iter != range.second; ++iter)
59  {
60  m_g4particle = iter->second;
61  if (m_g4particle->get_parent_id() == 0)
62  {
63  const CLHEP::Hep3Vector momentum_particle( m_g4particle->get_px()
64  , m_g4particle->get_py()
65  , m_g4particle->get_pz());
66 
67  m_pdg_id = m_g4particle->get_pid();
68  m_track_id = m_g4particle->get_track_id();
69  m_barcode = m_g4particle->get_barcode();
70  m_particle_eta = momentum_particle.eta();
71  m_particle_px = momentum_particle.x();
72  m_particle_py = momentum_particle.y();
73  m_particle_pz = momentum_particle.z();
74  m_particle_pt = momentum_particle.perp();
75  m_particle_p = momentum_particle.mag();
76  m_particle_phi = momentum_particle.phi();
77 
78  std::vector<double> delta_phi;
79  //Has to be done this way as truth eval functions are sPHENIX specific
80  PHG4HitContainer::ConstRange hit_range = hit_container->getHits();
81  double shortest_distance = DBL_MAX;
82  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first;
83  hit_iter != hit_range.second; ++hit_iter)
84  {
85  if (hit_iter->second->get_trkid() == m_g4particle->get_track_id())
86  {
87  double hit_distance = sqrt(pow(hit_iter->second->get_x(0), 2) + pow(hit_iter->second->get_y(0), 2) + pow(hit_iter->second->get_z(0), 2));
88  if (hit_distance < shortest_distance)
89  {
90  the_first_hit = new CLHEP::Hep3Vector( hit_iter->second->get_px(0)
91  , hit_iter->second->get_py(0)
92  , hit_iter->second->get_pz(0));
93  shortest_distance = hit_distance;
94  }
95  }
96  }
97  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first;
98  hit_iter != hit_range.second; ++hit_iter)
99  {
100  if (hit_iter->second->get_trkid() == m_g4particle->get_track_id())
101  {
102 
103  momentum_hit = new CLHEP::Hep3Vector( hit_iter->second->get_px(0)
104  , hit_iter->second->get_py(0)
105  , hit_iter->second->get_pz(0));
106 
107  const CLHEP::Hep3Vector momentum_first_hit = (*the_first_hit);
108  double hit_diff = momentum_hit->angle(momentum_first_hit);
109  //double hit_diff = momentum_hit->angle(momentum_particle);
110  delta_phi.push_back(hit_diff);
111  }
112  }
113 
114  double sum = std::accumulate(delta_phi.begin(), delta_phi.end(), 0.0);
115  double mean = sum / delta_phi.size();
116  double square_sum = std::inner_product(delta_phi.begin(), delta_phi.end(), delta_phi.begin(), 0.0);
117  m_delta_phi = sqrt(square_sum / delta_phi.size());
118  m_std_dev = sqrt(square_sum / delta_phi.size() - mean * mean);
119 
120  m_tree->Fill();
121  delta_phi.clear();
122  }
123  }
124 
125  ++m_event_number;
126 
127  return Fun4AllReturnCodes::EVENT_OK;
128 }
129 
130 //____________________________________________________________________________..
131 int getVectors::End(PHCompositeNode *topNode)
132 {
133  m_outfile->Write();
134  m_outfile->Close();
135  delete m_outfile;
136 
137  assert(topNode);
138 
139  return Fun4AllReturnCodes::EVENT_OK;
140 }
141 
142 void getVectors::initializeBranches()
143 {
144  m_outfile = new TFile(m_outfile_name.c_str(), "RECREATE");
145  delete m_tree;
146  m_tree = new TTree("dRICH", "dRICH");
147  m_tree->OptimizeBaskets();
148  m_tree->SetAutoSave(-5e6); // Save the output file every 5MB
149 
150  m_tree->Branch("EventNumber", &m_event_number, "EventNumber/I");
151  m_tree->Branch("particle_PDG_ID", &m_pdg_id, "particle_PDG_ID/I");
152  m_tree->Branch("particle_track_ID", &m_track_id, "particle_track_ID/I");
153  m_tree->Branch("particle_barcode", &m_barcode, "particle_barcode/I");
154  m_tree->Branch("particle_eta", &m_particle_eta, "particle_eta/F");
155  m_tree->Branch("particle_px", &m_particle_px, "particle_px/F");
156  m_tree->Branch("particle_py", &m_particle_py, "particle_py/F");
157  m_tree->Branch("particle_pz", &m_particle_pz, "particle_pz/F");
158  m_tree->Branch("particle_pt", &m_particle_pt, "particle_pt/F");
159  m_tree->Branch("particle_p", &m_particle_p, "particle_p/F");
160  m_tree->Branch("particle_phi", &m_particle_phi, "particle_phi/F");
161  m_tree->Branch("delta_phi", &m_delta_phi, "delta_phi/F");
162  m_tree->Branch("std_dev", &m_std_dev, "std_dev/F");
163 }