Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTruthVertexing.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTruthVertexing.cc
1 #include "PHTruthVertexing.h"
2 
7 #include <g4main/PHG4VtxPoint.h>
8 
10 
11 #include <phool/PHRandomSeed.h>
12 #include <phool/getClass.h>
13 #include <phool/phool.h>
14 
15 // gsl
16 #include <gsl/gsl_randist.h>
17 #include <gsl/gsl_rng.h>
18 
19 #include <iostream> // for operator<<, basic_ostream
20 #include <set> // for _Rb_tree_iterator, set
21 #include <utility> // for pair
22 
23 class PHCompositeNode;
24 
25 using namespace std;
26 
28  : PHInitVertexing(name)
29  , _g4truth_container(nullptr)
30  , _vertex_error({0.0005, 0.0005, 0.0005})
31  , _embed_only(false)
32 {
33 
34  m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
35 
36  unsigned int seed = PHRandomSeed(); // fixed seed is handled in this funtcion
37  // cout << Name() << " random seed: " << seed << endl;
38  gsl_rng_set(m_RandomGenerator, seed);
39 }
40 
41 
42 
44  gsl_rng_free(m_RandomGenerator);
45 }
46 
48 {
49  int ret = PHInitVertexing::Setup(topNode);
50  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
51 
52  ret = GetNodes(topNode);
53  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
54 
56 }
57 
59 {
60  /*
61  if(Verbosity() > 1)
62  {
63  std::cout << "embed only: " << std::boolalpha << _embed_only << std::endl;
64  }
65  */
66 
67  // we just copy all vertices from the truth container to the SvtxVertexMap
68 
69  auto vrange = _g4truth_container->GetPrimaryVtxRange();
70  for (auto iter = vrange.first; iter != vrange.second; ++iter) // process all primary vertexes
71  {
72  const int point_id = iter->first;
73  int gembed = _g4truth_container->isEmbededVtx(point_id);
74 
75  std::vector<float> pos;
76  pos.clear();
77  pos.assign(3, 0.0);
78 
79  pos[0] = iter->second->get_x();
80  pos[1] = iter->second->get_y();
81  pos[2] = iter->second->get_z();
82 
83  if(Verbosity() > 1)
84  {
85  cout << " PHTruthVertexing: gembed: " << gembed << " nominal position: " << " vx " << pos[0] << " vy " << pos[1] << " vz " << pos[2] << endl;
86  }
87 
88  // We take all primary vertices and copy them to the vertex map, regardless of track embedding ID
89 
90  pos[0] += _vertex_error[0] * gsl_ran_ugaussian(m_RandomGenerator);
91  pos[1] += _vertex_error[1] * gsl_ran_ugaussian(m_RandomGenerator);
92  pos[2] += _vertex_error[2] * gsl_ran_ugaussian(m_RandomGenerator);
93 
94 
95  if (Verbosity() > 0)
96  {
97  cout << __LINE__ << " PHTruthVertexing::Process: point_id " << point_id << " gembed " << gembed << " {" << pos[0]
98  << ", " << pos[1] << ", " << pos[2] << "} +- {"
99  << _vertex_error[0] << ", " << _vertex_error[1] << ", "
100  << _vertex_error[2] << "}" << endl;
101  }
102 
104 
105  vertex->set_x(pos[0]);
106  vertex->set_y(pos[1]);
107  vertex->set_z(pos[2]);
108 
109  for (int j = 0; j < 3; ++j)
110  {
111  for (int i = j; i < 3; ++i)
112  {
113  vertex->set_error(i, j,
114  (i == j ? _vertex_error[i] * _vertex_error[i] : 0));
115  }
116  }
117 
118  vertex->set_id(0);
119  vertex->set_t(0);
120  vertex->set_chisq(0);
121  vertex->set_ndof(1);
122  _vertex_map->insert(vertex);
123 
124  }
125 
126  if (Verbosity() > 0)
128 
130  { assignTracksVertices(topNode); }
131 
133 }
134 
136 {
137  SvtxTrackMap * trackMap = findNode::getClass<SvtxTrackMap>(topNode, _track_map_name.c_str());
138  if(!trackMap)
139  {
140  std::cout << PHWHERE
141  << "Can't find requested track map. Exiting"
142  << std::endl;
143  return;
144  }
145 
146  for(SvtxTrackMap::Iter iter = trackMap->begin();
147  iter != trackMap->end();
148  ++iter)
149  {
150  auto trackKey = iter->first;
151  auto track = iter->second;
152 
153  if(Verbosity() > 3)
154  track->identify();
155 
156  const double trackZ = track->get_z();
157 
158  double dz = 9999.;
159  int trackVertexId = 9999;
160 
161  for(SvtxVertexMap::Iter viter = _vertex_map->begin();
162  viter != _vertex_map->end();
163  ++viter)
164  {
165  auto vertexKey = viter->first;
166  auto vertex = viter->second;
167  if(Verbosity() > 3)
168  vertex->identify();
169 
170  const double vertexZ = vertex->get_z();
171 
172  if( fabs(trackZ - vertexZ) < dz )
173  {
174  dz = fabs(trackZ - vertexZ);
175  trackVertexId = vertexKey;
176  }
177 
178  }
179 
180  track->set_vertex_id(trackVertexId);
181  auto vertex = _vertex_map->find(trackVertexId)->second;
182  vertex->insert_track(trackKey);
183 
184  }
185 
186 
187 }
188 
190 {
191  _g4truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
192  if (!_g4truth_container)
193  {
194  cerr << PHWHERE << " ERROR: Can't find node G4TruthInfo" << endl;
196  }
198 }
199 
201 {
203 }