17 const std::string& name, ISvcLocator* pSvcLocator)
39 ATH_MSG_ERROR(
"Use either track finder or track reader, not both.");
40 return StatusCode::FAILURE;
57 ATH_MSG_INFO(
"Applying the following cuts during GNN-based track reconstruction: ");
66 ATH_MSG_INFO(
"Applying the following eta dependant track cuts after GNN-based track reconstruction: ");
67 std::vector <double> etaBins, minPT, maxz0, maxd0;
68 std::vector <int> minClusters, minPixelHits, maxHoles;
85 return StatusCode::SUCCESS;
94 ATH_CHECK(outputEdgeScores.
record(std::make_unique<std::vector<std::vector<float>>>()));
97 uint32_t runNumber = ctx.eventID().run_number();
98 uint32_t eventNumber = ctx.eventID().event_number();
104 std::vector<const Trk::PrepRawData*> allClusters =
getClustersInEvent(ctx, eventNumber);
107 std::vector<std::vector<uint32_t> > TT;
108 std::vector<std::vector<uint32_t> > clusterTracks;
109 std::unordered_map<int, std::unordered_map<int, float>> edgeMap;
123 ATH_MSG_ERROR(
"Both GNNTrackFinder and GNNTrackReader are not set");
124 return StatusCode::FAILURE;
127 ATH_MSG_DEBUG(
"Event " << eventNumber <<
" obtained " << TT.size() <<
" Tracks");
131 int trackCounter = -1;
132 std::vector<int> status_codes;
134 for (
auto& trackIndices : TT) {
139 std::pair<std::vector<const Trk::SpacePoint*>, std::vector<uint32_t> > trackInfos =
getSpacePoints(trackIndices, spacePoints);
140 std::vector<const Trk::SpacePoint*> trackCandidate = trackInfos.first;
141 std::vector<uint32_t> sortedID = trackInfos.second;
146 getClusters (clusterTracks, allClusters, trackCounter)
156 auto [fitSuccess, passTrackCut, track] =
doFitAndCut(ctx, trackCandidate, clusters, trackCounter);
158 if (not fitSuccess) {
162 if (passTrackCut <= 0) {
163 outputTracks->push_back(track.release());
166 outputEdgeScores->push_back(edges);
170 status_codes.push_back(passTrackCut);
175 ATH_MSG_INFO(
"Event " << eventNumber <<
" has " << status_codes.size() <<
" tracks found, "
176 << std::count(status_codes.begin(), status_codes.end(), 0)
177 <<
" tracks remains after applying track cuts");
179 ATH_MSG_INFO(
"Event " << eventNumber <<
" has " << status_codes.size() <<
" tracks found, all tracks are kept");
182 return StatusCode::SUCCESS;
186 const std::vector<uint32_t>& sortedID,
187 const std::unordered_map<
int, std::unordered_map<int, float>>& edgeMap
190 std::vector<float> edges;
191 edges.reserve(sortedID.size() - 1);
193 if (sortedID.size() < 2) {
194 ATH_MSG_WARNING(
"Not enough SP in this track, returning empty edge list!");
197 if (edgeMap.empty()) {
202 for (std::size_t i = 0; i < sortedID.size() - 1; ++i) {
203 uint32_t src = sortedID[i];
204 uint32_t dst = sortedID[i+1];
206 auto srcMap = edgeMap.find(src);
207 if (srcMap == edgeMap.end()) {
208 edges.push_back(-1.f);
212 auto dstMap = srcMap->second.find(dst);
213 if (dstMap == srcMap->second.end()) {
214 edges.push_back(-1.f);
218 float score = dstMap->second;
219 edges.push_back(score);
232 const Trk::Perigee* origPerigee = track.perigeeParameters();
233 double pt = origPerigee->
pT();
235 double eta = std::abs(origPerigee->
eta());
237 double d0 = std::abs(origPerigee->parameters()[
Trk::d0]);
239 double z0 = std::abs(origPerigee->parameters()[
Trk::z0]);
246 int nClusters = nPixels + nStrips;
248 ATH_MSG_DEBUG(
"track params: (pt,eta,d0,z0)=(" << pt <<
"," <<
eta <<
"," << d0 <<
"," << z0
249 <<
") Nclusters:" << nClusters <<
" pixel "<< nPixels
250 <<
" strips " << nStrips);
253 if (nClusters < m_etaDependentCutsSvc->getMinSiHitsAtEta(
eta))
257 if (nPixels < m_etaDependentCutsSvc->getMinPixelHitsAtEta(
eta))
261 if (pt < m_etaDependentCutsSvc->getMinPtAtEta(
eta))
280 const EventContext& ctx,
283 std::vector<const Trk::SpacePoint*> spacePoints;
285 int npixsp(0), nstrip(0), n_overlap(0);
287 spacePoints.push_back(
sp);
291 spacePoints.push_back(
sp);
295 spacePoints.push_back(
sp);
299 <<
" pixel space points, " << nstrip
300 <<
" strips space points, " << n_overlap
301 <<
" overlapping spacepoints, " << spacePoints.size()
308 const EventContext& ctx,
311 std::vector<const Trk::SpacePoint*> spacePoints;
312 if (not containerKey.
empty()){
316 if (container.isValid()){
318 auto spc = container->begin();
319 auto spce = container->end();
320 for(; spc != spce; ++spc){
322 auto sp = spCollection->
begin();
323 auto spe = spCollection->
end();
324 for(;
sp != spe; ++
sp) {
326 spacePoints.push_back(spacePoint);
336 const EventContext& ctx,
340 std::vector<const Trk::SpacePoint*> spacePoints;
342 if (not containerKey.
empty()){
348 spacePoints.push_back(
sp);
357 const EventContext& ctx,
360 std::vector<const Trk::PrepRawData*> allClusters;
372 if (!strip_container.
isValid()) {
377 auto pixcollection = pixcontainer->begin();
378 auto pixcollectionEnd = pixcontainer->end();
379 for (; pixcollection != pixcollectionEnd; ++pixcollection) {
380 if ((*pixcollection)->empty()) {
384 auto const* clusterCollection = (*pixcollection);
385 auto thisCluster = clusterCollection->begin();
386 auto clusterEnd = clusterCollection->end();
387 for (; thisCluster != clusterEnd; ++thisCluster) {
389 allClusters.push_back(cl);
393 auto strip_collection = strip_container->begin();
394 auto strip_collectionEnd = strip_container->end();
395 for (; strip_collection != strip_collectionEnd; ++strip_collection) {
396 if ((*strip_collection)->empty()) {
400 auto const* clusterCollection = (*strip_collection);
401 auto thisCluster = clusterCollection->begin();
402 auto clusterEnd = clusterCollection->end();
403 for (; thisCluster != clusterEnd; ++thisCluster) {
405 allClusters.push_back(cl);
409 ATH_MSG_DEBUG(
"Event " << eventNumber <<
" has " << allClusters.size()
417 const std::vector<uint32_t>& trackIndices,
418 const std::vector<const Trk::SpacePoint*>& allSpacePoints
421 std::vector<const Trk::SpacePoint*> trackCandidate;
422 std::vector<uint32_t> sorted_idx;
423 trackCandidate.reserve(trackIndices.size());
424 sorted_idx.reserve(trackIndices.size());
426 std::vector<std::tuple<double, const Trk::SpacePoint*, uint32_t> > distanceSortedSPs;
430 for (
auto&
id : trackIndices) {
432 if (
id > allSpacePoints.size()) {
434 <<
id <<
" out of range: " << allSpacePoints.size());
442 distanceSortedSPs.push_back(
444 pow(
sp->globalPosition().x(), 2) + pow(
sp->globalPosition().y(), 2),
453 std::sort(distanceSortedSPs.begin(), distanceSortedSPs.end());
456 for (
size_t i = 0; i < distanceSortedSPs.size(); i++) {
457 trackCandidate.push_back(std::get<1>(distanceSortedSPs[i]));
460 sorted_idx.push_back(std::get<2>(distanceSortedSPs[i]));
464 return std::make_pair(trackCandidate, sorted_idx);
467std::vector<const Trk::PrepRawData*> InDet::SiSPGNNTrackMaker :: spacePointsToClusters (
468 const std::vector<const Trk::SpacePoint*>& spacePoints
470 std::vector<const Trk::PrepRawData*> clusters;
472 clusters.push_back(
sp->clusterList().first);
473 if (
sp->clusterList().second !=
nullptr) {
474 clusters.push_back(
sp->clusterList().second);
480std::vector<const Trk::PrepRawData*> InDet::SiSPGNNTrackMaker :: getClusters (
481 const std::vector<std::vector<uint32_t>>& clusterTracks,
482 const std::vector<const Trk::PrepRawData*>& allClusters,
485 std::vector<uint32_t> clusterIndices = clusterTracks[trackNumber];
486 std::vector<const Trk::PrepRawData*> clusters;
487 clusters.reserve(clusterIndices.size());
488 for (uint32_t
id : clusterIndices) {
489 if (
id > allClusters.size()) {
493 clusters.push_back(allClusters[
id]);
532 const EventContext& ctx,
533 std::vector<const Trk::SpacePoint*>& spacePoints,
534 std::vector<const Trk::PrepRawData*>& clusters,
539 int nPIX(0), nStrip(0);
546 ATH_MSG_DEBUG(
"Track " << trackCounter <<
" has " << spacePoints.size()
547 <<
" space points, " << clusters.size()
548 <<
" clusters, " << nPIX <<
" pixel clusters, "
549 << nStrip <<
" strip clusters");
552 if (not
prefitCheck(nPIX, nStrip, clusters.size(), spacePoints.size())) {
553 ATH_MSG_DEBUG(
"Track " << trackCounter <<
" does not pass prefit cuts, skipping");
554 return std::make_tuple(
false, 999,
nullptr);
559 if (trkParameters ==
nullptr) {
561 return std::make_tuple(
false, 999,
nullptr);
564 std::unique_ptr<Trk::Track> track =
fitTrack(ctx, clusters, *trkParameters, trackCounter);
566 if (track ==
nullptr || track->perigeeParameters() ==
nullptr) {
568 <<
" fails the chi2 fit, skipping");
569 return std::make_tuple(
false, 999,
nullptr);
573 if (track->perigeeParameters()->pT() <
m_pTmin) {
575 <<
"with pt = " << track->perigeeParameters()->pT()
576 <<
" has pT too low, skipping track!");
577 return std::make_tuple(
false, 999,
nullptr);
581 if (std::abs(track->perigeeParameters()->eta()) >
m_etamax) {
583 << std::abs(track->perigeeParameters()->eta())
584 <<
" has eta too high, skipping track!");
585 return std::make_tuple(
false, 999,
nullptr);
590 ctx, *track,
false );
594 return std::make_tuple(
true, passTrackCut, std::move(track));
599 const EventContext& ctx,
600 std::vector<const Trk::PrepRawData*> clusters,
608 bool keepOnTrying =
true;
609 std::unique_ptr<Trk::Track> track;
610 while (keepOnTrying) {
611 track =
m_trackFitter->fit(ctx, clusters, initial_params,
false, matEffects);
613 if (track ==
nullptr || track->perigeeParameters() ==
nullptr) {
616 keepOnTrying =
false;
619 keepOnTrying =
false;
623 if (track ==
nullptr || track->perigeeParameters() ==
nullptr) {
625 <<
" fails the first chi2 fit, skipping");
630 if (track->perigeeParameters()->pT() <
m_pTmin) {
632 <<
" fails the first chi2 fit, skipping");
639 while (keepOnTrying) {
640 track =
m_trackFitter->fit(ctx, clusters, origPerigee,
true, matEffects);
643 if (track ==
nullptr || track->trackSummary() ==
nullptr || track->outliersOnTrack()->size()>=3) {
657 <<
" fails the second chi2 fit, skipping");
667MsgStream& InDet::operator <<
676std::ostream& InDet::operator <<
699 out<<
"| Location of output tracks | "
701 out<<
"|----------------------------------------------------------------"
702 <<
"----------------------------------------------------|"
Scalar eta() const
pseudorapidity method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_WARNING(x)
bool msgLvl(const MSG::Level lvl) const
An algorithm that can be simultaneously executed in multiple threads.
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
InDet::SiSPGNNTrackMaker is an algorithm that uses the GNN-based track finding tool to reconstruct tr...
SG::WriteHandleKey< std::vector< std::vector< float > > > m_outputEdgeScoresKey
ToolHandle< ISeedFitter > m_seedFitter
MsgStream & dump(MsgStream &out) const
std::unique_ptr< Trk::Track > fitTrack(const EventContext &ctx, std::vector< const Trk::PrepRawData * > clusters, const Trk::TrackParameters &initial_params, int trackCounter) const
UnsignedIntegerProperty m_minStripClusters
ToolHandle< IGNNTrackReaderTool > m_gnnTrackReader
std::pair< std::vector< const Trk::SpacePoint * >, std::vector< uint32_t > > getSpacePoints(const std::vector< uint32_t > &trackIndices, const std::vector< const Trk::SpacePoint * > &allSpacePoints) const
SG::ReadHandleKey< InDet::PixelClusterContainer > m_ClusterPixelKey
BooleanProperty m_doRecoverFailedFits
SG::ReadHandleKey< InDet::SCT_ClusterContainer > m_ClusterStripKey
UnsignedIntegerProperty m_minPixelClusters
BooleanProperty m_doRecoTrackCuts
BooleanProperty m_areInputClusters
int passEtaDepCuts(const Trk::Track &track) const
SiSPGNNTrackMaker(const std::string &name, ISvcLocator *pSvcLocator)
ServiceHandle< IInDetEtaDependentCutsSvc > m_etaDependentCutsSvc
std::vector< const Trk::PrepRawData * > getClustersInEvent(const EventContext &ctx, int eventNumber) const
std::vector< float > getEdgeScores(const std::vector< uint32_t > &sortedID, const std::unordered_map< int, std::unordered_map< int, float > > &edgeMap) const
SG::WriteHandleKey< TrackCollection > m_outputTracksKey
std::vector< const Trk::PrepRawData * > spacePointsToClusters(const std::vector< const Trk::SpacePoint * > &spacePoints) const
BooleanProperty m_saveEdgeScore
ToolHandle< Trk::ITrackFitter > m_trackFitter
Track Fitter.
SG::ReadHandleKey< SpacePointContainer > m_SpacePointsSCTKey
MsgStream & dumpevent(MsgStream &out) const
virtual StatusCode initialize() override
ToolHandle< IGNNTrackFinder > m_gnnTrackFinder
GNN-based track finding tool that produces track candidates.
bool prefitCheck(unsigned int nPix, unsigned int nStrip, unsigned int nClusters, unsigned int nSpacePoints) const
std::vector< const Trk::SpacePoint * > getSpacePointsInEvent(const EventContext &ctx, int eventNumber) const
UnsignedIntegerProperty m_minClusters
virtual StatusCode execute(const EventContext &ctx) const override
ToolHandle< Trk::IExtendedTrackSummaryTool > m_trackSummaryTool
std::vector< const Trk::PrepRawData * > getClusters(const std::vector< std::vector< uint32_t > > &clusterTracks, const std::vector< const Trk::PrepRawData * > &allClusters, int trackNumber) const
SG::ReadHandleKey< SpacePointOverlapCollection > m_SpacePointsOverlapKey
SG::ReadHandleKey< SpacePointContainer > m_SpacePointsPixelKey
std::tuple< bool, int, std::unique_ptr< Trk::Track > > doFitAndCut(const EventContext &ctx, std::vector< const Trk::SpacePoint * > &spacePoints, std::vector< const Trk::PrepRawData * > &clusters, int &trackCounter) const
Fits a track and applies quality cuts to determine if it should be kept.
MsgStream & dumptools(MsgStream &out) const
Property holding a SG store/key/clid from which a ReadHandle is made.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
bool empty() const
Test if the key is blank.
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
double eta() const
Access method for pseudorapidity - from momentum.
double pT() const
Access method for transverse momentum.
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
ParticleHypothesis
Enumeration for Particle hypothesis respecting the interaction with material.
ParametersBase< TrackParametersDim, Charged > TrackParameters
@ numberOfSCTHits
number of SCT holes
@ numberOfPixelHits
number of pixel layers on track with absence of hits
@ numberOfSCTHoles
number of Holes in both sides of a SCT module
@ numberOfPixelHoles
number of pixels which have a ganged ambiguity.
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.