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