221 {
222 ACTS_DEBUG("TrackFindingGNNAlg::execute() - begin");
223
224 std::optional<Athena::Chrono>
timer;
226
227 Acts::GeometryContext gctx =
229 Acts::MagneticFieldContext mctx =
232
234
235
238 const auto &pixelSPContainer = *pixelSPHandle.cptr();
239
242 const auto &stripSPContainer = *stripSPHandle.cptr();
243
244 auto stripSPOVHandle =
247 const auto &stripSPOVContainer = *stripSPOVHandle.cptr();
248
249 constexpr std::size_t nFeatures = 12;
250 std::size_t nSP = pixelSPContainer.size() + stripSPContainer.size() +
251 stripSPOVContainer.size();
252
253 ACTS_DEBUG("Number spacepoints: "
254 << nSP << " (" << "pixel: " << pixelSPContainer.size() << ", "
255 << "strip: " << stripSPContainer.size() << ", "
256 << "strip overlap: " << stripSPOVContainer.size() << ")");
257
258
260
261 std::vector<std::uint64_t> moduleIds;
262 moduleIds.reserve(nSP);
263 std::vector<const xAOD::SpacePoint *> allSPPtrs;
264 allSPPtrs.reserve(nSP);
265 std::vector<Acts::GeometryIdentifier> geoIds, sortedGeoIds(nSP);
266 geoIds.reserve(nSP);
267
268 std::size_t skipped = 0;
269 for (const auto &spc :
270 {pixelSPContainer, stripSPContainer, stripSPOVContainer}) {
271 for (
auto sp : spc) {
272 auto cl1 =
sp->measurements().front();
273 auto geoIdCl1 =
276
277 if (
sp->measurements().size() == 2) {
278 auto cl2 =
sp->measurements().at(1);
280
281 auto overlapFlag =
286
287 if (overlapFlag == 2 || overlapFlag == 3) {
288 skipped++;
289 ACTS_VERBOSE("Skip phi overlap spacepoint (flag=" << overlapFlag
290 << ")");
291 continue;
292 }
293 }
294
295 geoIds.push_back(geoIdCl1);
296 moduleIds.push_back(atlasIdCl1.get_compact());
297 allSPPtrs.push_back(
sp);
298 }
299 }
300
301 ACTS_DEBUG("Skipped " << skipped << " SPs because of phi overlap");
302 nSP = allSPPtrs.size();
303 ACTS_DEBUG("Keep " << nSP << " SPs for feature creation");
304
306
307 std::vector<std::size_t> idxs(nSP);
308 std::iota(idxs.begin(), idxs.end(), 0);
309
310 std::ranges::sort(
311 idxs, [&](
auto a,
auto b) {
return moduleIds.at(
a) < moduleIds.at(b); });
312 std::ranges::sort(moduleIds);
313
314 std::vector<float> features(nFeatures * nSP);
315 std::vector<boost::container::static_vector<Acts::SourceLink, 2>> sourceLinks(
316 nSP);
317 std::vector<int>
id(nSP);
318
319 for (
auto k = 0ul;
k < nSP;
k++) {
322
323 std::span<float>
f(features.data() + k * nFeatures, nFeatures);
324 const auto &
sp = *allSPPtrs.at(i);
325
326 using namespace Acts::VectorHelpers;
327 using namespace Acts::AngleHelpers;
328
329 Acts::Vector3 spp{
sp.x(),
sp.y(),
sp.z()};
330
331 if (
sp.measurements().size() == 1) {
332 for (
auto j = 0ul;
j < nFeatures;
j += 4) {
333 f[
j + 0] =
perp(spp) / 1000.f;
334 f[
j + 1] =
phi(spp) / std::numbers::pi_v<float>;
335 f[
j + 2] =
sp.z() / 1000.f;
337 }
338 } else {
340 f[
j + 0] =
perp(spp) / 1000.f;
341 f[
j + 1] =
phi(spp) / std::numbers::pi_v<float>;
342 f[
j + 2] =
sp.z() / 1000.f;
344
345 for (
auto m :
sp.measurements()) {
347 auto gp =
cl->globalPosition();
349 f[
j + 0] =
perp(gp) / 1000.f;
350 f[
j + 1] =
phi(gp) / std::numbers::pi_v<float>;
351 f[
j + 2] = gp.z() / 1000.f;
353 }
354 }
355
358 }
359
360 sortedGeoIds.at(k) = geoIds.at(i);
361 }
362
365
366 m_gpuInstanceCount->acquire();
369 m_gpuInstanceCount->release();
370
371 ACTS_DEBUG(
"Have " <<
candidates.size() <<
" candidates after GNN");
372
373
374 auto candidateSelector = [&](
const std::vector<int> &
c) {
375 bool tooFewMeasurements = std::accumulate(
c.begin(),
c.end(), 0ul, [&](
auto sum,
auto spi) {
376 return sum + allSPPtrs.at(spi)->measurements().size();
378 bool noPixelHits = !std::ranges::any_of(c, [&](auto spi) { return allSPPtrs.at(spi)->measurements().size() == 1; });
379 return tooFewMeasurements || noPixelHits;
380 };
381
386
389
390 Acts::VectorTrackContainer trackBackend;
391 Acts::VectorMultiTrajectory trackStateBackend;
392 constexpr std::size_t nTracksExpected = 3000;
393 trackBackend.reserve(nTracksExpected);
394 trackStateBackend.reserve(nTracksExpected * 30);
396
397
398 ActsTrk::SeedContainer seedContainer;
399
400 auto makeSeedFromCandidate = [&](const std::vector<int> &cand) -> std::optional<boost::container::small_vector<const xAOD::SpacePoint*, 3>> {
401
402 boost::container::small_vector<const xAOD::SpacePoint*, 3> picked;
403 if (cand.empty()) return std::nullopt;
405 Acts::Vector3
v{
sp->x(),
sp->y(),
sp->z()};
407 };
409 picked.push_back(last);
410 for (std::size_t i = 1;
i < cand.size() && picked.size() < 3; ++i) {
413 picked.push_back(
sp);
415 }
416 }
417 if (picked.size() < 3) return std::nullopt;
418 return picked;
419 };
420
421 auto retrieveSurface = [&](
const ActsTrk::Seed&
seed,
bool useTopSp) ->
const Acts::Surface& {
425 if (!surface) {
426 throw std::runtime_error("retrieveSurface: no Acts surface for GeometryIdentifier " + std::to_string(geoId.value()));
427 }
428 return *surface;
429 };
430
432 return Acts::fastHypot(
sp->x(),
sp->y(),
sp->z());
433 };
434
435 for (const auto &cand : candidates) {
436 auto pickedOpt = makeSeedFromCandidate(cand);
437 if (!pickedOpt.has_value()) continue;
438
439 auto picked = *pickedOpt;
442 return R_of(a) < R_of(b);
443 });
445 ActsTrk::SpacePointRange(picked.data(), picked.size()), 0.f, 0.f);
446
448 seed, true, gctx, mctx, retrieveSurface);
449 if (!initialParamsOpt.has_value()) continue;
450
451 boost::container::small_vector<const xAOD::SpacePoint*, 16> sortedSP;
452 sortedSP.reserve(cand.size());
453 for (int spi : cand) sortedSP.push_back(allSPPtrs.at(spi));
454 std::sort(sortedSP.begin(), sortedSP.end(),
456 return R_of(a) < R_of(b);
457 });
458
459 std::vector<const xAOD::UncalibratedMeasurement*> measList;
460 measList.reserve(sortedSP.size() * 2);
463 measList.push_back(m);
464 }
465 }
466
467 auto fitted =
m_fitterTool->fit(measList, *initialParamsOpt, gctx, mctx, cctx);
468 if (fitted) {
469 for (auto track : *fitted) {
470 auto newTrack = tracks.makeTrack();
471 newTrack.copyFrom(track);
472 }
473 }
474 }
475
476 ACTS_DEBUG(
"After track fit: " << tracks.size() <<
" / " <<
candidates.size()
477 << " successfull");
478
479
480 if (
candidates.size() == 1 && tracks.size() == 1) {
481 const auto &
t = *tracks.begin();
482 ACTS_DEBUG(
"Single particle case: " <<
candidates.front().size() <<
" -> "
484 << " measurements");
485 }
486
489
490 Acts::VectorTrackContainer selTrackBackend;
491 selTrackBackend.reserve(trackBackend.size());
493
495 for (auto track : tracks) {
497 auto newTrack = selectedTracks.makeTrack();
498
499
500 newTrack.copyFrom(track);
501 }
502 }
503
504 ACTS_DEBUG(
"GNN cand: " <<
candidates.size() <<
", fitted: " << tracks.size()
505 << ", selected: " << selectedTracks.size());
506
507
508 Acts::ConstVectorTrackContainer constTrackBackend(std::move(selTrackBackend));
509 Acts::ConstVectorMultiTrajectory constTrackStateBackend(std::move(trackStateBackend));
510 std::unique_ptr<ActsTrk::TrackContainer> constTracksContainer
511 = std::make_unique<ActsTrk::TrackContainer>(std::move(constTrackBackend), std::move(constTrackStateBackend) );
512
515 ATH_CHECK(trackContainerHandle.
record(std::move(constTracksContainer)));
516
517 return StatusCode::SUCCESS;
518}
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
ToolHandle< ActsTrk::IFitterTool > m_fitterTool
Acts::TrackSelector::EtaBinnedConfig m_trackSelectorConfig
ToolHandle< ActsTrk::ITrackParamsEstimationTool > m_paramEstimationTool
std::unique_ptr< ActsPlugins::GnnPipeline > m_gnnPipeline
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodStripSpacePointOverlapContainerKey
SG::WriteHandleKey< ActsTrk::TrackContainer > m_trackContainerKey
PublicToolHandle< ITrackingGeometryTool > m_trackingGeometryTool
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodPixelSpacePointContainerKey
Gaudi::Property< unsigned int > m_minCandidateMeasurements
ToolHandle< ActsTrk::IExtrapolationTool > m_extrapolationTool
Gaudi::Property< double > m_minDeltaR
ServiceHandle< IChronoStatSvc > m_chronoSvc
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodStripSpacePointContainerKey
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)
Acts::CalibrationContext getCalibrationContext(const EventContext &ctx)
The Acts::Calibration context is piped through the Acts fitters to (re)calibrate the Acts::SourceLink...
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)