28 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
38 , m_svtxEvaluator(svtxEvaluator)
39 , m_truthInfo(nullptr)
41 , m_svtxEvalStack(nullptr)
42 , m_actsTrackKeyMap(nullptr)
43 , m_actsFitResults(nullptr)
44 , m_tGeometry(nullptr)
57 std::cout <<
"Starting ActsEvaluator::Init" << std::endl;
64 std::cout <<
"Finished ActsEvaluator::Init" << std::endl;
74 std::cout <<
"Starting ActsEvaluator at event " <<
m_eventNr
89 std::cout <<
"Finished ActsEvaluator::process_event" << std::endl;
98 std::cout <<
"Evaluating Acts track fits" << std::endl;
102 std::map<const unsigned int, Trajectory>::iterator trajIter;
111 unsigned int trackKey = trajIter->first;
115 std::cout <<
"Starting trajectory " << iTraj
116 <<
" with trackKey corresponding to track seed "
117 << trackKey << std::endl;
119 const auto &[trackTips, mj] = traj.trajectory();
123 if (trackTips.empty())
126 std::cout <<
"TrackTips empty in ActsEvaluator" << std::endl;
137 std::map<const size_t, const unsigned int> trackKeyMap;
142 for(
const size_t &trackTip : trackTips)
145 std::cout <<
"beginning trackTip " << trackTip
146 <<
" corresponding to track " << iTrack
147 <<
" in trajectroy " << iTraj << std::endl;
154 std::map<const size_t, const unsigned int>::iterator ckfiter
155 = trackKeyMap.find(trackTip);
156 if(ckfiter == trackKeyMap.end())
159 std::cout <<
"Track result flagged as bad, continuing"
163 trackKey = ckfiter->second;
167 std::cout<<
"Evaluating track key " << trackKey
168 <<
" for track tip " << trackTip << std::endl;
173 if(vertexId == UINT_MAX)
183 std::cout <<
"Analyzing SvtxTrack "<< trackKey << std::endl;
185 std::cout <<
"TruthParticle : " << g4particle->
get_px()
186 <<
", " << g4particle->
get_py() <<
", "
187 << g4particle->
get_pz() <<
", "<< g4particle->
get_e()
195 Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
199 if(traj.hasTrackParameters(trackTip))
201 std::cout <<
"Fitted params : "
203 << std::endl << traj.trackParameters(trackTip).momentum()
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;
231 std::cout <<
"Finished track " << iTrack <<std::endl;
239 std::cout <<
"Analyzed " << iTrack <<
" tracks in trajectory number "
240 << iTraj << std::endl;
247 std::cout <<
"Finished evaluating track fits" << std::endl;
271 const size_t &trackTip,
276 std::cout <<
"Begin visit track states" << std::endl;
278 const auto &[trackTips, mj] = traj.trajectory();
280 mj.visitBackwards(trackTip, [&](
const auto &state) {
282 auto typeFlags = state.typeFlags();
289 auto geoID = state.referenceSurface().geometryId();
295 std::cout <<
"Cluster volume : layer : sensitive " << geoID.volume()
296 <<
" : " << geoID.layer() <<
" : "
297 << geoID.sensitive() << std::endl;
299 auto meas = std::get<Measurement>(*state.uncalibrated());
313 auto cov = meas.covariance();
327 float gx = globalTruthPos(0);
328 float gy = globalTruthPos(1);
329 float gz = globalTruthPos(2);
332 const float r = sqrt(gx * gx + gy * gy + gz * gz);
335 auto vecResult = meas.referenceObject().globalToLocal(
344 m_t_r.push_back(sqrt(gx * gx + gy * gy));
353 float truthTHETA = 0;
363 truthLOC0 = truthLocVec.x();
364 truthLOC1 = truthLocVec.y();
372 truthPHI =
phi(globalTruthUnitDir);
373 truthTHETA =
theta(globalTruthUnitDir);
383 m_t_eT.push_back(truthTIME);
386 bool predicted =
false;
387 if (state.hasPredicted())
393 auto covariance = state.predictedCovariance();
396 auto H = meas.projector();
397 auto resCov = cov +
H * covariance *
H.transpose();
411 m_dim_hit.push_back(state.calibratedSize());
437 sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
439 sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
441 sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
443 sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
445 sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
447 sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
452 sqrt(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)));
455 sqrt(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)));
458 sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)));
460 (
parameters[Acts::eBoundTheta] - truthTHETA) /
461 sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)));
464 sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)));
467 sqrt(covariance(Acts::eBoundTime, Acts::eBoundTime)));
470 Acts::detail::transformBoundToFreeParameters(state.referenceSurface(),
481 m_pT_prt.push_back(
p * std::hypot(freeParams[Acts::eFreeDir0],
482 freeParams[Acts::eFreeDir1]));
530 bool filtered =
false;
531 if (state.hasFiltered())
537 state.referenceSurface().getSharedPtr(),
539 state.filteredCovariance());
541 auto covariance = state.filteredCovariance();
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()));
607 m_chi2.push_back(state.chi2());
647 bool smoothed =
false;
648 if (state.hasSmoothed())
654 state.referenceSurface().getSharedPtr(),
656 state.smoothedCovariance());
658 auto covariance = state.smoothedCovariance();
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()));
762 m_prt.push_back(predicted);
763 m_flt.push_back(filtered);
764 m_smt.push_back(smoothed);
772 std::cout <<
"Finished track states" << std::endl;
794 gx = truth_cluster->getX();
795 gy = truth_cluster->getY();
796 gz = truth_cluster->getZ();
837 auto iter = surfMap.find(hitsetkey);
838 if(iter != surfMap.end())
855 auto surfvec = iter->second;
856 return surfvec.at(surfkey);
875 std::cout <<
"Filling proto track seed quantities" << std::endl;
878 track->
get_z() * 10);
895 auto key = *clusIter;
900 cluster->getLocalY() * 10);
904 cluster->getY() * 10,
905 cluster->getZ() * 10);
907 m_SLx.push_back(globalPos(0));
908 m_SLy.push_back(globalPos(1));
909 m_SLz.push_back(globalPos(2));
918 float gx = globalTruthPos(0);
919 float gy = globalTruthPos(1);
920 float gz = globalTruthPos(2);
923 const float r = sqrt(gx * gx + gy * gy + gz * gz);
926 auto surf =
getSurface(key, cluster->getSubSurfKey());
953 std::cout <<
"Filled proto track" << std::endl;
958 const size_t &trackTip,
964 std::cout <<
"Filling fitted track parameters" << std::endl;
967 if (traj.hasTrackParameters(trackTip))
970 const auto &boundParam = traj.trackParameters(trackTip);
971 const auto ¶meter = boundParam.parameters();
972 const auto &covariance = *boundParam.covariance();
990 m_px_fit = boundParam.momentum()(0);
991 m_py_fit = boundParam.momentum()(1);
992 m_pz_fit = boundParam.momentum()(2);
1024 std::cout <<
"Finished fitted track params" << std::endl;
1040 if(param.covariance())
1041 cov = param.covariance().value();
1044 for(
int i = 0; i < 3; ++i)
1046 for(
int j = 0; j < 3; ++j)
1048 posCov(i,j) = cov(i,j);
1053 float phi = atan2(
r(1),
r(0));
1057 rot(0,0) = cos(phi);
1058 rot(0,1) = -sin(phi);
1060 rot(1,0) = sin(phi);
1061 rot(1,1) = cos(phi);
1067 rot_T = rot.transpose();
1088 const auto vtx = trutheval->
get_vertex(part);
1093 std::cout <<
"truth vertex : (" <<
m_t_vx <<
", " <<
m_t_vy
1094 <<
", " <<
m_t_vz <<
")" << std::endl;
1126 m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
1129 std::cout <<
PHWHERE <<
"Acts surface maps not on node tree."
1134 m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode,
"SvtxVertexMap");
1137 std::cout <<
PHWHERE <<
"SvtxVertexMap not found, cannot continue!"
1142 m_truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
1146 std::cout <<
PHWHERE <<
"PHG4TruthInfoContainer not found, cannot continue!"
1152 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
1156 std::cout <<
PHWHERE <<
"No Acts Tracking geometry on node tree. Bailing"
1163 std::map<
const size_t,
1164 const unsigned int>>>
1165 (topNode,
"ActsTrackKeys");
1170 std::cout <<
PHWHERE <<
"No acts CKF track map on node tree."
1172 <<
"If you are analyzing the CKF, your results will be incorrect."
1177 m_actsFitResults = findNode::getClass<std::map<const unsigned int, Trajectory>>
1178 (topNode,
"ActsTrajectories");
1182 std::cout <<
PHWHERE <<
"No Acts fit results on node tree. Bailing."
1187 m_trackMap = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
1191 std::cout <<
PHWHERE <<
"No SvtxTrackMap on node tree. Bailing."
1200 std::cout <<
PHWHERE <<
"No Acts proto tracks on node tree. Bailing."
1205 m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
1208 std::cout <<
PHWHERE <<
"No clusters, bailing"
1381 m_trackTree =
new TTree(
"tracktree",
"A tree with Acts KF track information");