219 {
220 ACTS_DEBUG("TrackFindingGNNAlg::execute() - begin");
221
222 std::optional<Athena::Chrono>
timer;
224
225 const Acts::GeometryContext gctx =
m_ctxProvider.getGeometryContext(ctx);
226 const Acts::MagneticFieldContext mctx =
m_ctxProvider.getMagneticFieldContext(ctx);
227 const Acts::CalibrationContext cctx =
m_ctxProvider.getCalibrationContext(ctx);
228
230
231
234 const auto &pixelSPContainer = *pixelSPHandle.cptr();
235
238 const auto &stripSPContainer = *stripSPHandle.cptr();
239
240 auto stripSPOVHandle =
243 const auto &stripSPOVContainer = *stripSPOVHandle.cptr();
244
245 constexpr std::size_t nFeatures = 12;
246 std::size_t nSP = pixelSPContainer.size() + stripSPContainer.size() +
247 stripSPOVContainer.size();
248
249 ACTS_DEBUG("Number spacepoints: "
250 << nSP << " (" << "pixel: " << pixelSPContainer.size() << ", "
251 << "strip: " << stripSPContainer.size() << ", "
252 << "strip overlap: " << stripSPOVContainer.size() << ")");
253
254
256
257 std::vector<std::uint64_t> moduleIds;
258 moduleIds.reserve(nSP);
259 std::vector<const xAOD::SpacePoint *> allSPPtrs;
260 allSPPtrs.reserve(nSP);
261 std::vector<Acts::GeometryIdentifier> geoIds, sortedGeoIds(nSP);
262 geoIds.reserve(nSP);
263
264 std::size_t skipped = 0;
265 for (const auto &spc :
266 {pixelSPContainer, stripSPContainer, stripSPOVContainer}) {
267 for (
auto sp : spc) {
268 auto cl1 =
sp->measurements().front();
269 auto geoIdCl1 =
272
273 if (
sp->measurements().size() == 2) {
274 auto cl2 =
sp->measurements().at(1);
276
277 auto overlapFlag =
282
283 if (overlapFlag == 2 || overlapFlag == 3) {
284 skipped++;
285 ACTS_VERBOSE("Skip phi overlap spacepoint (flag=" << overlapFlag
286 << ")");
287 continue;
288 }
289 }
290
291 geoIds.push_back(geoIdCl1);
292 moduleIds.push_back(atlasIdCl1.get_compact());
293 allSPPtrs.push_back(
sp);
294 }
295 }
296
297 ACTS_DEBUG("Skipped " << skipped << " SPs because of phi overlap");
298 nSP = allSPPtrs.size();
299 ACTS_DEBUG("Keep " << nSP << " SPs for feature creation");
300
302
303 std::vector<std::size_t> idxs(nSP);
304 std::iota(idxs.begin(), idxs.end(), 0);
305
306 std::ranges::sort(
307 idxs, [&](
auto a,
auto b) {
return moduleIds.at(
a) < moduleIds.at(b); });
308 std::ranges::sort(moduleIds);
309
310 std::vector<float> features(nFeatures * nSP);
311 std::vector<boost::container::static_vector<Acts::SourceLink, 2>> sourceLinks(
312 nSP);
313 std::vector<int>
id(nSP);
314
315 for (
auto k = 0ul;
k < nSP;
k++) {
318
319 std::span<float>
f(features.data() + k * nFeatures, nFeatures);
320 const auto &
sp = *allSPPtrs.at(i);
321
322 using namespace Acts::VectorHelpers;
323 using namespace Acts::AngleHelpers;
324
325 Acts::Vector3 spp{
sp.x(),
sp.y(),
sp.z()};
326
327 if (
sp.measurements().size() == 1) {
328 for (
auto j = 0ul;
j < nFeatures;
j += 4) {
329 f[
j + 0] =
perp(spp) / 1000.f;
330 f[
j + 1] =
phi(spp) / std::numbers::pi_v<float>;
331 f[
j + 2] =
sp.z() / 1000.f;
333 }
334 } else {
336 f[
j + 0] =
perp(spp) / 1000.f;
337 f[
j + 1] =
phi(spp) / std::numbers::pi_v<float>;
338 f[
j + 2] =
sp.z() / 1000.f;
340
341 for (
auto m :
sp.measurements()) {
343 auto gp =
cl->globalPosition();
345 f[
j + 0] =
perp(gp) / 1000.f;
346 f[
j + 1] =
phi(gp) / std::numbers::pi_v<float>;
347 f[
j + 2] = gp.z() / 1000.f;
349 }
350 }
351
354 }
355
356 sortedGeoIds.at(k) = geoIds.at(i);
357 }
358
361
362 m_gpuInstanceCount->acquire();
365 m_gpuInstanceCount->release();
366
367 ACTS_DEBUG(
"Have " <<
candidates.size() <<
" candidates after GNN");
368
369
370 auto candidateSelector = [&](
const std::vector<int> &
c) {
371 bool tooFewMeasurements = std::accumulate(
c.begin(),
c.end(), 0ul, [&](
auto sum,
auto spi) {
372 return sum + allSPPtrs.at(spi)->measurements().size();
374 bool noPixelHits = !std::ranges::any_of(c, [&](auto spi) { return allSPPtrs.at(spi)->measurements().size() == 1; });
375 return tooFewMeasurements || noPixelHits;
376 };
377
382
385
386 Acts::VectorTrackContainer trackBackend;
387 Acts::VectorMultiTrajectory trackStateBackend;
388 constexpr std::size_t nTracksExpected = 3000;
389 trackBackend.reserve(nTracksExpected);
390 trackStateBackend.reserve(nTracksExpected * 30);
392
393
394 ActsTrk::SeedContainer seedContainer;
395
396 auto makeSeedFromCandidate = [&](const std::vector<int> &cand) -> std::optional<boost::container::small_vector<const xAOD::SpacePoint*, 3>> {
397
398 boost::container::small_vector<const xAOD::SpacePoint*, 3> picked;
399 if (cand.empty()) return std::nullopt;
401 Acts::Vector3
v{
sp->x(),
sp->y(),
sp->z()};
403 };
405 picked.push_back(last);
406 for (std::size_t i = 1;
i < cand.size() && picked.size() < 3; ++i) {
409 picked.push_back(
sp);
411 }
412 }
413 if (picked.size() < 3) return std::nullopt;
414 return picked;
415 };
416
417 auto retrieveSurface = [&](
const ActsTrk::Seed&
seed,
bool useTopSp) ->
const Acts::Surface& {
421 if (!surface) {
422 throw std::runtime_error("retrieveSurface: no Acts surface for GeometryIdentifier " + std::to_string(geoId.value()));
423 }
424 return *surface;
425 };
426
428 return Acts::fastHypot(
sp->x(),
sp->y(),
sp->z());
429 };
430
431 for (const auto &cand : candidates) {
432 auto pickedOpt = makeSeedFromCandidate(cand);
433 if (!pickedOpt.has_value()) continue;
434
435 auto picked = *pickedOpt;
438 return R_of(a) < R_of(b);
439 });
441 ActsTrk::SpacePointRange(picked.data(), picked.size()), 0.f, 0.f);
442
444 seed, true, gctx, mctx, retrieveSurface);
445 if (!initialParamsOpt.has_value()) continue;
446
447 boost::container::small_vector<const xAOD::SpacePoint*, 16> sortedSP;
448 sortedSP.reserve(cand.size());
449 for (int spi : cand) sortedSP.push_back(allSPPtrs.at(spi));
450 std::sort(sortedSP.begin(), sortedSP.end(),
452 return R_of(a) < R_of(b);
453 });
454
455 std::vector<const xAOD::UncalibratedMeasurement*> measList;
456 measList.reserve(sortedSP.size() * 2);
459 measList.push_back(m);
460 }
461 }
462
463 auto fitted =
m_fitterTool->fit(measList, *initialParamsOpt, gctx, mctx, cctx);
464 if (fitted) {
465 for (auto track : *fitted) {
466 auto newTrack = tracks.makeTrack();
467 newTrack.copyFrom(track);
468 }
469 }
470 }
471
472 ACTS_DEBUG(
"After track fit: " << tracks.size() <<
" / " <<
candidates.size()
473 << " successfull");
474
475
476 if (
candidates.size() == 1 && tracks.size() == 1) {
477 const auto &
t = *tracks.begin();
478 ACTS_DEBUG(
"Single particle case: " <<
candidates.front().size() <<
" -> "
480 << " measurements");
481 }
482
485
486 Acts::VectorTrackContainer selTrackBackend;
487 selTrackBackend.reserve(trackBackend.size());
489
491 for (auto track : tracks) {
493 auto newTrack = selectedTracks.makeTrack();
494
495
496 newTrack.copyFrom(track);
497 }
498 }
499
500 ACTS_DEBUG(
"GNN cand: " <<
candidates.size() <<
", fitted: " << tracks.size()
501 << ", selected: " << selectedTracks.size());
502
503
504 Acts::ConstVectorTrackContainer constTrackBackend(std::move(selTrackBackend));
505 Acts::ConstVectorMultiTrajectory constTrackStateBackend(std::move(trackStateBackend));
506 std::unique_ptr<ActsTrk::TrackContainer> constTracksContainer
507 = std::make_unique<ActsTrk::TrackContainer>(std::move(constTrackBackend), std::move(constTrackStateBackend) );
508
511 ATH_CHECK(trackContainerHandle.
record(std::move(constTracksContainer)));
512
513 return StatusCode::SUCCESS;
514}
Scalar eta() const
pseudorapidity method
Scalar perp() const
perp method - perpendicular length
Scalar phi() const
phi method
#define ATH_CHECK
Evaluate an expression and check for errors.
const SCT_ID * m_stripIdHelper
Acts::TrackSelector::EtaBinnedConfig m_trackSelectorConfig
ToolHandle< ITrackParamsEstimationTool > m_paramEstimationTool
std::unique_ptr< ActsPlugins::GnnPipeline > m_gnnPipeline
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodStripSpacePointOverlapContainerKey
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodPixelSpacePointContainerKey
Gaudi::Property< unsigned int > m_minCandidateMeasurements
ServiceHandle< ActsTrk::ITrackingGeometrySvc > m_trackingGeometrySvc
SG::WriteHandleKey< TrackContainer > m_trackContainerKey
Gaudi::Property< double > m_minDeltaR
ContextUtility m_ctxProvider
Utility to fetch the geometry, magnetic field and calibration context in the event.
ServiceHandle< IChronoStatSvc > m_chronoSvc
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodStripSpacePointContainerKey
ToolHandle< IFitterTool > m_fitterTool
Gaudi::Property< int > m_cudaDeviceIndex
static Acts::SourceLink pack(const Ptr_t &measurement)
Pack the measurement type pointer to an Acts::SourceLink including the intermediate conversion into a...
unsigned long long value_type
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
Acts::TrackContainer< Acts::VectorTrackContainer, Acts::VectorMultiTrajectory > RecoTrackContainer
Acts::GeometryIdentifier getSurfaceGeometryIdOfMeasurement(const DetectorElementToActsGeometryIdMap &detector_element_to_geoid, const xAOD::UncalibratedMeasurement &measurement)
int compute_overlap_SP_flag(const int &eta_module_cl1, const int &phi_module_cl1, const int &eta_module_cl2, const int &phi_module_cl2)
SG::ReadCondHandle< T > makeHandle(const SG::ReadCondHandleKey< T > &key, const EventContext &ctx=Gaudi::Hive::currentContext())
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)
cl
print [x.__class__ for x in toList(dqregion.getSubRegions()) ]
timer(name, disabled=False)
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
DataModel_detail::iterator< DVL > remove_if(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end, Predicate pred)
Specialization of remove_if for DataVector/List.
StripCluster_v1 StripCluster
Define the version of the strip cluster class.
UncalibratedMeasurement_v1 UncalibratedMeasurement
Define the version of the uncalibrated measurement class.
Seed push_back(SpacePointRange spacePoints, float quality, float vertexZ)