Staggering + linearisation references + misc

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<std::vector<std::uint32_t>>;
 
@@ -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<bool> 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<TrackKinematics> 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<o2::track::TrackParCov> innerStates(long long entry);

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<float>::quiet_NaN(); ///< generated momentum at production [GeV]
+  float etaGen = std::numeric_limits<float>::quiet_NaN();
+  std::uint8_t origin = 0xFFU;                                   ///< ParticleOrigin; unknown for unlabelled tracks
+  float invPtResidual = std::numeric_limits<float>::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<std::uint64_t, std::uint64_t>;
 /// \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<const QATrack> tracks, const char* tree = "tracks");
+
 /// \a entries into \a output's `resolution` tree.
 void streamResolution(o2::utils::TreeStreamRedirector& output, std::span<const QAResolution> 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
 

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<Det>()), etaReach(PlacedEnvelopes<Det>::all()), comparisonSettings<Det>(), mKinematics ? mKinematics->getNEvents(0) : 0U);
+      streamRun(*mQAOutput, workflow::iterationCount(recoParam<Det>()), etaReach(PlacedEnvelopes<Det>::all()), comparisonSettings<Det>(), mKinematics ? mKinematics->getNEvents(0) : 0U, std::is_same_v<Det, detectors::TRKFT3> ? 0U : (std::uint64_t{1} << Detector<Det>::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<std::uint64_t, std::uint8_t> outward;
+    std::unordered_map<std::uint64_t, std::uint64_t> outwardMasks;
     std::unordered_map<std::uint64_t, std::vector<float>> legMomenta; // --species-oracle
     if (mHits) {
       std::unordered_map<int, std::vector<std::uint64_t>> byEvent;
@@ -3427,6 +3430,13 @@ class MeasurementSinkTask : public Task
           const auto& truth = particles.at(id);
           const auto leg = outwardLeg(crossings, {static_cast<double>(truth.vx), static_cast<double>(truth.vy), static_cast<double>(truth.vz)});
           outward[id] = outwardLayers(leg);
+          if constexpr (std::is_same_v<Det, detectors::TRKFT3>) {
+            if (outward[id] >= comparisonSettings<Det>().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<Det, detectors::TRKFT3>) {
+        if (const auto mask = outwardMasks.find(id); mask != outwardMasks.end()) {
+          extra.expectedLayerMask = mask->second;
+        }
+      } else {
+        extra.expectedLayerMask = (std::uint64_t{1} << Detector<Det>::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<Det>()), comparisonSettings<Det>()));
+    if (theirs != nullptr) {
+      const auto nativeKinematics = mO2Tracks->kinematics(static_cast<long long>(mTimeframe));
+      const auto nativeQA = trackingQA(rows, particles, *theirs, nativeKinematics, {}, extras, nullptr, 1, comparisonSettings<Det>());
+      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<std::uint32_t>(track));
       }
     }
-    const auto theirs = mO2Tracks->read(static_cast<long long>(mTimeframe));
+    const auto theirs = o2Clusters(frame);
     constexpr auto kMaxClusters = static_cast<std::size_t>(Detector<Det>::nLayers) + 1U;
     std::array<std::size_t, kMaxClusters> o2Total{};
     std::array<std::size_t, kMaxClusters> 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<long long>(mTimeframe));
