diff --git a/Detectors/GenericTracking/study/include/GenericTracking/Study/Alice/Tracks.h b/Detectors/GenericTracking/study/include/GenericTracking/Study/Alice/Tracks.h index 04a434bb60bc7..e5ecdc29ffa0b 100644 --- a/Detectors/GenericTracking/study/include/GenericTracking/Study/Alice/Tracks.h +++ b/Detectors/GenericTracking/study/include/GenericTracking/Study/Alice/Tracks.h @@ -11,7 +11,7 @@ /// \file Tracks.h /// \brief O2's ITS and MFT tracks, as the clusters they hold, their MC labels' -/// fake flags and, for the ITS, their states at the innermost cluster. +/// fake flags, fitted kinematics and ITS innermost-cluster states. /// /// All a comparison with our tracks needs, because both trackers read the /// same clusters in the same order (Clusters.h): a track is its list of O2 @@ -34,6 +34,8 @@ class TTree; namespace o2::gtrk::study { +struct TrackKinematics; + /// One timeframe's tracks, each as its cluster indices. using TrackClusters = std::vector>; @@ -56,6 +58,10 @@ class TrackFile /// Whether each track's own MC label is flagged fake, as O2's tracker set it. [[nodiscard]] std::vector fakeFlags(long long entry); + /// Fitted inward kinematics in the same order as read(). Native outputs + /// contain accepted tracks; fitValid checks finite, positive momentum. + [[nodiscard]] std::vector kinematics(long long entry); + /// Every ITS track's state at its innermost cluster. Throws /// std::runtime_error for the MFT, whose tracks are parameterised in z. [[nodiscard]] std::vector innerStates(long long entry); diff --git a/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h b/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h index ff6cb3ea1e004..57dcc235b28a2 100644 --- a/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h +++ b/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h @@ -25,7 +25,9 @@ /// (StageFunnel.h); its PDG codes, decay geometry, layers, /// the funnel's verdict and the other tracker's /// (`QAParticle`) -/// tracks: the fit's pT, eta and phi, its iteration, fake and clone +/// tracks: the fit's pT, eta and phi, its iteration, fake and clone; +/// fit status, generated pT/eta, origin and relative inverse-pT +/// residual for matched tracks, without requiring hit files /// /// Every efficiency follows from the two; for iteration i, `bit` = 1 << i: /// @@ -95,6 +97,7 @@ struct TrackKinematics { std::uint8_t iteration = 0U; float chi2 = 0.F; std::int32_t ndf = 0; + bool fitValid = true; ///< the refit's propagation status is kOk }; /// Iterations a particle's `foundBy` holds. @@ -147,7 +150,8 @@ struct ParticleExtras { /// own part, a FitLoss, kOther if the fit takes it; kQANotObserved /// otherwise. std::uint8_t raisedFit = kQANotObserved; - std::uint8_t outwardLayers = 0U; ///< the layers of its MC hits' outward leg (`outwardLayers`) + std::uint8_t outwardLayers = 0U; ///< the layers of its MC hits' outward leg (`outwardLayers`) + std::uint64_t expectedLayerMask = 0U; ///< full-length reference: detector layers, or ALICE3 outward MC traversal }; /// A findable particle of the `particles` tree. @@ -231,6 +235,14 @@ struct QATrack { bool clone = false; float chi2 = 0.F; ///< the fit's; TMath::Prob(chi2, ndf) is its probability std::int32_t ndf = 0; + bool fitValid = true; + std::uint64_t expectedLayerMask = 0U; ///< zero when a full-length reference is unavailable + std::uint64_t layerMask = 0U; ///< distinct reconstructed detector layers + std::uint64_t mcParticle = kNoTrackParticle; ///< source/event/track identity, local to this timeframe + float ptGen = std::numeric_limits::quiet_NaN(); ///< generated momentum at production [GeV] + float etaGen = std::numeric_limits::quiet_NaN(); + std::uint8_t origin = 0xFFU; ///< ParticleOrigin; unknown for unlabelled tracks + float invPtResidual = std::numeric_limits::quiet_NaN(); ///< (1/pT - 1/pTgen) * pTgen }; /// A good track against its particle's MC crossing of its innermost @@ -343,13 +355,18 @@ using SeededBy = std::unordered_map; /// \a qa's records into \a output's `particles` and `tracks` trees. void stream(o2::utils::TreeStreamRedirector& output, const TrackingQA& qa); +/// Write the same track schema for a comparison tracker, without duplicating +/// the generated-particle tree. Truth and clone rules are those of trackingQA(). +void streamTracks(o2::utils::TreeStreamRedirector& output, std::span tracks, const char* tree = "tracks"); + /// \a entries into \a output's `resolution` tree. void streamResolution(o2::utils::TreeStreamRedirector& output, std::span entries); /// The run's one entry of the `run` tree: what the macro's axes need, /// \a settings' findable rules (`minLayers`, `minPtGeV`, `minOutwardLayers`) -/// for its definitions, and the run's MC \a events. -void streamRun(o2::utils::TreeStreamRedirector& output, int iterations, EtaRange eta, const ComparisonSettings& settings, std::size_t events); +/// for its definitions, and the run's MC \a events. A nonzero fullLengthMask +/// requires all detector layers; zero denotes per-particle MC traversal masks. +void streamRun(o2::utils::TreeStreamRedirector& output, int iterations, EtaRange eta, const ComparisonSettings& settings, std::size_t events, std::uint64_t fullLengthMask = 0U); } // namespace o2::gtrk::study diff --git a/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx b/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx index 48088f1403cf7..3449b52350cd4 100644 --- a/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx +++ b/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx @@ -9,6 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +#include "DataFormatsITSMFT/ClusterID.h" #include "GenericTracking/Study/MCVertexReaderSpec.h" #include "GenericTracking/Study/Digest.h" #include "GenericTracking/Study/MeasurementSinkSpec.h" @@ -3237,7 +3238,7 @@ class MeasurementSinkTask : public Task const std::string name = mQAOutput->GetFile()->GetName(); if (mInterface.configured()) { // The MC events: the signal source's, as the truth is. - streamRun(*mQAOutput, workflow::iterationCount(recoParam()), etaReach(PlacedEnvelopes::all()), comparisonSettings(), mKinematics ? mKinematics->getNEvents(0) : 0U); + streamRun(*mQAOutput, workflow::iterationCount(recoParam()), etaReach(PlacedEnvelopes::all()), comparisonSettings(), mKinematics ? mKinematics->getNEvents(0) : 0U, std::is_same_v ? 0U : (std::uint64_t{1} << Detector::nLayers) - 1U); } mQAOutput->Close(); mQAOutput.reset(); @@ -3347,6 +3348,7 @@ class MeasurementSinkTask : public Task .iteration = view.iteration[track], .chi2 = view.chi2[track], .ndf = view.ndf[track], + .fitValid = view.status[track] == PropagationStatus::kOk, }; // A good track against its particle's crossing, with the hits. const auto particle = goodParticle(rows, ours[track]); @@ -3397,6 +3399,7 @@ class MeasurementSinkTask : public Task // Each signal particle's outward leg in its MC hits (TrackingQA.h), its // event's segments ordered by track once. std::unordered_map outward; + std::unordered_map outwardMasks; std::unordered_map> legMomenta; // --species-oracle if (mHits) { std::unordered_map> byEvent; @@ -3427,6 +3430,13 @@ class MeasurementSinkTask : public Task const auto& truth = particles.at(id); const auto leg = outwardLeg(crossings, {static_cast(truth.vx), static_cast(truth.vy), static_cast(truth.vz)}); outward[id] = outwardLayers(leg); + if constexpr (std::is_same_v) { + if (outward[id] >= comparisonSettings().minLayers) { + for (const auto& crossing : leg) { + outwardMasks[id] |= std::uint64_t{1} << crossing.layer; + } + } + } if (mSpeciesOracle && (std::abs(truth.pdg) == 11 || (std::abs(truth.pdg) == 211 && id % 16U == 0U))) { legMomenta[id] = momentaAlongLeg(leg); } @@ -3441,6 +3451,13 @@ class MeasurementSinkTask : public Task } const auto& helix = label->second.helix; auto& extra = extras[id]; + if constexpr (std::is_same_v) { + if (const auto mask = outwardMasks.find(id); mask != outwardMasks.end()) { + extra.expectedLayerMask = mask->second; + } + } else { + extra.expectedLayerMask = (std::uint64_t{1} << Detector::nLayers) - 1U; + } if (const auto leg = outward.find(id); leg != outward.end()) { extra.outwardLayers = leg->second; } @@ -3479,6 +3496,11 @@ class MeasurementSinkTask : public Task } } stream(*mQAOutput, trackingQA(rows, particles, ours, kinematics, funnel != nullptr ? funnel->seededBy() : SeededBy{}, extras, theirs, workflow::iterationCount(recoParam()), comparisonSettings())); + if (theirs != nullptr) { + const auto nativeKinematics = mO2Tracks->kinematics(static_cast(mTimeframe)); + const auto nativeQA = trackingQA(rows, particles, *theirs, nativeKinematics, {}, extras, nullptr, 1, comparisonSettings()); + streamTracks(*mQAOutput, nativeQA.tracks, "tracksO2"); + } streamResolution(*mQAOutput, resolution); streamOracle(*mQAOutput, oracle); if (!legMomenta.empty()) { @@ -3638,7 +3660,7 @@ class MeasurementSinkTask : public Task tracksOfRow[view.clusters[track].index[i]].push_back(static_cast(track)); } } - const auto theirs = mO2Tracks->read(static_cast(mTimeframe)); + const auto theirs = o2Clusters(frame); constexpr auto kMaxClusters = static_cast(Detector::nLayers) + 1U; std::array o2Total{}; std::array o2Found{}; @@ -3842,12 +3864,36 @@ class MeasurementSinkTask : public Task return rowOfCluster; } + /// Decode O2's per-layer references to the converter's dense cluster order. + /// Monolithic references (stream zero) already are dense indices. + [[nodiscard]] TrackClusters o2Clusters(const Frame& frame) const + { + auto tracks = mO2Tracks->read(static_cast(mTimeframe)); + std::array::nLayers + 1U> offsets{}; + for (Frame::Index row = 0U; row < frame.size(); ++row) { + ++offsets[mInterface.geometry().sensor(frame.measurement(row).surface).layer + 1U]; + } + std::partial_sum(offsets.begin(), offsets.end(), offsets.begin()); + for (auto& track : tracks) { + for (auto& cluster : track) { + const auto layer = static_cast(o2::itsmft::clusID2Layer(static_cast(cluster))); + if (layer == 0U) { + continue; + } + const auto index = static_cast(o2::itsmft::clusID2Index(static_cast(cluster))); + cluster = layer + 1U < offsets.size() && index < offsets[layer + 1U] - offsets[layer] + ? offsets[layer] + index : std::numeric_limits::max(); + } + } + return tracks; + } + /// O2's tracks of this timeframe as our measurement rows: O2's tracks name - /// O2's cluster indices, which are our measurements' clusterIndex. + /// cluster references, decoded to our measurements' clusterIndex. [[nodiscard]] TrackRows o2TrackRows(const Frame& frame) const { const auto rowOfCluster = rowsOfClusters(frame); - const auto clusters = mO2Tracks->read(static_cast(mTimeframe)); + const auto clusters = o2Clusters(frame); TrackRows theirs; theirs.reserve(clusters.size()); for (const auto& track : clusters) { diff --git a/Detectors/GenericTracking/study/src/TrackingQA.cxx b/Detectors/GenericTracking/study/src/TrackingQA.cxx index 1568cfb075124..518958c82a8fa 100644 --- a/Detectors/GenericTracking/study/src/TrackingQA.cxx +++ b/Detectors/GenericTracking/study/src/TrackingQA.cxx @@ -221,8 +221,24 @@ TrackingQA trackingQA(std::span rows, const ParticleTable& truth .fake = particle == kNoTrackParticle, .chi2 = k.chi2, .ndf = k.ndf, + .fitValid = k.fitValid, + .mcParticle = particle, }; + for (const auto row : tracks[track]) { + record.layerMask |= std::uint64_t{1} << rows[row].layer; + } if (!record.fake) { + if (const auto extra = extras.find(particle); extra != extras.end()) { + record.expectedLayerMask = extra->second.expectedLayerMask; + } + if (const auto generated = truth.find(particle); generated != truth.end()) { + record.ptGen = generated->second.pt; + record.etaGen = generated->second.eta; + record.origin = static_cast(generated->second.origin); + if (record.fitValid && std::isfinite(k.pt) && k.pt > 0.F && std::isfinite(record.ptGen) && record.ptGen > 0.F) { + record.invPtResidual = static_cast((static_cast(record.ptGen) / static_cast(k.pt)) - 1.); + } + } auto& mask = foundBy[particle]; record.clone = mask != 0U; mask |= std::uint64_t{1} << k.iteration; @@ -279,8 +295,13 @@ void stream(o2::utils::TreeStreamRedirector& output, const TrackingQA& qa) const auto& e = p.extras; output << "particles" << "pt=" << p.pt << "eta=" << p.eta << "phi=" << p.phi << "origin=" << p.origin << "foundBy=" << p.foundBy << "seededBy=" << p.seededBy << "pdg=" << p.pdg << "motherPdg=" << p.motherPdg << "vr=" << e.vr << "vz=" << p.vz << "d0=" << e.d0 << "z0=" << e.z0 << "openingAngle=" << e.openingAngle << "motherMiss=" << e.motherMiss << "layers=" << p.layers << "lostAt=" << e.lostAt << "pairLoss=" << e.pairLoss << "roadLoss=" << e.roadLoss << "fitLoss=" << e.fitLoss << "fitRoadLayers=" << e.fitRoadLayers << "roadLosses=" << e.roadLosses << "fitLosses=" << e.fitLosses << "raisedFit=" << e.raisedFit << "foundByO2=" << p.foundByO2 << "findable=" << p.findable << "outwardLayers=" << e.outwardLayers << "\n"; } - for (const auto& t : qa.tracks) { - output << "tracks" << "pt=" << t.pt << "eta=" << t.eta << "phi=" << t.phi << "iteration=" << t.iteration << "fake=" << t.fake << "clone=" << t.clone << "chi2=" << t.chi2 << "ndf=" << t.ndf << "\n"; + streamTracks(output, qa.tracks); +} + +void streamTracks(o2::utils::TreeStreamRedirector& output, std::span tracks, const char* tree) +{ + for (const auto& t : tracks) { + output << tree << "pt=" << t.pt << "eta=" << t.eta << "phi=" << t.phi << "iteration=" << t.iteration << "fake=" << t.fake << "clone=" << t.clone << "chi2=" << t.chi2 << "ndf=" << t.ndf << "expectedLayerMask=" << t.expectedLayerMask << "layerMask=" << t.layerMask << "fitValid=" << t.fitValid << "mcParticle=" << t.mcParticle << "ptGen=" << t.ptGen << "etaGen=" << t.etaGen << "origin=" << t.origin << "dInvPtRel=" << t.invPtResidual << "\n"; } } @@ -306,6 +327,6 @@ void streamMomenta(o2::utils::TreeStreamRedirector& output, std::span #include +#include #include #include +#include namespace o2::gtrk::study { @@ -63,6 +66,34 @@ TrackClusters readEntry(TTree& tree, long long entry, const char* trackBranch, c return result; } +template +std::vector readKinematics(TTree& tree, long long entry, const char* branch) +{ + std::vector tracks; + auto* pointer = &tracks; + tree.SetBranchAddress(branch, &pointer); + tree.GetEntry(entry); + tree.ResetBranchAddresses(); + std::vector result; + result.reserve(tracks.size()); + for (const auto& track : tracks) { + const auto pt = static_cast(track.getPt()); + const auto eta = static_cast(track.getEta()); + const auto phi = static_cast(track.getPhi()); + float chi2 = 0.F; + int ndf = 0; + if constexpr (std::is_same_v) { + chi2 = track.getChi2(); + ndf = 2 * track.getNumberOfClusters() - 5; + } else { + chi2 = static_cast(track.getTrackChi2()); + ndf = track.getNDF(); + } + result.push_back({.pt = pt, .eta = eta, .phi = phi, .chi2 = chi2, .ndf = ndf, .fitValid = std::isfinite(pt) && pt > 0.F && std::isfinite(eta) && std::isfinite(phi)}); + } + return result; +} + } // namespace TrackFile::TrackFile(const std::string& path, Chips chips) : mChips(chips), mFile(TFile::Open(path.c_str())) @@ -107,6 +138,16 @@ std::vector TrackFile::fakeFlags(long long entry) return result; } +std::vector TrackFile::kinematics(long long entry) +{ + if (entry >= mTree->GetEntries()) { + return {}; + } + return mChips == Chips::kITS + ? readKinematics(*mTree, entry, "ITSTrack") + : readKinematics(*mTree, entry, "MFTTrack"); +} + std::vector TrackFile::innerStates(long long entry) { if (mChips != Chips::kITS) { diff --git a/Detectors/GenericTracking/study/tools/CompleteAnatomy.h b/Detectors/GenericTracking/study/tools/CompleteAnatomy.h index 921306a73b6b2..13389c6b4a4ef 100644 --- a/Detectors/GenericTracking/study/tools/CompleteAnatomy.h +++ b/Detectors/GenericTracking/study/tools/CompleteAnatomy.h @@ -4342,7 +4342,11 @@ class CompleteAnatomy // The road fit, its last step nominal and allowed. const RowChain chain{.frame = &frame, .rows = road.data()}; const auto count = static_cast(road.size()); - const auto plain = fitOutward(propagator, geometry, material, chain, count, parameters.roadChi2MaxByNdf, nominalBudget); + auto plainParameters = parameters; + plainParameters.kinkScale2 = 1.F; + plainParameters.maxKinks = nominalBudget.max; + plainParameters.kinkChi2Min = nominalBudget.gate; + const auto plain = fitRoad(plainParameters, propagator, geometry, material, chain, count); const auto allowed = fitRoad(parameters, propagator, geometry, material, chain, count); if (plain.status != RoadFitStatus::kAccepted && allowed.status == RoadFitStatus::kAccepted) { need |= kLastStep; diff --git a/Detectors/GenericTracking/study/tools/StudyWorkflow.h b/Detectors/GenericTracking/study/tools/StudyWorkflow.h index 3d870b820f4c5..e5584c393dba6 100644 --- a/Detectors/GenericTracking/study/tools/StudyWorkflow.h +++ b/Detectors/GenericTracking/study/tools/StudyWorkflow.h @@ -37,6 +37,7 @@ #include "GenericTracking/Study/MCVertexReaderSpec.h" #include "GenericTracking/Study/MeasurementSinkSpec.h" #include "GenericTracking/Workflow/ClusterConverterSpec.h" +#include "GenericTracking/Workflow/ITSMFTOptions.h" #include "GenericTracking/Workflow/MeasurementProducerSpec.h" #include "GenericTracking/Workflow/TrackerSpec.h" @@ -54,6 +55,7 @@ namespace o2::gtrk::study /// Every study workflow's options. inline void addStudyOptions(std::vector& options) { + o2::itsmft::DPLAlpideParamInitializer::addConfigOption(options); using o2::framework::ConfigParamSpec; using o2::framework::VariantType; options.push_back(ConfigParamSpec{"disable-mc", VariantType::Bool, false, {"do not propagate MC labels"}}); @@ -86,7 +88,7 @@ template if (!context.options().get("clusters-from-upstream")) { specs.insert(specs.end(), readers.begin(), readers.end()); } - specs.push_back(gw::getClusterConverterSpec(useMC)); + specs.push_back(gw::getClusterConverterSpec(useMC, gw::staggeredInput(context))); specs.push_back(gw::getMeasurementProducerSpec()); specs.push_back(getMeasurementSinkSpec(useMC, mcVertices, products)); if (products != TrackerProducts::kNone) { diff --git a/Detectors/GenericTracking/study/tools/TrackBench.h b/Detectors/GenericTracking/study/tools/TrackBench.h index 46bc3abc5bf73..999976c646eca 100644 --- a/Detectors/GenericTracking/study/tools/TrackBench.h +++ b/Detectors/GenericTracking/study/tools/TrackBench.h @@ -51,6 +51,7 @@ #include "DetectorsBase/MatLayerCylSet.h" #include +#include #include "FitInputs.h" #include @@ -139,6 +140,7 @@ struct Options { // NOLINT(clang-analyzer-optin.performance.Padding) std::string trace; std::string json; std::string tracksOut; + std::string fitResults; std::string clonesOut; ///< --clones-out std::string rowsOut; ///< --rows-out std::string fakesOut; ///< --fakes-out @@ -546,6 +548,57 @@ void writeTracks(const std::string& path, const typename TrackRefit::Tracks } } +/// Per-track fit output for comparisons, outside the measured tracking loop. +/// Particle IDs include source/event/track; the caller supplies the TF identity. +template +void writeFitResults(const std::string& path, const typename TrackRefit::Tracks& tracks, const TimeFrame& frame, const study::TrackingTruth* truth) +{ + std::ofstream out{path}; + if (!out) { + throw std::runtime_error("cannot write " + path); + } + out << "track,iteration,status,count,layerMask,chi2,ndf,kinks,good,particle,ptGen,etaGen,physicalPrimary"; + for (int i = 0; i < kNBoundParams; ++i) { + out << ",par" << i; + } + for (int i = 0; i < kBoundCovSize; ++i) { + out << ",cov" << i; + } + out << ",clusters\n" << std::setprecision(std::numeric_limits::max_digits10); + const auto view = tracks.view(); + const auto columns = frame.measurements(); + for (std::size_t track = 0; track < tracks.size(); ++track) { + const auto rows = std::span{view.clusters[track].index, view.count[track]}; + const auto good = truth ? study::goodParticle(truth->rows, rows) : study::kNoTrackParticle; + const auto particle = good != study::kNoTrackParticle ? good : truth ? study::majorityParticle(truth->rows, rows).particle : study::kNoTrackParticle; + const study::ParticleTruth* mc = nullptr; + std::uint64_t layerMask = 0U; + if (truth) { + const auto it = truth->particles.find(particle); + if (it != truth->particles.end()) { + mc = &it->second; + } + for (const auto row : rows) { + layerMask |= std::uint64_t{1} << truth->rows[row].layer; + } + } + out << track << ',' << static_cast(view.iteration[track]) << ',' << static_cast(view.status[track]) << ',' << rows.size() << ',' << layerMask + << ',' << view.chi2[track] << ',' << view.ndf[track] << ',' << static_cast(view.kinks[track]) << ',' << (good != study::kNoTrackParticle) << ',' << particle + << ',' << (mc ? mc->pt : 0.F) << ',' << (mc ? mc->eta : 0.F) << ',' << (mc && mc->origin == study::ParticleOrigin::kPhysicalPrimary); + for (const auto parameter : view.state[track].par) { + out << ',' << parameter; + } + for (const auto covariance : view.state[track].cov) { + out << ',' << covariance; + } + out << ','; + for (const auto row : rows) { + out << columns.clusterIndex[row] << ';'; + } + out << '\n'; + } +} + /// The checkpoint file of --save-state and --resume-state: a header, then /// each part raw -- the first iterations' settings, which --resume-state checks against its own, their summaries, the used mask and /// the tracks. The header's fingerprints hold what else the iterations read: @@ -1920,6 +1973,9 @@ int track(const Options& options) // NOLINT(readability-function-size): the benc if (!options.tracksOut.empty() && last) { writeTracks(options.tracksOut, deviceTracks, frame); } + if (!options.fitResults.empty() && last) { + writeFitResults(options.fitResults, deviceTracks, frame, truth ? &*truth : nullptr); + } std::cout << std::format("\n=== run {}{}: tracks, {}, {}: {} in {:.4f} s (plan {:.4f} s, upload {:.4f} s, iterations", repetition, repetition == 0U ? " (cold)" : "", label, options.backend, deviceTracks.size(), planSeconds + deviceSeconds, planSeconds, deviceRun.uploadSeconds); for (const double seconds : deviceRun.iterationSeconds) { std::cout << std::format(" {:.4f}", seconds); @@ -2042,6 +2098,9 @@ int track(const Options& options) // NOLINT(readability-function-size): the benc if (!options.tracksOut.empty() && last) { writeTracks(options.tracksOut, *result.tracks, frame); } + if (!options.fitResults.empty() && last) { + writeFitResults(options.fitResults, *result.tracks, frame, truth ? &*truth : nullptr); + } if (truth && last) { writeTruthStudies(options, *result.tracks, *truth, comparison, frame, tracking, *material, propagator, studies.followed, [&](std::size_t iteration) -> const IterationParameters& { return tracker.iteration(iteration); }, [&](std::uint64_t id, std::size_t iteration) { return funnel ? funnel->stages(id, iteration) : std::nullopt; }, [&](std::size_t iteration, std::span rows) { return funnel ? funnel->recordChi2(iteration, rows) : std::nullopt; }); } diff --git a/Detectors/GenericTracking/study/tools/checkMatLUT.cxx b/Detectors/GenericTracking/study/tools/checkMatLUT.cxx index ae63215338ee3..f1f78a5193216 100644 --- a/Detectors/GenericTracking/study/tools/checkMatLUT.cxx +++ b/Detectors/GenericTracking/study/tools/checkMatLUT.cxx @@ -534,7 +534,14 @@ void write(const std::string& path, const std::vector& steps) template int runFor(const Options& options) { - o2::base::GeometryManager::loadGeometry(options.geometry.empty() ? options.sim + "/sgn_geometry.root" : options.geometry); + // O2 exports the ideal geometry before applying alignment, then transports + // with the aligned geometry (O2MCApplicationBase::MisalignGeometry). + auto geometryPath = options.geometry; + if (geometryPath.empty()) { + const auto aligned = options.sim + "/sgn_geometry-aligned.root"; + geometryPath = std::filesystem::exists(aligned) ? aligned : options.sim + "/sgn_geometry.root"; + } + o2::base::GeometryManager::loadGeometry(geometryPath); const auto geometry = std::make_unique>(o2::gtrk::alice::placeChecked()); const std::unique_ptr lut{o2::base::MatLayerCylSet::loadFromFile(options.lut)}; if (!lut) { @@ -577,9 +584,9 @@ int run(int argc, char** argv) description.add_options() // ("help,h", "this message") // ("detector,d", bpo::value(&options.detector)->default_value(options.detector), "ITS, MFT or TRKFT3") // - ("sim", bpo::value(&options.sim), "the simulation's directory: sgn_Hits*.root, sgn_Kine.root, sgn_geometry.root") // + ("sim", bpo::value(&options.sim), "the simulation's directory: sgn_Hits*.root, sgn_Kine.root, sgn_geometry[-aligned].root") // ("lut", bpo::value(&options.lut), "the material LUT") // - ("geometry", bpo::value(&options.geometry), "the geometry, if not /sgn_geometry.root (the hits' own transport geometry)") // + ("geometry", bpo::value(&options.geometry), "override geometry; default /sgn_geometry-aligned.root if present, otherwise sgn_geometry.root") // ("output,o", bpo::value(&options.output), "write the steps to this ROOT file") // ("events,n", bpo::value(&options.events)->default_value(options.events), "the first MC events") // ("min-pt", bpo::value(&options.minPt)->default_value(options.minPt), "the particles' pT floor [GeV]") // diff --git a/Detectors/GenericTracking/study/tools/gtrk-its-study-workflow.cxx b/Detectors/GenericTracking/study/tools/gtrk-its-study-workflow.cxx index b249b9d84feb0..3593854de694a 100644 --- a/Detectors/GenericTracking/study/tools/gtrk-its-study-workflow.cxx +++ b/Detectors/GenericTracking/study/tools/gtrk-its-study-workflow.cxx @@ -40,7 +40,7 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) o2::conf::ConfigurableParam::updateFromString(context.options().get("configKeyValues")); const bool useMC = !context.options().get("disable-mc"); - auto specs = o2::gtrk::study::studyChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, false, true, false)}); + auto specs = o2::gtrk::study::studyChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, o2::gtrk::workflow::staggeredInput(context), true, false)}); const o2::raw::HBFUtilsInitializer hbfIni(context, specs); diff --git a/Detectors/GenericTracking/study/tools/gtrk-mft-study-workflow.cxx b/Detectors/GenericTracking/study/tools/gtrk-mft-study-workflow.cxx index 2a5e060e598d8..e1aebcc237470 100644 --- a/Detectors/GenericTracking/study/tools/gtrk-mft-study-workflow.cxx +++ b/Detectors/GenericTracking/study/tools/gtrk-mft-study-workflow.cxx @@ -40,7 +40,7 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) o2::conf::ConfigurableParam::updateFromString(context.options().get("configKeyValues")); const bool useMC = !context.options().get("disable-mc"); - auto specs = o2::gtrk::study::studyChain(context, {o2::itsmft::getMFTClusterReaderSpec(useMC, false, true, false)}); + auto specs = o2::gtrk::study::studyChain(context, {o2::itsmft::getMFTClusterReaderSpec(useMC, o2::gtrk::workflow::staggeredInput(context), true, false)}); const o2::raw::HBFUtilsInitializer hbfIni(context, specs); diff --git a/Detectors/GenericTracking/study/tools/trackBench.cxx b/Detectors/GenericTracking/study/tools/trackBench.cxx index 399304e0c0e79..97eb2e364239a 100644 --- a/Detectors/GenericTracking/study/tools/trackBench.cxx +++ b/Detectors/GenericTracking/study/tools/trackBench.cxx @@ -166,6 +166,7 @@ static int run(int argc, char** argv) ("fit-inputs", bpo::value(&options.fitInputs)->default_value(""), "write the last run's triplet and neighbour fit inputs here, for o2-gtrk-fit-bench") // ("trace", bpo::value(&options.trace)->default_value(""), "write the last run's zones and chunks here, as a Chrome trace") // ("json", bpo::value(&options.json)->default_value(""), "append one JSON line per run here") // + ("fit-results", bpo::value(&options.fitResults)->default_value(""), "write the last run's fit parameters, covariance, cluster identities and MC labels as CSV, outside timing") // ("tracks-out", bpo::value(&options.tracksOut)->default_value(""), "write the last run's tracks here, one line each: iteration, clusters, chi2") // ("follow-steps", bpo::value(&options.followSteps)->default_value(""), "write the follow's steps on correct roads against the truth here, replayed after each iteration") // diff --git a/Detectors/GenericTracking/tools/src/importGeometry.cxx b/Detectors/GenericTracking/tools/src/importGeometry.cxx index e7323dcbda0dc..a357acedbd0c4 100644 --- a/Detectors/GenericTracking/tools/src/importGeometry.cxx +++ b/Detectors/GenericTracking/tools/src/importGeometry.cxx @@ -95,7 +95,7 @@ double checkTRKPixels(std::span sensors) const int subDet = disc ? ft3.getSubDetID(id) : geometry.getSubDetID(id); const int layer = disc ? ft3.getLayer(id) : geometry.getLayer(id); // within the sub-detector const auto& spec = Geometry::layerSpec(sensors[chip].layer); - const int nRows = subDet == 0 ? o2::trk::constants::VD::petal::layer::nRows[static_cast(layer)] : o2::trk::constants::moduleMLOT::chip::nRows; + const int nRows = subDet == 0 ? Segmentation::getNRowsVD(layer) : o2::trk::constants::moduleMLOT::chip::nRows; const int nCols = subDet == 0 ? o2::trk::constants::VD::petal::layer::nCols : o2::trk::constants::moduleMLOT::chip::nCols; for (const auto [row, col] : std::array, 5>{{{0, 0}, {nRows - 1, 0}, {0, nCols - 1}, {nRows - 1, nCols - 1}, {nRows / 2, nCols / 2}}}) { float xRow = 0.F; diff --git a/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh b/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh index 0af400d69b538..fb3fed5735999 100644 --- a/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh +++ b/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh @@ -641,7 +641,7 @@ class DeviceTrackerImpl final : public DeviceTracker mStream.synchronize(); } - [[nodiscard]] StageStatus findVertices(std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result) final + [[nodiscard]] StageStatus findVertices(std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result, std::vector* claims) final { requireNoIteration("findVertices"); if (pass >= mIterationCount) { @@ -649,6 +649,9 @@ class DeviceTrackerImpl final : public DeviceTracker } // Slice by slice, each its own frames' vertices; the lists in time order. result = {}; + if (claims != nullptr) { + claims->clear(); + } const auto& slices = mPlan->slices; StageStatus status = StageStatus::kComplete; for (std::size_t slice = 0; slice < slices.size() && status == StageStatus::kComplete; ++slice) { @@ -658,14 +661,20 @@ class DeviceTrackerImpl final : public DeviceTracker DeviceStageOutput triplets; status = tripletStages(pass, pass, slice, tracklets, permutation, triplets); VertexerResult found; + const auto claimBegin = claims == nullptr ? 0U : claims->size(); if (status == StageStatus::kComplete) { - status = vertex_glue::findVertices(*mRuntime, mArena, mStream, std::as_const(mColumns).view(), static_cast(mRows), constView(tracklets), constView(triplets), triplets.size(), mParameters.data() + pass, settings, beamSpot, ownedFrames(mHostFrame->grid().geometry(), slices[slice].slice.coreBlocks), found); + status = vertex_glue::findVertices(*mRuntime, mArena, mStream, std::as_const(mColumns).view(), static_cast(mRows), constView(tracklets), constView(triplets), triplets.size(), mParameters.data() + pass, settings, beamSpot, ownedFrames(mHostFrame->grid().geometry(), slices[slice].slice.coreBlocks), found, claims); } // A complete vertexer has synchronised after its last readback. if (status != StageStatus::kComplete) { mStream.synchronize(); } closeFrom(Scope::kSlice); + if (claims != nullptr) { + for (std::size_t c = claimBegin; c < claims->size(); ++c) { + (*claims)[c].vertex += static_cast(result.vertices.size()); + } + } appendVertices(result, found); } return status; @@ -1105,7 +1114,7 @@ class DeviceTrackerImpl final : public DeviceTracker if (!tracks.allocate(mArena, trackCount) || !proposals.allocate(mArena, followCount) || !kinkSteps.allocate(mArena, kinkStepCount) || !followRows.allocate(mArena, followCount) || !followCounts.allocate(mArena, followCount) || !endStarts.allocate(mArena, hostParameters.followEnd ? trackCount : 0U)) { return StageStatus::kOutOfMemory; } - const TrackRefitElement refit{mFrame.data(), mGeometry.data(), mMaterial.data(), mPropagator, records.clusters, records.count, records.outwardState, records.chi2, std::as_const(selected.selected).data(), tracks.view(), 0U, static_cast(mIteration), kinkBudget(hostParameters)}; + const TrackRefitElement refit{mFrame.data(), mGeometry.data(), mMaterial.data(), mPropagator, records.clusters, records.count, records.outwardState, records.chi2, std::as_const(selected.selected).data(), tracks.view(), 0U, static_cast(mIteration), kinkBudget(hostParameters), hostParameters.fitReference}; algorithms::forEach(*mRuntime, mStream, trackCount, refit); // The follow, the host's sequence diff --git a/Detectors/GenericTracking/tracking/GPU/cuda/VertexGlue.cuh b/Detectors/GenericTracking/tracking/GPU/cuda/VertexGlue.cuh index 42482ebf39af3..293ac87ebd241 100644 --- a/Detectors/GenericTracking/tracking/GPU/cuda/VertexGlue.cuh +++ b/Detectors/GenericTracking/tracking/GPU/cuda/VertexGlue.cuh @@ -680,14 +680,29 @@ struct RoundStorage { } // namespace } // namespace detail +/// Download only contributor identities; vertex finding stays on the device. +struct ClaimRows { + const Index* sortedTriplets; + const std::int32_t* claims; + TrackletConstView tracklets; + TripletConstView triplets; + TripletClaim* output; + + GPUdi() void operator()(Index i) const + { + output[i] = {.rows = tripletRows(tracklets, triplets, sortedTriplets[i]), .vertex = claims[i]}; + } +}; + /// The vertices of the pass's \a tripletCount triplets on \a measurements, /// \a rows of them, into \a result (`DeviceTracker::findVertices`); \a pass /// the pass's parameters on the device, for the candidates' errors; the /// vertices of \a owned's frames, a slice's, as `Vertexer::find` finds them. /// Its arrays on the base end, in the caller's scope; kOutOfMemory if the /// arena refuses. +/// With \a keptClaims, append the claimed triplets in the output vertex order. template -[[nodiscard]] StageStatus findVertices(DeviceRuntime& runtime, DeviceArena& arena, DeviceStream& stream, MeasurementView measurements, Index rows, TrackletConstView tracklets, TripletConstView triplets, Index tripletCount, const IterationParameters* pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, const OwnedFrames& owned, VertexerResult& result) +[[nodiscard]] StageStatus findVertices(DeviceRuntime& runtime, DeviceArena& arena, DeviceStream& stream, MeasurementView measurements, Index rows, TrackletConstView tracklets, TripletConstView triplets, Index tripletCount, const IterationParameters* pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, const OwnedFrames& owned, VertexerResult& result, std::vector* keptClaims) { using namespace detail; result = {}; @@ -1059,12 +1074,24 @@ template download(selectedScatterSum, scatterSum.data(), vertexSlots, stream); download(selectedCount, members.data(), vertexSlots, stream); } + std::vector downloadedClaims; + DeviceSpan claimRows; + if (keptClaims != nullptr) { + if (!claimRows.allocate(arena, count)) { + return StageStatus::kOutOfMemory; + } + algorithms::forEach(runtime, stream, count, ClaimRows{std::as_const(sortedTriplets).data(), std::as_const(claims).data(), tracklets, triplets, claimRows.data()}); + downloadedClaims.resize(count); + download(claimRows, downloadedClaims.data(), count, stream); + } stream.synchronize(); + std::vector vertexOfSlot(keptClaims == nullptr ? 0U : vertexSlots, VertexerResult::kNone); std::uint32_t mostInFrame = 0U; for (Index f = 0; f < frames; ++f) { struct Found { VertexSeed seed; Index contributors = 0U; + std::size_t slot = 0U; }; std::vector inFrame; for (std::size_t r = 0; r < found; ++r) { @@ -1080,15 +1107,24 @@ template .earliestPS = earliest[at], .latestPS = latest[at], }; - inFrame.push_back({.seed = vertexSeedOf(sums, beamSpot), .contributors = contributors[at]}); + inFrame.push_back({.seed = vertexSeedOf(sums, beamSpot), .contributors = contributors[at], .slot = at}); } mostInFrame = std::max(mostInFrame, static_cast(inFrame.size())); std::ranges::stable_sort(inFrame, {}, [](const Found& vertex) { return vertex.seed.z; }); for (const auto& vertex : inFrame) { + if (keptClaims != nullptr) { + vertexOfSlot[vertex.slot] = static_cast(result.vertices.size()); + } result.vertices.push_back(vertex.seed); result.contributors.push_back(vertex.contributors); } } + for (auto claim : downloadedClaims) { + if (claim.vertex >= 0) { + claim.vertex = vertexOfSlot[static_cast(claim.vertex)]; + keptClaims->push_back(claim); + } + } result.rounds = mostInFrame + 1U; return StageStatus::kComplete; } diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/CovarianceTransport.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/CovarianceTransport.h index 30a4ebcc86ee9..a1532ef255fe4 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/CovarianceTransport.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/CovarianceTransport.h @@ -182,17 +182,11 @@ GPUhdi() bool transportFactors(const SolenoidJacobian& jac, const StepNoise& noi /// The reference, not the filtered state. A Kalman filter that linearises /// about its own estimate is an *extended* filter, and its linearisation point /// carries whatever error the filter currently has; linearising about a -/// reference trajectory - O2 does this with the `linRef` overload of -/// `PropagatorImpl::PropagateToXBxByBz` - keeps F and Q independent of the -/// measurements. For a two-pass fit that is a correctness property rather than -/// a refinement: the outward and inward passes then see *identical* F and Q -/// for the same physical step, which is what lets them be combined. -/// -/// Unlike O2 we do not also transport the state through the reference. O2 -/// writes `state = ref + F (state - ref)` because its exact propagation is a -/// stepped field integration and doing it once for the reference is a real -/// saving. Ours is a closed form, so linearising the state would trade -/// accuracy for a speed-up we do not need. +/// reference trajectory allows the fitter to evaluate transport independently +/// of its measurement updates. The caller chooses the linearisation state; +/// this routine returns the exact propagated reference and its Jacobian. +/// The reference filter uses them for the affine prediction +/// f(reference) + F (state - reference). /// /// The derivation, for a uniform Bz. With `b = Bz * B2C` and `kp = (q/p) b`, /// the turn angle is `psi = kp * s` - **independent of theta**, because the diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Follow.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Follow.h index 7e562f5b7ce62..2144df0a489bf 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Follow.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Follow.h @@ -1429,7 +1429,7 @@ class FollowStage // Every lane holds the best so far, each passing row shared as // it comes. const auto one = follow::shared(group, candidate, from, branch.state); - if (++passing == 1U || follow::beats(one.row, one.chi2, {.row = child.row, .chi2 = child.chi2})) { + if (++passing == 1U || follow::beats(one.row, one.chi2, follow::Witness{.row = child.row, .chi2 = child.chi2})) { child = one; } continue; @@ -1513,7 +1513,7 @@ class FollowStage const bool confirmed = confirmation(branch, box, direction, taken, next, witness, WorkCounter::kFollowUpdates, group); const Scalar key = confirmed ? taken.chi2 + witness : taken.chi2; observe.lookaheadConfirmation(taken.row, taken.chi2, next, witness); - const bool better = best.child.row == ClusterList::kNoCluster || (confirmed != best.confirmed ? confirmed : follow::beats(taken.row, key, {.row = best.child.row, .chi2 = best.key})); + const bool better = best.child.row == ClusterList::kNoCluster || (confirmed != best.confirmed ? confirmed : follow::beats(taken.row, key, follow::Witness{.row = best.child.row, .chi2 = best.key})); if (better) { best.child = taken; best.confirmed = confirmed; @@ -1557,7 +1557,7 @@ class FollowStage continue; } const Scalar score = taken.chi2 + witness; - if (follow::beats(taken.row, score, {.row = child.row, .chi2 = bestScore})) { + if (follow::beats(taken.row, score, follow::Witness{.row = child.row, .chi2 = bestScore})) { child = taken; bestScore = score; onNext = next != ClusterList::kNoCluster; @@ -1909,7 +1909,7 @@ struct Refit { /// followed again with `followEarlier`. [[nodiscard]] GPUhdi() GPUAlwaysInline() bool accept(std::uint32_t track, const ClusterList& list, std::uint8_t n, std::uint8_t outward, std::uint8_t inward) const { - const auto fit = fitMeasurements(propagator, *geometry, *material, RowChain{.frame = frame, .rows = list.index}, static_cast(n), kinkBudget(*parameters)); + const auto fit = fitMeasurementsReference(propagator, *geometry, *material, RowChain{.frame = frame, .rows = list.index}, static_cast(n), parameters->fitReference, kinkBudget(*parameters)); if (!passesRoadChi2(*parameters, fit)) { return false; } @@ -2013,7 +2013,7 @@ struct TrimAmbiguous { rows.index[count - 1U] = ClusterList::kNoCluster; const auto shorter = static_cast(count - 1U); countWork(WorkCounter::kFollowTrimFits); - const auto fit = fitMeasurements(propagator, *geometry, *material, RowChain{.frame = frame, .rows = rows.index}, static_cast(shorter), kinkBudget(*parameters)); + const auto fit = fitMeasurementsReference(propagator, *geometry, *material, RowChain{.frame = frame, .rows = rows.index}, static_cast(shorter), parameters->fitReference, kinkBudget(*parameters)); if (!passesRoadChi2(*parameters, fit) || !keptAfterFollow(*parameters, shorter, fit.state)) { return; } diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/GPU/DeviceTracker.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/GPU/DeviceTracker.h index d2d13a614fa31..f2c29de310dd5 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/GPU/DeviceTracker.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/GPU/DeviceTracker.h @@ -235,10 +235,10 @@ class DeviceTracker /// (`OwnedFrames`), \a beamSpot giving the vertices' transverse spread. /// The vertices, sorted by /// time and each frame's by z, and their contributors into \a result, as - /// `VertexingPass::find` gives them; its counts too, but no per-triplet - /// claims and no beam estimate. Not with an + /// `VertexingPass::find` gives them; its counts too, and optionally their + /// claimed triplets in \a claims; no beam estimate. Not with an /// iteration open; the used mask and the tracks stay as they are. - [[nodiscard]] virtual StageStatus findVertices(std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result) = 0; + [[nodiscard]] virtual StageStatus findVertices(std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result, std::vector* claims = nullptr) = 0; /// Replace the open timeframe's plan with \a plan, of the same frame and /// the same run (std::invalid_argument before anything changes @@ -269,7 +269,7 @@ class DeviceTracker /// the pass, iteration \a pass, into \a result (findVertices). The /// timeframe stays open for trackTimeframe; \a frame and \a plan must /// outlive it. - [[nodiscard]] StageStatus findTimeframeVertices(const TimeFrame& frame, const TimeframePlan& plan, std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result, TimeframeReport& report) + [[nodiscard]] StageStatus findTimeframeVertices(const TimeFrame& frame, const TimeframePlan& plan, std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result, TimeframeReport& report, std::vector* claims = nullptr) { report = {}; auto start = Clock::now(); @@ -279,7 +279,7 @@ class DeviceTracker return status; } start = Clock::now(); - status = findVertices(pass, settings, beamSpot, result); + status = findVertices(pass, settings, beamSpot, result, claims); report.vertexSeconds = since(start); return status; } diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Iteration.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Iteration.h index e2d8d711b4f5f..f4e4a546426cb 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Iteration.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Iteration.h @@ -33,6 +33,7 @@ #include "GenericTracking/Base/LayerLinks.h" #include "GenericTracking/Base/LayerSpec.h" #include "GenericTracking/Tracking/Definitions.h" +#include "GenericTracking/Tracking/TrackFitter.h" #include "GenericTracking/Tracking/LinkHoles.h" #include "GenericTracking/Tracking/MathUtils.h" #include "GenericTracking/Tracking/Measurement.h" @@ -70,6 +71,7 @@ struct IterationSettings { int minClusters = 4; /// The kinks a fit pass may take; 0 for none. int maxKinks = 1; + FitReference fitReference = FitReference::kFilteredState; /// The triplet and neighbour cuts' kink allowance: the factor on theta0 at pT_min a triplet's middle /// point may have scattered by, beyond the nominal cuts, to be kept as a /// kinked triplet; 1 for none. @@ -294,6 +296,7 @@ struct IterationParameters { Scalar roadChi2MaxByNdf[(2 * maxChainLayers)-4]{}; /// The kinks a fit pass may take. std::uint8_t maxKinks = 0U; + FitReference fitReference = FitReference::kFilteredState; /// `IterationSettings::followEnd`, `innerTrim` /// and `crossingRule`. bool followEnd = false; @@ -568,6 +571,7 @@ template parameters.roadChi2MaxByNdf[ndf] = static_cast(chiSquareQuantile(ndf, settings.windows.nSigma)); } parameters.maxKinks = static_cast(std::max(settings.maxKinks, 0)); + parameters.fitReference = settings.fitReference; parameters.innerTrim = static_cast(std::clamp(settings.innerTrim, 0, static_cast(maxChainLayers))); parameters.crossingRule = settings.crossingRule; parameters.kinkChi2Min = static_cast(chiSquareQuantile(2, settings.windows.nSigma)); diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/RoadFit.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/RoadFit.h index fda09b03fb760..33504001cef66 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/RoadFit.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/RoadFit.h @@ -133,7 +133,7 @@ template template [[nodiscard]] GPUdi() bool clearsBeamLine(const Propagator& propagator, const IterationParameters& parameters, const TimeFrame& frame, const TrackingGeometry& geometry, const Material& material, const std::uint32_t* rows, std::uint32_t count, const BoundTrackState& outwardState, Scalar& chi2) { - const auto inward = fitInward(propagator, geometry, material, RowChain{.frame = &frame, .rows = rows}, static_cast(count), outwardState, kinkBudget(parameters)); + const auto inward = fitInwardReference(propagator, geometry, material, RowChain{.frame = &frame, .rows = rows}, static_cast(count), outwardState, parameters.fitReference, kinkBudget(parameters)); if (inward.status != PropagationStatus::kOk) { return false; } diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackFitter.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackFitter.h index d9e5446e389d5..a1c6c4d6ceffb 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackFitter.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackFitter.h @@ -10,12 +10,13 @@ // or submit itself to any jurisdiction. /// \file TrackFitter.h -/// \brief Independent outward and inward Kalman filters over a chain, each step linearised about its filtered state. +/// \brief Independent outward and inward Kalman filters over a chain. #ifndef GTRK_TRACKING_TRACKFITTER_H_ #define GTRK_TRACKING_TRACKFITTER_H_ #include "GenericTracking/Base/LayerSpec.h" +#include "GenericTracking/Base/LayerLinks.h" #include "GenericTracking/Tracking/CovarianceTransport.h" #include "GenericTracking/Tracking/Definitions.h" #include "GenericTracking/Tracking/Material.h" @@ -44,6 +45,18 @@ struct KinkBudget { Scalar gate = 0.F; ///< `chiSquareQuantile(2, nSigma)` }; +/// Linearisation of the road fits and final/follower refits. +/// The outward reference starts from the whole-chain circle seed. In +/// bidirectional mode the inward reference starts from the outward fitted +/// parameters. Each reference is transported independently and position-anchored +/// at every measurement, without Kalman updates to its direction or momentum. +/// Neither mode changes the reset covariance. +enum class FitReference : std::uint8_t { + kFilteredState = 0, + kCircleOutward = 1, + kIndependentBoth = 2 +}; + /// An iteration's kink budget, from its parameters' `maxKinks` and /// `kinkChi2Min`: the refit's and the extension's. template @@ -183,54 +196,84 @@ GPUdi() PropagationStatus transportStep(const Prop& prop, const StepMaterial& st return PropagationStatus::kOk; } -/// One step of a filter pass: \a state carried to \a measurement's plane -/// and updated with it, linearised about itself -- the filtered state the step -/// starts from, as ACTS's Kalman filter (a fixed whole-chain reference, never -/// corrected by the measurements, falls short of the outer planes of a low-momentum track near -/// its turning point). The covariance is \a factors, the pass's -/// (CovarianceFactors): `state.cov` is neither read nor written. \a Scaled, -/// the step's scattering theta0 is multiplied by \a noiseScale, both halves: -/// the kink allowance's last step (`allowanceStep`); the nominal step is its -/// own instantiation. \a Rule is the propagation's choice of crossing -/// (`CrossingRule`), and with it every pass below. Never inlined on the -/// device: one copy for each instantiation, not one in every fit of every -/// stage variant. -template -GPUdnii() PropagationStatus filterStep(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, [[maybe_unused]] Scalar noiseScale = 1.F) +/// One filter step, with the pass covariance held in factors. Normally the +/// filtered state supplies the linearisation. With Reference, F and material +/// are evaluated on an independent trajectory and the mean is predicted as +/// f(reference) + F (state - reference). The propagated reference, including +/// mean energy loss, is optionally saved for the next step. After a successful +/// update its position is anchored to the measurement, as in native ITS; its +/// direction and momentum keep their propagated values. Fits accept an +/// intersection outside the sensor's active area. +/// Scaled multiplies both halves' scattering theta0 by noiseScale. +template +GPUdnii() PropagationStatus filterStep(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, [[maybe_unused]] Scalar noiseScale = 1.F, [[maybe_unused]] const BoundTrackState* reference = nullptr, [[maybe_unused]] BoundTrackState* nextReference = nullptr) { - const auto& from = geometry.planeOf(state); + BoundTrackState prediction; + if constexpr (Reference) { + prediction = *reference; + } + auto& transportedState = Reference ? prediction : state; + const auto& from = geometry.planeOf(transportedState); const auto& to = geometry.planeOf(measurement); Scalar start[3]{}; - surfaceToGlobal(from, state.loc0(), state.loc1(), start); - const Scalar leaving = state.theta(); + surfaceToGlobal(from, transportedState.loc0(), transportedState.loc1(), start); + const Scalar leaving = transportedState.theta(); // An intersection just outside the sensor's active area is still the point // to update at: the cluster is there, so the track was. Bounds are for // searches, not fits (`jacobianAnalytic` takes kOutOfBounds). SolenoidJacobian jacobian; - const auto jacobianStatus = jacobianAnalytic(prop, state, from, to, direction, jacobian, &state); + const auto jacobianStatus = jacobianAnalytic(prop, transportedState, from, to, direction, jacobian, &transportedState); if (jacobianStatus != PropagationStatus::kOk) { return jacobianStatus; } Scalar end[3]{}; - surfaceToGlobal(to, state.loc0(), state.loc1(), end); + surfaceToGlobal(to, transportedState.loc0(), transportedState.loc1(), end); const auto step = stepMaterial(material, start, end); - const Scalar momentum = 1.F / o2::gpu::GPUCommonMath::Abs(state.qOverP()); + const Scalar momentum = 1.F / o2::gpu::GPUCommonMath::Abs(transportedState.qOverP()); const Scalar betaP = betaMomentum(momentum, prop.massGeV); // The nearer half before the transport (StepMaterial), at the direction on // the surface the step leaves; the farther half after it. const auto before = scatteringNoise(scaledTheta0(highlandTheta0(step.nearerX2X0, betaP), noiseScale), leaving); - const auto transported = transportStep(prop, step, betaP, before, direction, jacobian, state, factors, noiseScale); + const auto transported = transportStep(prop, step, betaP, before, direction, jacobian, transportedState, factors, noiseScale); if (transported != PropagationStatus::kOk) { return transported; } + if constexpr (Reference) { + if (nextReference != nullptr) { + *nextReference = prediction; + } + Scalar delta[kNBoundParams]; + for (int i = 0; i < kNBoundParams; ++i) { + delta[i] = state.par[i] - reference->par[i]; + } + o2::math_utils::detail::bringToPMPi(delta[kPhi]); + const Scalar* rows[3] = {jacobian.loc0, jacobian.loc1, jacobian.phi}; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < kNBoundParams; ++j) { + prediction.par[i] += rows[i][j] * delta[j]; + } + } + // These rows are the identity, as in the existing energy-loss transport. + prediction.par[kTheta] += delta[kTheta]; + prediction.par[kQOverP] += delta[kQOverP]; + state = prediction; + } const auto update = updateWithMeasurement(state, factors, measurement); if (!updateUsable(update, factors)) { return PropagationStatus::kNotFinite; } + if constexpr (Reference) { + if (nextReference != nullptr) { + // Re-anchor the reference to the observed position, without feeding the + // filtered direction, momentum or covariance into it (ITS ShiftRefToCluster). + nextReference->par[kLoc0] = measurement.loc0; + nextReference->par[kLoc1] = measurement.loc1; + } + } chi2 += update.chi2; return PropagationStatus::kOk; } @@ -263,16 +306,16 @@ GPUhdi() GPUAlwaysInline() void freeDirection(BoundTrackState& state, Covariance /// \a followed -- a measurement comes after this one in the pass -- and its /// increment is beyond the gate, taken again from the surface before with the /// direction freed, and counted in \a kinks. -template -GPUdi() PropagationStatus kinkStep(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, const KinkBudget& budget, bool followed, int& kinks) +template +GPUdi() PropagationStatus kinkStep(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, const KinkBudget& budget, bool followed, int& kinks, const BoundTrackState* reference = nullptr, BoundTrackState* nextReference = nullptr) { if (kinks >= budget.max || !followed) { - return filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + return filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); } const BoundTrackState keptState = state; const CovarianceFactors keptFactors = factors; const Scalar keptChi2 = chi2; - auto status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + auto status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); if (status != PropagationStatus::kOk || !(chi2 - keptChi2 > budget.gate)) { return status; } @@ -280,7 +323,7 @@ GPUdi() PropagationStatus kinkStep(const Prop& prop, const TrackingGeometry factors = keptFactors; chi2 = keptChi2; freeDirection(state, factors); - status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); kinks += status == PropagationStatus::kOk ? 1 : 0; return status; } @@ -292,24 +335,24 @@ GPUdi() PropagationStatus kinkStep(const Prop& prop, const TrackingGeometry /// \a noiseScale, as the extension's kink step, and counted in \a kinks; /// \a refused if the increment is still beyond the gate. It frees no /// parameter. -template -GPUdi() PropagationStatus allowanceStep(const Propagator& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, const KinkBudget& budget, Scalar noiseScale, int& kinks, bool& refused) +template +GPUdi() PropagationStatus allowanceStep(const Propagator& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, const KinkBudget& budget, Scalar noiseScale, int& kinks, bool& refused, const BoundTrackState* reference = nullptr, BoundTrackState* nextReference = nullptr) { refused = false; if (kinks >= budget.max) { - return filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + return filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); } const BoundTrackState keptState = state; const CovarianceFactors keptFactors = factors; const Scalar keptChi2 = chi2; - auto status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + auto status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); if (status != PropagationStatus::kOk || !(chi2 - keptChi2 > budget.gate)) { return status; } state = keptState; factors = keptFactors; chi2 = keptChi2; - status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, noiseScale); + status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, noiseScale, reference, nextReference); if (status == PropagationStatus::kOk) { ++kinks; refused = !(chi2 - keptChi2 <= budget.gate); @@ -344,7 +387,7 @@ GPUhdi() GPUAlwaysInline() MeasurementUpdateResult startOnMeasurement(BoundTrack /// A filter pass from \a seed over the chain. Its covariance is carried as /// factors from the first update to the last (CovarianceFactors) and written /// to \a result's `cov` once, at the end. -template +template GPUdi() PropagationStatus runPass(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, PropagationDirection direction, BoundTrackState seed, BoundTrackState& result, Scalar& chi2, std::int16_t& failedMeasurement, const KinkBudget& passBudget, int& kinks, std::int16_t& kinkMeasurement, const Observer& observe) { const auto budget = kinkBudgetFor(passBudget, n); @@ -359,17 +402,27 @@ GPUdi() PropagationStatus runPass(const Prop& prop, const TrackingGeometry& return PropagationStatus::kNotFinite; } chi2 = firstUpdate.chi2; + BoundTrackState reference; + if constexpr (Reference) { + // Begin each pass from its own seed: circle outward, outward fit inward. + // startOnMeasurement has already anchored the seed's position. + reference = state; + } observe(direction, first, chi2); for (int step = 1; step < n; ++step) { const int index = first + step * increment; const Scalar before = chi2; const int kinksBefore = kinks; const auto& stepProp = stepPropagator(prop, geometry, measurements, index - increment, index); - const auto status = kinkStep(stepProp, geometry, material, measurements[index], direction, state, factors, chi2, budget, step < n - 1, kinks); + BoundTrackState nextReference; + const auto status = kinkStep(stepProp, geometry, material, measurements[index], direction, state, factors, chi2, budget, step < n - 1, kinks, &reference, &nextReference); if (status != PropagationStatus::kOk) { failedMeasurement = static_cast(index); return status; } + if constexpr (Reference) { + reference = nextReference; + } if (kinks != kinksBefore) { kinkMeasurement = static_cast(index); } @@ -382,7 +435,7 @@ GPUdi() PropagationStatus runPass(const Prop& prop, const TrackingGeometry& /// The inward pass of `fitMeasurements`: its seed, the outward pass's end /// state with a diffuse covariance, then the filter from there in. -template +template GPUdi() void inwardPass(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, Scalar qOverPSigma, const KinkBudget& kinkBudget, const Observer& observe, TrackFitResult& fit) { // From the outward pass's end state, its covariance made diffuse with the @@ -393,7 +446,7 @@ GPUdi() void inwardPass(const Prop& prop, const TrackingGeometry& geometry, setDiffuseCovariance(outerSeed, qOverPSigma); int kinks = 0; - const auto status = runPass(prop, geometry, material, measurements, n, PropagationDirection::kBackward, outerSeed, fit.state, fit.chi2, fit.failedMeasurement, kinkBudget, kinks, fit.kinkMeasurement, observe); + const auto status = runPass(prop, geometry, material, measurements, n, PropagationDirection::kBackward, outerSeed, fit.state, fit.chi2, fit.failedMeasurement, kinkBudget, kinks, fit.kinkMeasurement, observe); if (status != PropagationStatus::kOk) { fit.status = status; fit.failedPass = FitPass::kInward; @@ -441,6 +494,39 @@ GPUdi() TrackFitResult fitMeasurements(const Prop& prop, const TrackingGeometry< return fit; } +/// Refit using the whole-chain circle/longitudinal seed as an independent +/// linearisation trajectory outward. In bidirectional mode, an independent +/// reference starts inward from the outward fitted parameters. Both references +/// receive mean energy loss and position anchoring at each measurement; their +/// direction and momentum are not Kalman-updated within a pass. Outward-only +/// mode keeps the existing inward filter. No posterior covariance is recycled +/// as a prior. +/// Failed reference intersections are reported as failed fits, without fallback. +template +GPUdnii() TrackFitResult fitMeasurementsReference(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, FitReference referenceMode, const KinkBudget& kinkBudget = {}, const Observer& observe = {}) +{ + if (referenceMode == FitReference::kFilteredState) { + return fitMeasurements(prop, geometry, material, measurements, n, kinkBudget, observe); + } + TrackFitResult fit{.ndf = static_cast((2 * n) - kNBoundParams)}; + const auto seed = makeChainSeed(prop, geometry, measurements, n); + int outwardKinks = 0; + std::int16_t outwardKinkMeasurement = -1; + const auto status = detail::runPass(prop, geometry, material, measurements, n, PropagationDirection::kForward, seed.state, fit.outwardState, fit.outwardChi2, fit.failedMeasurement, kinkBudget, outwardKinks, outwardKinkMeasurement, observe); + fit.outwardKinks = static_cast(outwardKinks); + if (status != PropagationStatus::kOk) { + fit.status = status; + fit.failedPass = FitPass::kOutward; + return fit; + } + if (referenceMode == FitReference::kIndependentBoth) { + detail::inwardPass(prop, geometry, material, measurements, n, seed.qOverPSigma, kinkBudget, observe, fit); + } else { + detail::inwardPass(prop, geometry, material, measurements, n, seed.qOverPSigma, kinkBudget, observe, fit); + } + return fit; +} + /// `fitMeasurements` given its outward pass's end state, \a outwardState /// -- as the road fit keeps it for an accepted /// road (`RoadFitResult::state`): only the inward pass runs. The outward @@ -455,6 +541,22 @@ GPUdi() TrackFitResult fitInward(const Prop& prop, const TrackingGeometry& return fit; } +/// Finish a road fit without repeating its outward KF. Bidirectional mode +/// starts an independent inward reference from the outward fitted parameters; +/// outward-only and ordinary modes use the existing inward filter. +template +GPUdi() TrackFitResult fitInwardReference(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, const BoundTrackState& outwardState, FitReference mode, const KinkBudget& kinkBudget = {}, const Observer& observe = {}) +{ + if (mode != FitReference::kIndependentBoth) { + return fitInward(prop, geometry, material, measurements, n, outwardState, kinkBudget, observe); + } + TrackFitResult fit{.ndf = static_cast((2 * n) - kNBoundParams)}; + fit.outwardState = outwardState; + const auto seed = makeChainSeed(prop, geometry, measurements, n); + detail::inwardPass(prop, geometry, material, measurements, n, seed.qOverPSigma, kinkBudget, observe, fit); + return fit; +} + /// The road fit's candidate pass: the outward filter /// alone, as in `fitMeasurements`. struct RoadFitResult { @@ -480,7 +582,7 @@ struct RoadFitResult { /// beyond the gate, it frees no degree of freedom, and the result's kinks /// count it. A compile-time choice: without it the last step is a `kinkStep` /// like the others. -template +template GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, const Scalar* chi2MaxByNdf, const KinkBudget& kinkBudget = {}, [[maybe_unused]] Scalar lastStepScale = 1.F) { const auto budget = detail::kinkBudgetFor(kinkBudget, n); @@ -498,6 +600,10 @@ GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& return result; } result.chi2 = firstUpdate.chi2; + BoundTrackState reference; + if constexpr (Reference) { + reference = state; + } const int ndf = (2 * n) - kNBoundParams; int kinks = 0; // The allowance's last step, which frees no degree of freedom; counted by its @@ -515,11 +621,12 @@ GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& } countWork(WorkCounter::kRoadSteps); const auto& stepProp = stepPropagator(prop, geometry, measurements, index - 1, index); + BoundTrackState nextReference; if constexpr (LastStepAllowance) { if (index == n - 1) { int taken = kinks; bool refused = false; - const auto status = detail::allowanceStep(stepProp, geometry, material, measurements[index], PropagationDirection::kForward, state, factors, result.chi2, budget, lastStepScale, taken, refused); + const auto status = detail::allowanceStep(stepProp, geometry, material, measurements[index], PropagationDirection::kForward, state, factors, result.chi2, budget, lastStepScale, taken, refused, &reference, &nextReference); if (status != PropagationStatus::kOk) { countWork(WorkCounter::kRoadPropagation); result.status = RoadFitStatus::kPropagationFailed; @@ -531,15 +638,21 @@ GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& return result; } allowed = taken - kinks; + if constexpr (Reference) { + reference = nextReference; + } continue; } } - const auto status = detail::kinkStep(stepProp, geometry, material, measurements[index], PropagationDirection::kForward, state, factors, result.chi2, budget, index < n - 1, kinks); + const auto status = detail::kinkStep(stepProp, geometry, material, measurements[index], PropagationDirection::kForward, state, factors, result.chi2, budget, index < n - 1, kinks, &reference, &nextReference); if (status != PropagationStatus::kOk) { countWork(WorkCounter::kRoadPropagation); result.status = RoadFitStatus::kPropagationFailed; return result; } + if constexpr (Reference) { + reference = nextReference; + } } result.kinks = static_cast(kinks + allowed); covarianceFromFactors(factors, state.cov); @@ -559,6 +672,9 @@ GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& template [[nodiscard]] GPUdi() RoadFitResult fitRoad(const Parameters& parameters, const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n) { + if (parameters.fitReference != FitReference::kFilteredState) { + return fitOutward(prop, geometry, material, measurements, n, parameters.roadChi2MaxByNdf, kinkBudget(parameters), o2::gpu::CAMath::Sqrt(parameters.kinkScale2)); + } if constexpr (LastStepAllowance) { return fitOutward(prop, geometry, material, measurements, n, parameters.roadChi2MaxByNdf, kinkBudget(parameters), o2::gpu::CAMath::Sqrt(parameters.kinkScale2)); } else { diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackRefit.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackRefit.h index 32853986a3a14..98a7778d17b21 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackRefit.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackRefit.h @@ -111,7 +111,8 @@ GPUhdi() GPUAlwaysInline() void writeFit(const Tracks& tracks, std::uint32_t tra /// One element per selected record: the inward pass on its rows from the road /// fit's outward state (`fitInward`), and the track's time from their boxes; -/// the pass may take the iteration's kinks. Its frame, geometry and material +/// the pass may take the iteration's kinks. Bidirectional reference mode starts +/// its inward reference from that outward state. Its frame, geometry and material /// are pointers, which may be a device's. template class TrackRefitElement @@ -120,7 +121,7 @@ class TrackRefitElement static constexpr std::size_t kMaxClusters = maxChainLayers; using View = typename BoundedColumns::template RecordT>::View; - TrackRefitElement(const TimeFrame* frame, const TrackingGeometry* geometry, const Material* material, const Propagator& propagator, const ClusterList* clusters, const std::uint8_t* count, const BoundTrackState* outward, const Scalar* outwardChi2, const std::uint32_t* selected, View tracks, std::uint32_t offset, std::uint8_t iteration, KinkBudget kinks) + TrackRefitElement(const TimeFrame* frame, const TrackingGeometry* geometry, const Material* material, const Propagator& propagator, const ClusterList* clusters, const std::uint8_t* count, const BoundTrackState* outward, const Scalar* outwardChi2, const std::uint32_t* selected, View tracks, std::uint32_t offset, std::uint8_t iteration, KinkBudget kinks, FitReference reference = FitReference::kFilteredState) : mFrame(frame), mGeometry(geometry), mMaterial(material), @@ -133,7 +134,8 @@ class TrackRefitElement mTracks(tracks), mOffset(offset), mIteration(iteration), - mKinks(kinks) + mKinks(kinks), + mReference(reference) { } @@ -144,7 +146,8 @@ class TrackRefitElement const auto& rows = mClusters[record]; const std::uint8_t count = mCount[record]; // The road fit ran the outward pass already and kept its end state. - const auto fit = fitInward(mPropagator, *mGeometry, *mMaterial, RowChain{.frame = mFrame, .rows = rows.index}, static_cast(count), mOutward[record], mKinks); + const auto chain = RowChain{.frame = mFrame, .rows = rows.index}; + const auto fit = fitInwardReference(mPropagator, *mGeometry, *mMaterial, chain, static_cast(count), mOutward[record], mReference, mKinks); writeFit(mTracks, track, fit, rows, count, trackTimeBox(mFrame->measurements(), rows, count)); mTracks.outwardChi2[track] = mOutwardChi2[record]; mTracks.iteration[track] = mIteration; @@ -166,6 +169,7 @@ class TrackRefitElement std::uint32_t mOffset = 0U; std::uint8_t mIteration = 0U; KinkBudget mKinks; + FitReference mReference; }; /// The timeframe's tracks: each iteration's appended after the ones before, @@ -200,7 +204,7 @@ class TrackRefit const auto offset = static_cast(mTracks.size()); mTracks.resize(offset + selected.size()); const auto records = selection.records(); - const TrackRefitElement element{&frame, &geometry, &material, propagator, records.clusters, records.count, records.outwardState, records.chi2, selected.data(), mTracks.view(), offset, iteration, kinkBudget(parameters)}; + const TrackRefitElement element{&frame, &geometry, &material, propagator, records.clusters, records.count, records.outwardState, records.chi2, selected.data(), mTracks.view(), offset, iteration, kinkBudget(parameters), parameters.fitReference}; forEach(static_cast(selected.size()), element); return StageStatus::kComplete; }); diff --git a/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx b/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx index f67c63d8544b1..a2944f3788c9f 100644 --- a/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx +++ b/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx @@ -172,6 +172,7 @@ struct Options { /// never both in one iteration. double kinkSigma = 1.; bool followEnd = false; + FitReference reference = FitReference::kFilteredState; }; std::string describe(const Options& options) { return std::string{options.source == VertexSource::kBeamLine ? "beam line" : "vertex list"} + ", width " + std::to_string(options.width) + ", shared " + std::to_string(options.sharedInnermost) + (options.extension ? ", follow" : ", no follow") + (options.kinkSigma > 1. ? ", kink allowance " + std::to_string(options.kinkSigma) + " sigma" : "") + (options.followEnd ? ", end again" : ""); } @@ -182,6 +183,7 @@ std::unique_ptr makeHost(const Options& options) static const UniformMaterial none{}; auto settings = iterations(options.source); for (auto& iteration : settings) { + iteration.fitReference = options.reference; iteration.sharedInnermost = options.sharedInnermost; iteration.follow = options.extension; } @@ -232,7 +234,44 @@ bool sameTracks(const View& host, std::size_t hostRows, const View& device, std: std::ranges::sort(result); return result; }; - return hostRows == deviceRows && rowsOf(host, hostRows) == rowsOf(device, deviceRows); + if (hostRows != deviceRows || rowsOf(host, hostRows) != rowsOf(device, deviceRows)) { + return false; + } + // Check states explicitly: a non-arithmetic memlayout column has no scalar + // maxAbsDifference in StageComparison, so matching identities is not a + // numerical check of its parameters and covariance. + const auto order = [](const View& view, std::size_t rows) { + std::vector, std::size_t>> sorted; + for (std::size_t i = 0; i < rows; ++i) { + sorted.emplace_back(std::vector{view.clusters[i].index, view.clusters[i].index + view.count[i]}, i); + } + std::ranges::sort(sorted); + return sorted; + }; + const auto lhs = order(host, hostRows); + const auto rhs = order(device, deviceRows); + const auto closeState = [](const BoundTrackState& a, const BoundTrackState& b) { + constexpr Scalar tolerance = 1.e-4F; + if (a.surface != b.surface || std::abs(a.surfacePhi - b.surfacePhi) > tolerance) { + return false; + } + for (int i = 0; i < kNBoundParams; ++i) { + Scalar delta = a.par[i] - b.par[i]; + if (i == kPhi) { + o2::math_utils::detail::bringToPMPi(delta); + } + if (!(std::abs(delta) <= tolerance)) { + return false; + } + } + return std::ranges::equal(a.cov, b.cov, [=](Scalar x, Scalar y) { return std::abs(x - y) <= tolerance; }); + }; + for (std::size_t i = 0; i < lhs.size(); ++i) { + if (!closeState(host.state[lhs[i].second], device.state[rhs[i].second]) || !closeState(host.outwardState[lhs[i].second], device.outwardState[rhs[i].second])) { + return false; + } + } + return true; #endif } @@ -397,22 +436,25 @@ BOOST_AUTO_TEST_CASE(EverySliceAndIterationMatchesTheHostTracker) auto* device = runtime(); const auto tracker = kernels().makeToyBarrel8Tracker(*device, block(kArenaBytes)); BOOST_TEST(tracker->arenaBytes() == kArenaBytes); - for (const auto& options : configurations()) { - BOOST_TEST_CONTEXT(describe(options)) - { - const auto totals = compare(*tracker, options); - BOOST_TEST(totals.triplets > 0U); - BOOST_TEST(totals.neighbours > 0U); - BOOST_TEST(totals.chained > 0U); - BOOST_TEST(totals.roads > 0U); - BOOST_TEST(totals.records > 0U); - BOOST_TEST(totals.gathered == totals.records); - BOOST_TEST(totals.tracks == kTracks); - BOOST_TEST(totals.maskedSlices > 0U); - if (options.kinkSigma > 1.) { - BOOST_TEST(totals.extended > 0U); + for (const auto reference : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + for (auto options : configurations()) { + options.reference = reference; + BOOST_TEST_CONTEXT(describe(options) << ", reference " << static_cast(reference)) + { + const auto totals = compare(*tracker, options); + BOOST_TEST(totals.triplets > 0U); + BOOST_TEST(totals.neighbours > 0U); + BOOST_TEST(totals.chained > 0U); + BOOST_TEST(totals.roads > 0U); + BOOST_TEST(totals.records > 0U); + BOOST_TEST(totals.gathered == totals.records); + BOOST_TEST(totals.tracks == kTracks); + BOOST_TEST(totals.maskedSlices > 0U); + if (options.kinkSigma > 1.) { + BOOST_TEST(totals.extended > 0U); + } + BOOST_TEST_MESSAGE(device->backend() << " " << describe(options) << ": " << totals.slices << " slices, " << totals.tracklets << " tracklets, " << totals.triplets << " triplets, " << totals.neighbours << " neighbours, " << totals.roads << " roads, " << totals.records << " records, " << totals.rounds << " rounds, " << totals.tracks << " tracks, " << totals.extended << " extended"); } - BOOST_TEST_MESSAGE(device->backend() << " " << describe(options) << ": " << totals.slices << " slices, " << totals.tracklets << " tracklets, " << totals.triplets << " triplets, " << totals.neighbours << " neighbours, " << totals.roads << " roads, " << totals.records << " records, " << totals.rounds << " rounds, " << totals.tracks << " tracks, " << totals.extended << " extended"); } } } @@ -764,7 +806,7 @@ BOOST_AUTO_TEST_CASE(TheVertexingPassGivesTheHostsListAndTracks) // The host: the pass's list, and the tracking from it. VertexingPass pass{passSettings(), VertexerSettings{}, Fixture::kBeamConstraint, scattering, pool, 2U, 2}; - const auto expected = pass.find(*fixture.frame, 13, 0); + const auto expected = pass.find(*fixture.frame, 13, 0, true); BOOST_REQUIRE(expected.status == StageStatus::kComplete); BOOST_REQUIRE(expected.found.vertices.size() == kCollisionsNS.size()); const auto host = makeHost({.source = VertexSource::kVertexList}); @@ -786,7 +828,8 @@ BOOST_AUTO_TEST_CASE(TheVertexingPassGivesTheHostsListAndTracks) const auto empty = planTimeframe(*fixture.frame, parameters, windows, host->beamLine(), {}, 13, 0); BOOST_REQUIRE((tracker->beginTimeframe(*fixture.frame, empty) == StageStatus::kComplete)); VertexerResult found; - BOOST_REQUIRE((tracker->findVertices(host->iterationCount(), VertexerSettings{}, Fixture::kBeamConstraint, found) == StageStatus::kComplete)); + std::vector claims; + BOOST_REQUIRE((tracker->findVertices(host->iterationCount(), VertexerSettings{}, Fixture::kBeamConstraint, found, &claims) == StageStatus::kComplete)); BOOST_REQUIRE(found.vertices.size() == expected.found.vertices.size()); BOOST_TEST(found.candidates == expected.found.candidates); for (std::size_t v = 0; v < found.vertices.size(); ++v) { @@ -799,6 +842,18 @@ BOOST_AUTO_TEST_CASE(TheVertexingPassGivesTheHostsListAndTracks) } } + // Contributor identities must survive both vertex sorting and slice offsets. + auto expectedClaims = expected.claims; + std::erase_if(expectedClaims, [](const TripletClaim& claim) { return claim.vertex < 0; }); + const auto order = [](const TripletClaim& claim) { return std::pair{claim.vertex, claim.rows}; }; + std::ranges::sort(expectedClaims, {}, order); + std::ranges::sort(claims, {}, order); + BOOST_REQUIRE(claims.size() == expectedClaims.size()); + for (std::size_t c = 0; c < claims.size(); ++c) { + BOOST_TEST(claims[c].vertex == expectedClaims[c].vertex); + BOOST_TEST(claims[c].rows == expectedClaims[c].rows, boost::test_tools::per_element()); + } + // The tracking's plan with the device's list replaces the empty one. const auto trackingPlan = planTimeframe(*fixture.frame, parameters, windows, host->beamLine(), found.vertices, 13, 0); BOOST_REQUIRE((tracker->replacePlan(trackingPlan) == StageStatus::kComplete)); diff --git a/Detectors/GenericTracking/tracking/test/testRoadFit.cxx b/Detectors/GenericTracking/tracking/test/testRoadFit.cxx index 6631e5776d53b..afe14afb233cd 100644 --- a/Detectors/GenericTracking/tracking/test/testRoadFit.cxx +++ b/Detectors/GenericTracking/tracking/test/testRoadFit.cxx @@ -220,23 +220,28 @@ BOOST_AUTO_TEST_CASE(SmearedRoadsHaveTheirDegreesOfFreedom) hit.point.z = global[2]; } const Fixture fixture{hits, 100.}; - const auto cut = study::roadChi2Cut(parameters(), kLayers); - double sum = 0.; - std::size_t accepted = 0U; - for (std::size_t track = 0; track < kTracks; ++track) { - const auto rows = trackRows(fixture, track); - const RowChain chain{.frame = fixture.frame.get(), .rows = rows.data()}; - const auto fit = fitOutward(kPropagator, trackingGeometry(), kNoMaterial, chain, static_cast(kLayers), kNoCut.data()); - BOOST_REQUIRE((fit.status == RoadFitStatus::kAccepted)); - sum += static_cast(fit.chi2); - accepted += fit.chi2 <= cut ? 1U : 0U; + for (const auto mode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + const auto cut = study::roadChi2Cut(parameters(), kLayers); + auto iteration = parameters(); + iteration.fitReference = mode; + std::fill(std::begin(iteration.roadChi2MaxByNdf), std::end(iteration.roadChi2MaxByNdf), std::numeric_limits::infinity()); + double sum = 0.; + std::size_t accepted = 0U; + for (std::size_t track = 0; track < kTracks; ++track) { + const auto rows = trackRows(fixture, track); + const RowChain chain{.frame = fixture.frame.get(), .rows = rows.data()}; + const auto fit = fitRoad(iteration, kPropagator, trackingGeometry(), kNoMaterial, chain, static_cast(kLayers)); + BOOST_REQUIRE((fit.status == RoadFitStatus::kAccepted)); + sum += static_cast(fit.chi2); + accepted += fit.chi2 <= cut ? 1U : 0U; + } + const double mean = sum / static_cast(kTracks); + BOOST_TEST_MESSAGE("mean chi2 " << mean << " (11), within the cut " << accepted << "/" << kTracks); + // The mean's error over 400 fits is sqrt(2 * 11 / 400) = 0.23. + BOOST_TEST(mean > 10.F); + BOOST_TEST(mean < 12.F); + BOOST_TEST(accepted >= 0.99 * kTracks); } - const double mean = sum / static_cast(kTracks); - BOOST_TEST_MESSAGE("mean chi2 " << mean << " (11), within the cut " << accepted << "/" << kTracks); - // The mean's error over 400 fits is sqrt(2 * 11 / 400) = 0.23. - BOOST_TEST(mean > 10.F); - BOOST_TEST(mean < 12.F); - BOOST_TEST(accepted >= 0.99 * kTracks); } /// A helix whose direction turns by \a turn on the sensor of layer \a kink -- @@ -572,78 +577,81 @@ BOOST_AUTO_TEST_CASE(SelectedTracksAreRefittedWithBothPasses) // fit, the selection, the refit. constexpr std::size_t kTracks = 6U; const Fixture fixture{trackHits(kTracks)}; - const auto iteration = parameters(); - const RoadsRun roads{fixture}; - const auto& run = roads.result; - const auto& count = roads.count; - const SerialForEach serial{}; - const Stage stage{iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, *run.tracklets, *run.triplets, count.walk(*run.triplets, *run.neighbours)}; - RecordOutput records; - records.setMemoryPool(std::make_shared()); - records.prepare(count.roadCount()); - BOOST_REQUIRE(runSerial(records, stage) == StageStatus::kComplete); - GlobalSelection selection; - selection.setMemoryPool(std::make_shared()); - BOOST_REQUIRE(selection.append(records, SerialForEach{}) == StageStatus::kComplete); - std::vector mask(fixture.frame->size(), kFree); - BOOST_REQUIRE(selection.select(serial, [](std::span indices, const auto& less) { std::ranges::sort(indices, less); }, 0U, mask) == StageStatus::kComplete); - BOOST_REQUIRE(selection.selected().size() == kTracks); - - TrackRefit refit; - refit.setMemoryPool(std::make_shared()); - BOOST_REQUIRE(refit.run(serial, selection, iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, 0U) == StageStatus::kComplete); - BOOST_REQUIRE(refit.size() == kTracks); - const auto tracks = refit.tracks().view(); - for (std::uint32_t track = 0; track < kTracks; ++track) { - const auto& rows = tracks.clusters[track]; - const auto index = fixture.frame->measurement(rows.index[0]).clusterIndex / kLayers; - const auto h = helix(index); - BOOST_TEST_CONTEXT("track " << index) - { - // Both passes, as fitMeasurements gives them on the same rows. - BOOST_REQUIRE(tracks.status[track] == PropagationStatus::kOk); - BOOST_TEST(tracks.count[track] == kLayers); - std::array copied{}; - const RowChain chain{.frame = fixture.frame.get(), .rows = rows.index}; - for (std::size_t i = 0; i < kLayers; ++i) { - copied[i] = chain[static_cast(i)]; + for (const auto mode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + auto iteration = parameters(); + iteration.fitReference = mode; + const RoadsRun roads{fixture}; + const auto& run = roads.result; + const auto& count = roads.count; + const SerialForEach serial{}; + const Stage stage{iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, *run.tracklets, *run.triplets, count.walk(*run.triplets, *run.neighbours)}; + RecordOutput records; + records.setMemoryPool(std::make_shared()); + records.prepare(count.roadCount()); + BOOST_REQUIRE(runSerial(records, stage) == StageStatus::kComplete); + GlobalSelection selection; + selection.setMemoryPool(std::make_shared()); + BOOST_REQUIRE(selection.append(records, SerialForEach{}) == StageStatus::kComplete); + std::vector mask(fixture.frame->size(), kFree); + BOOST_REQUIRE(selection.select(serial, [](std::span indices, const auto& less) { std::ranges::sort(indices, less); }, 0U, mask) == StageStatus::kComplete); + BOOST_REQUIRE(selection.selected().size() == kTracks); + + TrackRefit refit; + refit.setMemoryPool(std::make_shared()); + BOOST_REQUIRE(refit.run(serial, selection, iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, 0U) == StageStatus::kComplete); + BOOST_REQUIRE(refit.size() == kTracks); + const auto tracks = refit.tracks().view(); + for (std::uint32_t track = 0; track < kTracks; ++track) { + const auto& rows = tracks.clusters[track]; + const auto index = fixture.frame->measurement(rows.index[0]).clusterIndex / kLayers; + const auto h = helix(index); + BOOST_TEST_CONTEXT("track " << index) + { + // Both passes, as fitMeasurements gives them on the same rows. + BOOST_REQUIRE(tracks.status[track] == PropagationStatus::kOk); + BOOST_TEST(tracks.count[track] == kLayers); + std::array copied{}; + const RowChain chain{.frame = fixture.frame.get(), .rows = rows.index}; + for (std::size_t i = 0; i < kLayers; ++i) { + copied[i] = chain[static_cast(i)]; + } + const auto full = fitMeasurementsReference(kPropagator, trackingGeometry(), kNoMaterial, copied.data(), static_cast(kLayers), mode); + BOOST_TEST(tracks.chi2[track] == full.chi2); + + BOOST_TEST(tracks.ndf[track] == (2 * static_cast(kLayers)) - 5); + // On the innermost surface, and the generating helix. + const auto& state = tracks.state[track]; + const auto first = fixture.frame->measurement(rows.index[0]); + BOOST_TEST(state.surface == first.surface); + BOOST_TEST(std::abs(state.loc0() - first.loc0) < 1.e-4F); + BOOST_TEST(std::abs(state.loc1() - first.loc1) < 1.e-4F); + BOOST_TEST(state.theta() == static_cast(std::atan2(1., h.slope)), boost::test_tools::tolerance(1.e-4F)); + const double momentum = h.ptGeV * std::hypot(1., h.slope); + BOOST_TEST(std::abs(state.qOverP()) == static_cast(1. / momentum), boost::test_tools::tolerance(1.e-3F)); } - const auto full = fitMeasurements(kPropagator, trackingGeometry(), kNoMaterial, copied.data(), static_cast(kLayers)); - BOOST_TEST(tracks.chi2[track] == full.chi2); - - BOOST_TEST(tracks.ndf[track] == (2 * static_cast(kLayers)) - 5); - // On the innermost surface, and the generating helix. - const auto& state = tracks.state[track]; - const auto first = fixture.frame->measurement(rows.index[0]); - BOOST_TEST(state.surface == first.surface); - BOOST_TEST(std::abs(state.loc0() - first.loc0) < 1.e-4F); - BOOST_TEST(std::abs(state.loc1() - first.loc1) < 1.e-4F); - BOOST_TEST(state.theta() == static_cast(std::atan2(1., h.slope)), boost::test_tools::tolerance(1.e-4F)); - const double momentum = h.ptGeV * std::hypot(1., h.slope); - BOOST_TEST(std::abs(state.qOverP()) == static_cast(1. / momentum), boost::test_tools::tolerance(1.e-3F)); } - } - for (const int threads : {1, 2, 4}) { - detail::TaskArena cpuArena{threads}; - CpuStage::template RecordT> cpu{cpuArena}; - cpu.setMemoryPool(std::make_shared()); - TrackRefit parallel; - parallel.setMemoryPool(std::make_shared()); - BOOST_REQUIRE(parallel.run([&](std::uint32_t elements, const auto& body) { cpu.forEachElement(elements, 1U, body); }, selection, iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, 0U) == StageStatus::kComplete); - // Bit for bit. - const auto other = parallel.tracks().view(); - const auto bits = [](Scalar value) { return std::bit_cast(value); }; - for (std::uint32_t track = 0; track < kTracks; ++track) { - const auto& a = other.state[track]; - const auto& b = tracks.state[track]; - BOOST_TEST(std::ranges::equal(a.par, b.par, {}, bits, bits)); - BOOST_TEST(std::ranges::equal(a.cov, b.cov, {}, bits, bits)); - BOOST_TEST(other.timePS[track] == tracks.timePS[track]); - BOOST_TEST(other.widthPS[track] == tracks.widthPS[track]); - BOOST_TEST(a.surface == b.surface); - BOOST_TEST(bits(other.chi2[track]) == bits(tracks.chi2[track])); - BOOST_TEST(std::ranges::equal(other.clusters[track].index, tracks.clusters[track].index)); + for (const int threads : {1, 2, 4}) { + detail::TaskArena cpuArena{threads}; + CpuStage::template RecordT> cpu{cpuArena}; + cpu.setMemoryPool(std::make_shared()); + TrackRefit parallel; + parallel.setMemoryPool(std::make_shared()); + BOOST_REQUIRE(parallel.run([&](std::uint32_t elements, const auto& body) { cpu.forEachElement(elements, 1U, body); }, selection, iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, 0U) == StageStatus::kComplete); + // Bit for bit. + const auto other = parallel.tracks().view(); + const auto bits = [](Scalar value) { return std::bit_cast(value); }; + for (std::uint32_t track = 0; track < kTracks; ++track) { + const auto& a = other.state[track]; + const auto& b = tracks.state[track]; + BOOST_TEST(std::ranges::equal(a.par, b.par, {}, bits, bits)); + BOOST_TEST(std::ranges::equal(a.cov, b.cov, {}, bits, bits)); + BOOST_TEST(other.timePS[track] == tracks.timePS[track]); + BOOST_TEST(other.widthPS[track] == tracks.widthPS[track]); + BOOST_TEST(a.surface == b.surface); + BOOST_TEST(bits(other.chi2[track]) == bits(tracks.chi2[track])); + BOOST_TEST(std::ranges::equal(other.clusters[track].index, tracks.clusters[track].index)); + } } } } diff --git a/Detectors/GenericTracking/tracking/test/testTrackFitter.cxx b/Detectors/GenericTracking/tracking/test/testTrackFitter.cxx index 1d49207d7c256..3569abfaf4c50 100644 --- a/Detectors/GenericTracking/tracking/test/testTrackFitter.cxx +++ b/Detectors/GenericTracking/tracking/test/testTrackFitter.cxx @@ -216,19 +216,24 @@ BOOST_AUTO_TEST_CASE(ThreePointSeedRemainsFiniteAsCurvatureDiverges) BOOST_AUTO_TEST_CASE(IndependentPassesFitTheSameNoiselessChain) { - const SyntheticChain chain; - const UniformMaterial material; - const auto fit = fitMeasurements(chain.propagator, chain.tracking, material, chain.measurements.data(), static_cast(chain.measurements.size())); - - BOOST_TEST(static_cast(fit.status) == static_cast(PropagationStatus::kOk)); - BOOST_TEST(fit.ndf == 1); - BOOST_TEST(fit.outwardState.surface == chain.truth[2].surface); - BOOST_TEST(fit.state.surface == chain.truth[0].surface); - BOOST_TEST(std::abs(fit.state.phi() - chain.truth[0].phi()) < 1e-5F); - BOOST_TEST(std::abs(fit.state.theta() - chain.truth[0].theta()) < 1e-5F); - BOOST_TEST(std::abs(fit.state.qOverP() - chain.truth[0].qOverP()) < 1e-3F); - BOOST_TEST(std::abs(fit.chi2) < 1e-4F); - BOOST_TEST(std::abs(fit.outwardChi2) < 1e-4F); + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) + { + const SyntheticChain chain; + const UniformMaterial material; + const auto fit = fitMeasurementsReference(chain.propagator, chain.tracking, material, chain.measurements.data(), static_cast(chain.measurements.size()), referenceMode); + + BOOST_TEST(static_cast(fit.status) == static_cast(PropagationStatus::kOk)); + BOOST_TEST(fit.ndf == 1); + BOOST_TEST(fit.outwardState.surface == chain.truth[2].surface); + BOOST_TEST(fit.state.surface == chain.truth[0].surface); + BOOST_TEST(std::abs(fit.state.phi() - chain.truth[0].phi()) < 1e-5F); + BOOST_TEST(std::abs(fit.state.theta() - chain.truth[0].theta()) < 1e-5F); + BOOST_TEST(std::abs(fit.state.qOverP() - chain.truth[0].qOverP()) < 1e-3F); + BOOST_TEST(std::abs(fit.chi2) < 1e-4F); + BOOST_TEST(std::abs(fit.outwardChi2) < 1e-4F); + } + } } /// A primary from a vertex in the diamond's tail crosses the far end of barrel @@ -239,105 +244,115 @@ BOOST_AUTO_TEST_CASE(IndependentPassesFitTheSameNoiselessChain) /// last. BOOST_AUTO_TEST_CASE(LayerZeroAtShallowIncidenceFitsBothWays) { - using Mixed = detectors::ToyBarrel4Disc8; - const auto detector = toyGeometry(); - const TrackingGeometry tracking{detector}; - const Propagator propagator; - const UniformMaterial material; - - SensorID start = firstSensorOfLayer(0); - for (SensorID sensor = start; sensor < firstSensorOfLayer(1); ++sensor) { - if (tracking.surface(sensor).origin[2] > tracking.surface(start).origin[2]) { - start = sensor; - } - } - const auto& first = tracking.surface(start); - std::array truth; - truth[0].par[kPhi] = std::atan2(first.origin[1], first.origin[0]); - truth[0].par[kTheta] = 0.048F; - truth[0].par[kQOverP] = 0.05F; - truth[0].surface = start; - const double cosIncidence = std::sin(0.048) * (std::cos(static_cast(truth[0].phi())) * static_cast(first.normal[0]) + std::sin(static_cast(truth[0].phi())) * static_cast(first.normal[1])); - BOOST_TEST_REQUIRE(std::abs(cosIncidence) < 0.05); - BOOST_TEST_REQUIRE(std::abs(cosIncidence) > static_cast(propagator.minAbsCosIncidence)); - - constexpr std::array kDiscs{8, 9, 10, 11}; - for (std::size_t i = 0; i < kDiscs.size(); ++i) { - const auto& from = tracking.surface(truth[i].surface); - for (SensorID sensor = firstSensorOfLayer(kDiscs[i]); sensor < firstSensorOfLayer(static_cast(kDiscs[i] + 1U)); ++sensor) { - BoundTrackState probe = truth[i]; - if (propagate(propagator, probe, from, tracking.surface(sensor), PropagationDirection::kForward).status == PropagationStatus::kOk) { - truth[i + 1U] = probe; - break; + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) + { + using Mixed = detectors::ToyBarrel4Disc8; + const auto detector = toyGeometry(); + const TrackingGeometry tracking{detector}; + const Propagator propagator; + const UniformMaterial material; + + SensorID start = firstSensorOfLayer(0); + for (SensorID sensor = start; sensor < firstSensorOfLayer(1); ++sensor) { + if (tracking.surface(sensor).origin[2] > tracking.surface(start).origin[2]) { + start = sensor; + } + } + const auto& first = tracking.surface(start); + std::array truth; + truth[0].par[kPhi] = std::atan2(first.origin[1], first.origin[0]); + truth[0].par[kTheta] = 0.048F; + truth[0].par[kQOverP] = 0.05F; + truth[0].surface = start; + const double cosIncidence = std::sin(0.048) * (std::cos(static_cast(truth[0].phi())) * static_cast(first.normal[0]) + std::sin(static_cast(truth[0].phi())) * static_cast(first.normal[1])); + BOOST_TEST_REQUIRE(std::abs(cosIncidence) < 0.05); + BOOST_TEST_REQUIRE(std::abs(cosIncidence) > static_cast(propagator.minAbsCosIncidence)); + + constexpr std::array kDiscs{8, 9, 10, 11}; + for (std::size_t i = 0; i < kDiscs.size(); ++i) { + const auto& from = tracking.surface(truth[i].surface); + for (SensorID sensor = firstSensorOfLayer(kDiscs[i]); sensor < firstSensorOfLayer(static_cast(kDiscs[i] + 1U)); ++sensor) { + BoundTrackState probe = truth[i]; + if (propagate(propagator, probe, from, tracking.surface(sensor), PropagationDirection::kForward).status == PropagationStatus::kOk) { + truth[i + 1U] = probe; + break; + } + } + BOOST_TEST_REQUIRE(truth[i + 1U].surface != truth[i].surface); + } + std::array measurements; + for (std::size_t i = 0; i < measurements.size(); ++i) { + measurements[i] = measurementFromState(truth[i], static_cast(i)); } - } - BOOST_TEST_REQUIRE(truth[i + 1U].surface != truth[i].surface); - } - std::array measurements; - for (std::size_t i = 0; i < measurements.size(); ++i) { - measurements[i] = measurementFromState(truth[i], static_cast(i)); - } - const auto fit = fitMeasurements(propagator, tracking, material, measurements.data(), static_cast(measurements.size())); - BOOST_TEST(static_cast(fit.status) == static_cast(PropagationStatus::kOk)); - BOOST_TEST(fit.state.surface == start); - BOOST_TEST(std::abs(fit.state.theta() - truth[0].theta()) < 1e-5F); + const auto fit = fitMeasurementsReference(propagator, tracking, material, measurements.data(), static_cast(measurements.size()), referenceMode); + BOOST_TEST(static_cast(fit.status) == static_cast(PropagationStatus::kOk)); + BOOST_TEST(fit.state.surface == start); + BOOST_TEST(std::abs(fit.state.theta() - truth[0].theta()) < 1e-5F); - const Propagator former{.minAbsCosIncidence = 0.05F}; - const auto failed = fitMeasurements(former, tracking, material, measurements.data(), static_cast(measurements.size())); - BOOST_TEST(static_cast(failed.status) == static_cast(PropagationStatus::kGrazing)); - BOOST_TEST(static_cast(failed.failedPass) == static_cast(FitPass::kInward)); - BOOST_TEST(failed.failedMeasurement == 0); + const Propagator former{.minAbsCosIncidence = 0.05F}; + const auto failed = fitMeasurementsReference(former, tracking, material, measurements.data(), static_cast(measurements.size()), referenceMode); + BOOST_TEST(static_cast(failed.status) == static_cast(PropagationStatus::kGrazing)); + BOOST_TEST(static_cast(failed.failedPass) == static_cast(FitPass::kInward)); + BOOST_TEST(failed.failedMeasurement == 0); + } + } } BOOST_AUTO_TEST_CASE(SmearedChainsGiveUnbiasedParametersAndUnitPulls) { - // A fit can be self-consistent and still wrong in - // exactly one way that matters here: if the inward pass inherited the - // outward pass's covariance it would consume every measurement twice, the - // estimate would be too confident by root two, and the pull widths would - // come out near 0.707 while every other test in this file still passed. - // That number is asserted against explicitly. - const SyntheticChain chain; - const UniformMaterial material; - - constexpr int kTrials = 3000; - const Scalar sigmaLoc = 1.5811388e-3F; // sqrt(2.5e-6) cm - - std::mt19937 rng(20260919); - std::normal_distribution gauss(0., 1.); - - std::array, kNBoundParams> pull; - - for (int trial = 0; trial < kTrials; ++trial) { - auto smeared = chain.measurements; - for (auto& m : smeared) { - m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - } - - const auto fit = fitMeasurements(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size())); - BOOST_REQUIRE(fit.status == PropagationStatus::kOk); - - collectPulls(trial, fit.state, chain.truth[0], pull); - } + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) + { + // A fit can be self-consistent and still wrong in + // exactly one way that matters here: if the inward pass inherited the + // outward pass's covariance it would consume every measurement twice, the + // estimate would be too confident by root two, and the pull widths would + // come out near 0.707 while every other test in this file still passed. + // That number is asserted against explicitly. + const SyntheticChain chain; + const UniformMaterial material; + + constexpr int kTrials = 3000; + const Scalar sigmaLoc = 1.5811388e-3F; // sqrt(2.5e-6) cm + + std::mt19937 rng(20260919); + std::normal_distribution gauss(0., 1.); + + std::array, kNBoundParams> pull; + + for (int trial = 0; trial < kTrials; ++trial) { + auto smeared = chain.measurements; + for (auto& m : smeared) { + m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + } - const char* names[kNBoundParams] = {"loc0", "loc1", "phi", "theta", "qOverP"}; - for (int i = 0; i < kNBoundParams; ++i) { - const auto& p = pull[static_cast(i)]; - const auto n = static_cast(p.size()); + const auto fit = fitMeasurementsReference(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size()), referenceMode); + BOOST_REQUIRE(fit.status == PropagationStatus::kOk); - const auto [meanPull, width] = moments(p); + collectPulls(trial, fit.state, chain.truth[0], pull); + } - BOOST_TEST_CONTEXT("parameter " << names[i] << " pull mean " << meanPull << " width " << width) - { - // The mean of kTrials unit-width pulls has standard error 1/sqrt(N), - // which is 0.018 here; four of those is a comfortable band that still - // catches any real bias. - BOOST_TEST(std::abs(meanPull) < 4. / std::sqrt(n)); - // The width itself is known to 1/sqrt(2N) = 1.3 %, so 10 % is generous - // -- and 0.707 is 29 % low, far outside it. - BOOST_TEST(std::abs(width - 1.) < 0.10); + const char* names[kNBoundParams] = {"loc0", "loc1", "phi", "theta", "qOverP"}; + for (int i = 0; i < kNBoundParams; ++i) { + const auto& p = pull[static_cast(i)]; + const auto n = static_cast(p.size()); + + const auto [meanPull, width] = moments(p); + + BOOST_TEST_CONTEXT("parameter " << names[i] << " pull mean " << meanPull << " width " << width) + { + // The mean of kTrials unit-width pulls has standard error 1/sqrt(N), + // which is 0.018 here; four of those is a comfortable band that still + // catches any real bias. + BOOST_TEST(std::abs(meanPull) < 4. / std::sqrt(n)); + // The width itself is known to 1/sqrt(2N) = 1.3 %, so 10 % is generous + // -- and 0.707 is 29 % low, far outside it. + BOOST_TEST(std::abs(width - 1.) < 0.10); + } + } } } } @@ -547,64 +562,69 @@ BOOST_AUTO_TEST_CASE(ChainSeedOfAStraightTrackHasNoCurvature) /// is. BOOST_AUTO_TEST_CASE(LowMomentumChainsFitWithUnitPulls) { - constexpr Scalar kQOverP = 2.5F; - constexpr Scalar kX2X0PerCm = 1.5e-3F; - constexpr std::size_t kLayers = 4; - BarrelChain chain(kQOverP); - BOOST_REQUIRE(chain.built); - const UniformMaterial material{.x2X0PerCm = kX2X0PerCm}; - - constexpr int kTrials = 1000; - const Scalar sigmaLoc = 1.5811388e-3F; // sqrt(2.5e-6) cm - - std::mt19937 rng(20260923); - std::normal_distribution gauss(0., 1.); - std::array, kNBoundParams> inwardPull; - std::array, kNBoundParams> outwardPull; - double inwardChi2PerDof = 0.; - double outwardChi2PerDof = 0.; - - for (int trial = 0; trial < kTrials; ++trial) { - // A kink can send the track through the 0.5 mm gap between two modules; - // draw again rather than fit a chain with a missing layer. - for (int attempt = 0; attempt < 10; ++attempt) { - chain.generate(kQOverP, kX2X0PerCm, &rng); - if (chain.built) { - break; + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) + { + constexpr Scalar kQOverP = 2.5F; + constexpr Scalar kX2X0PerCm = 1.5e-3F; + constexpr std::size_t kLayers = 4; + BarrelChain chain(kQOverP); + BOOST_REQUIRE(chain.built); + const UniformMaterial material{.x2X0PerCm = kX2X0PerCm}; + + constexpr int kTrials = 1000; + const Scalar sigmaLoc = 1.5811388e-3F; // sqrt(2.5e-6) cm + + std::mt19937 rng(20260923); + std::normal_distribution gauss(0., 1.); + std::array, kNBoundParams> inwardPull; + std::array, kNBoundParams> outwardPull; + double inwardChi2PerDof = 0.; + double outwardChi2PerDof = 0.; + + for (int trial = 0; trial < kTrials; ++trial) { + // A kink can send the track through the 0.5 mm gap between two modules; + // draw again rather than fit a chain with a missing layer. + for (int attempt = 0; attempt < 10; ++attempt) { + chain.generate(kQOverP, kX2X0PerCm, &rng); + if (chain.built) { + break; + } + } + BOOST_REQUIRE(chain.built); + auto smeared = chain.measurements; + for (auto& m : smeared) { + m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + } + const auto fit = fitMeasurementsReference(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size()), referenceMode); + BOOST_TEST_CONTEXT("trial " << trial << ", propagation status " << static_cast(fit.status)) { BOOST_REQUIRE(fit.status == PropagationStatus::kOk); } + collectPulls(trial, fit.state, chain.truth[0], inwardPull); + collectPulls(trial, fit.outwardState, chain.truth[kLayers - 1], outwardPull); + inwardChi2PerDof += static_cast(fit.chi2) / fit.ndf / kTrials; + outwardChi2PerDof += static_cast(fit.outwardChi2) / fit.ndf / kTrials; } - } - BOOST_REQUIRE(chain.built); - auto smeared = chain.measurements; - for (auto& m : smeared) { - m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - } - const auto fit = fitMeasurements(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size())); - BOOST_TEST_CONTEXT("trial " << trial << ", propagation status " << static_cast(fit.status)) { BOOST_REQUIRE(fit.status == PropagationStatus::kOk); } - collectPulls(trial, fit.state, chain.truth[0], inwardPull); - collectPulls(trial, fit.outwardState, chain.truth[kLayers - 1], outwardPull); - inwardChi2PerDof += static_cast(fit.chi2) / fit.ndf / kTrials; - outwardChi2PerDof += static_cast(fit.outwardChi2) / fit.ndf / kTrials; - } - const char* names[kNBoundParams] = {"loc0", "loc1", "phi", "theta", "qOverP"}; - for (const auto* pass : {&inwardPull, &outwardPull}) { - for (int i = 0; i < kNBoundParams; ++i) { - const auto& p = (*pass)[static_cast(i)]; - const auto n = static_cast(p.size()); - const auto [mean, width] = moments(p); - BOOST_TEST_CONTEXT((pass == &inwardPull ? "inward" : "outward") << " parameter " << names[i] << " pull mean " << mean << " width " << width) + const char* names[kNBoundParams] = {"loc0", "loc1", "phi", "theta", "qOverP"}; + for (const auto* pass : {&inwardPull, &outwardPull}) { + for (int i = 0; i < kNBoundParams; ++i) { + const auto& p = (*pass)[static_cast(i)]; + const auto n = static_cast(p.size()); + const auto [mean, width] = moments(p); + BOOST_TEST_CONTEXT((pass == &inwardPull ? "inward" : "outward") << " parameter " << names[i] << " pull mean " << mean << " width " << width) + { + BOOST_TEST(std::abs(mean) < 4. / std::sqrt(n)); + BOOST_TEST(std::abs(width - 1.) < 0.10); + } + } + } + BOOST_TEST_CONTEXT("inward chi2 per spatial dof " << inwardChi2PerDof << ", outward " << outwardChi2PerDof) { - BOOST_TEST(std::abs(mean) < 4. / std::sqrt(n)); - BOOST_TEST(std::abs(width - 1.) < 0.10); + BOOST_TEST(std::abs(inwardChi2PerDof - 1.) < 0.15); + BOOST_TEST(std::abs(outwardChi2PerDof - 1.) < 0.15); } } } - BOOST_TEST_CONTEXT("inward chi2 per spatial dof " << inwardChi2PerDof << ", outward " << outwardChi2PerDof) - { - BOOST_TEST(std::abs(inwardChi2PerDof - 1.) < 0.15); - BOOST_TEST(std::abs(outwardChi2PerDof - 1.) < 0.15); - } } /// The same chains losing energy: the chain seed is a circle through the @@ -617,61 +637,66 @@ BOOST_AUTO_TEST_CASE(LowMomentumChainsFitWithUnitPulls) /// the fit's model of where the loss is taken. BOOST_AUTO_TEST_CASE(LowMomentumChainsLosingEnergyFitWithUnitPulls) { - constexpr Scalar kQOverP = 2.5F; - constexpr Scalar kX2X0PerCm = 1.5e-3F; - constexpr Scalar kRho = 0.3F; // about 2.5 % of the momentum over the chain - constexpr std::size_t kLayers = 4; - BarrelChain chain(kQOverP); - BOOST_REQUIRE(chain.built); - const UniformMaterial material{.x2X0PerCm = kX2X0PerCm, .rho = kRho}; - - constexpr int kTrials = 1000; - const Scalar sigmaLoc = 1.5811388e-3F; - std::mt19937 rng(20261021); - std::normal_distribution gauss(0., 1.); - std::array, kNBoundParams> inwardPull; - double inwardChi2PerDof = 0.; - double outwardChi2PerDof = 0.; - double lostShare = 0.; - for (int trial = 0; trial < kTrials; ++trial) { - for (int attempt = 0; attempt < 10; ++attempt) { - chain.generate(kQOverP, kX2X0PerCm, &rng, kRho); - if (chain.built) { - break; - } - } - BOOST_REQUIRE(chain.built); - lostShare += (1. - static_cast(kQOverP / chain.truth[kLayers - 1].qOverP())) / kTrials; - auto smeared = chain.measurements; - for (auto& m : smeared) { - m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - } - const auto fit = fitMeasurements(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size())); - BOOST_TEST_CONTEXT("trial " << trial << ", propagation status " << static_cast(fit.status)) { BOOST_REQUIRE(fit.status == PropagationStatus::kOk); } - collectPulls(trial, fit.state, chain.truth[0], inwardPull); - inwardChi2PerDof += static_cast(fit.chi2) / fit.ndf / kTrials; - outwardChi2PerDof += static_cast(fit.outwardChi2) / fit.ndf / kTrials; - } - BOOST_TEST_MESSAGE("momentum lost over the chain " << lostShare << ", chi2 per dof inward " << inwardChi2PerDof << " outward " << outwardChi2PerDof); - BOOST_TEST(lostShare > 0.01); - for (int i = 0; i < kNBoundParams; ++i) { - const auto& p = inwardPull[static_cast(i)]; - const bool held = i != o2::gtrk::kQOverP; - const auto n = static_cast(p.size()); - const auto [mean, width] = moments(p); - BOOST_TEST_MESSAGE("inward parameter " << i << " pull mean " << mean << " width " << width); - BOOST_TEST_CONTEXT("inward parameter " << i << " pull mean " << mean << " width " << width) + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) { - BOOST_TEST((!held || std::abs(mean) < 4. / std::sqrt(n))); - BOOST_TEST((!held || std::abs(width - 1.) < 0.10)); + constexpr Scalar kQOverP = 2.5F; + constexpr Scalar kX2X0PerCm = 1.5e-3F; + constexpr Scalar kRho = 0.3F; // about 2.5 % of the momentum over the chain + constexpr std::size_t kLayers = 4; + BarrelChain chain(kQOverP); + BOOST_REQUIRE(chain.built); + const UniformMaterial material{.x2X0PerCm = kX2X0PerCm, .rho = kRho}; + + constexpr int kTrials = 1000; + const Scalar sigmaLoc = 1.5811388e-3F; + std::mt19937 rng(20261021); + std::normal_distribution gauss(0., 1.); + std::array, kNBoundParams> inwardPull; + double inwardChi2PerDof = 0.; + double outwardChi2PerDof = 0.; + double lostShare = 0.; + for (int trial = 0; trial < kTrials; ++trial) { + for (int attempt = 0; attempt < 10; ++attempt) { + chain.generate(kQOverP, kX2X0PerCm, &rng, kRho); + if (chain.built) { + break; + } + } + BOOST_REQUIRE(chain.built); + lostShare += (1. - static_cast(kQOverP / chain.truth[kLayers - 1].qOverP())) / kTrials; + auto smeared = chain.measurements; + for (auto& m : smeared) { + m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + } + const auto fit = fitMeasurementsReference(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size()), referenceMode); + BOOST_TEST_CONTEXT("trial " << trial << ", propagation status " << static_cast(fit.status)) { BOOST_REQUIRE(fit.status == PropagationStatus::kOk); } + collectPulls(trial, fit.state, chain.truth[0], inwardPull); + inwardChi2PerDof += static_cast(fit.chi2) / fit.ndf / kTrials; + outwardChi2PerDof += static_cast(fit.outwardChi2) / fit.ndf / kTrials; + } + BOOST_TEST_MESSAGE("momentum lost over the chain " << lostShare << ", chi2 per dof inward " << inwardChi2PerDof << " outward " << outwardChi2PerDof); + BOOST_TEST(lostShare > 0.01); + for (int i = 0; i < kNBoundParams; ++i) { + const auto& p = inwardPull[static_cast(i)]; + const bool held = i != o2::gtrk::kQOverP; + const auto n = static_cast(p.size()); + const auto [mean, width] = moments(p); + BOOST_TEST_MESSAGE("inward parameter " << i << " pull mean " << mean << " width " << width); + BOOST_TEST_CONTEXT("inward parameter " << i << " pull mean " << mean << " width " << width) + { + BOOST_TEST((!held || std::abs(mean) < 4. / std::sqrt(n))); + BOOST_TEST((!held || std::abs(width - 1.) < 0.10)); + } + } + BOOST_TEST_CONTEXT("inward chi2 per spatial dof " << inwardChi2PerDof << ", outward " << outwardChi2PerDof) + { + BOOST_TEST(std::abs(inwardChi2PerDof - 1.) < 0.15); + BOOST_TEST(std::abs(outwardChi2PerDof - 1.) < 0.15); + } } } - BOOST_TEST_CONTEXT("inward chi2 per spatial dof " << inwardChi2PerDof << ", outward " << outwardChi2PerDof) - { - BOOST_TEST(std::abs(inwardChi2PerDof - 1.) < 0.15); - BOOST_TEST(std::abs(outwardChi2PerDof - 1.) < 0.15); - } } BOOST_AUTO_TEST_CASE(PerigeePhiIsTheDirectionAtTheTransverseClosestApproach) diff --git a/Detectors/GenericTracking/workflow/CMakeLists.txt b/Detectors/GenericTracking/workflow/CMakeLists.txt index 67af4f0273036..fcb9eb0b4d33c 100644 --- a/Detectors/GenericTracking/workflow/CMakeLists.txt +++ b/Detectors/GenericTracking/workflow/CMakeLists.txt @@ -15,11 +15,11 @@ # unit, ALICE 3's (TRK.cxx) with the upgrades, and the tracking device # (TrackerDevice.h), one unit per detector, each instantiating the tracker. set(sources src/Clusters.cxx src/Placement.cxx src/RecoParam.cxx src/ITS.cxx src/MFT.cxx src/TrackerITS.cxx src/TrackerMFT.cxx src/AdaptorITS.cxx - src/AdaptorMFT.cxx src/GPURecoTracking.cxx) + src/AdaptorMFT.cxx src/GPURecoTracking.cxx src/ITSROFSelection.cxx) set(trackerUnits src/TrackerITS.cxx src/TrackerMFT.cxx) set(publicLibraries O2::GTrkTracking O2::Framework O2::CommonUtils O2::DataFormatsCalibration O2::DetectorsBase O2::DataFormatsITSMFT O2::ITSBase O2::MFTBase O2::DetectorsRaw ROOT::Core) -set(privateLibraries O2::ITSMFTBase O2::ITSMFTTracking O2::DataFormatsITS O2::DataFormatsMFT O2::SimulationDataFormat O2::CommonConstants +set(privateLibraries O2::ITStracking O2::DataFormatsTRD O2::ITSMFTBase O2::ITSMFTTracking O2::DataFormatsITS O2::DataFormatsMFT O2::SimulationDataFormat O2::CommonConstants O2::MathUtils O2::ReconstructionDataFormats ROOT::Geom ROOT::RIO ROOT::Tree) if(ENABLE_UPGRADES) list(APPEND sources src/TRK.cxx src/TrackerTRKFT3.cxx) diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/AdaptorSpec.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/AdaptorSpec.h index e3a0602f8bc90..2340881c6f6d4 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/AdaptorSpec.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/AdaptorSpec.h @@ -25,8 +25,10 @@ namespace o2::gtrk::workflow /// The adaptor of \a Det (the ITS or the MFT); with \a useMC the track labels /// from the clusters' (and, for the ITS, the vertices' outputs, without truth). +/// With staggered input, subscribe to each layer and publish composed O2 cluster +/// references; GTRK measurement indices remain dense in layer-stream order. template -o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC); +o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC, bool staggered = false); } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ClusterConverterSpec.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ClusterConverterSpec.h index e3d64259b0687..19b30553bca53 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ClusterConverterSpec.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ClusterConverterSpec.h @@ -24,10 +24,9 @@ namespace o2::gtrk::workflow { /// Defined in the detector's own unit (ITS.cxx, MFT.cxx, TRK.cxx). The ITS's -/// and MFT's readout frames must not be staggered: only subspecification 0 is -/// read. +/// and MFT's staggered input has one stream per layer, concatenated in layer order. template -o2::framework::DataProcessorSpec getClusterConverterSpec(bool useMC); +o2::framework::DataProcessorSpec getClusterConverterSpec(bool useMC, bool staggered = false); } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/Conditions.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/Conditions.h index 244febffe34d3..667d32a25b5bf 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/Conditions.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/Conditions.h @@ -157,24 +157,20 @@ class Conditions return false; } - /// The readout frames of the run's `DPLAlpideParam`, as O2's ITS tracker - /// takes them (`ITSTrackingInterface::updateTimeDependentParams`): frame `k` - /// at `k * length + delay + bias` BC, one per layer, every layer the same. - /// Throws std::runtime_error for staggered frames: the converter reads one - /// ROF stream. + /// Frame k starts at k * length + delay + bias BC, independently per layer. [[nodiscard]] std::vector frameTiming() const { const auto& param = AlpideParam::Instance(); - const FrameTiming result{.lengthBC = param.getROFLengthInBC(0), .firstStartBC = param.getROFDelayInBC(0) + param.getROFBiasInBC(0)}; - for (int layer = 1; layer < AlpideParam::getNLayers(); ++layer) { - if (param.getROFLengthInBC(layer) != result.lengthBC || param.getROFDelayInBC(layer) + param.getROFBiasInBC(layer) != result.firstStartBC) { - throw std::runtime_error(std::format("{}'s readout frames are staggered (layer {}): not read yet", Det::name, layer)); + std::vector result; + result.reserve(Detector::nLayers); + for (int layer = 0; layer < Detector::nLayers; ++layer) { + const auto length = param.getROFLengthInBC(layer); + if (length <= 0) { + throw std::runtime_error(std::format("{} layer {}: AlpideParam gives readout frames of {} BC", Det::name, layer, length)); } + result.push_back({.lengthBC = length, .firstStartBC = param.getROFDelayInBC(layer) + param.getROFBiasInBC(layer)}); } - if (result.lengthBC <= 0) { - throw std::runtime_error(std::format("{}'s AlpideParam gives readout frames of {} BC", Det::name, result.lengthBC)); - } - return std::vector(static_cast(Detector::nLayers), result); + return result; } /// The run's orbits per timeframe, its GRPECS's, as O2's ITS tracking diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h index 271c8b4b5130b..bf8c736a25fa3 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h @@ -102,6 +102,14 @@ inline constexpr std::string_view kVertexSeeds{"VTX"}; /// Each of those vertices' contributors, a std::uint32_t per vertex. inline constexpr std::string_view kVertexContributors{"VTXCONTR"}; +/// MC contributor identities: TripletClaim records with measurement rows and +/// indices into kVertexSeeds. Empty when MC propagation is disabled. +inline constexpr std::string_view kVertexClaims{"VTXCLAIMS"}; + +/// ITS clock-ROF acceptance after cluster/random and vertex multiplicity cuts. +/// One uint8_t per clock ROF; empty for a dropped/empty timeframe. +inline constexpr std::string_view kSelectedROFs{"ROFSEL"}; + /// Detector \a Det's name in lower case: "its", "mft", "trk". template [[nodiscard]] std::string detectorName() diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GPURecoTracking.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GPURecoTracking.h index 507e7d3ab1e2a..812b0d3d18ff2 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GPURecoTracking.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GPURecoTracking.h @@ -73,7 +73,7 @@ class GPURecoTracking /// The tracking of \a detector, "its" or "mft"; std::invalid_argument for /// another. -[[nodiscard]] std::unique_ptr makeGPURecoTracking(std::string_view detector); +[[nodiscard]] std::unique_ptr makeGPURecoTracking(std::string_view detector, bool useMC = false); } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GridConfig.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GridConfig.h index 1132b2e8cbe6d..91154eb704a8c 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GridConfig.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GridConfig.h @@ -108,6 +108,9 @@ namespace o2::gtrk::workflow template inline void checkRecoParam(const RecoSettings& param) { + if (param.fitReference < 0 || param.fitReference > 2) { + throw std::invalid_argument("GTRKReco.fitReference must be 0 (filtered state), 1 (circle outward), or 2 (circle outward, outward-fit reference inward)"); + } const auto margins = param.frameMarginBC.size(); if (margins != 0U && margins != static_cast(Detector::nLayers)) { throw std::invalid_argument(std::format("{}.frameMarginBC has {} entries for {}'s {} layers", RecoParamOf::type::Instance().getName(), margins, Det::name, Detector::nLayers)); @@ -253,6 +256,7 @@ inline void checkRecoParam(const RecoSettings& param) .vertexSource = entry(param.iterationVertexSource, static_cast(defaults.vertexSource)) == 0 ? VertexSource::kBeamLine : VertexSource::kVertexList, .minClusters = entry(param.iterationMinClusters, defaults.minClusters), .maxKinks = entry(param.iterationMaxKinks, defaults.maxKinks), + .fitReference = static_cast(param.fitReference), .kinkSigma = kinkSigma, .graphKinkSigma = entryDouble(param.iterationGraphKinkSigma, defaults.graphKinkSigma), .skipLinks = flag(param.iterationSkipLinks, defaults.skipLinks), diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSMFTOptions.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSMFTOptions.h new file mode 100644 index 0000000000000..4412af16c03cc --- /dev/null +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSMFTOptions.h @@ -0,0 +1,39 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef GTRK_WORKFLOW_ITSMFTOPTIONS_H_ +#define GTRK_WORKFLOW_ITSMFTOPTIONS_H_ + +#include "GenericTracking/Detectors/ITS.h" +#include "GenericTracking/Detectors/MFT.h" +#include "DataFormatsITSMFT/DPLAlpideParamInitializer.h" +#include "Framework/ConfigContext.h" + +#include + +namespace o2::gtrk::workflow +{ + +template +[[nodiscard]] bool staggeredInput(const o2::framework::ConfigContext& context) +{ + if constexpr (std::same_as) { + return o2::itsmft::DPLAlpideParamInitializer::isITSStaggeringEnabled(context); + } else if constexpr (std::same_as) { + return o2::itsmft::DPLAlpideParamInitializer::isMFTStaggeringEnabled(context); + } else { + return false; + } +} + +} // namespace o2::gtrk::workflow + +#endif // GTRK_WORKFLOW_ITSMFTOPTIONS_H_ diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSROFSelection.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSROFSelection.h new file mode 100644 index 0000000000000..f706a9ca7f799 --- /dev/null +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSROFSelection.h @@ -0,0 +1,68 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef GTRK_WORKFLOW_ITSROFSELECTION_H_ +#define GTRK_WORKFLOW_ITSROFSELECTION_H_ + +#include "GenericTracking/Base/FrameTiming.h" +#include "GenericTracking/Tracking/Measurement.h" +#include "GenericTracking/Tracking/PairCut.h" +#include "Framework/InputSpec.h" + +#include +#include +#include +#include + +namespace o2::framework +{ +class ProcessingContext; +} + +namespace o2::gtrk::workflow +{ + +/// ITS's native fastMultConfig selection: cluster/random selection before +/// vertexing, vertex-contributor selection before tracking. It consumes raw +/// cluster ROFs so that expanded measurement time margins cannot change a +/// cluster's ROF identity. Both host and device receive the selected rows. +class ITSROFSelection +{ + public: + static void addInputs(std::vector& inputs); + ITSROFSelection(o2::framework::ProcessingContext& pc, std::span timing, int orbitsPerTF); + ~ITSROFSelection(); + ITSROFSelection(const ITSROFSelection&) = delete; + ITSROFSelection& operator=(const ITSROFSelection&) = delete; + + /// Vertices remain in the output, as in native ITS; the cut disables their + /// overlapping ROFs for tracking. Returns whether a vertex cut is active. + bool selectWithVertices(std::span seeds, std::span contributors); + [[nodiscard]] std::vector clockMask() const; + + struct Measurements { + MeasurementStore columns; + std::vector offsets; + std::vector originalRows; ///< selected row -> original input measurement row + }; + + /// Preserve bucket ordering and original measurement identities. No copy + /// when every row passes. originalRows composes mappings after vertex cuts. + [[nodiscard]] std::optional filter(MeasurementView source, std::span offsets, std::span originalRows = {}) const; + + private: + struct State; + std::unique_ptr mState; +}; + +} // namespace o2::gtrk::workflow + +#endif // GTRK_WORKFLOW_ITSROFSELECTION_H_ diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h index 560bc3cd8da32..d4468644f8d33 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h @@ -64,6 +64,12 @@ struct RecoSettings { /// scattering, in the windows at pT_min and in the fit, and the fit's energy /// loss. The charged pion's by default. double massHypothesisGeV = o2::constants::physics::MassPionCharged; + /// Linearisation: 0 uses the evolving filtered state both ways (default). + /// 1 uses an independent circle reference outward and the filtered state inward. + /// 2 uses the circle reference outward and an independent reference seeded + /// from the outward fitted parameters inward. References are position-anchored + /// to clusters within each pass. Selected tracks reuse the road's outward fit. + int fitReference = 0; /// The uniform material per layer [X/X0] that sizes the windows and the /// grid, and the scattering when no material LUT is given; with a LUT the /// tracking takes its budgets. @@ -269,6 +275,10 @@ struct RecoSettings { /// empty products and is counted, instead of failing the run. Off, the /// default, as ITS's DropTFUponFailure: the run fails. bool dropTFUponFailure = false; + + /// ITS ROF selection uses native fastMultConfig. Trigger input: 0 none, + /// 1 ITS physics triggers, 2 TRD tracklet triggers. Other detectors ignore it. + int selectWithTriggers = 0; }; /// Detector \a N's settings (`o2::detectors::DetID`): `GTRKReco`. diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackerSpec.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackerSpec.h index dd6e3840b9b51..3458ab1a4f082 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackerSpec.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackerSpec.h @@ -37,14 +37,14 @@ using TrackRecord = typename TrackRefit::Tracks::Record; /// The device; with \a withVertices the native vertex seeds (kVertices) are /// an input. template -o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices); +o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices, bool useMC = false); class GPURecoTracking; /// \a Det's tracking inside the GPU reconstruction (GPURecoTracking.h), the /// same timeframe as the device's. template -std::unique_ptr makeGPURecoTrackingOf(); +std::unique_ptr makeGPURecoTrackingOf(bool useMC); } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h index fcb9c756acb4b..509e4d616e9f3 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h @@ -289,6 +289,7 @@ class TrackingInterface std::span seeds; const VertexerResult* found = nullptr; ///< the own vertexer's list and contributors, host or device; none without it const typename Vertexing::Result* pass = nullptr; ///< the own vertexer's on the host, with its claims; none without it or on a GPU + std::span claims; ///< the own vertexer's claimed triplets when requested double seconds = 0.; double planSeconds = 0.; ///< the initial device plan, before upload and vertex finding bool dropped = false; ///< out of memory, dropped (dropTFUponFailure): no tracking @@ -308,19 +309,19 @@ class TrackingInterface if (!beginDeviceTimeframe()) { return {.seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(), .planSeconds = planSeconds, .dropped = true}; } - const auto status = mDevice->findTimeframeVertices(frame, mDevicePlan, mDeviceIterations.size() - 1U, mVertexing->settings(), beamSpotConstraint(mBeamLine, param), mDeviceFound, mDeviceReport); + const auto status = mDevice->findTimeframeVertices(frame, mDevicePlan, mDeviceIterations.size() - 1U, mVertexing->settings(), beamSpotConstraint(mBeamLine, param), mDeviceFound, mDeviceReport, keepClaims ? &mDeviceClaims : nullptr); const double seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); if (status != StageStatus::kComplete) { return {.seconds = seconds, .planSeconds = planSeconds, .dropped = dropOrFail(status, "the vertexing pass")}; } - return {.seeds = mDeviceFound.vertices, .found = &mDeviceFound, .seconds = seconds, .planSeconds = planSeconds}; + return {.seeds = mDeviceFound.vertices, .found = &mDeviceFound, .claims = keepClaims ? std::span{mDeviceClaims} : std::span{}, .seconds = seconds, .planSeconds = planSeconds}; } mPass.emplace(mVertexing->find(frame, frame.grid().geometry().blocksPerSlice(param.sliceWidthPS), param.sliceHaloPS, keepClaims)); const double seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); if (mPass->status != StageStatus::kComplete) { return {.pass = &*mPass, .seconds = seconds, .dropped = dropOrFail(mPass->status, "the vertexing pass")}; } - return {.seeds = mPass->found.vertices, .found = &mPass->found, .pass = &*mPass, .seconds = seconds}; + return {.seeds = mPass->found.vertices, .found = &mPass->found, .pass = &*mPass, .claims = mPass->claims, .seconds = seconds}; } /// The tracker's run on \a frame from \a seeds. Valid until endTimeframe. @@ -363,6 +364,7 @@ class TrackingInterface mDevicePlan = {}; mDeviceTrackingPlan = {}; mDeviceFound = {}; + mDeviceClaims.clear(); mDeviceTimeframe = false; mVertexed = false; } @@ -495,9 +497,9 @@ class TrackingInterface [[nodiscard]] Tracked trackOnDevice(const Frame& frame, std::span seeds, std::int32_t blocksPerSlice, std::chrono::steady_clock::time_point start) { const auto& param = recoParam(); - const bool replaced = mVertexing != nullptr; + const bool replaced = mVertexed; const std::size_t iterations = mTracker->iterationCount(); - mDeviceTrackingPlan = replaced ? planTimeframe(frame, mDeviceIterations, mDeviceWindows, mTracker->beamLine(), seeds, blocksPerSlice, param.sliceHaloPS) : mTracker->planTimeframe(frame, seeds, blocksPerSlice, param.sliceHaloPS); + mDeviceTrackingPlan = mVertexing != nullptr ? planTimeframe(frame, mDeviceIterations, mDeviceWindows, mTracker->beamLine(), seeds, blocksPerSlice, param.sliceHaloPS) : mTracker->planTimeframe(frame, seeds, blocksPerSlice, param.sliceHaloPS); Tracked tracked; tracked.planSeconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); if (!beginDeviceTimeframe()) { @@ -553,6 +555,7 @@ class TrackingInterface /// tracks downloaded into the pool, the calls' report. TimeframePlan mDevicePlan; TimeframePlan mDeviceTrackingPlan; + std::vector mDeviceClaims; ///< optional MC contributor identities VertexerResult mDeviceFound; std::shared_ptr mDevicePool; typename TrackRefit::Tracks mDeviceTracks; diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/VertexLabel.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/VertexLabel.h new file mode 100644 index 0000000000000..6dd3fea663b70 --- /dev/null +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/VertexLabel.h @@ -0,0 +1,84 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef GTRK_WORKFLOW_VERTEXLABEL_H_ +#define GTRK_WORKFLOW_VERTEXLABEL_H_ + +#include "GenericTracking/Tracking/VertexStage.h" +#include "GenericTracking/Workflow/TrackLabel.h" + +#include +#include +#include +#include +#include +#include + +namespace o2::gtrk::workflow +{ + +struct VertexTruth { + o2::MCCompLabel label; + float purity = 0.F; +}; + +/// One vote per claimed, cluster-disjoint triplet, from the particle common +/// to its three measurements. Vote for its (source, event), ignoring track ID. +/// Mixed/noise/QED triplets cast no vote but remain in the purity denominator. +/// The largest vote wins (smallest source/event on ties); without a strict +/// majority its label is fake. No valid votes gives an unset label, purity 0. +/// RowLabels maps measurement rows to their original clusters' MC labels. +/// Claims must cover all contributors, before any adapter vertex reordering. +template +[[nodiscard]] std::vector vertexTruth(std::span claims, std::span contributors, const RowLabels& rowLabels) +{ + using Collision = std::pair; // source, event + std::vector> votes(contributors.size()); + std::vector counts(contributors.size(), 0U); + for (const auto& claim : claims) { + if (claim.vertex < 0) { + continue; + } + const auto v = static_cast(claim.vertex); + if (v >= contributors.size()) { + throw std::runtime_error("vertex MC claim refers to a missing vertex"); + } + ++counts[v]; + std::array rows; + for (std::size_t i = 0; i < rows.size(); ++i) { + rows[i] = rowParticles(rowLabels(claim.rows[i])); + } + constexpr std::array order{0U, 1U, 2U}; + if (const auto particle = goodParticle(rows, order); particle != kNoTrackParticle) { + const auto label = labelOfParticle(particle); + ++votes[v][{label.getSourceID(), label.getEventID()}]; + } + } + std::vector result(contributors.size()); + for (std::size_t v = 0; v < result.size(); ++v) { + if (counts[v] != contributors[v]) { + throw std::runtime_error("vertex MC claims do not match its contributor count"); + } + std::uint32_t best = 0U; + for (const auto& [collision, count] : votes[v]) { + if (count > best) { + best = count; + result[v].label = {o2::MCCompLabel::maxTrackID(), collision.second, collision.first, count <= contributors[v] / 2U}; + } + } + result[v].purity = contributors[v] == 0U ? 0.F : static_cast(best) / static_cast(contributors[v]); + } + return result; +} + +} // namespace o2::gtrk::workflow + +#endif // GTRK_WORKFLOW_VERTEXLABEL_H_ diff --git a/Detectors/GenericTracking/workflow/src/AdaptorDevice.h b/Detectors/GenericTracking/workflow/src/AdaptorDevice.h index b0452a5533e2a..12ff2ce1d7f56 100644 --- a/Detectors/GenericTracking/workflow/src/AdaptorDevice.h +++ b/Detectors/GenericTracking/workflow/src/AdaptorDevice.h @@ -23,6 +23,9 @@ #ifndef GTRK_WORKFLOW_ADAPTORDEVICE_H_ #define GTRK_WORKFLOW_ADAPTORDEVICE_H_ +#include "AdaptorHelpers.h" +#include "GenericTracking/Workflow/VertexLabel.h" + #include "GenericTracking/Alice/Clusters.h" #include "GenericTracking/Detectors/ITS.h" #include "GenericTracking/Detectors/MFT.h" @@ -80,101 +83,6 @@ namespace o2::gtrk::workflow namespace adaptor { -/// A particle's label, the fake flag unless \a good; O2's noise label for no -/// particle. -[[nodiscard]] inline o2::MCCompLabel labelOf(std::uint64_t particle, bool good) { return particle == kNoTrackParticle ? o2::MCCompLabel(true) : labelOfParticle(particle, !good); } - -/// A timeframe's clusters as the adaptor needs them: each one's readout frame -/// (its record in CLUSTERSROF) and that frame's start in BC from the -/// timeframe's start, without the delay and bias. -struct ClusterFrames { - std::vector record; ///< per cluster - std::vector recordBC; ///< per record -}; - -[[nodiscard]] inline ClusterFrames clusterFrames(std::span rofs, std::size_t clusters, std::uint32_t tfFirstOrbit) -{ - ClusterFrames result; - result.record.assign(clusters, std::numeric_limits::max()); - result.recordBC.reserve(rofs.size()); - for (std::size_t r = 0; r < rofs.size(); ++r) { - const auto& ir = rofs[r].getBCData(); - result.recordBC.push_back((static_cast(ir.orbit - tfFirstOrbit) * o2::constants::lhc::LHCMaxBunches) + ir.bc); - for (int i = 0; i < rofs[r].getNEntries(); ++i) { - result.record[static_cast(rofs[r].getFirstEntry()) + static_cast(i)] = static_cast(r); - } - } - return result; -} - -/// A track's time as ITS takes it (`TrackerTraits`): the intersection of its -/// clusters' frames `[start + delay + bias, + length)` in BC from the -/// timeframe's start; where they do not intersect -- a track across frames, -/// which only a margin lets the tracker make -- the record's box, the frames -/// widened by the margin. [lower, upper) in BC. -struct BCRange { - std::int64_t lower = 0; - std::int64_t upper = 0; -}; - -template -[[nodiscard]] BCRange trackTime(const Record& record, std::span clusterIndex, const ClusterFrames& frames, const FrameTiming& timing) -{ - BCRange nominal{.lower = std::numeric_limits::min(), .upper = std::numeric_limits::max()}; - for (std::uint8_t i = 0; i < record.count; ++i) { - const auto start = frames.recordBC[frames.record[clusterIndex[record.clusters.index[i]]]] + timing.firstStartBC; - nominal.lower = std::max(nominal.lower, start); - nominal.upper = std::min(nominal.upper, start + timing.lengthBC); - } - if (nominal.lower < nominal.upper) { - return nominal; - } - return {.lower = bcOfPS(record.timePS), .upper = bcOfPS(static_cast(record.timePS) + record.widthPS)}; -} - -/// The earliest frame (its CLUSTERSROF record) among a track's clusters. -template -[[nodiscard]] std::uint32_t earliestRecord(const Record& record, std::span clusterIndex, const ClusterFrames& frames) -{ - std::uint32_t earliest = std::numeric_limits::max(); - for (std::uint8_t i = 0; i < record.count; ++i) { - earliest = std::min(earliest, frames.record[clusterIndex[record.clusters.index[i]]]); - } - return earliest; -} - -/// Each cluster's pixel count, the converter's decoding (`localClusters`). -[[nodiscard]] inline std::vector pixelCounts(std::span compact, std::span patterns, const o2::itsmft::TopologyDictionary& dictionary) -{ - const auto local = alice::localClusters(compact, patterns, dictionary); - std::vector result(local.size()); - std::ranges::transform(local, result.begin(), [](const auto& cluster) { return cluster.nPixels; }); - return result; -} - -/// The track's particle by the shared rule (TrackLabel.h), its rows' -/// clusters' labels, innermost first; \a carried gets, per row, whether it -/// carries that particle. -template -[[nodiscard]] o2::MCCompLabel trackLabel(const Record& record, std::span clusterIndex, const o2::dataformats::MCTruthContainer& labels, std::span carried) -{ - constexpr std::size_t kRows = std::extent_v>; - std::array rows{}; - std::array track{}; - for (std::uint8_t i = 0; i < record.count; ++i) { - const auto all = labels.getLabels(clusterIndex[record.clusters.index[i]]); - rows[i] = rowParticles(std::span{all.data(), all.size()}); - track[i] = i; - } - const std::span order{track.data(), record.count}; - const auto majority = majorityParticle(rows, order); - const bool good = majority.particle != kNoTrackParticle && majority.carried == record.count; - for (std::uint8_t i = 0; i < record.count; ++i) { - carried[i] = majority.particle != kNoTrackParticle && rows[i].carries(majority.particle); - } - return labelOf(majority.particle, good); -} - /// The description of O2's tracker metadata of \a Det, META/TRACKER. template [[nodiscard]] o2::header::DataDescription metaDescription() @@ -200,6 +108,9 @@ void publishConfiguration(o2::framework::ProcessingContext& pc, bool& published) TMap md; md.SetOwnerKeyValue(); md.Add(new TObjString(name.c_str()), new TObjString(o2::conf::ConfigurableParam::asJSON(name).c_str())); + if constexpr (std::same_as) { + md.Add(new TObjString("fastMultConfig"), new TObjString(o2::conf::ConfigurableParam::asJSON("fastMultConfig").c_str())); + } pc.outputs().snapshot(o2::framework::Output{"META", metaDescription(), 0}, md); } @@ -209,10 +120,10 @@ class AdaptorTask : public o2::framework::Task { public: using Record = TrackRecord; - using Labels = o2::dataformats::MCTruthContainer; static constexpr bool kITS = std::same_as; - AdaptorTask(bool useMC, std::shared_ptr> conditions) : mUseMC(useMC), mConditions(std::move(conditions)) {} + AdaptorTask(bool useMC, bool staggered, std::shared_ptr> conditions) + : mUseMC(useMC), mStaggered(staggered), mConditions(std::move(conditions)) {} void init(o2::framework::InitContext& ic) final { mConditions->init(ic); } @@ -237,18 +148,13 @@ class AdaptorTask : public o2::framework::Task const auto orbit = pc.services().get().firstTForbit; const auto records = pc.inputs().get>("tracks"); const auto rows = pc.inputs().get>("clusterIndex"); - const auto rofs = pc.inputs().get>("rofs"); - const auto compact = pc.inputs().get>("compclusters"); - const auto patterns = pc.inputs().get>("patterns"); - std::unique_ptr labels; - if (mUseMC) { - labels = pc.inputs().get("labels"); - } + const auto input = readITSMFTInput(pc, mUseMC, mStaggered ? Detector::nLayers : 1); const std::span clusterIndex{rows.data(), rows.size()}; - const std::span frameRecords{rofs.data(), rofs.size()}; - const auto frames = adaptor::clusterFrames(frameRecords, compact.size(), orbit); - const auto pixels = compact.empty() ? std::vector{} : adaptor::pixelCounts({compact.data(), compact.size()}, {patterns.data(), patterns.size()}, *mDictionary); - const auto timing = mConditions->frameTiming().front(); + const auto frames = adaptor::clusterFrames(input, orbit); + const auto pixels = adaptor::pixelCounts(input, *mDictionary); + const auto timings = mConditions->frameTiming(); + const auto clock = mStaggered ? adaptor::clockLayer(timings) : 0U; + const auto frameRecords = input.streams()[clock].rofs; const double bz = nominalBzKiloGauss(); // The fitted tracks, in the tracker's order. @@ -259,21 +165,22 @@ class AdaptorTask : public o2::framework::Task fitted.push_back(track); } } - const auto layerOf = [&](std::uint32_t row) { return mGeometry->surface(compact[clusterIndex[row]].getSensorID()).layer; }; + const auto layerOf = [&](std::uint32_t row) { return mGeometry->surface(input.cluster(clusterIndex[row]).getSensorID()).layer; }; if constexpr (kITS) { - runITS(pc, records, fitted, clusterIndex, frameRecords, frames, pixels, labels.get(), timing, bz, orbit, layerOf); + runITS(pc, records, fitted, clusterIndex, frameRecords, frames, pixels, input, timings, clock, bz, orbit, layerOf); } else { - runMFT(pc, records, fitted, clusterIndex, frameRecords, frames, pixels, labels.get(), bz, orbit, layerOf); + runMFT(pc, records, fitted, clusterIndex, frameRecords, frames, pixels, input, timings, clock, bz, orbit, layerOf); } } private: template - void runITS(o2::framework::ProcessingContext& pc, gsl::span records, std::span fitted, std::span clusterIndex, std::span rofs, const adaptor::ClusterFrames& frames, std::span pixels, const Labels* labels, const FrameTiming& timing, double bz, std::uint32_t orbit, const LayerOf& layerOf) + void runITS(o2::framework::ProcessingContext& pc, gsl::span records, std::span fitted, std::span clusterIndex, std::span rofs, const adaptor::ClusterFrames& frames, std::span pixels, const ITSMFTInput& input, std::span timings, std::size_t clock, double bz, std::uint32_t orbit, const LayerOf& layerOf) { using o2::framework::Output; using AlpideParam = o2::itsmft::DPLAlpideParam; const auto& origin = detectorOrigin(); + const auto selectedROFs = pc.inputs().get>("selectedROFs"); auto& irFrames = pc.outputs().make>(Output{origin, "IRFRAMES", 0}); auto& clusterIds = pc.outputs().make>(Output{origin, "TRACKCLSID", 0}); auto& tracks = pc.outputs().make>(Output{origin, "TRACKS", 0}); @@ -281,14 +188,14 @@ class AdaptorTask : public o2::framework::Task auto& vertices = pc.outputs().make>(Output{origin, "VERTICES", 0}); auto& vertexROFs = pc.outputs().make>(Output{origin, "VERTICESROF", 0}); - const auto length = static_cast(timing.lengthBC); - const auto delay = static_cast(AlpideParam::Instance().getROFDelayInBC(0)); + const auto length = static_cast(timings[clock].lengthBC); + const auto delay = static_cast(AlpideParam::Instance().getROFDelayInBC(static_cast(clock))); using BC = o2::itsmft::tracking::ROFTimingLayer::BCType; // ITS's frame of a time: its lower edge's, 0 before the first. const o2::itsmft::tracking::ROFTimingLayer rofTiming{ .mROFLength = static_cast(length), .mROFDelay = static_cast(delay), - .mROFBias = static_cast(AlpideParam::Instance().getROFBiasInBC(0)), + .mROFBias = static_cast(AlpideParam::Instance().getROFBiasInBC(static_cast(clock))), }; // A vertex's frame as a track's: the lower edge of its time with the half // width at most half a frame. Its box's own lower edge lies a margin @@ -328,19 +235,19 @@ class AdaptorTask : public o2::framework::Task out.track.getParamOut() = toTrackParCov(record.outwardState, mGeometry->planeOf(record.outwardState), outerAlpha, bz); out.track.setChi2(record.chi2); out.track.setNumberOfClusters(record.count); - if (labels != nullptr) { - out.label = adaptor::trackLabel(record, clusterIndex, *labels, std::span{carried.data(), record.count}); + if (mUseMC) { + out.label = adaptor::trackLabel(record, clusterIndex, input, std::span{carried.data(), record.count}); } std::uint32_t pattern = 0U; for (std::uint8_t i = 0; i < record.count; ++i) { const auto layer = layerOf(record.clusters.index[i]); pattern |= 1U << layer; - pattern |= labels != nullptr && carried[i] ? 1U << (16U + layer) : 0U; + pattern |= mUseMC && carried[i] ? 1U << (16U + layer) : 0U; out.track.setClusterSize(layer, static_cast(pixels[clusterIndex[record.clusters.index[i]]])); } out.track.setPattern(pattern); out.track.setSharedClusters(innermost[clusterIndex[record.clusters.index[0]]] > 1U); - const auto time = adaptor::trackTime(record, clusterIndex, frames, timing); + const auto time = adaptor::trackTime(record, clusterIndex, frames, timings, layerOf); const float half = std::min(static_cast(time.upper - time.lower) / 2.F, static_cast(length) / 2.F); const float centre = static_cast(time.lower + time.upper) / 2.F; out.track.getTimeStamp() = o2::its::TimeStamp(centre, half); @@ -349,19 +256,45 @@ class AdaptorTask : public o2::framework::Task // ITS's order: by the time's lower edge, the tracker's order within. std::ranges::stable_sort(converted, {}, [](const Converted& c) { return c.track.getTimeStamp().getTimeStamp() - c.track.getTimeStamp().getTimeStampError(); }); + // Keep each vertex's MC truth beside it through the output ordering. + struct ConvertedVertex { + o2::its::Vertex vertex; + VertexTruth truth; + }; + std::vector convertedVertices; // The own vertexer's vertices, at the beam line's x and y. const auto shipped = pc.inputs().get>("measurementGeometry"); const auto seeds = pc.inputs().get>("vertices"); const auto contributors = pc.inputs().get>("contributors"); + if (contributors.size() != seeds.size()) { + throw std::runtime_error("vertex seeds and contributor counts differ in size"); + } + std::vector truth; + if (mUseMC) { + const auto claims = pc.inputs().get>("vertexClaims"); + truth = vertexTruth({claims.data(), claims.size()}, {contributors.data(), contributors.size()}, [&](std::uint32_t row) { + if (row >= clusterIndex.size()) { + throw std::runtime_error("vertex MC claim refers to a missing measurement"); + } + const auto labels = input.labels(clusterIndex[row]); + return std::span{labels.data(), labels.size()}; + }); + } + convertedVertices.reserve(seeds.size()); for (std::size_t v = 0; v < seeds.size(); ++v) { const auto& seed = seeds[v]; const std::array cov{seed.varTransverse, 0.F, seed.varTransverse, 0.F, 0.F, seed.varZ}; - auto& vertex = vertices.emplace_back(o2::math_utils::Point3D(static_cast(shipped.beamLine.x), static_cast(shipped.beamLine.y), seed.z), cov, static_cast(std::min(contributors[v], std::numeric_limits::max())), 0.F); + auto& convertedVertex = convertedVertices.emplace_back(); + convertedVertex.truth = mUseMC ? truth[v] : VertexTruth{}; + auto& vertex = convertedVertex.vertex; + vertex = o2::its::Vertex(o2::math_utils::Point3D(static_cast(shipped.beamLine.x), static_cast(shipped.beamLine.y), seed.z), cov, static_cast(std::min(contributors[v], std::numeric_limits::max())), 0.F); const auto lower = bcOfPS(seed.timePS); vertex.setTimeStamp(o2::its::TimeEstBC(static_cast(lower), static_cast(bcOfPS(static_cast(seed.timePS) + seed.widthPS) - lower))); } // By frame, as the frame records count them, then by time and contributors. - std::ranges::stable_sort(vertices, [&](const o2::its::Vertex& a, const o2::its::Vertex& b) { + std::ranges::stable_sort(convertedVertices, [&](const ConvertedVertex& ca, const ConvertedVertex& cb) { + const auto& a = ca.vertex; + const auto& b = cb.vertex; const auto frameA = vertexFrame(a); const auto frameB = vertexFrame(b); if (frameA != frameB) { @@ -370,6 +303,11 @@ class AdaptorTask : public o2::framework::Task return a.getTimeStamp().lower() != b.getTimeStamp().lower() ? a.getTimeStamp().lower() < b.getTimeStamp().lower() : a.getNContributors() > b.getNContributors(); }); + vertices.reserve(convertedVertices.size()); + for (const auto& convertedVertex : convertedVertices) { + vertices.push_back(convertedVertex.vertex); + } + // The clock's frames of the timeframe, as ITS's (`setBCData`), and more // if a track or vertex lies beyond them. std::int64_t frameCount = (o2::constants::lhc::LHCMaxBunches / length) * mConditions->orbitsPerTimeframe(); @@ -379,16 +317,8 @@ class AdaptorTask : public o2::framework::Task for (const auto& vertex : vertices) { frameCount = std::max(frameCount, vertexFrame(vertex) + 1); } - for (auto* frameRecords : {&trackROFs, &vertexROFs}) { - frameRecords->resize(static_cast(frameCount)); - for (std::int64_t i = 0; i < frameCount; ++i) { - auto& rof = (*frameRecords)[static_cast(i)]; - rof.setBCData(o2::InteractionRecord(static_cast(((i * length) % o2::constants::lhc::LHCMaxBunches) + delay), orbit + static_cast(i * length / o2::constants::lhc::LHCMaxBunches))); - rof.setROFrame(static_cast(i)); - rof.setNEntries(0); - rof.setFirstEntry(-1); - } - } + adaptor::fillClockRecords(trackROFs, static_cast(frameCount), static_cast(length), static_cast(delay), orbit); + vertexROFs = trackROFs; // As ITS: the standard vertexing mode on the input's frames, with the vertexer. if (recoParam().vertexer) { for (std::size_t i = 0; i < std::min(rofs.size(), vertexROFs.size()); ++i) { @@ -402,7 +332,7 @@ class AdaptorTask : public o2::framework::Task const auto& record = records[c.record]; c.track.setFirstClusterEntry(static_cast(clusterIds.size())); for (std::uint8_t i = record.count; i-- > 0;) { // outermost first - clusterIds.push_back(static_cast(clusterIndex[record.clusters.index[i]])); + clusterIds.push_back(input.externalIndex(clusterIndex[record.clusters.index[i]])); } tracks.push_back(c.track); ++entries[static_cast(c.frame)]; @@ -412,9 +342,10 @@ class AdaptorTask : public o2::framework::Task trackROFs[i].setFirstEntry(entries[i]); trackROFs[i].setNEntries(entries[i + 1] - entries[i]); trackROFs[i].setFlags(vertexROFs[i].getFlags()); - // Every frame is processed: no multiplicity selection here. - auto& irFrame = irFrames.emplace_back(trackROFs[i].getBCData(), trackROFs[i].getBCData() + (length - 1)); - irFrame.info = static_cast(trackROFs[i].getNEntries()); + if (i < selectedROFs.size() && selectedROFs[i] != 0U) { + auto& irFrame = irFrames.emplace_back(trackROFs[i].getBCData(), trackROFs[i].getBCData() + (length - 1)); + irFrame.info = static_cast(trackROFs[i].getNEntries()); + } } std::ranges::fill(entries, 0); for (const auto& vertex : vertices) { @@ -432,15 +363,20 @@ class AdaptorTask : public o2::framework::Task for (const auto& c : converted) { trackLabels.push_back(c.label); } - // No vertex truth: the products do not say which clusters contributed. - pc.outputs().make>(Output{origin, "VERTICESMCTR", 0}, vertices.size()); - pc.outputs().make>(Output{origin, "VERTICESMCPUR", 0}, vertices.size(), 0.F); + auto& vertexLabels = pc.outputs().make>(Output{origin, "VERTICESMCTR", 0}); + auto& vertexPurities = pc.outputs().make>(Output{origin, "VERTICESMCPUR", 0}); + vertexLabels.reserve(convertedVertices.size()); + vertexPurities.reserve(convertedVertices.size()); + for (const auto& convertedVertex : convertedVertices) { + vertexLabels.push_back(convertedVertex.truth.label); + vertexPurities.push_back(convertedVertex.truth.purity); + } } LOGP(info, "{}: the timeframe at orbit {}: {} tracks ({} failed refits left out, {} in a frame turned by pi), {} vertices, {} readout frames", Det::name, orbit, tracks.size(), records.size() - fitted.size(), turned, vertices.size(), trackROFs.size()); } template - void runMFT(o2::framework::ProcessingContext& pc, gsl::span records, std::span fitted, std::span clusterIndex, std::span rofs, const adaptor::ClusterFrames& frames, std::span pixels, const Labels* labels, double bz, std::uint32_t orbit, const LayerOf& layerOf) + void runMFT(o2::framework::ProcessingContext& pc, gsl::span records, std::span fitted, std::span clusterIndex, std::span rofs, const adaptor::ClusterFrames& frames, std::span pixels, const ITSMFTInput& input, std::span timings, std::size_t clock, double bz, std::uint32_t orbit, const LayerOf& layerOf) { using o2::framework::Output; const auto& origin = detectorOrigin(); @@ -448,6 +384,13 @@ class AdaptorTask : public o2::framework::Task auto& clusterIds = pc.outputs().make>(Output{origin, "TRACKCLSID", 0}); auto& tracks = pc.outputs().make>(Output{origin, "TRACKS", 0}); + using AlpideParam = o2::itsmft::DPLAlpideParam; + using BC = o2::itsmft::tracking::ROFTimingLayer::BCType; + const auto length = timings[clock].lengthBC; + const auto delay = AlpideParam::Instance().getROFDelayInBC(static_cast(clock)); + const o2::itsmft::tracking::ROFTimingLayer rofTiming{ + .mROFLength = static_cast(length), .mROFDelay = static_cast(delay), .mROFBias = static_cast(AlpideParam::Instance().getROFBiasInBC(static_cast(clock)))}; + struct Converted { o2::mft::TrackMFT track; o2::MCCompLabel label; @@ -467,21 +410,37 @@ class AdaptorTask : public o2::framework::Task out.track.setOutParam(toTrackParCovFwd(record.outwardState, mGeometry->planeOf(record.outwardState), bz, record.outwardChi2)); out.track.setCA(true); out.track.setNumberOfPoints(record.count); - if (labels != nullptr) { - out.label = adaptor::trackLabel(record, clusterIndex, *labels, std::span{carried.data(), record.count}); + if (mUseMC) { + out.label = adaptor::trackLabel(record, clusterIndex, input, std::span{carried.data(), record.count}); } for (std::uint8_t i = 0; i < record.count; ++i) { out.track.setClusterSize(layerOf(record.clusters.index[i]), static_cast(pixels[clusterIndex[record.clusters.index[i]]])); } - // TrackMFT has no time: a track across frames is timed by the earlier. - out.frame = adaptor::earliestRecord(record, clusterIndex, frames); - for (std::uint8_t i = 0; i < record.count; ++i) { - if (frames.record[clusterIndex[record.clusters.index[i]]] != out.frame) { - ++crossing; - break; + // Preserve the monolithic MFT's earliest-record convention. Per-layer + // record numbers have different clocks: use the interval's lower edge + // on the fastest layer's clock instead. + if (mStaggered) { + const auto time = adaptor::trackTime(record, clusterIndex, frames, timings, layerOf); + const float half = std::min(static_cast(time.upper - time.lower) / 2.F, static_cast(length) / 2.F); + const float centre = static_cast(time.lower + time.upper) / 2.F; + out.frame = rofTiming.getROF(o2::its::TimeStamp(centre, half)); + } else { + out.frame = adaptor::earliestRecord(record, clusterIndex, frames); + for (std::uint8_t i = 0; i < record.count; ++i) { + if (frames.record[clusterIndex[record.clusters.index[i]]] != out.frame) { + ++crossing; + break; + } } } } + if (mStaggered) { + std::size_t count = static_cast(o2::constants::lhc::LHCMaxBunches / length) * mConditions->orbitsPerTimeframe(); + for (const auto& track : converted) { + count = std::max(count, static_cast(track.frame) + 1U); + } + adaptor::fillClockRecords(trackROFs, count, length, delay, orbit); + } std::ranges::stable_sort(converted, {}, &Converted::frame); std::vector entries(trackROFs.size() + 1U, 0); @@ -490,7 +449,7 @@ class AdaptorTask : public o2::framework::Task const auto& record = records[c.record]; c.track.setExternalClusterIndexOffset(static_cast(clusterIds.size())); for (std::uint8_t i = 0; i < record.count; ++i) { // innermost first - clusterIds.push_back(static_cast(clusterIndex[record.clusters.index[i]])); + clusterIds.push_back(input.externalIndex(clusterIndex[record.clusters.index[i]])); } tracks.push_back(c.track); ++entries[c.frame]; @@ -507,10 +466,12 @@ class AdaptorTask : public o2::framework::Task trackLabels.push_back(c.label); } } - LOGP(info, "{}: the timeframe at orbit {}: {} tracks ({} failed refits left out, {} across readout frames, timed by the earlier), {} readout frames", Det::name, orbit, tracks.size(), records.size() - fitted.size(), crossing, trackROFs.size()); + LOGP(info, "{}: the timeframe at orbit {}: {} tracks ({} failed refits left out), {} readout frames; {}", Det::name, orbit, tracks.size(), records.size() - fitted.size(), trackROFs.size(), + mStaggered ? std::format("clock layer {}", clock) : std::format("{} across readout frames, timed by the earlier", crossing)); } bool mUseMC = false; + bool mStaggered = false; bool mConfigurationPublished = false; std::shared_ptr> mConditions; ///< the run's geometry, field, frames and orbits per timeframe std::unique_ptr> mGeometry; ///< the planes the tracker's states are bound to @@ -520,16 +481,21 @@ class AdaptorTask : public o2::framework::Task } // namespace adaptor template -o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC) +o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC, bool staggered) { using namespace o2::framework; constexpr bool kITS = std::same_as; const auto origin = detectorOrigin(); - std::vector inputs{InputSpec{"tracks", gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, InputSpec{"clusterIndex", gDataOrigin, description(kMeasurements), static_cast(MeasurementColumn::clusterIndex), Lifetime::Timeframe}, InputSpec{"rofs", origin, "CLUSTERSROF", 0, Lifetime::Timeframe}, InputSpec{"compclusters", origin, "COMPCLUSTERS", 0, Lifetime::Timeframe}, InputSpec{"patterns", origin, "PATTERNS", 0, Lifetime::Timeframe}, InputSpec{"dictionary", origin, "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec(origin.template as() + "/Calib/ClusterDictionary")}}; + std::vector inputs{ + InputSpec{"tracks", gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, + InputSpec{"clusterIndex", gDataOrigin, description(kMeasurements), static_cast(MeasurementColumn::clusterIndex), Lifetime::Timeframe}, + InputSpec{"dictionary", origin, "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec(origin.template as() + "/Calib/ClusterDictionary")}}; + addITSMFTClusterInputs(inputs, origin, useMC, staggered ? Detector::nLayers : 1); std::vector outputs{{origin, "TRACKS", 0, Lifetime::Timeframe}, {origin, "TRACKCLSID", 0, Lifetime::Timeframe}, {"META", adaptor::metaDescription(), 0, Lifetime::Sporadic}}; if constexpr (kITS) { inputs.emplace_back("measurementGeometry", gDataOrigin, description(kMeasurementGeometry), 0, Lifetime::Timeframe); inputs.emplace_back("vertices", gDataOrigin, description(kVertexSeeds), 0, Lifetime::Timeframe); + inputs.emplace_back("selectedROFs", gDataOrigin, description(kSelectedROFs), 0, Lifetime::Timeframe); inputs.emplace_back("contributors", gDataOrigin, description(kVertexContributors), 0, Lifetime::Timeframe); outputs.emplace_back(origin, "ITSTrackROF", 0, Lifetime::Timeframe); outputs.emplace_back(origin, "VERTICES", 0, Lifetime::Timeframe); @@ -539,9 +505,9 @@ o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC) outputs.emplace_back(origin, "MFTTrackROF", 0, Lifetime::Timeframe); } if (useMC) { - inputs.emplace_back("labels", origin, "CLUSTERSMCTR", 0, Lifetime::Timeframe); outputs.emplace_back(origin, "TRACKSMCTR", 0, Lifetime::Timeframe); if constexpr (kITS) { + inputs.emplace_back("vertexClaims", gDataOrigin, description(kVertexClaims), 0, Lifetime::Timeframe); outputs.emplace_back(origin, "VERTICESMCTR", 0, Lifetime::Timeframe); outputs.emplace_back(origin, "VERTICESMCPUR", 0, Lifetime::Timeframe); } @@ -551,7 +517,7 @@ o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC) .name = deviceName("adaptor"), .inputs = inputs, .outputs = outputs, - .algorithm = AlgorithmSpec{adaptFromTask>(useMC, conditions)}, + .algorithm = AlgorithmSpec{adaptFromTask>(useMC, staggered, conditions)}, }; } diff --git a/Detectors/GenericTracking/workflow/src/AdaptorHelpers.h b/Detectors/GenericTracking/workflow/src/AdaptorHelpers.h new file mode 100644 index 0000000000000..136fb86eb1e81 --- /dev/null +++ b/Detectors/GenericTracking/workflow/src/AdaptorHelpers.h @@ -0,0 +1,199 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef GTRK_WORKFLOW_ADAPTORHELPERS_H_ +#define GTRK_WORKFLOW_ADAPTORHELPERS_H_ + +#include "ITSMFTInput.h" +#include "GenericTracking/Alice/Clusters.h" +#include "GenericTracking/Base/FrameTiming.h" +#include "GenericTracking/Workflow/TrackLabel.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace o2::gtrk::workflow::adaptor +{ + +/// A particle's label, the fake flag unless \a good; O2's noise label for no +/// particle. +[[nodiscard]] inline o2::MCCompLabel labelOf(std::uint64_t particle, bool good) { return particle == kNoTrackParticle ? o2::MCCompLabel(true) : labelOfParticle(particle, !good); } + +/// A timeframe's clusters as the adaptor needs them: each one's readout frame +/// (its record in CLUSTERSROF) and that frame's start in BC from the +/// timeframe's start, without the delay and bias. +struct ClusterFrames { + std::vector record; ///< per cluster + std::vector recordBC; ///< per record +}; + +[[nodiscard]] inline ClusterFrames clusterFrames(std::span rofs, std::size_t clusters, std::uint32_t tfFirstOrbit) +{ + ClusterFrames result; + result.record.assign(clusters, std::numeric_limits::max()); + result.recordBC.reserve(rofs.size()); + std::size_t next = 0U; + for (std::size_t r = 0; r < rofs.size(); ++r) { + if (rofs[r].getFirstEntry() < 0 || static_cast(rofs[r].getFirstEntry()) != next || + rofs[r].getNEntries() < 0 || static_cast(rofs[r].getNEntries()) > clusters - next) { + throw std::runtime_error("cluster ROFs do not cover the input in order"); + } + next += static_cast(rofs[r].getNEntries()); + const auto& ir = rofs[r].getBCData(); + result.recordBC.push_back(((static_cast(ir.orbit) - tfFirstOrbit) * o2::constants::lhc::LHCMaxBunches) + ir.bc); + for (int i = 0; i < rofs[r].getNEntries(); ++i) { + result.record[static_cast(rofs[r].getFirstEntry()) + static_cast(i)] = static_cast(r); + } + } + if (next != clusters) { + throw std::runtime_error("cluster ROFs do not cover the input"); + } + return result; +} + +/// A track's time as ITS takes it (`TrackerTraits`): the intersection of its +/// clusters' frames `[start + delay + bias, + length)` in BC from the +/// timeframe's start; where they do not intersect -- a track across frames, +/// which only a margin lets the tracker make -- the record's box, the frames +/// widened by the margin. [lower, upper) in BC. +struct BCRange { + std::int64_t lower = 0; + std::int64_t upper = 0; +}; + +template +[[nodiscard]] BCRange trackTime(const Record& record, std::span clusterIndex, const ClusterFrames& frames, std::span timings, const LayerOf& layerOf) +{ + BCRange nominal{.lower = std::numeric_limits::min(), .upper = std::numeric_limits::max()}; + for (std::uint8_t i = 0; i < record.count; ++i) { + const auto row = record.clusters.index[i]; + const auto& timing = timings[layerOf(row)]; + const auto start = frames.recordBC[frames.record[clusterIndex[row]]] + timing.firstStartBC; + nominal.lower = std::max(nominal.lower, start); + nominal.upper = std::min(nominal.upper, start + timing.lengthBC); + } + if (nominal.lower < nominal.upper) { + return nominal; + } + return {.lower = bcOfPS(record.timePS), .upper = bcOfPS(static_cast(record.timePS) + record.widthPS)}; +} + +/// The earliest frame (its CLUSTERSROF record) among a track's clusters. +template +[[nodiscard]] std::uint32_t earliestRecord(const Record& record, std::span clusterIndex, const ClusterFrames& frames) +{ + std::uint32_t earliest = std::numeric_limits::max(); + for (std::uint8_t i = 0; i < record.count; ++i) { + earliest = std::min(earliest, frames.record[clusterIndex[record.clusters.index[i]]]); + } + return earliest; +} + +/// Each cluster's pixel count, the converter's decoding (`localClusters`). +[[nodiscard]] inline std::vector pixelCounts(std::span compact, std::span patterns, const o2::itsmft::TopologyDictionary& dictionary) +{ + const auto local = alice::localClusters(compact, patterns, dictionary); + std::vector result(local.size()); + std::ranges::transform(local, result.begin(), [](const auto& cluster) { return cluster.nPixels; }); + return result; +} + +/// The track's particle by the shared rule (TrackLabel.h), its rows' +/// clusters' labels, innermost first; \a carried gets, per row, whether it +/// carries that particle. +template +[[nodiscard]] o2::MCCompLabel trackLabel(const Record& record, std::span clusterIndex, const ITSMFTInput& input, std::span carried) +{ + constexpr std::size_t kRows = std::extent_v>; + std::array rows{}; + std::array track{}; + for (std::uint8_t i = 0; i < record.count; ++i) { + const auto all = input.labels(clusterIndex[record.clusters.index[i]]); + rows[i] = rowParticles(std::span{all.data(), all.size()}); + track[i] = i; + } + const std::span order{track.data(), record.count}; + const auto majority = majorityParticle(rows, order); + const bool good = majority.particle != kNoTrackParticle && majority.carried == record.count; + for (std::uint8_t i = 0; i < record.count; ++i) { + carried[i] = majority.particle != kNoTrackParticle && rows[i].carries(majority.particle); + } + return labelOf(majority.particle, good); +} + +/// Flatten ROF identities in the same layer order as the converter's clusters. +[[nodiscard]] inline ClusterFrames clusterFrames(const ITSMFTInput& input, std::uint32_t orbit) +{ + ClusterFrames result; + result.record.reserve(input.size()); + for (const auto& stream : input.streams()) { + auto frames = clusterFrames(stream.rofs, stream.compact.size(), orbit); + const auto offset = static_cast(result.recordBC.size()); + for (const auto record : frames.record) { + result.record.push_back(offset + record); + } + result.recordBC.insert(result.recordBC.end(), frames.recordBC.begin(), frames.recordBC.end()); + } + return result; +} + +[[nodiscard]] inline std::vector pixelCounts(const ITSMFTInput& input, const o2::itsmft::TopologyDictionary& dictionary) +{ + std::vector result; + result.reserve(input.size()); + for (const auto& stream : input.streams()) { + if (!stream.compact.empty()) { + const auto pixels = pixelCounts(stream.compact, stream.patterns, dictionary); + result.insert(result.end(), pixels.begin(), pixels.end()); + } + } + return result; +} + +/// O2's output clock is the layer with the most frames per timeframe. Equal +/// lengths retain the first layer, independently of the input's occupancy. +[[nodiscard]] inline std::size_t clockLayer(std::span timings) +{ + if (timings.empty() || std::ranges::any_of(timings, [](const auto& timing) { return timing.lengthBC <= 0; })) { + throw std::invalid_argument("output clock needs positive per-layer frame lengths"); + } + return static_cast(std::ranges::min_element(timings, {}, &FrameTiming::lengthBC) - timings.begin()); +} + +/// Frame records carry the delay, whereas frame intervals also carry the bias. +/// Normalise a delayed record across an orbit boundary. +template +inline void fillClockRecords(Records& result, std::size_t count, int lengthBC, int delayBC, std::uint32_t orbit) +{ + if (lengthBC <= 0 || delayBC < 0) { + throw std::invalid_argument("invalid output clock record"); + } + result.resize(count); + for (std::size_t i = 0; i < count; ++i) { + const auto bc = static_cast(i) * lengthBC + delayBC; + auto& rof = result[i]; + rof.setBCData(o2::InteractionRecord(static_cast(bc % o2::constants::lhc::LHCMaxBunches), + orbit + static_cast(bc / o2::constants::lhc::LHCMaxBunches))); + rof.setROFrame(static_cast(i)); + rof.setNEntries(0); + rof.setFirstEntry(-1); + } +} + +} // namespace o2::gtrk::workflow::adaptor + +#endif // GTRK_WORKFLOW_ADAPTORHELPERS_H_ diff --git a/Detectors/GenericTracking/workflow/src/AdaptorITS.cxx b/Detectors/GenericTracking/workflow/src/AdaptorITS.cxx index 4dd2b5194eace..6db81098acd95 100644 --- a/Detectors/GenericTracking/workflow/src/AdaptorITS.cxx +++ b/Detectors/GenericTracking/workflow/src/AdaptorITS.cxx @@ -14,6 +14,6 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getAdaptorSpec(bool useMC); +template o2::framework::DataProcessorSpec workflow::getAdaptorSpec(bool useMC, bool staggered); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/AdaptorMFT.cxx b/Detectors/GenericTracking/workflow/src/AdaptorMFT.cxx index 658d1bd9ae11e..29bc019b55c55 100644 --- a/Detectors/GenericTracking/workflow/src/AdaptorMFT.cxx +++ b/Detectors/GenericTracking/workflow/src/AdaptorMFT.cxx @@ -14,6 +14,6 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getAdaptorSpec(bool useMC); +template o2::framework::DataProcessorSpec workflow::getAdaptorSpec(bool useMC, bool staggered); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/ClusterConverterITSMFT.h b/Detectors/GenericTracking/workflow/src/ClusterConverterITSMFT.h index a9604105ec6a9..2354bcee31aa0 100644 --- a/Detectors/GenericTracking/workflow/src/ClusterConverterITSMFT.h +++ b/Detectors/GenericTracking/workflow/src/ClusterConverterITSMFT.h @@ -16,6 +16,8 @@ #ifndef GTRK_WORKFLOW_CLUSTERCONVERTERITSMFT_H_ #define GTRK_WORKFLOW_CLUSTERCONVERTERITSMFT_H_ +#include "ITSMFTInput.h" + #include "GenericTracking/Alice/Clusters.h" #include "GenericTracking/Clustering/Cluster.h" #include "GenericTracking/Workflow/Conditions.h" @@ -46,15 +48,18 @@ namespace o2::gtrk::workflow /// O2's compact clusters of the ITS or the MFT, as the O2 cluster reader /// publishes them, turned into GTRK/CLUSTERS for the measurement producer /// and GTRK/CLUSMCTR for the MC truth. The order is O2's, so a -/// measurement's cluster index is O2's cluster index. +/// measurement's cluster index addresses streams concatenated in layer order. template class ClusterConverterTask : public o2::framework::Task { public: - explicit ClusterConverterTask(bool useMC) : mUseMC(useMC) {} + ClusterConverterTask(bool useMC, bool staggered) : mUseMC(useMC), mStreams(staggered ? Detector::nLayers : 1) {} void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final { + if (matcher == o2::framework::ConcreteDataMatcher(Conditions::origin(), "ALPIDEPARAM", 0)) { + return; // the CCDB fetch updates the singleton + } if (matcher == o2::framework::ConcreteDataMatcher(Conditions::origin(), "CLUSDICT", 0)) { mDictionary = static_cast(obj); LOGP(info, "{} topology dictionary from the CCDB: {} patterns", Det::name, mDictionary->getSize()); @@ -64,52 +69,61 @@ class ClusterConverterTask : public o2::framework::Task void run(o2::framework::ProcessingContext& pc) final { pc.inputs().get("dictionary"); // delivers a new one through finaliseCCDB - const auto compact = pc.inputs().get>("compclusters"); - const auto patterns = pc.inputs().get>("patterns"); - const auto rofs = pc.inputs().get>("rofs"); - const auto& timing = pc.services().get(); - - std::vector clusters; - if (!rofs.empty()) { - const auto local = alice::localClusters(std::span{compact.data(), compact.size()}, std::span{patterns.data(), patterns.size()}, *mDictionary); - const std::span frames{rofs.data(), rofs.size()}; - clusters = alice::toClusters(local, frames, alice::readoutFrames(frames), timing.firstTForbit); + using AlpideParam = o2::itsmft::DPLAlpideParam::value>; + pc.inputs().get("alpideParam"); + const auto input = readITSMFTInput(pc, mUseMC, mStreams); + const auto& param = AlpideParam::Instance(); + const auto orbit = pc.services().get().firstTForbit; + + std::vector timings; + timings.reserve(Detector::nLayers); + for (int layer = 0; layer < Detector::nLayers; ++layer) { + timings.push_back({.lengthBC = param.getROFLengthInBC(layer), .firstStartBC = param.getROFDelayInBC(layer) + param.getROFBiasInBC(layer)}); + if (mStreams == 1 && timings.back() != timings.front()) { + throw std::runtime_error("staggered ITS/MFT timing needs per-layer cluster input (--enable-its-staggering or --enable-mft-staggering)"); + } } - LOGP(info, "{}: {} clusters in {} readout frames for the timeframe at orbit {}", Det::name, clusters.size(), rofs.size(), timing.firstTForbit); + const auto clusters = convertITSMFTInput(input, *mDictionary, timings, orbit); + std::size_t rofCount = 0U; + for (const auto& stream : input.streams()) { + rofCount += stream.rofs.size(); + } + LOGP(info, "{}: {} clusters in {} readout frames from {} streams for the timeframe at orbit {}", Det::name, clusters.size(), rofCount, mStreams, orbit); pc.outputs().snapshot(o2::framework::Output{gDataOrigin, description(kClusters), 0}, clusters); if (mUseMC) { - const auto labels = pc.inputs().get*>("labels"); - if (labels->getIndexedSize() != clusters.size() && !clusters.empty()) { - LOGP(fatal, "{}: {} labelled clusters for {} clusters", Det::name, labels->getIndexedSize(), clusters.size()); - } auto& flat = pc.outputs().make>(o2::framework::Output{gDataOrigin, description(kClustersMCTruth), 0}); - labels->flatten_to(flat); + concatenateLabels(input).flatten_to(flat); } } private: bool mUseMC; + int mStreams; const o2::itsmft::TopologyDictionary* mDictionary = nullptr; ///< the CCDB cache's, from the run's ITS/MFT Calib/ClusterDictionary }; /// The converter's device for the ITS or the MFT. template -o2::framework::DataProcessorSpec getITSMFTClusterConverterSpec(bool useMC) +o2::framework::DataProcessorSpec getITSMFTClusterConverterSpec(bool useMC, bool staggered) { using namespace o2::framework; const auto origin = Conditions::origin(); - std::vector inputs{InputSpec{"compclusters", origin, "COMPCLUSTERS", 0, Lifetime::Timeframe}, InputSpec{"patterns", origin, "PATTERNS", 0, Lifetime::Timeframe}, InputSpec{"rofs", origin, "CLUSTERSROF", 0, Lifetime::Timeframe}, InputSpec{"dictionary", origin, "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec(origin.template as() + "/Calib/ClusterDictionary")}}; + const int streams = staggered ? Detector::nLayers : 1; + std::vector inputs; + addITSMFTClusterInputs(inputs, origin, useMC, streams); + const auto detector = origin.template as(); + inputs.emplace_back("dictionary", origin, "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec(detector + "/Calib/ClusterDictionary")); + inputs.emplace_back("alpideParam", origin, "ALPIDEPARAM", 0, Lifetime::Condition, ccdbParamSpec(detector + "/Config/AlpideParam")); std::vector outputs{OutputSpec{gDataOrigin, description(kClusters), 0, Lifetime::Timeframe}}; if (useMC) { - inputs.emplace_back("labels", origin, "CLUSTERSMCTR", 0, Lifetime::Timeframe); outputs.emplace_back(gDataOrigin, description(kClustersMCTruth), 0, Lifetime::Timeframe); } return DataProcessorSpec{ .name = deviceName("cluster-converter"), .inputs = inputs, .outputs = outputs, - .algorithm = AlgorithmSpec{adaptFromTask>(useMC)}, + .algorithm = AlgorithmSpec{adaptFromTask>(useMC, staggered)}, .options = Options{}, }; } diff --git a/Detectors/GenericTracking/workflow/src/GPURecoTracking.cxx b/Detectors/GenericTracking/workflow/src/GPURecoTracking.cxx index d52520edb214e..a3a619f20e273 100644 --- a/Detectors/GenericTracking/workflow/src/GPURecoTracking.cxx +++ b/Detectors/GenericTracking/workflow/src/GPURecoTracking.cxx @@ -23,13 +23,13 @@ namespace o2::gtrk::workflow { -std::unique_ptr makeGPURecoTracking(std::string_view detector) +std::unique_ptr makeGPURecoTracking(std::string_view detector, bool useMC) { if (detector == "its") { - return makeGPURecoTrackingOf(); + return makeGPURecoTrackingOf(useMC); } if (detector == "mft") { - return makeGPURecoTrackingOf(); + return makeGPURecoTrackingOf(useMC); } throw std::invalid_argument("GenericTracking: no tracking inside the GPU reconstruction for " + std::string{detector}); } diff --git a/Detectors/GenericTracking/workflow/src/ITS.cxx b/Detectors/GenericTracking/workflow/src/ITS.cxx index 46aad4dc52184..9032eb7887b35 100644 --- a/Detectors/GenericTracking/workflow/src/ITS.cxx +++ b/Detectors/GenericTracking/workflow/src/ITS.cxx @@ -36,9 +36,9 @@ std::vector alice::placeChips() } template <> -o2::framework::DataProcessorSpec workflow::getClusterConverterSpec(bool useMC) +o2::framework::DataProcessorSpec workflow::getClusterConverterSpec(bool useMC, bool staggered) { - return workflow::getITSMFTClusterConverterSpec(useMC); + return workflow::getITSMFTClusterConverterSpec(useMC, staggered); } template o2::framework::DataProcessorSpec workflow::getMeasurementProducerSpec(); diff --git a/Detectors/GenericTracking/workflow/src/ITSMFTInput.h b/Detectors/GenericTracking/workflow/src/ITSMFTInput.h new file mode 100644 index 0000000000000..0e5db921d16b1 --- /dev/null +++ b/Detectors/GenericTracking/workflow/src/ITSMFTInput.h @@ -0,0 +1,181 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef GTRK_WORKFLOW_ITSMFTINPUT_H_ +#define GTRK_WORKFLOW_ITSMFTINPUT_H_ + +#include "GenericTracking/Alice/Clusters.h" +#include "GenericTracking/Base/FrameTiming.h" +#include "DataFormatsITSMFT/ClusterID.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "Framework/InputSpec.h" +#include "Framework/InputRecord.h" +#include "Framework/ProcessingContext.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2::gtrk::workflow +{ + +inline std::string clusterInputBinding(const char* name, int stream, int streams) +{ + return streams == 1 ? name : std::format("{}_{}", name, stream); +} + +/// The same subscription and ordering for the converter and the output adaptor. +inline void addITSMFTClusterInputs(std::vector& inputs, o2::header::DataOrigin origin, bool useMC, int streams) +{ + using o2::framework::Lifetime; + for (int stream = 0; stream < streams; ++stream) { + inputs.emplace_back(clusterInputBinding("compclusters", stream, streams), origin, "COMPCLUSTERS", stream, Lifetime::Timeframe); + inputs.emplace_back(clusterInputBinding("patterns", stream, streams), origin, "PATTERNS", stream, Lifetime::Timeframe); + inputs.emplace_back(clusterInputBinding("rofs", stream, streams), origin, "CLUSTERSROF", stream, Lifetime::Timeframe); + if (useMC) { + inputs.emplace_back(clusterInputBinding("labels", stream, streams), origin, "CLUSTERSMCTR", stream, Lifetime::Timeframe); + } + } +} + +/// GTRK indices address the streams concatenated in layer order. O2 references +/// retain the input stream and its local index, including when a stream is empty. +class ITSMFTInput +{ + public: + using Labels = o2::dataformats::MCTruthContainer; + struct Stream { + std::span compact; + std::span patterns; + std::span rofs; + std::unique_ptr labels; + }; + struct Source { + std::size_t stream; + std::uint32_t index; + }; + + explicit ITSMFTInput(std::vector streams) : mStreams(std::move(streams)) + { + if (mStreams.empty() || mStreams.size() > o2::itsmft::MaxClusLayers) { + throw std::invalid_argument("ITS/MFT cluster input: invalid stream count"); + } + for (const auto& stream : mStreams) { + if (stream.compact.size() > static_cast(o2::itsmft::ClusIndexMask) + 1U || + stream.compact.size() > std::numeric_limits::max() - mOffsets.back()) { + throw std::invalid_argument("ITS/MFT cluster input: too many clusters"); + } + if (stream.labels && stream.labels->getIndexedSize() != stream.compact.size() && !stream.compact.empty()) { + throw std::invalid_argument("ITS/MFT cluster input: labels do not cover the clusters"); + } + mOffsets.push_back(mOffsets.back() + static_cast(stream.compact.size())); + } + } + + [[nodiscard]] std::span streams() const { return mStreams; } + [[nodiscard]] std::uint32_t size() const { return mOffsets.back(); } + [[nodiscard]] std::uint32_t offset(std::size_t stream) const { return mOffsets.at(stream); } + + [[nodiscard]] Source source(std::uint32_t index) const + { + if (index >= size()) { + throw std::out_of_range("ITS/MFT cluster input: cluster index out of range"); + } + const auto stream = static_cast(std::upper_bound(mOffsets.begin(), mOffsets.end(), index) - mOffsets.begin() - 1); + return {stream, index - mOffsets[stream]}; + } + + [[nodiscard]] int externalIndex(std::uint32_t index) const + { + const auto from = source(index); + return o2::itsmft::composeClusID(static_cast(from.stream), static_cast(from.index)); + } + + [[nodiscard]] const o2::itsmft::CompClusterExt& cluster(std::uint32_t index) const + { + const auto from = source(index); + return mStreams[from.stream].compact[from.index]; + } + + [[nodiscard]] auto labels(std::uint32_t index) const + { + const auto from = source(index); + return mStreams[from.stream].labels->getLabels(from.index); + } + + private: + std::vector mStreams; + std::vector mOffsets{0U}; +}; + +/// Decode streams separately, then concatenate clusters and their labels in +/// the same order. The run's frame lengths also cover single-ROF inputs. +inline std::vector convertITSMFTInput(const ITSMFTInput& input, const o2::itsmft::TopologyDictionary& dictionary, + std::span timings, std::uint32_t orbit) +{ + if (timings.size() < input.streams().size()) { + throw std::invalid_argument("ITS/MFT input: missing per-layer timing"); + } + std::vector result; + result.reserve(input.size()); + for (std::size_t stream = 0; stream < input.streams().size(); ++stream) { + const auto& source = input.streams()[stream]; + if (source.rofs.size() > 1U && alice::readoutFrames(source.rofs).lengthBC != timings[stream].lengthBC) { + throw std::runtime_error(std::format("ITS/MFT stream {}: ROF spacing differs from the run's frame length {} BC", stream, timings[stream].lengthBC)); + } + const auto local = alice::localClusters(source.compact, source.patterns, dictionary); + const auto clusters = alice::toClusters(local, source.rofs, {.lengthBC = timings[stream].lengthBC}, orbit); + result.insert(result.end(), clusters.begin(), clusters.end()); + } + return result; +} + +inline ITSMFTInput::Labels concatenateLabels(const ITSMFTInput& input) +{ + ITSMFTInput::Labels labels; + for (const auto& stream : input.streams()) { + if (!stream.compact.empty()) { + labels.mergeAtBack(*stream.labels); + } + } + return labels; +} + +inline ITSMFTInput readITSMFTInput(o2::framework::ProcessingContext& pc, bool useMC, int streams) +{ + std::vector result; + result.reserve(streams); + for (int stream = 0; stream < streams; ++stream) { + const auto compact = pc.inputs().get>(clusterInputBinding("compclusters", stream, streams).c_str()); + const auto patterns = pc.inputs().get>(clusterInputBinding("patterns", stream, streams).c_str()); + const auto rofs = pc.inputs().get>(clusterInputBinding("rofs", stream, streams).c_str()); + result.push_back({.compact = {compact.data(), compact.size()}, .patterns = {patterns.data(), patterns.size()}, .rofs = {rofs.data(), rofs.size()}, + .labels = useMC ? pc.inputs().get(clusterInputBinding("labels", stream, streams).c_str()) : nullptr}); + } + return ITSMFTInput{std::move(result)}; +} + +} // namespace o2::gtrk::workflow + +#endif // GTRK_WORKFLOW_ITSMFTINPUT_H_ diff --git a/Detectors/GenericTracking/workflow/src/ITSROFSelection.cxx b/Detectors/GenericTracking/workflow/src/ITSROFSelection.cxx new file mode 100644 index 0000000000000..6db18eb4f2f21 --- /dev/null +++ b/Detectors/GenericTracking/workflow/src/ITSROFSelection.cxx @@ -0,0 +1,247 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "GenericTracking/Workflow/ITSROFSelection.h" +#include "GenericTracking/Workflow/RecoParam.h" +#include "DataFormatsITSMFT/DPLAlpideParam.h" +#include "DataFormatsITSMFT/PhysTrigger.h" +#include "DataFormatsTRD/TriggerRecord.h" +#include "Framework/DataRefUtils.h" +#include "Framework/InputRecord.h" +#include "Framework/InputRecordWalker.h" +#include "Framework/ProcessingContext.h" +#include "Framework/TimingInfo.h" +#include "ITStracking/FastMultEst.h" + +#include +#include +#include +#include + +namespace o2::gtrk::workflow +{ + +namespace +{ +using Estimator = o2::its::FastMultEst; +using Overlap = Estimator::ROFOverlapTableN; +using Mask = Estimator::ROFMaskTableN; +using Timing = o2::itsmft::tracking::ROFTimingLayer; +using Alpide = o2::itsmft::DPLAlpideParam; +constexpr int kLayers = Estimator::NLayers; + +int triggerType() +{ + const auto type = recoParam().selectWithTriggers; + if (type < 0 || type > 2) { + throw std::invalid_argument("GTRKRecoITS.selectWithTriggers must be 0 (none), 1 (ITS physics) or 2 (TRD)"); + } + return type; +} +} // namespace + +struct ITSROFSelection::State { + struct ClusterROF { + int layer; + int rof; + }; + Overlap overlap; + Mask mask; + std::vector clusterROFs; +}; + +void ITSROFSelection::addInputs(std::vector& inputs) +{ + using namespace o2::framework; + inputs.emplace_back("its_selection_clusters", ConcreteDataTypeMatcher{"ITS", "COMPCLUSTERS"}, Lifetime::Timeframe); + inputs.emplace_back("its_selection_rofs", ConcreteDataTypeMatcher{"ITS", "CLUSTERSROF"}, Lifetime::Timeframe); + if (triggerType() == 1) { + inputs.emplace_back("its_selection_triggers", "ITS", "PHYSTRIG", 0, Lifetime::Timeframe); + } else if (triggerType() == 2) { + inputs.emplace_back("its_selection_triggers", "TRD", "TRKTRGRD", 0, Lifetime::Timeframe); + } + if (o2::its::FastMultEstConfig::Instance().isVtxMultCutRequested() && !recoParam().vertexer) { + throw std::invalid_argument("fastMultConfig vertex cuts require GTRKRecoITS.vertexer=true (upstream VertexSeed inputs have no contributor counts)"); + } +} + +ITSROFSelection::ITSROFSelection(o2::framework::ProcessingContext& pc, std::span timing, int orbitsPerTF) + : mState(std::make_unique()) +{ + using namespace o2::framework; + std::array, kLayers> clusters{}; + std::array, kLayers> rofs{}; + std::array seen{}; + for (const auto& ref : InputRecordWalker{pc.inputs()}) { + const auto* header = DataRefUtils::getHeader(ref); + if (header->dataOrigin != o2::header::DataOrigin{"ITS"}) { + continue; + } + const auto layer = header->subSpecification; + if (header->dataDescription == o2::header::DataDescription{"COMPCLUSTERS"}) { + if (layer >= kLayers) { + throw std::runtime_error("ITS selection: invalid cluster stream"); + } + clusters[layer] = pc.inputs().get>(ref); + } else if (header->dataDescription == o2::header::DataDescription{"CLUSTERSROF"}) { + if (layer >= kLayers) { + throw std::runtime_error("ITS selection: invalid ROF stream"); + } + rofs[layer] = pc.inputs().get>(ref); + seen[layer] = true; + } + } + const auto streams = std::ranges::count(seen, true); + if (streams != 1 && streams != kLayers) { + throw std::runtime_error("ITS selection requires one monolithic stream or seven layer streams"); + } + const bool staggered = streams == kLayers; + if (timing.size() != kLayers || !seen[0]) { + throw std::runtime_error("ITS selection: missing layer timing or stream zero"); + } + auto& state = *mState; + const auto& alpide = Alpide::Instance(); + const auto& config = recoParam(); + for (int layer = 0; layer < kLayers; ++layer) { + const int nROFs = std::max(orbitsPerTF * (o2::constants::lhc::LHCMaxBunches / timing[layer].lengthBC), static_cast(rofs[staggered ? layer : 0].size())); + state.overlap.defineLayer(layer, {static_cast(nROFs), static_cast(timing[layer].lengthBC), static_cast(alpide.getROFDelayInBC(layer)), static_cast(alpide.getROFBiasInBC(layer)), static_cast(config.frameMarginBC.empty() ? 0 : config.frameMarginBC[layer])}); + } + state.overlap.init(); + state.mask = Mask{state.overlap}; + static const o2::itsmft::ChipMappingITS mapping; + for (int stream = 0; stream < streams; ++stream) { + const auto offset = state.clusterROFs.size(); + state.clusterROFs.resize(offset + clusters[stream].size(), {-1, -1}); + std::size_t next = 0; + for (std::size_t r = 0; r < rofs[stream].size(); ++r) { + const auto& record = rofs[stream][r]; + if (record.getFirstEntry() < 0 || static_cast(record.getFirstEntry()) != next || record.getNEntries() < 0 || static_cast(record.getNEntries()) > clusters[stream].size() - next) { + throw std::runtime_error("ITS selection: cluster ROFs do not cover the stream in order"); + } + for (int i = 0; i < record.getNEntries(); ++i, ++next) { + const int layer = staggered ? stream : mapping.getLayer(clusters[stream][next].getSensorID()); + state.clusterROFs[offset + next] = {layer, static_cast(r)}; + } + } + if (next != clusters[stream].size()) { + throw std::runtime_error("ITS selection: cluster ROFs do not cover the stream"); + } + } + const auto firstOrbit = pc.services().get().firstTForbit; + const o2::InteractionRecord firstIR{0, firstOrbit}; + std::vector triggers; + if (triggerType() == 1) { + const auto input = pc.inputs().get>("its_selection_triggers"); + triggers.assign(input.begin(), input.end()); + } else if (triggerType() == 2) { + auto lastIR = firstIR; + const auto input = pc.inputs().get>("its_selection_triggers"); + for (const auto& trigger : input) { + if (trigger.getBCData() >= lastIR && trigger.getNumberOfTracklets()) { + lastIR = trigger.getBCData(); + triggers.push_back({.ir = lastIR, .data = 0}); + } + } + } + // The native estimator indexes the selection array with the trigger ROF; + // discard triggers outside this input before handing it that array. + const auto& mult = o2::its::FastMultEstConfig::Instance(); + const int layer = mult.isMultCutRequested() ? std::clamp(mult.cutMultClusLayer, 0, kLayers - 1) : state.overlap.getView().getClock(); + const auto count = rofs[staggered ? layer : 0].size(); + std::erase_if(triggers, [&](const auto& trigger) { + const auto bc = trigger.ir.differenceInBC(firstIR); + return bc < 0 || bc > std::numeric_limits::max() || static_cast(state.overlap.getLayer(layer).getROF(static_cast(bc))) >= count; + }); + Estimator estimator; + const auto selected = estimator.selectROFs(rofs, clusters, triggers, firstOrbit, staggered, state.overlap.getView(), state.mask); + LOGP(info, "ITS ROF selection: {} of {} ROFs on selection layer {}", selected, count, layer); +} + +ITSROFSelection::~ITSROFSelection() = default; + +bool ITSROFSelection::selectWithVertices(std::span seeds, std::span contributors) +{ + if (!o2::its::FastMultEstConfig::Instance().isVtxMultCutRequested()) { + return false; + } + if (seeds.size() != contributors.size()) { + throw std::runtime_error("ITS vertex selection: contributor counts do not match the vertices"); + } + std::vector vertices; + vertices.reserve(seeds.size()); + for (std::size_t v = 0; v < seeds.size(); ++v) { + auto& vertex = vertices.emplace_back(); + vertex.setNContributors(static_cast(std::min(contributors[v], std::numeric_limits::max()))); + const auto start = bcOfPS(seeds[v].timePS); + const auto end = bcOfPS(static_cast(seeds[v].timePS) + seeds[v].widthPS); + vertex.setTimeStamp({static_cast(start), static_cast(end - start)}); + } + Estimator{}.selectROFsWithVertices(vertices, mState->overlap.getView(), mState->mask); + return true; +} + +std::vector ITSROFSelection::clockMask() const +{ + const auto clock = mState->overlap.getView().getClock(); + std::vector result(mState->overlap.getLayer(clock).mNROFsTF); + const auto mask = mState->mask.getView(); + for (std::size_t r = 0; r < result.size(); ++r) { + result[r] = mask.isROFEnabled(clock, static_cast(r)); + } + return result; +} + +std::optional ITSROFSelection::filter(MeasurementView source, std::span offsets, std::span originalRows) const +{ + const auto mask = mState->mask.getView(); + const auto accepted = [&](std::uint32_t row) { + const auto cluster = source.clusterIndex[row]; + if (cluster >= mState->clusterROFs.size()) { + throw std::runtime_error("ITS selection: measurement cluster index out of range"); + } + const auto [layer, rof] = mState->clusterROFs[cluster]; + return mask.isROFEnabled(layer, rof); + }; + std::uint32_t firstRejected = 0U; + while (firstRejected < offsets.back() && accepted(firstRejected)) { + ++firstRejected; + } + if (firstRejected == offsets.back()) { + return std::nullopt; + } + std::vector rows; + for (std::uint32_t r = 0; r < firstRejected; ++r) { + rows.push_back(r); + } + for (std::uint32_t r = firstRejected + 1U; r < offsets.back(); ++r) { + if (accepted(r)) { + rows.push_back(r); + } + } + Measurements selected; + resizeMeasurementColumns(selected.columns, rows.size()); + selected.originalRows.reserve(rows.size()); + for (std::size_t i = 0; i < rows.size(); ++i) { + Measurement measurement; + measurement.apply(source[rows[i]], [](auto& to, const auto& from) -> auto& { to = from; return to; }); + writeMeasurement(selected.columns, i, measurement); + selected.originalRows.push_back(originalRows.empty() ? rows[i] : originalRows[rows[i]]); + } + selected.offsets.reserve(offsets.size()); + auto at = rows.begin(); + for (const auto offset : offsets) { + at = std::lower_bound(at, rows.end(), offset); + selected.offsets.push_back(static_cast(at - rows.begin())); + } + return selected; +} + +} // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/src/MFT.cxx b/Detectors/GenericTracking/workflow/src/MFT.cxx index 70176703632a4..6f191c05cca3d 100644 --- a/Detectors/GenericTracking/workflow/src/MFT.cxx +++ b/Detectors/GenericTracking/workflow/src/MFT.cxx @@ -36,9 +36,9 @@ std::vector alice::placeChips() } template <> -o2::framework::DataProcessorSpec workflow::getClusterConverterSpec(bool useMC) +o2::framework::DataProcessorSpec workflow::getClusterConverterSpec(bool useMC, bool staggered) { - return workflow::getITSMFTClusterConverterSpec(useMC); + return workflow::getITSMFTClusterConverterSpec(useMC, staggered); } template o2::framework::DataProcessorSpec workflow::getMeasurementProducerSpec(); diff --git a/Detectors/GenericTracking/workflow/src/RecoWorkflow.h b/Detectors/GenericTracking/workflow/src/RecoWorkflow.h index 78edc26f25a51..48ccfbf5bbd17 100644 --- a/Detectors/GenericTracking/workflow/src/RecoWorkflow.h +++ b/Detectors/GenericTracking/workflow/src/RecoWorkflow.h @@ -17,6 +17,7 @@ #define GTRK_WORKFLOW_RECOWORKFLOW_H_ #include "GenericTracking/Workflow/ClusterConverterSpec.h" +#include "GenericTracking/Workflow/ITSMFTOptions.h" #include "GenericTracking/Workflow/MeasurementProducerSpec.h" #include "GenericTracking/Workflow/TrackerSpec.h" @@ -44,6 +45,7 @@ inline void addRecoOptions(std::vector& options, options.push_back(ConfigParamSpec{"clusters-from-upstream", VariantType::Bool, false, {"the clusters from the workflow before this one in the pipe, not O2's readers"}}); options.push_back(ConfigParamSpec{"vertices-from-upstream", VariantType::Bool, false, {"the native vertex seeds (GTRK/VERTICES) from the workflow before this one in the pipe"}}); if (o2Format) { + o2::itsmft::DPLAlpideParamInitializer::addConfigOption(options); options.push_back(ConfigParamSpec{"disable-adaptor", VariantType::Bool, false, {"only the native tracks and vertices (GTRK/...), not O2's formats nor their files"}}); options.push_back(ConfigParamSpec{"disable-root-output", VariantType::Bool, false, {"O2's formats without O2's writers"}}); options.push_back(ConfigParamSpec{"tracking-in-gpu-reco", @@ -79,14 +81,14 @@ template if (!context.options().get("clusters-from-upstream")) { specs.insert(specs.end(), readers.begin(), readers.end()); } - specs.push_back(getClusterConverterSpec(useMC)); + specs.push_back(getClusterConverterSpec(useMC, staggeredInput(context))); specs.push_back(getMeasurementProducerSpec()); const bool inGPUReco = trackingInGPUReco(context); if (inGPUReco && context.options().get("vertices-from-upstream")) { throw std::invalid_argument("--tracking-in-gpu-reco takes no native vertex seeds (--vertices-from-upstream)"); } if (!inGPUReco) { - specs.push_back(getTrackerSpec(context.options().get("vertices-from-upstream"))); + specs.push_back(getTrackerSpec(context.options().get("vertices-from-upstream"), useMC)); } return specs; } diff --git a/Detectors/GenericTracking/workflow/src/TRK.cxx b/Detectors/GenericTracking/workflow/src/TRK.cxx index a6aeb78294e95..fcebdf2f14e03 100644 --- a/Detectors/GenericTracking/workflow/src/TRK.cxx +++ b/Detectors/GenericTracking/workflow/src/TRK.cxx @@ -56,9 +56,11 @@ #include #include #include +#include #include #include #include +#include #include namespace o2::gtrk::alice @@ -69,11 +71,10 @@ namespace // TRKFT3.h restates O2's TRK sensors; here they are held to O2's, to the bit. namespace specs = o2::trk::constants; +// The four-petal VD widths are derived from these radii and gaps. O2 now +// computes widths at run time from the loaded geometry's petal count. static_assert(detectors::trkft3::kVDRadius == specs::VD::petal::layer::radii); static_assert(detectors::trkft3::kVDPetalGap == specs::VD::petal::layer::gaps); -static_assert(exactlyEqual(detectors::trkft3::vdWidth(0), specs::VD::petal::layer::width[0])); -static_assert(exactlyEqual(detectors::trkft3::vdWidth(1), specs::VD::petal::layer::width[1])); -static_assert(exactlyEqual(detectors::trkft3::vdWidth(2), specs::VD::petal::layer::width[2])); static_assert(exactlyEqual(detectors::trkft3::kVDLength, specs::VD::petal::layer::length)); static_assert(exactlyEqual(detectors::trkft3::kVDThickness, specs::VD::silicon::thickness)); static_assert(exactlyEqual(detectors::trkft3::kVDPitch, specs::VD::petal::layer::pitchX) && exactlyEqual(detectors::trkft3::kVDPitch, specs::VD::petal::layer::pitchZ)); @@ -247,12 +248,24 @@ class TrkClusterConverterTask : public Task [[nodiscard]] static o2::header::DataOrigin originOf(int layer) { return layer < kBarrelLayers ? o2::header::DataOrigin{"TRK"} : o2::header::DataOrigin{"FT3"}; } [[nodiscard]] static o2::header::DataHeader::SubSpecificationType subspecOf(int layer) { return static_cast(layer < kBarrelLayers ? layer : layer - kBarrelLayers); } - TrkClusterConverterTask(bool useMC, TrkDictionaryConditions dictionaries) : mUseMC(useMC), mDictionaries(dictionaries) {} + TrkClusterConverterTask(bool useMC, TrkDictionaryConditions dictionaries, std::shared_ptr geometryRequest) + : mUseMC(useMC), mDictionaries(dictionaries), mGeometryRequest(std::move(geometryRequest)) {} - void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final { mDictionaries.finaliseCCDB(matcher, obj); } + void init(InitContext&) final { o2::base::GRPGeomHelper::instance().setRequest(mGeometryRequest); } + + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final + { + if (!o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + mDictionaries.finaliseCCDB(matcher, obj); + } + } void run(ProcessingContext& pc) final { + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + // Dictionary coordinates use the VD segmentation's petal count, which + // GeometryTGeo derives from the loaded geometry in this device's process. + o2::trk::GeometryTGeo::Instance(); TrkDictionaryConditions::checkUpdates(pc); const auto& timing = pc.services().get(); std::vector clusters; @@ -300,12 +313,13 @@ class TrkClusterConverterTask : public Task private: bool mUseMC; TrkDictionaryConditions mDictionaries; + std::shared_ptr mGeometryRequest; }; } // namespace template <> -DataProcessorSpec getClusterConverterSpec(bool useMC) +DataProcessorSpec getClusterConverterSpec(bool useMC, bool /*staggered*/) { std::vector inputs; std::vector outputs{OutputSpec{gDataOrigin, description(kClusters), 0, Lifetime::Timeframe}}; @@ -323,11 +337,12 @@ DataProcessorSpec getClusterConverterSpec(bool useMC) outputs.emplace_back(gDataOrigin, description(kClustersMCTruth), 0, Lifetime::Timeframe); } const TrkDictionaryConditions dictionaries(inputs); + auto geometryRequest = std::make_shared(false, false, false, false, false, o2::base::GRPGeomRequest::Aligned, inputs); return DataProcessorSpec{ .name = deviceName("cluster-converter"), .inputs = inputs, .outputs = outputs, - .algorithm = AlgorithmSpec{adaptFromTask(useMC, dictionaries)}, + .algorithm = AlgorithmSpec{adaptFromTask(useMC, dictionaries, geometryRequest)}, .options = Options{}, }; } diff --git a/Detectors/GenericTracking/workflow/src/TrackerDevice.h b/Detectors/GenericTracking/workflow/src/TrackerDevice.h index 2594ee157889a..9d334b7f45d9d 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerDevice.h +++ b/Detectors/GenericTracking/workflow/src/TrackerDevice.h @@ -31,6 +31,7 @@ #include "GenericTracking/Workflow/DataSpecs.h" #include "GenericTracking/Workflow/GPURecoTracking.h" #include "GenericTracking/Workflow/GridConfig.h" +#include "GenericTracking/Workflow/ITSROFSelection.h" #include "GenericTracking/Workflow/RecoParam.h" #include "GenericTracking/Workflow/TrackerSpec.h" #include "GenericTracking/Workflow/TrackingInterface.h" @@ -60,7 +61,7 @@ namespace o2::gtrk::workflow { /// The tracking device's products made in \a pc's outputs, empty, published on -/// every way out: (tracks, vertices, contributors), DPL's vectors. +/// every way out: (tracks, vertices, contributors, claims, selectedROFs), DPL's vectors. template [[nodiscard]] auto makeTrackerProducts(o2::framework::ProcessingContext& pc) { @@ -70,7 +71,10 @@ template auto& tracks = pc.outputs().make>>(Output{gDataOrigin, description(kTracks), 0}); auto& vertices = pc.outputs().make>(Output{gDataOrigin, description(kVertexSeeds), 0}); auto& contributors = pc.outputs().make>(Output{gDataOrigin, description(kVertexContributors), 0}); - return std::tie(tracks, vertices, contributors); + static_assert(o2::framework::is_messageable::value); + auto& claims = pc.outputs().make>(Output{gDataOrigin, description(kVertexClaims), 0}); + auto& selectedROFs = pc.outputs().make>(Output{gDataOrigin, description(kSelectedROFs), 0}); + return std::tie(tracks, vertices, contributors, claims, selectedROFs); } /// The products' output specs. @@ -78,7 +82,7 @@ template [[nodiscard]] std::vector trackerOutputs() { using o2::framework::Lifetime; - return {{gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexSeeds), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexContributors), 0, Lifetime::Timeframe}}; + return {{gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexSeeds), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexContributors), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexClaims), 0, Lifetime::Timeframe}, {gDataOrigin, description(kSelectedROFs), 0, Lifetime::Timeframe}}; } /// One timeframe of \a Det's tracking, on \a interface's host or device @@ -86,11 +90,15 @@ template /// tracks, the products. The tracking device's run, and the GPU /// reconstruction's (GPURecoTracking.h). template -void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface& interface, RunConditions& conditions, bool withVertices) +void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface& interface, RunConditions& conditions, bool withVertices, bool useMC) { conditions.checkUpdates(pc); const auto orbit = pc.services().get().firstTForbit; - auto [tracks, vertices, contributors] = makeTrackerProducts(pc); + auto [tracks, vertices, contributors, claims, selectedROFs] = makeTrackerProducts(pc); + // Keep copied inputs alive until device work is synchronized on every + // exit, including an allocation failure during upload or vertexing. + std::optional beforeVertices, afterVertices; + std::vector trackingSeeds; const auto ended = interface.endTimeframeOnExit(); const auto input = readTimeframe(pc, withVertices); @@ -104,25 +112,81 @@ void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface return; } - const TimeFrame frame{measurementView(input.columns), interface.grid(), input.offsets}; - const auto found = interface.vertices(frame, input.seeds, false); + std::unique_ptr selection; + auto measurements = measurementView(input.columns); + auto offsets = input.offsets; + std::span originalRows; + if constexpr (std::same_as) { + selection = std::make_unique(pc, conditions.frameTiming(), conditions.orbitsPerTimeframe()); + beforeVertices = selection->filter(measurements, offsets); + if (beforeVertices) { + measurements = measurementView(beforeVertices->columns); + offsets = beforeVertices->offsets; + originalRows = beforeVertices->originalRows; + } + if (offsets.back() == 0U) { + const auto mask = selection->clockMask(); + selectedROFs.assign(mask.begin(), mask.end()); + return; + } + } + const TimeFrame vertexFrame{measurements, interface.grid(), offsets}; + const auto found = interface.vertices(vertexFrame, input.seeds, useMC); if (found.dropped) { return; } - const auto tracked = interface.track(frame, found.seeds); + // Snapshot before a vertex-multiplicity cut can end the vertexing pass + // and rebuild the tracking input. Publish nothing if tracking then fails. + std::vector foundVertices; + std::vector foundContributors; + std::vector foundClaims; + trackingSeeds.assign(found.seeds.begin(), found.seeds.end()); + if (found.found != nullptr) { + foundVertices = found.found->vertices; + foundContributors = found.found->contributors; + for (auto claim : found.claims) { + if (claim.vertex >= 0) { + if (!originalRows.empty()) { + for (auto& row : claim.rows) { + row = originalRows[row]; + } + } + foundClaims.push_back(claim); + } + } + } + if (selection && selection->selectWithVertices(foundVertices, foundContributors)) { + afterVertices = selection->filter(measurements, offsets, originalRows); + if (afterVertices) { + interface.endTimeframe(); + measurements = measurementView(afterVertices->columns); + offsets = afterVertices->offsets; + originalRows = afterVertices->originalRows; + } + } + const TimeFrame frame{measurements, interface.grid(), offsets}; + const auto tracked = interface.track(frame, trackingSeeds); if (tracked.dropped) { return; } - // The vertices with the tracks only: a dropped timeframe publishes nothing. - if (found.found != nullptr) { - vertices.assign(found.found->vertices.begin(), found.found->vertices.end()); - contributors.assign(found.found->contributors.begin(), found.found->contributors.end()); + vertices.assign(foundVertices.begin(), foundVertices.end()); + contributors.assign(foundContributors.begin(), foundContributors.end()); + claims.assign(foundClaims.begin(), foundClaims.end()); + if (selection) { + const auto mask = selection->clockMask(); + selectedROFs.assign(mask.begin(), mask.end()); } const auto& result = *tracked.result.tracks; tracks.reserve(result.size()); std::size_t failed = 0U; for (std::size_t track = 0; track < result.size(); ++track) { - tracks.push_back(result.read(track)); + auto record = result.read(track); + if (!originalRows.empty()) { + for (std::uint8_t c = 0; c < record.count; ++c) { + record.clusters.index[c] = originalRows[record.clusters.index[c]]; + } + } + tracks.push_back(record); failed += tracks.back().status == PropagationStatus::kOk ? 0U : 1U; } LOGP(info, "{}: the timeframe at orbit {}: {} tracks ({} with a failed refit), {} vertices, {:.3f} s tracking, {:.3f} s vertexing", Det::name, orbit, tracks.size(), failed, vertices.size(), tracked.seconds, found.seconds); @@ -153,8 +217,9 @@ class TrackerTask : public o2::framework::Task public: using RunConditions = Conditions; - TrackerTask(bool withVertices, std::shared_ptr conditions) + TrackerTask(bool withVertices, bool useMC, std::shared_ptr conditions) : mWithVertices(withVertices), + mUseMC(useMC), mInterface(withVertices), mConditions(std::move(conditions)) { @@ -176,27 +241,31 @@ class TrackerTask : public o2::framework::Task void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final { mConditions->finaliseCCDB(matcher, obj); } - void run(o2::framework::ProcessingContext& pc) final { trackTimeframe(pc, mInterface, *mConditions, mWithVertices); } + void run(o2::framework::ProcessingContext& pc) final { trackTimeframe(pc, mInterface, *mConditions, mWithVertices, mUseMC); } void endOfStream(o2::framework::EndOfStreamContext& /*ec*/) final { reportDropped(mInterface); } private: bool mWithVertices = false; + bool mUseMC = false; TrackingInterface mInterface; std::shared_ptr mConditions; ///< the run's field, geometry, MeanVertex and material LUT (Conditions.h) }; template -o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices) +o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices, bool useMC) { using namespace o2::framework; auto inputs = measurementInputs(withVertices); + if constexpr (std::same_as) { + ITSROFSelection::addInputs(inputs); + } auto conditions = std::make_shared>(inputs); return DataProcessorSpec{ .name = deviceName("tracker"), .inputs = inputs, .outputs = trackerOutputs(), - .algorithm = AlgorithmSpec{adaptFromTask>(withVertices, conditions)}, + .algorithm = AlgorithmSpec{adaptFromTask>(withVertices, useMC, conditions)}, .options = Options{{"gtrk-gpu-backend", VariantType::String, "", {"the GPU backend to track on, cuda or hip; the host if empty"}}, {"gtrk-gpu-device", VariantType::Int, 0, {"the GPU backend's device index"}}, {"gtrk-gpu-arena-mb", VariantType::Int, 0, {"the device tracker's arena in O2's GPU pool [MB]; required with a GPU backend"}}}, }; } @@ -212,11 +281,14 @@ class GPURecoTrackingOf final : public GPURecoTracking /// No native seeds: the vertex list, where an iteration needs one, is the /// own vertexer's. - GPURecoTrackingOf() : mInterface(false) {} + explicit GPURecoTrackingOf(bool useMC) : mUseMC(useMC), mInterface(false) {} void appendInputs(std::vector& inputs) final { auto own = measurementInputs(false); + if constexpr (std::same_as) { + ITSROFSelection::addInputs(own); + } mConditions = std::make_unique(own); inputs.insert(inputs.end(), own.begin(), own.end()); } @@ -245,7 +317,7 @@ class GPURecoTrackingOf final : public GPURecoTracking // Given back after the timeframe's end, which detaches the tracker. const auto borrowed = mInterface.context().borrow(tag()); mInterface.useBlock(borrowed->block()); - trackTimeframe(pc, mInterface, *mConditions, false); + trackTimeframe(pc, mInterface, *mConditions, false, mUseMC); LOGP(info, "{}: the GPU reconstruction's pool lent {:.1f} MB, the arena's peak {:.1f} MB", Det::name, static_cast(borrowed->block().bytes) / 1.e6, static_cast(mInterface.lastArenaPeakBytes()) / 1.e6); } @@ -266,14 +338,15 @@ class GPURecoTrackingOf final : public GPURecoTracking return packed; } + bool mUseMC = false; TrackingInterface mInterface; std::unique_ptr mConditions; ///< the detector's own conditions; the run-wide ones are the GPU reconstruction's }; template -std::unique_ptr makeGPURecoTrackingOf() +std::unique_ptr makeGPURecoTrackingOf(bool useMC) { - return std::make_unique>(); + return std::make_unique>(useMC); } } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/src/TrackerITS.cxx b/Detectors/GenericTracking/workflow/src/TrackerITS.cxx index ac0802b711b4d..033287c486157 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerITS.cxx +++ b/Detectors/GenericTracking/workflow/src/TrackerITS.cxx @@ -21,7 +21,7 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices); -template std::unique_ptr workflow::makeGPURecoTrackingOf(); +template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices, bool useMC); +template std::unique_ptr workflow::makeGPURecoTrackingOf(bool useMC); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/TrackerMFT.cxx b/Detectors/GenericTracking/workflow/src/TrackerMFT.cxx index 6167ddfda6944..0a9df8abdb37e 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerMFT.cxx +++ b/Detectors/GenericTracking/workflow/src/TrackerMFT.cxx @@ -21,7 +21,7 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices); -template std::unique_ptr workflow::makeGPURecoTrackingOf(); +template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices, bool useMC); +template std::unique_ptr workflow::makeGPURecoTrackingOf(bool useMC); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/TrackerTRKFT3.cxx b/Detectors/GenericTracking/workflow/src/TrackerTRKFT3.cxx index f91794c10d631..0d6946d216083 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerTRKFT3.cxx +++ b/Detectors/GenericTracking/workflow/src/TrackerTRKFT3.cxx @@ -20,6 +20,6 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices); +template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices, bool useMC); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx b/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx index 36e47dec18fa9..9c3b51cf3f6b5 100644 --- a/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx +++ b/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx @@ -45,15 +45,21 @@ void customize(std::vector& policies) } // NOLINTNEXTLINE(misc-use-internal-linkage) -void customize(std::vector& workflowOptions) { o2::gtrk::workflow::addRecoOptions(workflowOptions, true); } +void customize(std::vector& workflowOptions) +{ + o2::gtrk::workflow::addRecoOptions(workflowOptions, true); + workflowOptions.push_back({"select-with-triggers", VariantType::String, "none", {"use triggers to prescale processed ROFs: phys, trd, none"}}); +} #include "GenericTracking/Detectors/ITS.h" #include "GenericTracking/Workflow/AdaptorSpec.h" +#include "GenericTracking/Workflow/RecoParam.h" #include "GlobalTrackingWorkflowWriters/IRFrameWriterSpec.h" #include "ITSMFTWorkflow/ClusterReaderSpec.h" #include "ITSWorkflow/TrackWriterSpec.h" +#include #include #include "Framework/runDataProcessing.h" // IWYU pragma: keep @@ -62,13 +68,22 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) { o2::conf::ConfigurableParam::updateFromString(context.options().get("configKeyValues")); + const auto trigger = context.options().get("select-with-triggers"); + if (trigger != "none") { + if (trigger != "phys" && trigger != "trd") { + throw std::invalid_argument("--select-with-triggers expects phys, trd or none"); + } + o2::conf::ConfigurableParam::updateFromString("GTRKRecoITS.selectWithTriggers=" + std::string(trigger == "phys" ? "1" : "2")); + } + namespace gw = o2::gtrk::workflow; using Det = o2::gtrk::detectors::ITS; const bool useMC = !context.options().get("disable-mc"); - auto specs = gw::trackingChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, false, true, false)}); + const bool staggered = gw::staggeredInput(context); + auto specs = gw::trackingChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, staggered, true, o2::gtrk::recoParam().selectWithTriggers == 1)}); // The adaptor reads the tracks: after the GPU reconstruction when the tracking runs inside it. if (!context.options().get("disable-adaptor") && !gw::trackingInGPUReco(context)) { - specs.push_back(gw::getAdaptorSpec(useMC)); + specs.push_back(gw::getAdaptorSpec(useMC, staggered)); if (!context.options().get("disable-root-output")) { specs.push_back(o2::its::getTrackWriterSpec(useMC)); specs.push_back(o2::globaltracking::getIRFrameWriterSpec("irfr:ITS/IRFRAMES/0", "o2_its_irframe.root", "irframe-writer-its")); diff --git a/Detectors/GenericTracking/workflow/src/gtrk-mft-reco-workflow.cxx b/Detectors/GenericTracking/workflow/src/gtrk-mft-reco-workflow.cxx index 105754b5fbc30..9b2bb5e1a02a4 100644 --- a/Detectors/GenericTracking/workflow/src/gtrk-mft-reco-workflow.cxx +++ b/Detectors/GenericTracking/workflow/src/gtrk-mft-reco-workflow.cxx @@ -63,10 +63,11 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) namespace gw = o2::gtrk::workflow; using Det = o2::gtrk::detectors::MFT; const bool useMC = !context.options().get("disable-mc"); - auto specs = gw::trackingChain(context, {o2::itsmft::getMFTClusterReaderSpec(useMC, false, true, false)}); + const bool staggered = gw::staggeredInput(context); + auto specs = gw::trackingChain(context, {o2::itsmft::getMFTClusterReaderSpec(useMC, staggered, true, false)}); // The adaptor reads the tracks: after the GPU reconstruction when the tracking runs inside it. if (!context.options().get("disable-adaptor") && !gw::trackingInGPUReco(context)) { - specs.push_back(gw::getAdaptorSpec(useMC)); + specs.push_back(gw::getAdaptorSpec(useMC, staggered)); if (!context.options().get("disable-root-output")) { specs.push_back(o2::mft::getTrackWriterSpec(useMC)); } diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h index d283d6eed2d45..d79d86e98975c 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h @@ -24,6 +24,7 @@ #include "ITStracking/FastMultEstConfig.h" #include "ITSMFTTracking/ROFLookupTables.h" #include +#include #include namespace o2::its @@ -62,10 +63,14 @@ struct FastMultEst { if (!multEstConf.isPassingVtxMultCut(vertex.getNContributors())) { const auto& timestamp{vertex.getTimeStamp()}; for (int layer = 0; layer < NLayers; ++layer) { - uint32_t startROF = sel.getLayer(layer).getROF(timestamp.lower()); - uint32_t endROF = sel.getLayer(layer).getROF(timestamp.upper()); - for (uint32_t rof = startROF; rof <= endROF; ++rof) { - sel.setROFsEnabled(layer, rof, 0); + const auto& timing = sel.getLayer(layer); + if (timing.mNROFsTF == 0) { + continue; + } + const auto startROF = timing.getROF(timestamp.lower()); + const auto endROF = std::min(timing.getROF(timestamp.upper()), timing.mNROFsTF - 1U); + if (startROF <= endROF && startROF < timing.mNROFsTF) { + sel.setROFsEnabled(layer, startROF, endROF - startROF + 1U, 0); } } } diff --git a/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx b/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx index 2168c4b8d3308..8f872361d35d8 100644 --- a/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx +++ b/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx @@ -315,7 +315,10 @@ void Clusterer::ClustererThread::finishChipSingleHitFast(uint32_t hit, ChipPixel if (labelsClusPtr) { // MC labels were requested int nlab = 0; - fetchMCLabels(curChipData->getStartID() + hit, labelsDigPtr, nlab); + // Squashing can remove earlier digits on the chip, so the surviving + // pixel's source index need not be contiguous with the chip's start. + const auto digitID = parent->mSquashingDepth ? curChipData->getOrderedPixId(hit) : curChipData->getStartID() + hit; + fetchMCLabels(digitID, labelsDigPtr, nlab); auto cnt = compClusPtr->size(); for (int i = 0; i < nlab; i++) { labelsClusPtr->addElement(cnt, labelsBuff[i]); diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/SegmentationChip.h b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/SegmentationChip.h index fac5198966fd6..117a2a2c2c590 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/SegmentationChip.h +++ b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/SegmentationChip.h @@ -72,6 +72,20 @@ class SegmentationChip static constexpr std::array radiiVD = constants::VD::petal::layer::radii; + static inline int sNPetalsVD = -1; // not set + static void setNumberOfPetalsVD(int n) { sNPetalsVD = n; } + static float getWidthVD(int layer) // get the width of the petal according to the simulated geometry (supported: kIRIS4 (4 petals), kIRIS4a (3 petals), kIRISFullCyl) + { + if (sNPetalsVD < 1) { + LOGP(fatal, "SegmentationChip: VD number of petals not set. GeometryTGeo::Build() must be called before using the VD segmentation"); + } + if (sNPetalsVD == 1) { // FullCyl: 360°, no gap + return 2 * M_PI * constants::VD::petal::layer::radii[layer]; + } + return constants::VD::petal::layer::radii[layer] * 2 * M_PI / sNPetalsVD - constants::VD::petal::layer::gaps[layer]; + } + static int getNRowsVD(int layer) { return static_cast(getWidthVD(layer) / PitchRowVD); } + /// Transformation from Geant detector centered local coordinates (cm) to /// Pixel cell numbers iRow and iCol. /// Returns kTRUE if point x,z is inside sensitive volume, kFALSE otherwise. @@ -113,7 +127,7 @@ class SegmentationChip if (subDetID == 0) { pitchRow = PitchRowVD; pitchCol = PitchColVD; - maxWidth = constants::VD::petal::layer::width[layer]; + maxWidth = getWidthVD(layer); maxLength = constants::VD::petal::layer::length; // TODO: change this to use the layer and disk } else if (subDetID == 1 || subDetID == 2) { @@ -128,11 +142,11 @@ class SegmentationChip }; // Check local coordinates (cm) validity. - static constexpr bool isValidLoc(float x, float z, int subDetID, int layer) noexcept + static bool isValidLoc(float x, float z, int subDetID, int layer) noexcept { float maxWidth(0), maxLength(0); if (subDetID == 0) { - maxWidth = constants::VD::petal::layer::width[layer]; + maxWidth = getWidthVD(layer); maxLength = constants::VD::petal::layer::length; // TODO: change this to use the layer and disk } else if (subDetID == 1 || subDetID == 2) { // ML/OT @@ -143,12 +157,12 @@ class SegmentationChip } // Check detector coordinates validity. - static constexpr bool isValidDet(int row, int col, int subDetID, int layer) noexcept + static bool isValidDet(int row, int col, int subDetID, int layer) noexcept { // Check if the row and column are within the valid range int nRows(0), nCols(0); if (subDetID == 0) { - nRows = constants::VD::petal::layer::nRows[layer]; + nRows = getNRowsVD(layer); nCols = constants::VD::petal::layer::nCols; // TODO: change this to use the layer and disk } else if (subDetID == 1 || subDetID == 2) { @@ -171,7 +185,7 @@ class SegmentationChip /// \param int subDetID Sub-detector ID (0 for VD, 1 for ML/OT) /// \param int layer Layer number (0 to 2 for VD, 0 to 7 for ML/OT) /// \param int disk Disk number (0 to 5 for VD) - static constexpr bool detectorToLocal(int iRow, int iCol, float& xRow, float& zCol, int subDetID, int layer, int disk) noexcept + static bool detectorToLocal(int iRow, int iCol, float& xRow, float& zCol, int subDetID, int layer, int disk) noexcept { if (!isValidDet(iRow, iCol, subDetID, layer)) { LOGP(debug, "Detector coordinates not valid: iRow = {}, iCol = {}", iRow, iCol); @@ -194,7 +208,7 @@ class SegmentationChip /// xRow = half chip width - iRow(center) * pitch /// zCol = iCol * pitch - half chip lenght if (subDetID == 0) { - xRow = 0.5 * (constants::VD::petal::layer::width[layer] - PitchRowVD) - (row * PitchRowVD); + xRow = 0.5 * (getWidthVD(layer) - PitchRowVD) - (row * PitchRowVD); zCol = col * PitchColVD + 0.5 * (PitchColVD - constants::VD::petal::layer::length); } else if (subDetID == 1 || subDetID == 2) { // ML/OT xRow = 0.5 * (constants::moduleMLOT::chip::width - constants::moduleMLOT::chip::passiveEdgeReadOut - PitchRowMLOT) - (row * PitchRowMLOT); @@ -251,9 +265,9 @@ class SegmentationChip /// Print segmentation info static void Print() noexcept { - LOG(info) << "Number of rows:\nVD L0: " << constants::VD::petal::layer::nRows[0] - << "\nVD L1: " << constants::VD::petal::layer::nRows[1] - << "\nVD L2: " << constants::VD::petal::layer::nRows[2] + LOG(info) << "Number of rows:\nVD L0: " << getNRowsVD(0) + << "\nVD L1: " << getNRowsVD(1) + << "\nVD L2: " << getNRowsVD(2) << "\nML/OT chip: " << constants::moduleMLOT::chip::nRows; LOG(info) << "Number of cols:\nVD: " << constants::VD::petal::layer::nCols diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/Specs.h b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/Specs.h index 9ac1f789d43ba..866178f820231 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/Specs.h +++ b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/Specs.h @@ -48,15 +48,13 @@ constexpr int nLayers{3}; // number of layers in each VD petal constexpr int nDisks{6}; // number of disks in each VD petal namespace layer { -constexpr double pitchX{10 * mu}; // pitch of the row -constexpr double pitchZ{10 * mu}; // pitch of the column -constexpr double totalThickness{silicon::thickness + metalstack::thickness}; // total thickness of the chip -constexpr std::array gaps{1.63 * mm, 1.2 * mm, 1.2 * mm}; // gaps between two consecutive petals -constexpr std::array radii{0.5 * cm, 1.2 * cm, 2.5 * cm}; // radius of layer in cm -constexpr std::array width{radii[0] * 2 * M_PI / 4 - gaps[0], radii[1] * 2 * M_PI / 4 - gaps[1], radii[2] * 2 * M_PI / 4 - gaps[2]}; // width of the quarter of layer in cm -constexpr double length{50 * cm}; // length of the layer -constexpr int nCols{static_cast(length / pitchZ)}; // number of columns in the chip -constexpr std::array nRows{static_cast(width[0] / pitchX), static_cast(width[1] / pitchX), static_cast(width[2] / pitchX)}; // number of rows in the chip. For the moment is different for each layer since a siner segmentation in repetitive units is stil to be implemented +constexpr double pitchX{10 * mu}; // pitch of the row +constexpr double pitchZ{10 * mu}; // pitch of the column +constexpr double totalThickness{silicon::thickness + metalstack::thickness}; // total thickness of the chip +constexpr std::array gaps{1.63 * mm, 1.2 * mm, 1.2 * mm}; // gaps between two consecutive petals +constexpr std::array radii{0.5 * cm, 1.2 * cm, 2.5 * cm}; // radius of layer in cm +constexpr double length{50 * cm}; // length of the layer +constexpr int nCols{static_cast(length / pitchZ)}; // number of columns in the chip } // namespace layer namespace disk @@ -223,6 +221,14 @@ constexpr double responseYShift{5 * mu}; /// center of the epitaxial layer constexpr double thickness{20 * mu}; } // namespace alice3resp +namespace alpideresp /// parameters for the alpide chip response +{ +constexpr double pitchX{26.88 * mu}; +constexpr double pitchZ{29.24 * mu}; +constexpr double responseYShift{5 * mu}; /// center of the epitaxial layer +constexpr double thickness{45 * mu}; +} // namespace alpideresp + namespace MLOTDisks { constexpr int nLayers{12}; // number of disks in the ML and OT diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/src/GeometryTGeo.cxx index ddfc844cc964d..699b65dd31df6 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/src/GeometryTGeo.cxx +++ b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/src/GeometryTGeo.cxx @@ -84,6 +84,9 @@ void GeometryTGeo::Build(int loadTrans) mNumberOfLayersMLOT = extractNumberOfLayersMLOT(); mNumberOfPetalsVD = extractNumberOfPetalsVD(); + SegmentationChip::setNumberOfPetalsVD(mNumberOfPetalsVD); + LOGP(info, "VD segmentation set for {} petals", mNumberOfPetalsVD); + mNumberOfActivePartsVD = extractNumberOfActivePartsVD(); mNumberOfLayersVD = extractNumberOfLayersVD(); mNumberOfDisksVD = extractNumberOfDisksVD(); diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/common/macros/CheckTopologyDictionary.C b/Detectors/Upgrades/ALICE3/TRKFT3/common/macros/CheckTopologyDictionary.C index f97a586373308..04e4a5ab686da 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/common/macros/CheckTopologyDictionary.C +++ b/Detectors/Upgrades/ALICE3/TRKFT3/common/macros/CheckTopologyDictionary.C @@ -39,6 +39,7 @@ #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/MCTruthContainer.h" #include "TRKBase/GeometryTGeo.h" +#include "TRKBase/SegmentationChip.h" #include "TRKBase/Specs.h" #include "TRKReconstruction/BuildTopologyDictionary.h" #include @@ -190,7 +191,7 @@ void forEachCluster(const std::string& directory, F&& f) const auto* pattern = candidate.pattern; // Left out: at the matrix's edge, or one of several pieces of a crossing. const bool vd = cluster.subDetID == 0; - const int nRows = vd ? o2::trk::constants::VD::petal::layer::nRows[cluster.layer] : o2::trk::constants::moduleMLOT::chip::nRows; + const int nRows = vd ? o2::trk::SegmentationChip::getNRowsVD(cluster.layer) : o2::trk::constants::moduleMLOT::chip::nRows; const int nCols = vd ? o2::trk::constants::VD::petal::layer::nCols : o2::trk::constants::moduleMLOT::chip::nCols; if (cluster.row == 0 || cluster.col == 0 || cluster.row + pattern[0] >= nRows || cluster.col + pattern[1] >= nCols || nClusters[candidate.key] > static_cast(candidate.segments.size())) { diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/src/Clusterer.cxx index 03e7e5285a31f..8e37df7bfe10f 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/src/Clusterer.cxx @@ -181,7 +181,7 @@ void Clusterer::ClustererThread::initChip(gsl::span digits, if (geom) { if (geom->getSubDetID(chipID) == 0) { // VD const int layer = geom->getLayer(chipID); - size = constants::VD::petal::layer::nRows[layer] + 2; + size = SegmentationChip::getNRowsVD(layer) + 2; } } diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/ChipDigitsContainer.h b/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/ChipDigitsContainer.h index 0fcbfd2591a0e..5176e8cbbcfe7 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/ChipDigitsContainer.h +++ b/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/ChipDigitsContainer.h @@ -96,7 +96,7 @@ void o2::trkft3::ChipDigitsContainer::addNoise(UInt_t rofMin, UInt_t rofMax, con int maxCols = 0; if (subDetID == 0) { - maxRows = o2::trk::constants::VD::petal::layer::nRows[layer]; + maxRows = o2::trk::SegmentationChip::getNRowsVD(layer); maxCols = o2::trk::constants::VD::petal::layer::nCols; } else { maxRows = o2::trk::constants::moduleMLOT::chip::nRows; diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/Digitizer.h b/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/Digitizer.h index d0cf085d971a0..65629cf7de6fc 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/Digitizer.h @@ -130,7 +130,7 @@ class Digitizer int getNRows(int subDetID, int layer) { if (subDetID == 0) { // VD - return o2::trk::constants::VD::petal::layer::nRows[layer]; + return o2::trk::SegmentationChip::getNRowsVD(layer); } else if (subDetID == 1 || subDetID == 2) { // ML/OT return o2::trk::constants::moduleMLOT::chip::nRows; } diff --git a/GPU/Workflow/src/GPUWorkflowGenericTracking.cxx b/GPU/Workflow/src/GPUWorkflowGenericTracking.cxx index ec91e0ee2ea50..32ada64ae827c 100644 --- a/GPU/Workflow/src/GPUWorkflowGenericTracking.cxx +++ b/GPU/Workflow/src/GPUWorkflowGenericTracking.cxx @@ -22,10 +22,10 @@ namespace o2::gpu void GPURecoWorkflowSpec::createGenericTracking() { if (mSpecConfig.runGenericTrackingITS) { - mGenericTracking.push_back(o2::gtrk::workflow::makeGPURecoTracking("its")); + mGenericTracking.push_back(o2::gtrk::workflow::makeGPURecoTracking("its", mSpecConfig.processMC)); } if (mSpecConfig.runGenericTrackingMFT) { - mGenericTracking.push_back(o2::gtrk::workflow::makeGPURecoTracking("mft")); + mGenericTracking.push_back(o2::gtrk::workflow::makeGPURecoTracking("mft", mSpecConfig.processMC)); } }