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