EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ActsEvaluator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ActsEvaluator.cc
1 #include "ActsEvaluator.h"
2 
4 #include <fun4all/Fun4AllReturnCodes.h>
5 #include <fun4all/Fun4AllServer.h>
6 #include <g4main/PHG4Hit.h>
7 #include <g4main/PHG4Particle.h>
8 #include <g4main/PHG4TruthInfoContainer.h>
9 #include <phool/PHCompositeNode.h>
10 #include <phool/getClass.h>
11 
12 #include <trackbase/TrkrClusterContainer.h>
14 #include <trackbase/TrkrClusterv2.h>
15 #include <trackbase_historic/SvtxTrack.h>
16 #include <trackbase_historic/SvtxTrackMap.h>
17 
18 #include <g4eval/SvtxClusterEval.h>
19 #include <g4eval/SvtxEvalStack.h>
20 #include <g4eval/SvtxEvaluator.h>
21 #include <g4eval/SvtxTrackEval.h>
22 
23 #include <trackbase_historic/SvtxVertexMap.h>
24 #include <trackbase_historic/SvtxVertex.h>
25 
26 #include <g4main/PHG4VtxPoint.h>
27 
28 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
29 #include <Acts/Utilities/Definitions.hpp>
30 #include <Acts/Utilities/Units.hpp>
31 
32 #include <TFile.h>
33 #include <TTree.h>
34 
35 ActsEvaluator::ActsEvaluator(const std::string &name,
36  SvtxEvaluator *svtxEvaluator)
37  : SubsysReco(name)
38  , m_svtxEvaluator(svtxEvaluator)
39  , m_truthInfo(nullptr)
40  , m_trackMap(nullptr)
41  , m_svtxEvalStack(nullptr)
42  , m_actsTrackKeyMap(nullptr)
43  , m_actsFitResults(nullptr)
44  , m_tGeometry(nullptr)
45  , m_evalCKF(false)
46 {
47 }
48 
49 ActsEvaluator::~ActsEvaluator()
50 {
51 }
52 
53 int ActsEvaluator::Init(PHCompositeNode */*topNode*/)
54 {
55  if (Verbosity() > 1)
56  {
57  std::cout << "Starting ActsEvaluator::Init" << std::endl;
58  }
59 
60  initializeTree();
61 
62  if (Verbosity() > 1)
63  {
64  std::cout << "Finished ActsEvaluator::Init" << std::endl;
65  }
66 
67  return Fun4AllReturnCodes::EVENT_OK;
68 }
69 
70 int ActsEvaluator::process_event(PHCompositeNode *topNode)
71 {
72  if (Verbosity() > 1)
73  {
74  std::cout << "Starting ActsEvaluator at event " << m_eventNr
75  << std::endl;
76  }
77 
78  if (getNodes(topNode) == Fun4AllReturnCodes::ABORTEVENT)
79  return Fun4AllReturnCodes::ABORTEVENT;
80 
81  m_svtxEvalStack = new SvtxEvalStack(topNode);
82  m_svtxEvalStack->next_event(topNode);
83 
84  evaluateTrackFits(topNode);
85 
86  m_eventNr++;
87 
88  if (Verbosity() > 1)
89  std::cout << "Finished ActsEvaluator::process_event" << std::endl;
90 
91  return Fun4AllReturnCodes::EVENT_OK;
92 }
93 
94 
95 void ActsEvaluator::evaluateTrackFits(PHCompositeNode *topNode)
96 {
97  if(Verbosity() > 5)
98  std::cout << "Evaluating Acts track fits" << std::endl;
99 
100  SvtxTrackEval *trackeval = m_svtxEvalStack->get_track_eval();
101 
102  std::map<const unsigned int, Trajectory>::iterator trajIter;
103  int iTraj = 0;
104  int iTrack = 0;
105 
106  for (trajIter = m_actsFitResults->begin();
107  trajIter != m_actsFitResults->end();
108  ++trajIter)
109  {
111  unsigned int trackKey = trajIter->first;
112  Trajectory traj = trajIter->second;
113 
114  if(Verbosity() > 2)
115  std::cout << "Starting trajectory " << iTraj
116  << " with trackKey corresponding to track seed "
117  << trackKey << std::endl;
118 
119  const auto &[trackTips, mj] = traj.trajectory();
120  m_trajNr = iTraj;
121 
123  if (trackTips.empty())
124  {
125  if (Verbosity() > 1)
126  std::cout << "TrackTips empty in ActsEvaluator" << std::endl;
127  continue;
128  }
129 
130  iTrack = 0;
131 
134 
135  SvtxTrack* actsProtoTrack = m_actsProtoTrackMap->find(trackKey)->second;
137  std::map<const size_t, const unsigned int> trackKeyMap;
138  if(m_evalCKF)
139  trackKeyMap = m_actsTrackKeyMap->find(trackKey)->second;
140 
142  for(const size_t &trackTip : trackTips)
143  {
144  if(Verbosity() > 2)
145  std::cout << "beginning trackTip " << trackTip
146  << " corresponding to track " << iTrack
147  << " in trajectroy " << iTraj << std::endl;
148 
152  if(m_evalCKF)
153  {
154  std::map<const size_t, const unsigned int>::iterator ckfiter
155  = trackKeyMap.find(trackTip);
156  if(ckfiter == trackKeyMap.end())
157  {
158  if(Verbosity() > 2)
159  std::cout << "Track result flagged as bad, continuing"
160  << std::endl;
161  continue;
162  }
163  trackKey = ckfiter->second;
164  }
165 
166  if(Verbosity() > 2)
167  std::cout<<"Evaluating track key " << trackKey
168  << " for track tip " << trackTip << std::endl;
169 
170  SvtxTrack *track = m_trackMap->find(trackKey)->second;
171  PHG4Particle *g4particle = trackeval->max_truth_particle_by_nclusters(track);
172  unsigned int vertexId = track->get_vertex_id();
173  if(vertexId == UINT_MAX)
174  vertexId = 0;
175  const SvtxVertex *svtxVertex = m_vertexMap->get(vertexId);
176  Acts::Vector3D vertex;
177  vertex(0) = svtxVertex->get_x() * Acts::UnitConstants::cm;
178  vertex(1) = svtxVertex->get_y() * Acts::UnitConstants::cm;
179  vertex(2) = svtxVertex->get_z() * Acts::UnitConstants::cm;
180 
181  if(Verbosity() > 1)
182  {
183  std::cout << "Analyzing SvtxTrack "<< trackKey << std::endl;
184 
185  std::cout << "TruthParticle : " << g4particle->get_px()
186  << ", " << g4particle->get_py() << ", "
187  << g4particle->get_pz() << ", "<< g4particle->get_e()
188  << std::endl;
189  }
190 
191  m_trackNr = iTrack;
192 
193 
194  auto trajState =
195  Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
196 
197  if(Verbosity() > 1)
198  {
199  if(traj.hasTrackParameters(trackTip))
200  {
201  std::cout << "Fitted params : "
202  << traj.trackParameters(trackTip).position(m_tGeometry->geoContext)
203  << std::endl << traj.trackParameters(trackTip).momentum()
204  << std::endl;
205  std::cout << "Track has " << trajState.nMeasurements
206  << " measurements and " << trajState.nHoles
207  << " holes and " << trajState.nOutliers
208  << " outliers and " << trajState.nStates
209  << " states " << std::endl;
210  }
211  }
212 
213  m_nMeasurements = trajState.nMeasurements;
214  m_nStates = trajState.nStates;
215  m_nOutliers = trajState.nOutliers;
216  m_nHoles = trajState.nHoles;
217  m_chi2_fit = trajState.chi2Sum;
218  m_ndf_fit = trajState.NDF;
219  m_quality = track->get_quality();
220 
221  fillG4Particle(g4particle);
222  fillProtoTrack(actsProtoTrack, topNode);
223  fillFittedTrackParams(traj, trackTip, vertex);
224  visitTrackStates(traj,trackTip, topNode);
225 
226  m_trackTree->Fill();
227 
229  clearTrackVariables();
230  if(Verbosity() > 1)
231  std::cout << "Finished track " << iTrack <<std::endl;
232 
233  iTrack++;
234  }
235 
236 
237  if(Verbosity() > 1)
238  {
239  std::cout << "Analyzed " << iTrack << " tracks in trajectory number "
240  << iTraj << std::endl;
241  }
242 
243  ++iTraj;
244  }
245 
246  if(Verbosity () > 5)
247  std::cout << "Finished evaluating track fits" << std::endl;
248 
249  return;
250 }
251 
252 
253 int ActsEvaluator::End(PHCompositeNode */*topNode*/)
254 {
255  m_trackFile->cd();
256  m_trackTree->Write();
257  m_trackFile->Close();
258 
259  return Fun4AllReturnCodes::EVENT_OK;
260 }
261 
262 int ActsEvaluator::ResetEvent(PHCompositeNode */*topNode*/)
263 {
264  m_trajNr = 0;
265  return Fun4AllReturnCodes::EVENT_OK;
266 }
267 
268 
269 
270 void ActsEvaluator::visitTrackStates(const Trajectory traj,
271  const size_t &trackTip,
272  PHCompositeNode *topNode)
273 {
274 
275  if(Verbosity() > 2)
276  std::cout << "Begin visit track states" << std::endl;
277 
278  const auto &[trackTips, mj] = traj.trajectory();
279 
280  mj.visitBackwards(trackTip, [&](const auto &state) {
282  auto typeFlags = state.typeFlags();
283  if (not typeFlags.test(Acts::TrackStateFlag::MeasurementFlag))
284  {
285  return true;
286  }
287 
289  auto geoID = state.referenceSurface().geometryId();
290  m_volumeID.push_back(geoID.volume());
291  m_layerID.push_back(geoID.layer());
292  m_moduleID.push_back(geoID.sensitive());
293 
294  if(Verbosity() > 3)
295  std::cout << "Cluster volume : layer : sensitive " << geoID.volume()
296  << " : " << geoID.layer() << " : "
297  << geoID.sensitive() << std::endl;
298 
299  auto meas = std::get<Measurement>(*state.uncalibrated());
300 
302  Acts::Vector2D local(meas.parameters()[Acts::eBoundLoc0],
303  meas.parameters()[Acts::eBoundLoc1]);
304 
308  Acts::Vector3D mom(1, 1, 1);
309  Acts::Vector3D global = meas.referenceObject().localToGlobal(m_tGeometry->geoContext,
310  local, mom);
311 
313  auto cov = meas.covariance();
314 
315  m_lx_hit.push_back(local.x());
316  m_ly_hit.push_back(local.y());
317  m_x_hit.push_back(global.x());
318  m_y_hit.push_back(global.y());
319  m_z_hit.push_back(global.z());
320 
324  TrkrDefs::cluskey cluskey = state.uncalibrated().cluskey();
325  float gt = -9999;
326  Acts::Vector3D globalTruthPos = getGlobalTruthHit(topNode, cluskey, gt);
327  float gx = globalTruthPos(0);
328  float gy = globalTruthPos(1);
329  float gz = globalTruthPos(2);
330 
332  const float r = sqrt(gx * gx + gy * gy + gz * gz);
333  Acts::Vector3D globalTruthUnitDir(gx / r, gy / r, gz / r);
334 
335  auto vecResult = meas.referenceObject().globalToLocal(
336  m_tGeometry->geoContext,
337  globalTruthPos,
338  globalTruthUnitDir);
339 
341  m_t_x.push_back(gx);
342  m_t_y.push_back(gy);
343  m_t_z.push_back(gz);
344  m_t_r.push_back(sqrt(gx * gx + gy * gy));
345  m_t_dx.push_back(gx / r);
346  m_t_dy.push_back(gy / r);
347  m_t_dz.push_back(gz / r);
348 
350  float truthLOC0 = 0;
351  float truthLOC1 = 0;
352  float truthPHI = 0;
353  float truthTHETA = 0;
354  float truthQOP = 0;
355  float truthTIME = 0;
356  float momentum = sqrt(m_t_px * m_t_px +
357  m_t_py * m_t_py +
358  m_t_pz * m_t_pz);
359 
360  if(vecResult.ok())
361  {
362  Acts::Vector2D truthLocVec = vecResult.value();
363  truthLOC0 = truthLocVec.x();
364  truthLOC1 = truthLocVec.y();
365  }
366  else
367  {
368  truthLOC0 = -9999.;
369  truthLOC1 = -9999.;
370  }
371 
372  truthPHI = phi(globalTruthUnitDir);
373  truthTHETA = theta(globalTruthUnitDir);
374  truthQOP =
375  m_t_charge / momentum;
376  truthTIME = gt;
377 
378  m_t_eLOC0.push_back(truthLOC0);
379  m_t_eLOC1.push_back(truthLOC1);
380  m_t_ePHI.push_back(truthPHI);
381  m_t_eTHETA.push_back(truthTHETA);
382  m_t_eQOP.push_back(truthQOP);
383  m_t_eT.push_back(truthTIME);
384 
386  bool predicted = false;
387  if (state.hasPredicted())
388  {
389  predicted = true;
390  m_nPredicted++;
391 
392  auto parameters = state.predicted();
393  auto covariance = state.predictedCovariance();
394 
396  auto H = meas.projector();
397  auto resCov = cov + H * covariance * H.transpose();
398  auto residual = meas.residual(parameters);
399  m_res_x_hit.push_back(residual(Acts::eBoundLoc0));
400  m_res_y_hit.push_back(residual(Acts::eBoundLoc1));
401  m_err_x_hit.push_back(
402  sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
403  m_err_y_hit.push_back(
404  sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
405  m_pull_x_hit.push_back(
406  residual(Acts::eBoundLoc0) /
407  sqrt(resCov(Acts::eBoundLoc0, Acts::eBoundLoc0)));
408  m_pull_y_hit.push_back(
409  residual(Acts::eBoundLoc1) /
410  sqrt(resCov(Acts::eBoundLoc1, Acts::eBoundLoc1)));
411  m_dim_hit.push_back(state.calibratedSize());
412 
414  m_eLOC0_prt.push_back(parameters[Acts::eBoundLoc0]);
415  m_eLOC1_prt.push_back(parameters[Acts::eBoundLoc1]);
416  m_ePHI_prt.push_back(parameters[Acts::eBoundPhi]);
417  m_eTHETA_prt.push_back(parameters[Acts::eBoundTheta]);
418  m_eQOP_prt.push_back(parameters[Acts::eBoundQOverP]);
419  m_eT_prt.push_back(parameters[Acts::eBoundTime]);
420 
422  m_res_eLOC0_prt.push_back(parameters[Acts::eBoundLoc0] -
423  truthLOC0);
424  m_res_eLOC1_prt.push_back(parameters[Acts::eBoundLoc1] -
425  truthLOC1);
426  m_res_ePHI_prt.push_back(parameters[Acts::eBoundPhi] -
427  truthPHI);
428  m_res_eTHETA_prt.push_back(
429  parameters[Acts::eBoundTheta] - truthTHETA);
430  m_res_eQOP_prt.push_back(parameters[Acts::eBoundQOverP] -
431  truthQOP);
432  m_res_eT_prt.push_back(parameters[Acts::eBoundTime] -
433  truthTIME);
434 
436  m_err_eLOC0_prt.push_back(
437  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
438  m_err_eLOC1_prt.push_back(
439  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
440  m_err_ePHI_prt.push_back(
441  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
442  m_err_eTHETA_prt.push_back(
443  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
444  m_err_eQOP_prt.push_back(
445  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
446  m_err_eT_prt.push_back(
447  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
448 
450  m_pull_eLOC0_prt.push_back(
451  (parameters[Acts::eBoundLoc0] - truthLOC0) /
452  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
453  m_pull_eLOC1_prt.push_back(
454  (parameters[Acts::eBoundLoc1] - truthLOC1) /
455  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
456  m_pull_ePHI_prt.push_back(
457  (parameters[Acts::eBoundPhi] - truthPHI) /
458  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
459  m_pull_eTHETA_prt.push_back(
460  (parameters[Acts::eBoundTheta] - truthTHETA) /
461  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
462  m_pull_eQOP_prt.push_back(
463  (parameters[Acts::eBoundQOverP] - truthQOP) /
464  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
465  m_pull_eT_prt.push_back(
466  (parameters[Acts::eBoundTime] - truthTIME) /
467  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
468 
469  Acts::FreeVector freeParams =
470  Acts::detail::transformBoundToFreeParameters(state.referenceSurface(),
471  m_tGeometry->geoContext,
472  parameters);
473 
474  m_x_prt.push_back(freeParams[Acts::eFreePos0]);
475  m_y_prt.push_back(freeParams[Acts::eFreePos1]);
476  m_z_prt.push_back(freeParams[Acts::eFreePos2]);
477  auto p = std::abs(1 / freeParams[Acts::eFreeQOverP]);
478  m_px_prt.push_back(p * freeParams[Acts::eFreeDir0]);
479  m_py_prt.push_back(p * freeParams[Acts::eFreeDir1]);
480  m_pz_prt.push_back(p * freeParams[Acts::eFreeDir2]);
481  m_pT_prt.push_back(p * std::hypot(freeParams[Acts::eFreeDir0],
482  freeParams[Acts::eFreeDir1]));
483  m_eta_prt.push_back(
484  Acts::VectorHelpers::eta(freeParams.segment<3>(Acts::eFreeDir0)));
485  }
486  else
487  {
489  m_res_x_hit.push_back(-9999);
490  m_res_y_hit.push_back(-9999);
491  m_err_x_hit.push_back(-9999);
492  m_err_y_hit.push_back(-9999);
493  m_pull_x_hit.push_back(-9999);
494  m_pull_y_hit.push_back(-9999);
495  m_dim_hit.push_back(-9999);
496  m_eLOC0_prt.push_back(-9999);
497  m_eLOC1_prt.push_back(-9999);
498  m_ePHI_prt.push_back(-9999);
499  m_eTHETA_prt.push_back(-9999);
500  m_eQOP_prt.push_back(-9999);
501  m_eT_prt.push_back(-9999);
502  m_res_eLOC0_prt.push_back(-9999);
503  m_res_eLOC1_prt.push_back(-9999);
504  m_res_ePHI_prt.push_back(-9999);
505  m_res_eTHETA_prt.push_back(-9999);
506  m_res_eQOP_prt.push_back(-9999);
507  m_res_eT_prt.push_back(-9999);
508  m_err_eLOC0_prt.push_back(-9999);
509  m_err_eLOC1_prt.push_back(-9999);
510  m_err_ePHI_prt.push_back(-9999);
511  m_err_eTHETA_prt.push_back(-9999);
512  m_err_eQOP_prt.push_back(-9999);
513  m_err_eT_prt.push_back(-9999);
514  m_pull_eLOC0_prt.push_back(-9999);
515  m_pull_eLOC1_prt.push_back(-9999);
516  m_pull_ePHI_prt.push_back(-9999);
517  m_pull_eTHETA_prt.push_back(-9999);
518  m_pull_eQOP_prt.push_back(-9999);
519  m_pull_eT_prt.push_back(-9999);
520  m_x_prt.push_back(-9999);
521  m_y_prt.push_back(-9999);
522  m_z_prt.push_back(-9999);
523  m_px_prt.push_back(-9999);
524  m_py_prt.push_back(-9999);
525  m_pz_prt.push_back(-9999);
526  m_pT_prt.push_back(-9999);
527  m_eta_prt.push_back(-9999);
528  }
529 
530  bool filtered = false;
531  if (state.hasFiltered())
532  {
533  filtered = true;
534  m_nFiltered++;
535 
536  Acts::BoundTrackParameters parameter(
537  state.referenceSurface().getSharedPtr(),
538  state.filtered(),
539  state.filteredCovariance());
540 
541  auto covariance = state.filteredCovariance();
542 
543  m_eLOC0_flt.push_back(parameter.parameters()[Acts::eBoundLoc0]);
544  m_eLOC1_flt.push_back(parameter.parameters()[Acts::eBoundLoc1]);
545  m_ePHI_flt.push_back(parameter.parameters()[Acts::eBoundPhi]);
546  m_eTHETA_flt.push_back(parameter.parameters()[Acts::eBoundTheta]);
547  m_eQOP_flt.push_back(parameter.parameters()[Acts::eBoundQOverP]);
548  m_eT_flt.push_back(parameter.parameters()[Acts::eBoundTime]);
549 
550  m_res_eLOC0_flt.push_back(parameter.parameters()[Acts::eBoundLoc0] -
551  truthLOC0);
552  m_res_eLOC1_flt.push_back(parameter.parameters()[Acts::eBoundLoc1] -
553  truthLOC1);
554  m_res_ePHI_flt.push_back(parameter.parameters()[Acts::eBoundPhi] -
555  truthPHI);
556  m_res_eTHETA_flt.push_back(parameter.parameters()[Acts::eBoundTheta] -
557  truthTHETA);
558  m_res_eQOP_flt.push_back(parameter.parameters()[Acts::eBoundQOverP] -
559  truthQOP);
560  m_res_eT_flt.push_back(parameter.parameters()[Acts::eBoundTime] -
561  truthTIME);
562 
564  m_err_eLOC0_flt.push_back(
565  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
566  m_err_eLOC1_flt.push_back(
567  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
568  m_err_ePHI_flt.push_back(
569  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
570  m_err_eTHETA_flt.push_back(
571  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
572  m_err_eQOP_flt.push_back(
573  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
574  m_err_eT_flt.push_back(
575  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
576 
578  m_pull_eLOC0_flt.push_back(
579  (parameter.parameters()[Acts::eBoundLoc0] - truthLOC0) /
580  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
581  m_pull_eLOC1_flt.push_back(
582  (parameter.parameters()[Acts::eBoundLoc1] - truthLOC1) /
583  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
584  m_pull_ePHI_flt.push_back(
585  (parameter.parameters()[Acts::eBoundPhi] - truthPHI) /
586  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
587  m_pull_eTHETA_flt.push_back(
588  (parameter.parameters()[Acts::eBoundTheta] - truthTHETA) /
589  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
590  m_pull_eQOP_flt.push_back(
591  (parameter.parameters()[Acts::eBoundQOverP] - truthQOP) /
592  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
593  m_pull_eT_flt.push_back(
594  (parameter.parameters()[Acts::eBoundTime] - truthTIME) /
595  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
596 
598  m_x_flt.push_back(parameter.position(m_tGeometry->geoContext).x());
599  m_y_flt.push_back(parameter.position(m_tGeometry->geoContext).y());
600  m_z_flt.push_back(parameter.position(m_tGeometry->geoContext).z());
601  m_px_flt.push_back(parameter.momentum().x());
602  m_py_flt.push_back(parameter.momentum().y());
603  m_pz_flt.push_back(parameter.momentum().z());
604  m_pT_flt.push_back(sqrt(parameter.momentum().x() * parameter.momentum().x()
605  + parameter.momentum().y() * parameter.momentum().y()));
606  m_eta_flt.push_back(eta(parameter.position(m_tGeometry->geoContext)));
607  m_chi2.push_back(state.chi2());
608 
609  }
610  else
611  {
613  m_eLOC0_flt.push_back(-9999);
614  m_eLOC1_flt.push_back(-9999);
615  m_ePHI_flt.push_back(-9999);
616  m_eTHETA_flt.push_back(-9999);
617  m_eQOP_flt.push_back(-9999);
618  m_eT_flt.push_back(-9999);
619  m_res_eLOC0_flt.push_back(-9999);
620  m_res_eLOC1_flt.push_back(-9999);
621  m_res_ePHI_flt.push_back(-9999);
622  m_res_eTHETA_flt.push_back(-9999);
623  m_res_eQOP_flt.push_back(-9999);
624  m_res_eT_flt.push_back(-9999);
625  m_err_eLOC0_flt.push_back(-9999);
626  m_err_eLOC1_flt.push_back(-9999);
627  m_err_ePHI_flt.push_back(-9999);
628  m_err_eTHETA_flt.push_back(-9999);
629  m_err_eQOP_flt.push_back(-9999);
630  m_err_eT_flt.push_back(-9999);
631  m_pull_eLOC0_flt.push_back(-9999);
632  m_pull_eLOC1_flt.push_back(-9999);
633  m_pull_ePHI_flt.push_back(-9999);
634  m_pull_eTHETA_flt.push_back(-9999);
635  m_pull_eQOP_flt.push_back(-9999);
636  m_pull_eT_flt.push_back(-9999);
637  m_x_flt.push_back(-9999);
638  m_y_flt.push_back(-9999);
639  m_z_flt.push_back(-9999);
640  m_py_flt.push_back(-9999);
641  m_pz_flt.push_back(-9999);
642  m_pT_flt.push_back(-9999);
643  m_eta_flt.push_back(-9999);
644  m_chi2.push_back(-9999);
645  }
646 
647  bool smoothed = false;
648  if (state.hasSmoothed())
649  {
650  smoothed = true;
651  m_nSmoothed++;
652 
653  Acts::BoundTrackParameters parameter(
654  state.referenceSurface().getSharedPtr(),
655  state.smoothed(),
656  state.smoothedCovariance());
657 
658  auto covariance = state.smoothedCovariance();
659 
660  m_eLOC0_smt.push_back(parameter.parameters()[Acts::eBoundLoc0]);
661  m_eLOC1_smt.push_back(parameter.parameters()[Acts::eBoundLoc1]);
662  m_ePHI_smt.push_back(parameter.parameters()[Acts::eBoundPhi]);
663  m_eTHETA_smt.push_back(parameter.parameters()[Acts::eBoundTheta]);
664  m_eQOP_smt.push_back(parameter.parameters()[Acts::eBoundQOverP]);
665  m_eT_smt.push_back(parameter.parameters()[Acts::eBoundTime]);
666 
667  m_res_eLOC0_smt.push_back(parameter.parameters()[Acts::eBoundLoc0] -
668  truthLOC0);
669  m_res_eLOC1_smt.push_back(parameter.parameters()[Acts::eBoundLoc1] -
670  truthLOC1);
671  m_res_ePHI_smt.push_back(parameter.parameters()[Acts::eBoundPhi] -
672  truthPHI);
673  m_res_eTHETA_smt.push_back(parameter.parameters()[Acts::eBoundTheta] -
674  truthTHETA);
675  m_res_eQOP_smt.push_back(parameter.parameters()[Acts::eBoundQOverP] -
676  truthQOP);
677  m_res_eT_smt.push_back(parameter.parameters()[Acts::eBoundTime] -
678  truthTIME);
679 
681  m_err_eLOC0_smt.push_back(
682  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
683  m_err_eLOC1_smt.push_back(
684  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
685  m_err_ePHI_smt.push_back(
686  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
687  m_err_eTHETA_smt.push_back(
688  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
689  m_err_eQOP_smt.push_back(
690  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
691  m_err_eT_smt.push_back(
692  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
693 
695  m_pull_eLOC0_smt.push_back(
696  (parameter.parameters()[Acts::eBoundLoc0] - truthLOC0) /
697  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
698  m_pull_eLOC1_smt.push_back(
699  (parameter.parameters()[Acts::eBoundLoc1] - truthLOC1) /
700  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
701  m_pull_ePHI_smt.push_back(
702  (parameter.parameters()[Acts::eBoundPhi] - truthPHI) /
703  sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
704  m_pull_eTHETA_smt.push_back(
705  (parameter.parameters()[Acts::eBoundTheta] - truthTHETA) /
706  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
707  m_pull_eQOP_smt.push_back(
708  (parameter.parameters()[Acts::eBoundQOverP] - truthQOP) /
709  sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
710  m_pull_eT_smt.push_back(
711  (parameter.parameters()[Acts::eBoundTime] - truthTIME) /
712  sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
713 
714  m_x_smt.push_back(parameter.position(m_tGeometry->geoContext).x());
715  m_y_smt.push_back(parameter.position(m_tGeometry->geoContext).y());
716  m_z_smt.push_back(parameter.position(m_tGeometry->geoContext).z());
717  m_px_smt.push_back(parameter.momentum().x());
718  m_py_smt.push_back(parameter.momentum().y());
719  m_pz_smt.push_back(parameter.momentum().z());
720  m_pT_smt.push_back(sqrt(parameter.momentum().x() * parameter.momentum().x()
721  + parameter.momentum().y() * parameter.momentum().y()));
722  m_eta_smt.push_back(eta(parameter.position(m_tGeometry->geoContext)));
723  }
724  else
725  {
727  m_eLOC0_smt.push_back(-9999);
728  m_eLOC1_smt.push_back(-9999);
729  m_ePHI_smt.push_back(-9999);
730  m_eTHETA_smt.push_back(-9999);
731  m_eQOP_smt.push_back(-9999);
732  m_eT_smt.push_back(-9999);
733  m_res_eLOC0_smt.push_back(-9999);
734  m_res_eLOC1_smt.push_back(-9999);
735  m_res_ePHI_smt.push_back(-9999);
736  m_res_eTHETA_smt.push_back(-9999);
737  m_res_eQOP_smt.push_back(-9999);
738  m_res_eT_smt.push_back(-9999);
739  m_err_eLOC0_smt.push_back(-9999);
740  m_err_eLOC1_smt.push_back(-9999);
741  m_err_ePHI_smt.push_back(-9999);
742  m_err_eTHETA_smt.push_back(-9999);
743  m_err_eQOP_smt.push_back(-9999);
744  m_err_eT_smt.push_back(-9999);
745  m_pull_eLOC0_smt.push_back(-9999);
746  m_pull_eLOC1_smt.push_back(-9999);
747  m_pull_ePHI_smt.push_back(-9999);
748  m_pull_eTHETA_smt.push_back(-9999);
749  m_pull_eQOP_smt.push_back(-9999);
750  m_pull_eT_smt.push_back(-9999);
751  m_x_smt.push_back(-9999);
752  m_y_smt.push_back(-9999);
753  m_z_smt.push_back(-9999);
754  m_px_smt.push_back(-9999);
755  m_py_smt.push_back(-9999);
756  m_pz_smt.push_back(-9999);
757  m_pT_smt.push_back(-9999);
758  m_eta_smt.push_back(-9999);
759  }
760 
762  m_prt.push_back(predicted);
763  m_flt.push_back(filtered);
764  m_smt.push_back(smoothed);
765 
766  return true;
767  }
768  );
769 
770 
771  if(Verbosity() > 2)
772  std::cout << "Finished track states" << std::endl;
773 
774  return;
775 }
776 
777 
778 
779 Acts::Vector3D ActsEvaluator::getGlobalTruthHit(PHCompositeNode */*topNode*/,
780  TrkrDefs::cluskey cluskey,
781  float &_gt)
782 {
783  SvtxClusterEval *clustereval = m_svtxEvalStack->get_cluster_eval();
784 
785  std::shared_ptr<TrkrCluster> truth_cluster = clustereval->max_truth_cluster_by_energy(cluskey);
786 
787  float gx = -9999;
788  float gy = -9999;
789  float gz = -9999;
790  float gt = -9999;
791 
792  if (truth_cluster)
793  {
794  gx = truth_cluster->getX();
795  gy = truth_cluster->getY();
796  gz = truth_cluster->getZ();
797  }
798 
800  gx *= Acts::UnitConstants::cm;
801  gy *= Acts::UnitConstants::cm;
802  gz *= Acts::UnitConstants::cm;
803 
804  Acts::Vector3D globalPos(gx, gy, gz);
805  _gt = gt;
806  return globalPos;
807 
808 }
809 
810 
811 //___________________________________________________________________________________
812 Surface ActsEvaluator::getSurface(TrkrDefs::cluskey cluskey, TrkrDefs::subsurfkey surfkey)
813 {
814  const auto trkrid = TrkrDefs::getTrkrId(cluskey);
815  const auto hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluskey);
816 
817  switch( trkrid )
818  {
819  case TrkrDefs::TrkrId::micromegasId: return getMMSurface( hitsetkey );
820  case TrkrDefs::TrkrId::tpcId: return getTpcSurface(hitsetkey, surfkey);
821  case TrkrDefs::TrkrId::mvtxId:
822  case TrkrDefs::TrkrId::inttId:
823  {
824  return getSiliconSurface(hitsetkey);
825  }
826  }
827 
828  // unreachable
829  return nullptr;
830 
831 }
832 
833 //___________________________________________________________________________________
834 Surface ActsEvaluator::getSiliconSurface(TrkrDefs::hitsetkey hitsetkey)
835 {
836  auto surfMap = m_surfMaps->siliconSurfaceMap;
837  auto iter = surfMap.find(hitsetkey);
838  if(iter != surfMap.end())
839  {
840  return iter->second;
841  }
842 
844  return nullptr;
845 
846 }
847 
848 //___________________________________________________________________________________
849 Surface ActsEvaluator::getTpcSurface(TrkrDefs::hitsetkey hitsetkey, TrkrDefs::subsurfkey surfkey)
850 {
851  unsigned int layer = TrkrDefs::getLayer(hitsetkey);
852  const auto iter = m_surfMaps->tpcSurfaceMap.find(layer);
853  if(iter != m_surfMaps->tpcSurfaceMap.end())
854  {
855  auto surfvec = iter->second;
856  return surfvec.at(surfkey);
857  }
858 
860  return nullptr;
861 }
862 
863 //___________________________________________________________________________________
864 Surface ActsEvaluator::getMMSurface(TrkrDefs::hitsetkey hitsetkey)
865 {
866  const auto iter = m_surfMaps->mmSurfaceMap.find( hitsetkey );
867  return (iter == m_surfMaps->mmSurfaceMap.end()) ? nullptr:iter->second;
868 }
869 
870 
871 void ActsEvaluator::fillProtoTrack(SvtxTrack* track, PHCompositeNode *topNode)
872 {
873 
874  if(Verbosity() > 2)
875  std::cout << "Filling proto track seed quantities" << std::endl;
876 
877  Acts::Vector3D position(track->get_x() * 10, track->get_y() * 10,
878  track->get_z() * 10);
879 
880  Acts::Vector3D momentum(track->get_px(),
881  track->get_py(),
882  track->get_pz());
883 
884  m_protoTrackPx = momentum(0);
885  m_protoTrackPy = momentum(1);
886  m_protoTrackPz = momentum(2);
887  m_protoTrackX = position(0);
888  m_protoTrackY = position(1);
889  m_protoTrackZ = position(2);
890 
891  for(SvtxTrack::ConstClusterKeyIter clusIter = track->begin_cluster_keys();
892  clusIter != track->end_cluster_keys();
893  ++clusIter)
894  {
895  auto key = *clusIter;
896  auto cluster = m_clusterContainer->findCluster(key);
897 
899  Acts::Vector2D loc(cluster->getLocalX() * 10,
900  cluster->getLocalY() * 10);
901 
902  Acts::Vector3D mom(0,0,0);
903  Acts::Vector3D globalPos(cluster->getX() * 10,
904  cluster->getY() * 10,
905  cluster->getZ() * 10);
906 
907  m_SLx.push_back(globalPos(0));
908  m_SLy.push_back(globalPos(1));
909  m_SLz.push_back(globalPos(2));
910  m_SL_lx.push_back(loc(0));
911  m_SL_ly.push_back(loc(1));
912 
914  float gt = -9999;
915 
916  Acts::Vector3D globalTruthPos = getGlobalTruthHit(topNode, key, gt);
917 
918  float gx = globalTruthPos(0);
919  float gy = globalTruthPos(1);
920  float gz = globalTruthPos(2);
921 
923  const float r = sqrt(gx * gx + gy * gy + gz * gz);
924  Acts::Vector3D globalTruthUnitDir(gx / r, gy / r, gz / r);
925 
926  auto surf = getSurface(key, cluster->getSubSurfKey());
927 
928  auto truthLocal = (*surf).globalToLocal(m_tGeometry->geoContext,
929  globalTruthPos,
930  globalTruthUnitDir);
931 
932  if(truthLocal.ok())
933  {
934  Acts::Vector2D truthLocalVec = truthLocal.value();
935 
936  m_t_SL_lx.push_back(truthLocalVec(0));
937  m_t_SL_ly.push_back(truthLocalVec(1));
938  m_t_SL_gx.push_back(gx);
939  m_t_SL_gy.push_back(gy);
940  m_t_SL_gz.push_back(gz);
941  }
942  else
943  {
944  m_t_SL_lx.push_back(-9999.);
945  m_t_SL_ly.push_back(-9999.);
946  m_t_SL_gx.push_back(-9999.);
947  m_t_SL_gy.push_back(-9999.);
948  m_t_SL_gz.push_back(-9999.);
949  }
950  }
951 
952  if(Verbosity() > 2)
953  std::cout << "Filled proto track" << std::endl;
954 
955 }
956 
957 void ActsEvaluator::fillFittedTrackParams(const Trajectory traj,
958  const size_t &trackTip,
959  const Acts::Vector3D vertex)
960 {
961  m_hasFittedParams = false;
962 
963  if(Verbosity() > 2)
964  std::cout << "Filling fitted track parameters" << std::endl;
965 
967  if (traj.hasTrackParameters(trackTip))
968  {
969  m_hasFittedParams = true;
970  const auto &boundParam = traj.trackParameters(trackTip);
971  const auto &parameter = boundParam.parameters();
972  const auto &covariance = *boundParam.covariance();
973  m_charge_fit = boundParam.charge();
974  m_eLOC0_fit = parameter[Acts::eBoundLoc0];
975  m_eLOC1_fit = parameter[Acts::eBoundLoc1];
976  m_ePHI_fit = parameter[Acts::eBoundPhi];
977  m_eTHETA_fit = parameter[Acts::eBoundTheta];
978  m_eQOP_fit = parameter[Acts::eBoundQOverP];
979  m_eT_fit = parameter[Acts::eBoundTime];
980  m_err_eLOC0_fit =
981  sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0));
982  m_err_eLOC1_fit =
983  sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1));
984  m_err_ePHI_fit = sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi));
985  m_err_eTHETA_fit =
986  sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta));
987  m_err_eQOP_fit = sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP));
988  m_err_eT_fit = sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime));
989 
990  m_px_fit = boundParam.momentum()(0);
991  m_py_fit = boundParam.momentum()(1);
992  m_pz_fit = boundParam.momentum()(2);
993  m_x_fit = boundParam.position(m_tGeometry->geoContext)(0);
994  m_y_fit = boundParam.position(m_tGeometry->geoContext)(1);
995  m_z_fit = boundParam.position(m_tGeometry->geoContext)(2);
996 
997  calculateDCA(boundParam, vertex);
998 
999  return;
1000  }
1001 
1003  m_eLOC0_fit = -9999;
1004  m_eLOC1_fit = -9999;
1005  m_ePHI_fit = -9999;
1006  m_eTHETA_fit = -9999;
1007  m_eQOP_fit = -9999;
1008  m_eT_fit = -9999;
1009  m_charge_fit = -9999;
1010  m_err_eLOC0_fit = -9999;
1011  m_err_eLOC1_fit = -9999;
1012  m_err_ePHI_fit = -9999;
1013  m_err_eTHETA_fit = -9999;
1014  m_err_eQOP_fit = -9999;
1015  m_err_eT_fit = -9999;
1016  m_px_fit = -9999;
1017  m_py_fit = -9999;
1018  m_pz_fit = -9999;
1019  m_x_fit = -9999;
1020  m_y_fit = -9999;
1021  m_z_fit = -9999;
1022 
1023  if(Verbosity() > 2)
1024  std::cout << "Finished fitted track params" << std::endl;
1025 
1026  return;
1027 }
1028 
1029 void ActsEvaluator::calculateDCA(const Acts::BoundTrackParameters param,
1030  const Acts::Vector3D vertex)
1031 {
1032 
1033  Acts::Vector3D pos = param.position(m_tGeometry->geoContext);
1034  Acts::Vector3D mom = param.momentum();
1035 
1037  pos -= vertex;
1038 
1039  Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
1040  if(param.covariance())
1041  cov = param.covariance().value();
1042 
1043  Acts::ActsSymMatrixD<3> posCov;
1044  for(int i = 0; i < 3; ++i)
1045  {
1046  for(int j = 0; j < 3; ++j)
1047  {
1048  posCov(i,j) = cov(i,j);
1049  }
1050  }
1051 
1052  Acts::Vector3D r = mom.cross(Acts::Vector3D(0.,0.,1.));
1053  float phi = atan2(r(1), r(0));
1054 
1055  Acts::RotationMatrix3D rot;
1056  Acts::RotationMatrix3D rot_T;
1057  rot(0,0) = cos(phi);
1058  rot(0,1) = -sin(phi);
1059  rot(0,2) = 0;
1060  rot(1,0) = sin(phi);
1061  rot(1,1) = cos(phi);
1062  rot(1,2) = 0;
1063  rot(2,0) = 0;
1064  rot(2,1) = 0;
1065  rot(2,2) = 1;
1066 
1067  rot_T = rot.transpose();
1068 
1069  Acts::Vector3D pos_R = rot * pos;
1070  Acts::ActsSymMatrixD<3> rotCov = rot * posCov * rot_T;
1071 
1072  m_dca3Dxy = pos_R(0);
1073  m_dca3Dz = pos_R(2);
1074  m_dca3DxyCov = rotCov(0,0);
1075  m_dca3DzCov = rotCov(2,2);
1076 
1077 }
1078 
1079 void ActsEvaluator::fillG4Particle(PHG4Particle *part)
1080 {
1081  SvtxTruthEval *trutheval = m_svtxEvalStack->get_truth_eval();
1082 
1083  if (part)
1084  {
1085  m_t_barcode = part->get_track_id();
1086  const auto pid = part->get_pid();
1087  m_t_charge = pid < 0 ? -1 : 1;
1088  const auto vtx = trutheval->get_vertex(part);
1089  m_t_vx = vtx->get_x() * Acts::UnitConstants::cm;
1090  m_t_vy = vtx->get_y() * Acts::UnitConstants::cm;
1091  m_t_vz = vtx->get_z() * Acts::UnitConstants::cm;
1092  if(Verbosity() > 1)
1093  std::cout << "truth vertex : (" << m_t_vx << ", " << m_t_vy
1094  << ", " << m_t_vz << ")" << std::endl;
1095  m_t_px = part->get_px();
1096  m_t_py = part->get_py();
1097  m_t_pz = part->get_pz();
1098  const double p = sqrt(m_t_px * m_t_px + m_t_py * m_t_py + m_t_pz * m_t_pz);
1099  m_t_theta = acos(m_t_pz / p);
1100  m_t_phi = atan(m_t_py / m_t_px);
1101  m_t_pT = sqrt(m_t_px * m_t_px + m_t_py * m_t_py);
1102  m_t_eta = atanh(m_t_pz / p);
1103 
1104  return;
1105  }
1106 
1108  m_t_barcode = -9999;
1109  m_t_charge = -9999;
1110  m_t_vx = -9999;
1111  m_t_vy = -9999;
1112  m_t_vz = -9999;
1113  m_t_px = -9999;
1114  m_t_py = -9999;
1115  m_t_pz = -9999;
1116  m_t_theta = -9999;
1117  m_t_phi = -9999;
1118  m_t_pT = -9999;
1119  m_t_eta = -9999;
1120 
1121  return;
1122 }
1123 
1124 int ActsEvaluator::getNodes(PHCompositeNode *topNode)
1125 {
1126  m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
1127  if(!m_surfMaps)
1128  {
1129  std::cout << PHWHERE << "Acts surface maps not on node tree."
1130  << std::endl;
1131  return Fun4AllReturnCodes::ABORTEVENT;
1132  }
1133 
1134  m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
1135  if(!m_vertexMap)
1136  {
1137  std::cout << PHWHERE << "SvtxVertexMap not found, cannot continue!"
1138  << std::endl;
1139  return Fun4AllReturnCodes::ABORTEVENT;
1140  }
1141 
1142  m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1143 
1144  if (!m_truthInfo)
1145  {
1146  std::cout << PHWHERE << "PHG4TruthInfoContainer not found, cannot continue!"
1147  << std::endl;
1148 
1149  return Fun4AllReturnCodes::ABORTEVENT;
1150  }
1151 
1152  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
1153 
1154  if (!m_tGeometry)
1155  {
1156  std::cout << PHWHERE << "No Acts Tracking geometry on node tree. Bailing"
1157  << std::endl;
1158 
1159  return Fun4AllReturnCodes::ABORTEVENT;
1160  }
1161 
1162  m_actsTrackKeyMap = findNode::getClass<std::map<const unsigned int,
1163  std::map<const size_t,
1164  const unsigned int>>>
1165  (topNode, "ActsTrackKeys");
1166 
1167  if (!m_actsTrackKeyMap)
1168  {
1169  if(Verbosity() > 1)
1170  std::cout << PHWHERE << "No acts CKF track map on node tree."
1171  << std::endl
1172  << "If you are analyzing the CKF, your results will be incorrect."
1173  << std::endl;
1174 
1175  }
1176 
1177  m_actsFitResults = findNode::getClass<std::map<const unsigned int, Trajectory>>
1178  (topNode, "ActsTrajectories");
1179 
1180  if (!m_actsFitResults)
1181  {
1182  std::cout << PHWHERE << "No Acts fit results on node tree. Bailing."
1183  << std::endl;
1184  return Fun4AllReturnCodes::ABORTEVENT;
1185  }
1186 
1187  m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
1188 
1189  if (!m_trackMap)
1190  {
1191  std::cout << PHWHERE << "No SvtxTrackMap on node tree. Bailing."
1192  << std::endl;
1193 
1194  return Fun4AllReturnCodes::ABORTEVENT;
1195  }
1196 
1197  m_actsProtoTrackMap = findNode::getClass<SvtxTrackMap>(topNode, "SeedTrackMap");
1198  if (!m_actsProtoTrackMap)
1199  {
1200  std::cout << PHWHERE << "No Acts proto tracks on node tree. Bailing."
1201  << std::endl;
1202  return Fun4AllReturnCodes::ABORTEVENT;
1203  }
1204 
1205  m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode,"TRKR_CLUSTER");
1206  if(!m_clusterContainer)
1207  {
1208  std::cout << PHWHERE << "No clusters, bailing"
1209  << std::endl;
1210  return Fun4AllReturnCodes::ABORTEVENT;
1211  }
1212 
1213  return Fun4AllReturnCodes::EVENT_OK;
1214 }
1215 
1216 void ActsEvaluator::clearTrackVariables()
1217 {
1218  m_t_x.clear();
1219  m_t_y.clear();
1220  m_t_z.clear();
1221  m_t_r.clear();
1222  m_t_dx.clear();
1223  m_t_dy.clear();
1224  m_t_dz.clear();
1225  m_t_eLOC0.clear();
1226  m_t_eLOC1.clear();
1227  m_t_ePHI.clear();
1228  m_t_eTHETA.clear();
1229  m_t_eQOP.clear();
1230  m_t_eT.clear();
1231 
1232  m_volumeID.clear();
1233  m_layerID.clear();
1234  m_moduleID.clear();
1235  m_lx_hit.clear();
1236  m_ly_hit.clear();
1237  m_x_hit.clear();
1238  m_y_hit.clear();
1239  m_z_hit.clear();
1240  m_res_x_hit.clear();
1241  m_res_y_hit.clear();
1242  m_err_x_hit.clear();
1243  m_err_y_hit.clear();
1244  m_pull_x_hit.clear();
1245  m_pull_y_hit.clear();
1246  m_dim_hit.clear();
1247 
1248  m_prt.clear();
1249  m_eLOC0_prt.clear();
1250  m_eLOC1_prt.clear();
1251  m_ePHI_prt.clear();
1252  m_eTHETA_prt.clear();
1253  m_eQOP_prt.clear();
1254  m_eT_prt.clear();
1255  m_res_eLOC0_prt.clear();
1256  m_res_eLOC1_prt.clear();
1257  m_res_ePHI_prt.clear();
1258  m_res_eTHETA_prt.clear();
1259  m_res_eQOP_prt.clear();
1260  m_res_eT_prt.clear();
1261  m_err_eLOC0_prt.clear();
1262  m_err_eLOC1_prt.clear();
1263  m_err_ePHI_prt.clear();
1264  m_err_eTHETA_prt.clear();
1265  m_err_eQOP_prt.clear();
1266  m_err_eT_prt.clear();
1267  m_pull_eLOC0_prt.clear();
1268  m_pull_eLOC1_prt.clear();
1269  m_pull_ePHI_prt.clear();
1270  m_pull_eTHETA_prt.clear();
1271  m_pull_eQOP_prt.clear();
1272  m_pull_eT_prt.clear();
1273  m_x_prt.clear();
1274  m_y_prt.clear();
1275  m_z_prt.clear();
1276  m_px_prt.clear();
1277  m_py_prt.clear();
1278  m_pz_prt.clear();
1279  m_eta_prt.clear();
1280  m_pT_prt.clear();
1281 
1282  m_flt.clear();
1283  m_eLOC0_flt.clear();
1284  m_eLOC1_flt.clear();
1285  m_ePHI_flt.clear();
1286  m_eTHETA_flt.clear();
1287  m_eQOP_flt.clear();
1288  m_eT_flt.clear();
1289  m_res_eLOC0_flt.clear();
1290  m_res_eLOC1_flt.clear();
1291  m_res_ePHI_flt.clear();
1292  m_res_eTHETA_flt.clear();
1293  m_res_eQOP_flt.clear();
1294  m_res_eT_flt.clear();
1295  m_err_eLOC0_flt.clear();
1296  m_err_eLOC1_flt.clear();
1297  m_err_ePHI_flt.clear();
1298  m_err_eTHETA_flt.clear();
1299  m_err_eQOP_flt.clear();
1300  m_err_eT_flt.clear();
1301  m_pull_eLOC0_flt.clear();
1302  m_pull_eLOC1_flt.clear();
1303  m_pull_ePHI_flt.clear();
1304  m_pull_eTHETA_flt.clear();
1305  m_pull_eQOP_flt.clear();
1306  m_pull_eT_flt.clear();
1307  m_x_flt.clear();
1308  m_y_flt.clear();
1309  m_z_flt.clear();
1310  m_px_flt.clear();
1311  m_py_flt.clear();
1312  m_pz_flt.clear();
1313  m_eta_flt.clear();
1314  m_pT_flt.clear();
1315  m_chi2.clear();
1316 
1317  m_smt.clear();
1318  m_eLOC0_smt.clear();
1319  m_eLOC1_smt.clear();
1320  m_ePHI_smt.clear();
1321  m_eTHETA_smt.clear();
1322  m_eQOP_smt.clear();
1323  m_eT_smt.clear();
1324  m_res_eLOC0_smt.clear();
1325  m_res_eLOC1_smt.clear();
1326  m_res_ePHI_smt.clear();
1327  m_res_eTHETA_smt.clear();
1328  m_res_eQOP_smt.clear();
1329  m_res_eT_smt.clear();
1330  m_err_eLOC0_smt.clear();
1331  m_err_eLOC1_smt.clear();
1332  m_err_ePHI_smt.clear();
1333  m_err_eTHETA_smt.clear();
1334  m_err_eQOP_smt.clear();
1335  m_err_eT_smt.clear();
1336  m_pull_eLOC0_smt.clear();
1337  m_pull_eLOC1_smt.clear();
1338  m_pull_ePHI_smt.clear();
1339  m_pull_eTHETA_smt.clear();
1340  m_pull_eQOP_smt.clear();
1341  m_pull_eT_smt.clear();
1342  m_x_smt.clear();
1343  m_y_smt.clear();
1344  m_z_smt.clear();
1345  m_px_smt.clear();
1346  m_py_smt.clear();
1347  m_pz_smt.clear();
1348  m_eta_smt.clear();
1349  m_pT_smt.clear();
1350 
1351  m_SLx.clear();
1352  m_SLy.clear();
1353  m_SLz.clear();
1354  m_SL_lx.clear();
1355  m_SL_ly.clear();
1356  m_t_SL_lx.clear();
1357  m_t_SL_ly.clear();
1358  m_t_SL_gx.clear();
1359  m_t_SL_gy.clear();
1360  m_t_SL_gz.clear();
1361 
1362  m_protoTrackPx = -9999.;
1363  m_protoTrackPy = -9999.;
1364  m_protoTrackPz = -9999.;
1365  m_protoTrackX = -9999.;
1366  m_protoTrackY = -9999.;
1367  m_protoTrackZ = -9999.;
1368  m_protoD0Cov = -9999.;
1369  m_protoZ0Cov = -9999.;
1370  m_protoPhiCov = -9999.;
1371  m_protoThetaCov= -9999.;
1372  m_protoQopCov = -9999.;
1373 
1374  return;
1375 }
1376 
1377 void ActsEvaluator::initializeTree()
1378 {
1379  m_trackFile = new TFile(Name().c_str(), "RECREATE");
1380 
1381  m_trackTree = new TTree("tracktree", "A tree with Acts KF track information");
1382 
1383  m_trackTree->Branch("event_nr", &m_eventNr);
1384  m_trackTree->Branch("traj_nr", &m_trajNr);
1385  m_trackTree->Branch("track_nr", &m_trackNr);
1386  m_trackTree->Branch("t_barcode", &m_t_barcode, "t_barcode/l");
1387  m_trackTree->Branch("t_charge", &m_t_charge);
1388  m_trackTree->Branch("t_time", &m_t_time);
1389  m_trackTree->Branch("t_vx", &m_t_vx);
1390  m_trackTree->Branch("t_vy", &m_t_vy);
1391  m_trackTree->Branch("t_vz", &m_t_vz);
1392  m_trackTree->Branch("t_px", &m_t_px);
1393  m_trackTree->Branch("t_py", &m_t_py);
1394  m_trackTree->Branch("t_pz", &m_t_pz);
1395  m_trackTree->Branch("t_theta", &m_t_theta);
1396  m_trackTree->Branch("t_phi", &m_t_phi);
1397  m_trackTree->Branch("t_eta", &m_t_eta);
1398  m_trackTree->Branch("t_pT", &m_t_pT);
1399 
1400  m_trackTree->Branch("t_x", &m_t_x);
1401  m_trackTree->Branch("t_y", &m_t_y);
1402  m_trackTree->Branch("t_z", &m_t_z);
1403  m_trackTree->Branch("t_r", &m_t_r);
1404  m_trackTree->Branch("t_dx", &m_t_dx);
1405  m_trackTree->Branch("t_dy", &m_t_dy);
1406  m_trackTree->Branch("t_dz", &m_t_dz);
1407  m_trackTree->Branch("t_eLOC0", &m_t_eLOC0);
1408  m_trackTree->Branch("t_eLOC1", &m_t_eLOC1);
1409  m_trackTree->Branch("t_ePHI", &m_t_ePHI);
1410  m_trackTree->Branch("t_eTHETA", &m_t_eTHETA);
1411  m_trackTree->Branch("t_eQOP", &m_t_eQOP);
1412  m_trackTree->Branch("t_eT", &m_t_eT);
1413 
1414 
1415  m_trackTree->Branch("hasFittedParams", &m_hasFittedParams);
1416  m_trackTree->Branch("chi2_fit", &m_chi2_fit);
1417  m_trackTree->Branch("quality", &m_quality);
1418  m_trackTree->Branch("ndf_fit", &m_ndf_fit);
1419  m_trackTree->Branch("eLOC0_fit", &m_eLOC0_fit);
1420  m_trackTree->Branch("eLOC1_fit", &m_eLOC1_fit);
1421  m_trackTree->Branch("ePHI_fit", &m_ePHI_fit);
1422  m_trackTree->Branch("eTHETA_fit", &m_eTHETA_fit);
1423  m_trackTree->Branch("eQOP_fit", &m_eQOP_fit);
1424  m_trackTree->Branch("eT_fit", &m_eT_fit);
1425  m_trackTree->Branch("charge_fit", &m_charge_fit);
1426  m_trackTree->Branch("err_eLOC0_fit", &m_err_eLOC0_fit);
1427  m_trackTree->Branch("err_eLOC1_fit", &m_err_eLOC1_fit);
1428  m_trackTree->Branch("err_ePHI_fit", &m_err_ePHI_fit);
1429  m_trackTree->Branch("err_eTHETA_fit", &m_err_eTHETA_fit);
1430  m_trackTree->Branch("err_eQOP_fit", &m_err_eQOP_fit);
1431  m_trackTree->Branch("err_eT_fit", &m_err_eT_fit);
1432  m_trackTree->Branch("g_px_fit", &m_px_fit);
1433  m_trackTree->Branch("g_py_fit", &m_py_fit);
1434  m_trackTree->Branch("g_pz_fit", &m_pz_fit);
1435  m_trackTree->Branch("g_x_fit" , &m_x_fit);
1436  m_trackTree->Branch("g_y_fit" , &m_y_fit);
1437  m_trackTree->Branch("g_z_fit" , &m_z_fit);
1438  m_trackTree->Branch("g_dca3Dxy_fit" , &m_dca3Dxy);
1439  m_trackTree->Branch("g_dca3Dz_fit" , &m_dca3Dz);
1440  m_trackTree->Branch("g_dca3Dxy_cov" , &m_dca3DxyCov);
1441  m_trackTree->Branch("g_dca3Dz_cov" , &m_dca3DzCov);
1442 
1443  m_trackTree->Branch("g_protoTrackPx", &m_protoTrackPx);
1444  m_trackTree->Branch("g_protoTrackPy", &m_protoTrackPy);
1445  m_trackTree->Branch("g_protoTrackPz", &m_protoTrackPz);
1446  m_trackTree->Branch("g_protoTrackX", &m_protoTrackX);
1447  m_trackTree->Branch("g_protoTrackY", &m_protoTrackY);
1448  m_trackTree->Branch("g_protoTrackZ", &m_protoTrackZ);
1449  m_trackTree->Branch("g_protoTrackD0Cov", &m_protoD0Cov);
1450  m_trackTree->Branch("g_protoTrackZ0Cov", &m_protoZ0Cov);
1451  m_trackTree->Branch("g_protoTrackPhiCov", &m_protoPhiCov);
1452  m_trackTree->Branch("g_protoTrackThetaCov", &m_protoThetaCov);
1453  m_trackTree->Branch("g_protoTrackQopCov", &m_protoQopCov);
1454  m_trackTree->Branch("g_SLx", &m_SLx);
1455  m_trackTree->Branch("g_SLy", &m_SLy);
1456  m_trackTree->Branch("g_SLz", &m_SLz);
1457  m_trackTree->Branch("l_SLx", &m_SL_lx);
1458  m_trackTree->Branch("l_SLy", &m_SL_ly);
1459  m_trackTree->Branch("t_SL_lx", &m_t_SL_lx);
1460  m_trackTree->Branch("t_SL_ly", &m_t_SL_ly);
1461  m_trackTree->Branch("t_SL_gx", &m_t_SL_gx);
1462  m_trackTree->Branch("t_SL_gy", &m_t_SL_gy);
1463  m_trackTree->Branch("t_SL_gz", &m_t_SL_gz);
1464 
1465  m_trackTree->Branch("nHoles", &m_nHoles);
1466  m_trackTree->Branch("nOutliers", &m_nOutliers);
1467  m_trackTree->Branch("nStates", &m_nStates);
1468  m_trackTree->Branch("nMeasurements", &m_nMeasurements);
1469  m_trackTree->Branch("volume_id", &m_volumeID);
1470  m_trackTree->Branch("layer_id", &m_layerID);
1471  m_trackTree->Branch("module_id", &m_moduleID);
1472  m_trackTree->Branch("l_x_hit", &m_lx_hit);
1473  m_trackTree->Branch("l_y_hit", &m_ly_hit);
1474  m_trackTree->Branch("g_x_hit", &m_x_hit);
1475  m_trackTree->Branch("g_y_hit", &m_y_hit);
1476  m_trackTree->Branch("g_z_hit", &m_z_hit);
1477  m_trackTree->Branch("res_x_hit", &m_res_x_hit);
1478  m_trackTree->Branch("res_y_hit", &m_res_y_hit);
1479  m_trackTree->Branch("err_x_hit", &m_err_x_hit);
1480  m_trackTree->Branch("err_y_hit", &m_err_y_hit);
1481  m_trackTree->Branch("pull_x_hit", &m_pull_x_hit);
1482  m_trackTree->Branch("pull_y_hit", &m_pull_y_hit);
1483  m_trackTree->Branch("dim_hit", &m_dim_hit);
1484 
1485  m_trackTree->Branch("nPredicted", &m_nPredicted);
1486  m_trackTree->Branch("predicted", &m_prt);
1487  m_trackTree->Branch("eLOC0_prt", &m_eLOC0_prt);
1488  m_trackTree->Branch("eLOC1_prt", &m_eLOC1_prt);
1489  m_trackTree->Branch("ePHI_prt", &m_ePHI_prt);
1490  m_trackTree->Branch("eTHETA_prt", &m_eTHETA_prt);
1491  m_trackTree->Branch("eQOP_prt", &m_eQOP_prt);
1492  m_trackTree->Branch("eT_prt", &m_eT_prt);
1493  m_trackTree->Branch("res_eLOC0_prt", &m_res_eLOC0_prt);
1494  m_trackTree->Branch("res_eLOC1_prt", &m_res_eLOC1_prt);
1495  m_trackTree->Branch("res_ePHI_prt", &m_res_ePHI_prt);
1496  m_trackTree->Branch("res_eTHETA_prt", &m_res_eTHETA_prt);
1497  m_trackTree->Branch("res_eQOP_prt", &m_res_eQOP_prt);
1498  m_trackTree->Branch("res_eT_prt", &m_res_eT_prt);
1499  m_trackTree->Branch("err_eLOC0_prt", &m_err_eLOC0_prt);
1500  m_trackTree->Branch("err_eLOC1_prt", &m_err_eLOC1_prt);
1501  m_trackTree->Branch("err_ePHI_prt", &m_err_ePHI_prt);
1502  m_trackTree->Branch("err_eTHETA_prt", &m_err_eTHETA_prt);
1503  m_trackTree->Branch("err_eQOP_prt", &m_err_eQOP_prt);
1504  m_trackTree->Branch("err_eT_prt", &m_err_eT_prt);
1505  m_trackTree->Branch("pull_eLOC0_prt", &m_pull_eLOC0_prt);
1506  m_trackTree->Branch("pull_eLOC1_prt", &m_pull_eLOC1_prt);
1507  m_trackTree->Branch("pull_ePHI_prt", &m_pull_ePHI_prt);
1508  m_trackTree->Branch("pull_eTHETA_prt", &m_pull_eTHETA_prt);
1509  m_trackTree->Branch("pull_eQOP_prt", &m_pull_eQOP_prt);
1510  m_trackTree->Branch("pull_eT_prt", &m_pull_eT_prt);
1511  m_trackTree->Branch("g_x_prt", &m_x_prt);
1512  m_trackTree->Branch("g_y_prt", &m_y_prt);
1513  m_trackTree->Branch("g_z_prt", &m_z_prt);
1514  m_trackTree->Branch("px_prt", &m_px_prt);
1515  m_trackTree->Branch("py_prt", &m_py_prt);
1516  m_trackTree->Branch("pz_prt", &m_pz_prt);
1517  m_trackTree->Branch("eta_prt", &m_eta_prt);
1518  m_trackTree->Branch("pT_prt", &m_pT_prt);
1519 
1520  m_trackTree->Branch("nFiltered", &m_nFiltered);
1521  m_trackTree->Branch("filtered", &m_flt);
1522  m_trackTree->Branch("eLOC0_flt", &m_eLOC0_flt);
1523  m_trackTree->Branch("eLOC1_flt", &m_eLOC1_flt);
1524  m_trackTree->Branch("ePHI_flt", &m_ePHI_flt);
1525  m_trackTree->Branch("eTHETA_flt", &m_eTHETA_flt);
1526  m_trackTree->Branch("eQOP_flt", &m_eQOP_flt);
1527  m_trackTree->Branch("eT_flt", &m_eT_flt);
1528  m_trackTree->Branch("res_eLOC0_flt", &m_res_eLOC0_flt);
1529  m_trackTree->Branch("res_eLOC1_flt", &m_res_eLOC1_flt);
1530  m_trackTree->Branch("res_ePHI_flt", &m_res_ePHI_flt);
1531  m_trackTree->Branch("res_eTHETA_flt", &m_res_eTHETA_flt);
1532  m_trackTree->Branch("res_eQOP_flt", &m_res_eQOP_flt);
1533  m_trackTree->Branch("res_eT_flt", &m_res_eT_flt);
1534  m_trackTree->Branch("err_eLOC0_flt", &m_err_eLOC0_flt);
1535  m_trackTree->Branch("err_eLOC1_flt", &m_err_eLOC1_flt);
1536  m_trackTree->Branch("err_ePHI_flt", &m_err_ePHI_flt);
1537  m_trackTree->Branch("err_eTHETA_flt", &m_err_eTHETA_flt);
1538  m_trackTree->Branch("err_eQOP_flt", &m_err_eQOP_flt);
1539  m_trackTree->Branch("err_eT_flt", &m_err_eT_flt);
1540  m_trackTree->Branch("pull_eLOC0_flt", &m_pull_eLOC0_flt);
1541  m_trackTree->Branch("pull_eLOC1_flt", &m_pull_eLOC1_flt);
1542  m_trackTree->Branch("pull_ePHI_flt", &m_pull_ePHI_flt);
1543  m_trackTree->Branch("pull_eTHETA_flt", &m_pull_eTHETA_flt);
1544  m_trackTree->Branch("pull_eQOP_flt", &m_pull_eQOP_flt);
1545  m_trackTree->Branch("pull_eT_flt", &m_pull_eT_flt);
1546  m_trackTree->Branch("g_x_flt", &m_x_flt);
1547  m_trackTree->Branch("g_y_flt", &m_y_flt);
1548  m_trackTree->Branch("g_z_flt", &m_z_flt);
1549  m_trackTree->Branch("px_flt", &m_px_flt);
1550  m_trackTree->Branch("py_flt", &m_py_flt);
1551  m_trackTree->Branch("pz_flt", &m_pz_flt);
1552  m_trackTree->Branch("eta_flt", &m_eta_flt);
1553  m_trackTree->Branch("pT_flt", &m_pT_flt);
1554  m_trackTree->Branch("chi2", &m_chi2);
1555 
1556  m_trackTree->Branch("nSmoothed", &m_nSmoothed);
1557  m_trackTree->Branch("smoothed", &m_smt);
1558  m_trackTree->Branch("eLOC0_smt", &m_eLOC0_smt);
1559  m_trackTree->Branch("eLOC1_smt", &m_eLOC1_smt);
1560  m_trackTree->Branch("ePHI_smt", &m_ePHI_smt);
1561  m_trackTree->Branch("eTHETA_smt", &m_eTHETA_smt);
1562  m_trackTree->Branch("eQOP_smt", &m_eQOP_smt);
1563  m_trackTree->Branch("eT_smt", &m_eT_smt);
1564  m_trackTree->Branch("res_eLOC0_smt", &m_res_eLOC0_smt);
1565  m_trackTree->Branch("res_eLOC1_smt", &m_res_eLOC1_smt);
1566  m_trackTree->Branch("res_ePHI_smt", &m_res_ePHI_smt);
1567  m_trackTree->Branch("res_eTHETA_smt", &m_res_eTHETA_smt);
1568  m_trackTree->Branch("res_eQOP_smt", &m_res_eQOP_smt);
1569  m_trackTree->Branch("res_eT_smt", &m_res_eT_smt);
1570  m_trackTree->Branch("err_eLOC0_smt", &m_err_eLOC0_smt);
1571  m_trackTree->Branch("err_eLOC1_smt", &m_err_eLOC1_smt);
1572  m_trackTree->Branch("err_ePHI_smt", &m_err_ePHI_smt);
1573  m_trackTree->Branch("err_eTHETA_smt", &m_err_eTHETA_smt);
1574  m_trackTree->Branch("err_eQOP_smt", &m_err_eQOP_smt);
1575  m_trackTree->Branch("err_eT_smt", &m_err_eT_smt);
1576  m_trackTree->Branch("pull_eLOC0_smt", &m_pull_eLOC0_smt);
1577  m_trackTree->Branch("pull_eLOC1_smt", &m_pull_eLOC1_smt);
1578  m_trackTree->Branch("pull_ePHI_smt", &m_pull_ePHI_smt);
1579  m_trackTree->Branch("pull_eTHETA_smt", &m_pull_eTHETA_smt);
1580  m_trackTree->Branch("pull_eQOP_smt", &m_pull_eQOP_smt);
1581  m_trackTree->Branch("pull_eT_smt", &m_pull_eT_smt);
1582  m_trackTree->Branch("g_x_smt", &m_x_smt);
1583  m_trackTree->Branch("g_y_smt", &m_y_smt);
1584  m_trackTree->Branch("g_z_smt", &m_z_smt);
1585  m_trackTree->Branch("px_smt", &m_px_smt);
1586  m_trackTree->Branch("py_smt", &m_py_smt);
1587  m_trackTree->Branch("pz_smt", &m_pz_smt);
1588  m_trackTree->Branch("eta_smt", &m_eta_smt);
1589  m_trackTree->Branch("pT_smt", &m_pT_smt);
1590 }