ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackFastSimEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackFastSimEval.cc
1 #include "TrackFastSimEval.h"
2 
6 #include <trackbase_historic/SvtxVertex.h> // for SvtxVertex
8 
9 #include <g4main/PHG4Hit.h>
11 #include <g4main/PHG4Particle.h>
13 #include <g4main/PHG4VtxPoint.h>
14 
16 #include <fun4all/PHTFileServer.h>
17 #include <fun4all/SubsysReco.h> // for SubsysReco
18 
19 #include <phool/getClass.h>
20 #include <phool/phool.h>
21 
22 #include <TH2.h>
23 #include <TSystem.h>
24 #include <TTree.h>
25 #include <TVector3.h>
26 
27 #include <cassert>
28 #include <cmath>
29 #include <iostream>
30 #include <map> // for _Rb_tree_const_ite...
31 #include <utility> // for pair
32 
33 #define LogError(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
34 #define LogWarning(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
35 
36 using namespace std;
37 
38 //----------------------------------------------------------------------------//
39 //-- Constructor:
40 //-- simple initialization
41 //----------------------------------------------------------------------------//
42 TrackFastSimEval::TrackFastSimEval(const string &name, const string &filename, const string &trackmapname)
43  : SubsysReco(name)
44  , _outfile_name(filename)
45  , _trackmapname(trackmapname)
46  , _event(0)
47  , _flags(NONE)
48  , _eval_tree_tracks(nullptr)
49  , _eval_tree_vertex(nullptr)
50  , _h2d_Delta_mom_vs_truth_mom(nullptr)
51  , _h2d_Delta_mom_vs_truth_eta(nullptr)
52  , _truth_container(nullptr)
53  , _trackmap(nullptr)
54  , _vertexmap(nullptr)
55 {
57 }
58 
59 //----------------------------------------------------------------------------//
60 //-- Init():
61 //-- Intialize all histograms, trees, and ntuples
62 //----------------------------------------------------------------------------//
64 {
65  cout << PHWHERE << " Openning file " << _outfile_name << endl;
66  PHTFileServer::get().open(_outfile_name, "RECREATE");
67 
68  // create TTree
69  _eval_tree_tracks = new TTree("tracks", "FastSim Eval => tracks");
70  _eval_tree_tracks->Branch("event", &event, "event/I");
71  _eval_tree_tracks->Branch("gtrackID", &gtrackID, "gtrackID/I");
72  _eval_tree_tracks->Branch("gflavor", &gflavor, "gflavor/I");
73  _eval_tree_tracks->Branch("gpx", &gpx, "gpx/F");
74  _eval_tree_tracks->Branch("gpy", &gpy, "gpy/F");
75  _eval_tree_tracks->Branch("gpz", &gpz, "gpz/F");
76  _eval_tree_tracks->Branch("gvx", &gvx, "gvx/F");
77  _eval_tree_tracks->Branch("gvy", &gvy, "gvy/F");
78  _eval_tree_tracks->Branch("gvz", &gvz, "gvz/F");
79  _eval_tree_tracks->Branch("gvt", &gvt, "gvt/F");
80  _eval_tree_tracks->Branch("trackID", &trackID, "trackID/I");
81  _eval_tree_tracks->Branch("charge", &charge, "charge/I");
82  _eval_tree_tracks->Branch("nhits", &nhits, "nhits/I");
83  _eval_tree_tracks->Branch("px", &px, "px/F");
84  _eval_tree_tracks->Branch("py", &py, "py/F");
85  _eval_tree_tracks->Branch("pz", &pz, "pz/F");
86  _eval_tree_tracks->Branch("pcax", &pcax, "pcax/F");
87  _eval_tree_tracks->Branch("pcay", &pcay, "pcay/F");
88  _eval_tree_tracks->Branch("pcaz", &pcaz, "pcaz/F");
89  _eval_tree_tracks->Branch("dca2d", &dca2d, "dca2d/F");
90  const string xyz[3] = {"x", "y", "z"};
91  for (map<string, int>::const_iterator iter = m_ProjectionNameMap.begin(); iter != m_ProjectionNameMap.end(); ++iter)
92  {
93  for (int i = 0; i < 3; i++)
94  {
95  string bname = iter->first + "_" + xyz[i];
96  string bdef = bname + "/F";
97  _eval_tree_tracks->Branch(bname.c_str(), &ref[i][iter->second], bdef.c_str());
98  }
99  for (int i = 0; i < 3; i++)
100  {
101  string bname = iter->first + "_p" + xyz[i];
102  string bdef = bname + "/F";
103  _eval_tree_tracks->Branch(bname.c_str(), &ref_p[i][iter->second], bdef.c_str());
104  }
105  for (int i = 0; i < 3; i++)
106  {
107  string bname = iter->first + "_proj_" + xyz[i];
108  string bdef = bname + "/F";
109  _eval_tree_tracks->Branch(bname.c_str(), &proj[i][iter->second], bdef.c_str());
110  }
111  for (int i = 0; i < 3; i++)
112  {
113  string bname = iter->first + "_proj_p" + xyz[i];
114  string bdef = bname + "/F";
115  _eval_tree_tracks->Branch(bname.c_str(), &proj_p[i][iter->second], bdef.c_str());
116  }
117  }
118  _h2d_Delta_mom_vs_truth_eta = new TH2D("_h2d_Delta_mom_vs_truth_eta",
119  "#frac{#Delta p}{truth p} vs. truth #eta", 54, -4.5, +4.5, 1000, -1,
120  1);
121 
122  _h2d_Delta_mom_vs_truth_mom = new TH2D("_h2d_Delta_mom_vs_truth_mom",
123  "#frac{#Delta p}{truth p} vs. truth p", 41, -0.5, 40.5, 1000, -1,
124  1);
125 
126  // create TTree - vertex
127  _eval_tree_vertex = new TTree("vertex", "FastSim Eval => vertces");
128  _eval_tree_vertex->Branch("event", &event, "event/I");
129  _eval_tree_vertex->Branch("gvx", &gvx, "gvx/F");
130  _eval_tree_vertex->Branch("gvy", &gvy, "gvy/F");
131  _eval_tree_vertex->Branch("gvz", &gvz, "gvz/F");
132  _eval_tree_vertex->Branch("gvt", &gvt, "gvt/F");
133  _eval_tree_vertex->Branch("vx", &vx, "vx/F");
134  _eval_tree_vertex->Branch("vy", &vy, "vy/F");
135  _eval_tree_vertex->Branch("vz", &vz, "vz/F");
136  _eval_tree_vertex->Branch("deltavx", &deltavx, "deltavx/F");
137  _eval_tree_vertex->Branch("deltavy", &deltavy, "deltavy/F");
138  _eval_tree_vertex->Branch("deltavz", &deltavz, "deltavz/F");
139  _eval_tree_vertex->Branch("gID", &gtrackID, "gID/I");
140  _eval_tree_vertex->Branch("ID", &trackID, "ID/I");
141  _eval_tree_vertex->Branch("ntracks", &ntracks, "ntracks/I");
142  _eval_tree_vertex->Branch("n_from_truth", &n_from_truth, "n_from_truth/I");
143 
145 }
146 
147 //----------------------------------------------------------------------------//
148 //-- process_event():
149 //-- Call user instructions for every event.
150 //-- This function contains the analysis structure.
151 //----------------------------------------------------------------------------//
153 {
154  _event++;
155  if (Verbosity() >= 2 and _event % 1000 == 0)
156  cout << PHWHERE << "Events processed: " << _event << endl;
157 
158  //std::cout << "Opening nodes" << std::endl;
159  GetNodes(topNode);
160 
161  //std::cout << "Filling trees" << std::endl;
162  fill_track_tree(topNode);
163  fill_vertex_tree(topNode);
164  //std::cout << "DONE" << std::endl;
165 
167 }
168 
169 //----------------------------------------------------------------------------//
170 //-- End():
171 //-- End method, wrap everything up
172 //----------------------------------------------------------------------------//
174 {
176 
177  _eval_tree_tracks->Write();
178  _eval_tree_vertex->Write();
179 
182 
183  //PHTFileServer::get().close();
184 
186 }
187 
188 //----------------------------------------------------------------------------//
189 //-- fill_tree():
190 //-- Fill the trees with truth, track fit, and cluster information
191 //----------------------------------------------------------------------------//
193 {
194  // Make sure to reset all the TTree variables before trying to set them.
195 
196  if (!_truth_container)
197  {
198  LogError("_truth_container not found!");
199  return;
200  }
201 
202  if (!_trackmap)
203  {
204  LogError("_trackmap not found!");
205  return;
206  }
207 
210  //std::cout << "A2" << std::endl;
211  for (PHG4TruthInfoContainer::ConstIterator truth_itr = range.first;
212  truth_itr != range.second; ++truth_itr)
213  {
214  reset_variables();
215  //std::cout << "A1" << std::endl;
216  event = _event;
217 
218  PHG4Particle *g4particle = truth_itr->second;
219  if (!g4particle)
220  {
221  continue;
222  }
223  //std::cout << "B1" << std::endl;
224 
225  SvtxTrack_FastSim *track = nullptr;
226 
227  //std::cout << "TRACKmap size " << _trackmap->size() << std::endl;
228  for (SvtxTrackMap::ConstIter track_itr = _trackmap->begin();
229  track_itr != _trackmap->end();
230  track_itr++)
231  {
232  //std::cout << "TRACK * " << track_itr->first << std::endl;
233  SvtxTrack_FastSim *temp = dynamic_cast<SvtxTrack_FastSim *>(track_itr->second);
234  if (!temp)
235  {
236  std::cout << "ERROR CASTING PARTICLE!" << std::endl;
237  continue;
238  }
239  //std::cout << " PARTICLE!" << std::endl;
240 
241  if ((temp->get_truth_track_id() - g4particle->get_track_id()) == 0)
242  {
243  track = temp;
244  }
245  }
246 
247  //std::cout << "B2" << std::endl;
248  gtrackID = g4particle->get_track_id();
249  gflavor = g4particle->get_pid();
250 
251  gpx = g4particle->get_px();
252  gpy = g4particle->get_py();
253  gpz = g4particle->get_pz();
254 
255  gvx = NAN;
256  gvy = NAN;
257  gvz = NAN;
258  gvt = NAN;
259  PHG4VtxPoint *vtx = _truth_container->GetVtx(g4particle->get_vtx_id());
260  if (vtx)
261  {
262  gvx = vtx->get_x();
263  gvy = vtx->get_y();
264  gvz = vtx->get_z();
265  gvt = vtx->get_t();
266  }
267 
268  if (track)
269  {
270  //std::cout << "C1" << std::endl;
271  trackID = track->get_id();
272  charge = track->get_charge();
273  nhits = track->size_clusters();
274 
275  px = track->get_px();
276  py = track->get_py();
277  pz = track->get_pz();
278  pcax = track->get_x();
279  pcay = track->get_y();
280  pcaz = track->get_z();
281  dca2d = track->get_dca2d();
282 
283  TVector3 truth_mom(gpx, gpy, gpz);
284  TVector3 reco_mom(px, py, pz);
285  //std::cout << "C2" << std::endl;
286 
287  _h2d_Delta_mom_vs_truth_mom->Fill(truth_mom.Mag(), (reco_mom.Mag() - truth_mom.Mag()) / truth_mom.Mag());
288  _h2d_Delta_mom_vs_truth_eta->Fill(truth_mom.Eta(), (reco_mom.Mag() - truth_mom.Mag()) / truth_mom.Mag());
289  // find projections
290  for (int k = 0; k < 3; k++)
291  {
292  for (int j = 0; j < nproj; j++)
293  {
294  proj[k][j] = -9999;
295  proj_p[k][j] = -9999;
296  ref[k][j] = -9999;
297  ref_p[k][j] = -9999;
298  }
299  }
300  for (SvtxTrack::ConstStateIter trkstates = track->begin_states();
301  trkstates != track->end_states();
302  ++trkstates)
303  {
304  // cout << "checking " << trkstates->second->get_name() << endl;
305  map<string, int>::const_iterator iter = m_ProjectionNameMap.find(trkstates->second->get_name());
306  if (iter != m_ProjectionNameMap.end())
307  {
308  // cout << "found " << trkstates->second->get_name() << endl;
309  // setting the projection (xyz and pxpypz)
310  proj[0][iter->second] = trkstates->second->get_x();
311  proj[1][iter->second] = trkstates->second->get_y();
312  proj[2][iter->second] = trkstates->second->get_z();
313  proj_p[0][iter->second] = trkstates->second->get_px();
314  proj_p[1][iter->second] = trkstates->second->get_py();
315  proj_p[2][iter->second] = trkstates->second->get_pz();
316 
317  string nodename = "G4HIT_" + trkstates->second->get_name();
318  PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode, nodename);
319  if (!hits)
320  {
321  cout << "could not find " << nodename << endl;
322  continue;
323  }
324  // cout << "number of hits: " << hits->size() << endl;
325  PHG4HitContainer::ConstRange hit_range = hits->getHits();
326  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
327  {
328  // cout << "checking hit id " << hit_iter->second->get_trkid() << " against " << track->get_truth_track_id() << endl;
329  if (hit_iter->second->get_trkid() - track->get_truth_track_id() == 0)
330  {
331  // cout << "found hit with id " << hit_iter->second->get_trkid() << endl;
332  if (iter->second >= nproj)
333  {
334  cout << "bad index: " << iter->second << endl;
335  gSystem->Exit(1);
336  }
337  ref[0][iter->second] = hit_iter->second->get_x(0);
338  ref[1][iter->second] = hit_iter->second->get_y(0);
339  ref[2][iter->second] = hit_iter->second->get_z(0);
340  ref_p[0][iter->second] = hit_iter->second->get_px(0);
341  ref_p[1][iter->second] = hit_iter->second->get_py(0);
342  ref_p[2][iter->second] = hit_iter->second->get_pz(0);
343  }
344  }
345  }
346  }
347  }
348  //std::cout << "B3" << std::endl;
349 
350  _eval_tree_tracks->Fill();
351  }
352  //std::cout << "A3" << std::endl;
353 
354  return;
355 }
356 
357 //----------------------------------------------------------------------------//
358 //-- fill_tree():
359 //-- Fill the trees with truth, track fit, and cluster information
360 //----------------------------------------------------------------------------//
362 {
363  if (!_truth_container)
364  {
365  LogError("_truth_container not found!");
366  return;
367  }
368 
369  if (!_trackmap)
370  {
371  LogError("_trackmap not found!");
372  return;
373  }
374 
375  if (!_vertexmap)
376  {
377  return;
378  }
379 
380  for (SvtxVertexMap::Iter iter = _vertexmap->begin();
381  iter != _vertexmap->end();
382  ++iter)
383  {
384  SvtxVertex *vertex = iter->second;
385 
386  // Make sure to reset all the TTree variables before trying to set them.
387  reset_variables();
388  //std::cout << "A1" << std::endl;
389  event = _event;
390 
391  if (!vertex)
392  {
393  continue;
394  }
395 
396  //std::cout << "C1" << std::endl;
397  trackID = vertex->get_id();
398  ntracks = vertex->size_tracks();
399 
400  vx = vertex->get_x();
401  vy = vertex->get_y();
402  vz = vertex->get_z();
403  deltavx = sqrt(vertex->get_error(1, 1));
404  deltavy = sqrt(vertex->get_error(2, 2));
405  deltavz = sqrt(vertex->get_error(3, 3));
406 
407  // best matched vertex
408  PHG4VtxPoint *best_vtx = nullptr;
409  int best_n_match = -1;
410  map<PHG4VtxPoint *, int> vertex_match_map;
411  for (auto iter = vertex->begin_tracks(); iter != vertex->end_tracks(); ++iter)
412  {
413  const auto &trackID = *iter;
414  const auto trackIter = _trackmap->find(trackID);
415 
416  if (trackIter == _trackmap->end()) continue;
417 
418  SvtxTrack_FastSim *temp = dynamic_cast<SvtxTrack_FastSim *>(trackIter->second);
419 
420  if (!temp) continue;
421 
422  const auto g4trackID = temp->get_truth_track_id();
423  const PHG4Particle *g4particle = _truth_container->GetParticle(g4trackID);
424  assert(g4particle);
425  PHG4VtxPoint *vtx = _truth_container->GetVtx(g4particle->get_vtx_id());
426 
427  int n_match = ++vertex_match_map[vtx];
428 
429  if (n_match > best_n_match)
430  {
431  best_n_match = n_match;
432  best_vtx = vtx;
433  }
434  }
435  if (best_vtx)
436  {
437  gvx = best_vtx->get_x();
438  gvy = best_vtx->get_y();
439  gvz = best_vtx->get_z();
440  gvt = best_vtx->get_t();
441 
442  n_from_truth = best_n_match;
443  gtrackID = best_vtx->get_id();
444  }
445  _eval_tree_vertex->Fill();
446  }
447  //std::cout << "B3" << std::endl;
448 
449  return;
450 }
451 
452 //----------------------------------------------------------------------------//
453 //-- reset_variables():
454 //-- Reset all the tree variables to their default values.
455 //-- Needs to be called at the start of every event
456 //----------------------------------------------------------------------------//
458 {
459  event = -9999;
460 
461  //-- truth
462  gtrackID = -9999;
463  gflavor = -9999;
464  gpx = NAN;
465  gpy = NAN;
466  gpz = NAN;
467 
468  gvx = NAN;
469  gvy = NAN;
470  gvz = NAN;
471  gvt = NAN;
472 
473  //-- reco
474  trackID = -9999;
475  charge = -9999;
476  nhits = -9999;
477  px = NAN;
478  py = NAN;
479  pz = NAN;
480  pcax = NAN;
481  pcay = NAN;
482  pcaz = NAN;
483  dca2d = NAN;
484 
485  vx = NAN;
486  vy = NAN;
487  vz = NAN;
488  deltavx = NAN;
489  deltavy = NAN;
490  deltavz = NAN;
491  ntracks = -9999;
492  n_from_truth = -9999;
493 }
494 
495 //----------------------------------------------------------------------------//
496 //-- GetNodes():
497 //-- Get all the all the required nodes off the node tree
498 //----------------------------------------------------------------------------//
500 {
501  //DST objects
502  //Truth container
503  _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
504  "G4TruthInfo");
505  if (!_truth_container && _event < 2)
506  {
507  cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree"
508  << endl;
510  }
511 
512  _trackmap = findNode::getClass<SvtxTrackMap>(topNode,
513  _trackmapname);
514  //std::cout << _trackmapname.c_str() << std::endl;
515  if (!_trackmap)
516  {
517  cout << PHWHERE << "SvtxTrackMap node with name "
518  << _trackmapname
519  << " not found on node tree"
520  << endl;
522  }
523 
524  _vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
525  if (!_vertexmap && Verbosity())
526  {
527  cout << PHWHERE << "SvtxTrackMap node with name SvtxVertexMap not found on node tree. Will not build the vertex eval tree"
528  << endl;
529  }
530 
532 }
533 
535 {
536  m_ProjectionNameMap.insert(make_pair(name, m_ProjectionNameMap.size()));
537  return;
538 }