+    std::array<std::uint32_t, Detector<Det>::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<std::size_t>(o2::itsmft::clusID2Layer(static_cast<int>(cluster)));
+        if (layer == 0U) {
+          continue;
+        }
+        const auto index = static_cast<std::uint32_t>(o2::itsmft::clusID2Index(static_cast<int>(cluster)));
+        cluster = layer + 1U < offsets.size() && index < offsets[layer + 1U] - offsets[layer]
+                    ? offsets[layer] + index : std::numeric_limits<std::uint32_t>::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<long long>(mTimeframe));
+    const auto clusters = o2Clusters(frame);
     TrackRows theirs;
     theirs.reserve(clusters.size());
     for (const auto& track : clusters) {

Detectors/GenericTracking/study/src/TrackingQA.cxx

@@ -221,8 +221,24 @@ TrackingQA trackingQA(std::span<const TruthRow> 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<std::uint8_t>(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<float>((static_cast<double>(record.ptGen) / static_cast<double>(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<const QATrack> 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<const QAMo
   }
 }
 
-void streamRun(o2::utils::TreeStreamRedirector& output, int iterations, EtaRange eta, const ComparisonSettings& settings, std::size_t events) { output << "run" << "iterations=" << iterations << "etaLow=" << eta.low << "etaHigh=" << eta.high << "minLayers=" << settings.minLayers << "minPtGeV=" << settings.minPtGeV << "minOutwardLayers=" << settings.minOutwardLayers << "events=" << static_cast<std::uint64_t>(events) << "\n"; }
+void streamRun(o2::utils::TreeStreamRedirector& output, int iterations, EtaRange eta, const ComparisonSettings& settings, std::size_t events, std::uint64_t fullLengthMask) { output << "run" << "iterations=" << iterations << "etaLow=" << eta.low << "etaHigh=" << eta.high << "minLayers=" << settings.minLayers << "minPtGeV=" << settings.minPtGeV << "minOutwardLayers=" << settings.minOutwardLayers << "events=" << static_cast<std::uint64_t>(events) << "fullLengthMask=" << fullLengthMask << "\n"; }
 
 } // namespace o2::gtrk::study

Detectors/GenericTracking/study/src/Tracks.cxx

@@ -10,6 +10,7 @@
 // or submit itself to any jurisdiction.
 
 #include "GenericTracking/Study/Alice/Tracks.h"
+#include "GenericTracking/Study/TrackingQA.h"
 
 #include "DataFormatsITS/TrackITS.h"
 #include "DataFormatsMFT/TrackMFT.h"
@@ -18,8 +19,10 @@
 #include <TFile.h>
 #include <TTree.h>
 
+#include <cmath>
 #include <format>
 #include <stdexcept>
+#include <type_traits>
 
 namespace o2::gtrk::study
 {
@@ -63,6 +66,34 @@ TrackClusters readEntry(TTree& tree, long long entry, const char* trackBranch, c
   return result;
 }
 
+template <typename Track>
+std::vector<TrackKinematics> readKinematics(TTree& tree, long long entry, const char* branch)
+{
+  std::vector<Track> tracks;
+  auto* pointer = &tracks;
+  tree.SetBranchAddress(branch, &pointer);
+  tree.GetEntry(entry);
+  tree.ResetBranchAddresses();
+  std::vector<TrackKinematics> result;
+  result.reserve(tracks.size());
+  for (const auto& track : tracks) {
+    const auto pt = static_cast<float>(track.getPt());
+    const auto eta = static_cast<float>(track.getEta());
+    const auto phi = static_cast<float>(track.getPhi());
+    float chi2 = 0.F;
+    int ndf = 0;
+    if constexpr (std::is_same_v<Track, o2::its::TrackITS>) {
+      chi2 = track.getChi2();
+      ndf = 2 * track.getNumberOfClusters() - 5;
+    } else {
+      chi2 = static_cast<float>(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<bool> TrackFile::fakeFlags(long long entry)
   return result;
 }
 
+std::vector<TrackKinematics> TrackFile::kinematics(long long entry)
+{
+  if (entry >= mTree->GetEntries()) {
+    return {};
+  }
+  return mChips == Chips::kITS
+           ? readKinematics<o2::its::TrackITS>(*mTree, entry, "ITSTrack")
+           : readKinematics<o2::mft::TrackMFT>(*mTree, entry, "MFTTrack");
+}
+
 std::vector<o2::track::TrackParCov> TrackFile::innerStates(long long entry)
 {
   if (mChips != Chips::kITS) {

Detectors/GenericTracking/study/tools/CompleteAnatomy.h

@@ -4342,7 +4342,11 @@ class CompleteAnatomy
       // The road fit, its last step nominal and allowed.
       const RowChain<Det> chain{.frame = &frame, .rows = road.data()};
       const auto count = static_cast<int>(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;

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<o2::framework::ConfigParamSpec>& 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 <DetectorLayout Det>
   if (!context.options().get<bool>("clusters-from-upstream")) {
     specs.insert(specs.end(), readers.begin(), readers.end());
   }
-  specs.push_back(gw::getClusterConverterSpec<Det>(useMC));
+  specs.push_back(gw::getClusterConverterSpec<Det>(useMC, gw::staggeredInput<Det>(context)));
   specs.push_back(gw::getMeasurementProducerSpec<Det>());
   specs.push_back(getMeasurementSinkSpec<Det>(useMC, mcVertices, products));
   if (products != TrackerProducts::kNone) {