36 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
38 #include <ActsExamples/EventData/Track.hpp>
39 #include <ActsExamples/Framework/AlgorithmContext.hpp>
47 , m_trajectories(nullptr)
53 std::cout <<
"Setup PHActsTrkFitter" << std::endl;
61 m_fitCfg.fit = ActsExamples::TrkrClusterFittingAlgorithm::makeFitterFunction(
65 m_fitCfg.dFit = ActsExamples::TrkrClusterFittingAlgorithm::makeFitterFunction(
74 h_fitTime =
new TH2F(
"h_fitTime",
";p_{T} [GeV];time [ms]",
75 80, 0, 40, 100000, 0, 1000);
79 h_rotTime =
new TH1F(
"h_rotTime",
";time [ms]",
86 std::cout <<
"Finish PHActsTrkFitter Setup" << std::endl;
93 PHTimer eventTimer(
"eventTimer");
103 std::cout <<
PHWHERE <<
"Events processed: " <<
m_event << std::endl;
104 std::cout <<
"Start PHActsTrkFitter::process_event" << std::endl;
128 std::cout <<
"PHActsTrkFitter total event time "
129 << eventTime << std::endl;
136 std::cout <<
"PHActsTrkFitter::process_event finished"
152 std::cout <<
"Reset PHActsTrkFitter" << std::endl;
178 std::cout <<
"The Acts track fitter had " <<
m_nBadFits
179 <<
" fits return an error" << std::endl;
181 std::cout <<
"Finished PHActsTrkFitter" << std::endl;
192 std::vector<unsigned int> badTracks;
199 PHTimer trackTimer(
"TrackTimer");
204 if(sourceLinks.size() == 0)
continue;
214 if( surfaces.empty() )
continue;
218 std::none_of( surfaces.begin(), surfaces.end(), [
this](
const auto&
surface )
231 auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
239 if(
m_fieldMap.find(
"3d") != std::string::npos)
256 Acts::PropagatorPlainOptions ppPlainOptions;
263 Acts::LoggerWrapper(*logger),
276 std::cout <<
"PHActsTrkFitter Acts fit time "
277 << fitTime << std::endl;
282 const FitResult& fitOutput = result.value();
287 .transverseMomentum(),
305 badTracks.push_back(trackKey);
314 std::cout <<
"PHActsTrkFitter total single track time "
315 << trackTime << std::endl;
320 for(
const auto& key : badTracks)
322 m_trackMap->erase(key);
356 auto iter = surfMap.find(hitsetkey);
357 if(iter != surfMap.end())
374 auto surfvec = iter->second;
375 return surfvec.at(surfkey);
399 auto key = *clusIter;
403 if(
Verbosity() > 0) std::cout <<
"Failed to get cluster with key " << key <<
" for track " << track->
get_id() << std::endl;
432 std::cout <<
"source link " << sl.cluskey() <<
", loc : "
433 << sl.location().transpose() << std::endl <<
", cov : " << sl.covariance().transpose() << std::endl
434 <<
" geo id " << sl.geoId() << std::endl;
435 std::cout <<
"Surface : " << std::endl;
437 std::cout << std::endl;
438 std::cout <<
"Cluster error " << cluster->getRPhiError() <<
" , " << cluster->getZError() << std::endl;
439 std::cout <<
"For key " << key <<
" with local pos " << std::endl
440 << cluster->getLocalX() <<
", " << cluster->getLocalY()
444 sourcelinks.push_back(sl);
456 std::vector<size_t> trackTips;
457 trackTips.push_back(fitOutput.
trackTip);
458 ActsExamples::IndexedParams indexedParams;
461 indexedParams.emplace(fitOutput.
trackTip,
468 std::cout <<
"Fitted parameters for track" << std::endl;
472 std::cout <<
"charge: "<<params.charge()<<std::endl;
473 std::cout <<
" momentum : " << params.momentum().transpose()
475 std::cout <<
"For trackTip == " << fitOutput.
trackTip << std::endl;
479 auto trajectory = std::make_unique<Trajectory>(fitOutput.
fittedStates,
480 trackTips, indexedParams,
488 PHTimer updateTrackTimer(
"UpdateTrackTimer");
489 updateTrackTimer.
stop();
494 updateTrackTimer.
stop();
498 std::cout <<
"PHActsTrkFitter update SvtxTrack time "
499 << updateTime << std::endl;
515 return m_fitCfg.dFit(sourceLinks, seed, kfOptions, surfSequence);
517 return m_fitCfg.fit(sourceLinks, seed, kfOptions);
526 std::cout <<
"Sorting " << sourceLinks.size() <<
" SLs" << std::endl;
528 for(
const auto& sl : sourceLinks)
531 { std::cout<<
"SL available on : " << sl.referenceSurface().geometryId()<<std::endl; }
540 siliconMMSls.push_back(sl);
541 surfaces.push_back(&sl.referenceSurface());
548 if(!surfaces.empty())
553 for(
const auto& surf : surfaces)
555 std::cout <<
"Surface vector : " << surf->geometryId() << std::endl;
565 for(
int i = 0; i < surfaces.size() - 1; i++)
568 auto thisVolume =
surface->geometryId().volume();
569 auto thisLayer =
surface->geometryId().layer();
571 auto nextSurface = surfaces.at(i+1);
572 auto nextVolume = nextSurface->geometryId().volume();
573 auto nextLayer = nextSurface->geometryId().layer();
576 if(nextVolume == thisVolume)
578 if(nextLayer < thisLayer)
582 <<
"Surface not in order... removing surface"
583 <<
surface->geometryId() << std::endl;
584 surfaces.erase(surfaces.begin() + i);
592 if(nextVolume < thisVolume)
596 <<
"Volume not in order... removing surface"
597 <<
surface->geometryId() << std::endl;
598 surfaces.erase(surfaces.begin() + i);
611 const auto &[trackTips, mj] = traj.trajectory();
613 auto &trackTip = trackTips.front();
617 std::cout <<
"Identify (proto) track before updating with acts results " << std::endl;
628 float pathlength = 0.0;
636 Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
638 const auto& params = traj.trackParameters(trackTip);
648 track->
set_px(params.momentum()(0));
649 track->
set_py(params.momentum()(1));
650 track->
set_pz(params.momentum()(2));
659 if(params.covariance())
665 for(
int i = 0; i < 6; i++)
667 for(
int j = 0; j < 6; j++)
671 params.covariance().value()(i,j));
679 PHTimer trackStateTimer(
"TrackStateTimer");
680 trackStateTimer.
stop();
687 trackStateTimer.
stop();
691 std::cout <<
"PHActsTrkFitter update SvtxTrackStates time "
692 << stateTime << std::endl;
699 std::cout <<
" Identify fitted track after updating track states:"
728 0., 0., 0.1, 0., 0., 0.,
729 0., 0., 0., 0.1, 0., 0.,
730 0., 0., 0., 0., 0.005 , 0.,
731 0., 0., 0., 0., 0., 1.;
737 0., 0., 0.05, 0., 0., 0.,
738 0., 0., 0., 0.05, 0., 0.,
739 0., 0., 0., 0., 0.00005 , 0.,
740 0., 0., 0., 0., 0., 1.;
748 std::cout <<
PHWHERE <<
" Processing proto track with id: "
749 << seed->
get_id() <<
" and position:"
750 << seed->
get_x() <<
", " << seed->
get_y() <<
", "
751 << seed->
get_z() << std::endl
752 <<
"momentum: " << seed->
get_px() <<
", " << seed->
get_py()
753 <<
", " << seed->
get_pz() << std::endl
767 std::cerr <<
"DST node is missing, quitting" << std::endl;
768 throw std::runtime_error(
"Failed to find DST node in PHActsTrkFitter::createNodes");
782 "SvtxSiliconMMTrackMap");
794 m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode,
"ActsTrajectories");
823 m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
826 std::cout <<
PHWHERE <<
"ActsSurfaceMaps not on node tree, bailing."
831 m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode,
"CORRECTED_TRKR_CLUSTER");
834 std::cout <<
" Using CORRECTED_TRKR_CLUSTER node " << std::endl;
838 std::cout <<
" CORRECTED_TRKR_CLUSTER node not found, using TRKR_CLUSTER" << std::endl;
839 m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
845 <<
"No trkr cluster container, exiting." << std::endl;
849 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
852 std::cout <<
"ActsTrackingGeometry not on node tree. Exiting."
862 std::cout <<
PHWHERE <<
"SvtxTrackMap not found on node tree. Exiting."