refactor: Rework initial qOverP sigma in Examples (#3422)

- Use relative `pt` resolution as parameter. Get initial `qOverP` sigma by combining relative `pt` resolution and `theta` sigma
- Apply to `ParticleSmearing`
- Wire to python
- Disable smearing for truth smeared seeding to allow track finding encounter measurements

.gitlab-ci.yml

@@ -207,7 +207,7 @@ linux_test_examples:
     - source build/this_acts_withdeps.sh
     - cd src
     - pip3 install -r Examples/Python/tests/requirements.txt
-    - pytest -rFsv -k "not exatrkx" -v
+    - pytest -rFsv -k "not exatrkx" -v -s
 
 linux_physmon:
   stage: test

CI/physmon/workflows/physmon_ckf_tracking.py

@@ -91,9 +91,19 @@ def run_ckf_tracking(truthSmearedSeeded, truthEstimatedSeeded, label):
             setup.trackingGeometry,
             setup.field,
             TruthSeedRanges(pt=(500 * u.MeV, None), nHits=(9, None)),
-            ParticleSmearingSigmas(
-                pRel=0.01
-            ),  # only used by SeedingAlgorithm.TruthSmeared
+            ParticleSmearingSigmas(  # only used by SeedingAlgorithm.TruthSmeared
+                # zero eveything so the CKF has a chance to find the measurements
+                d0=0,
+                d0PtA=0,
+                d0PtB=0,
+                z0=0,
+                z0PtA=0,
+                z0PtB=0,
+                t0=0,
+                phi=0,
+                theta=0,
+                ptRel=0,
+            ),
             SeedFinderConfigArg(
                 r=(33 * u.mm, 200 * u.mm),
                 deltaR=(1 * u.mm, 60 * u.mm),
@@ -125,9 +135,10 @@ def run_ckf_tracking(truthSmearedSeeded, truthEstimatedSeeded, label):
                 1 * u.mm,
                 1 * u.degree,
                 1 * u.degree,
-                0.1 / u.GeV,
+                0.1 * u.e / u.GeV,
                 1 * u.ns,
             ],
+            initialSigmaPtRel=0.01,
             initialVarInflation=[1.0] * 6,
             geoSelectionConfigFile=setup.geoSel,
             rnd=rnd,  # only used by SeedingAlgorithm.TruthSmeared

CI/physmon/workflows/physmon_vertexing.py

@@ -105,9 +105,10 @@ def run_vertexing(fitter, mu, events):
                 1 * u.mm,
                 1 * u.degree,
                 1 * u.degree,
-                0.1 / u.GeV,
+                0.1 * u.e / u.GeV,
                 1 * u.ns,
             ],
+            initialSigmaPtRel=0.1,
             initialVarInflation=[1.0] * 6,
             geoSelectionConfigFile=setup.geoSel,
         )

Core/src/Material/VolumeMaterialMapper.cpp

@@ -402,8 +402,6 @@ void Acts::VolumeMaterialMapper::mapMaterialTrack(
                                 << " at position = (" << mVolumes.position.x()
                                 << ", " << mVolumes.position.y() << ", "
                                 << mVolumes.position.z() << ")");
-
-    mappingVolumes.push_back(mVolumes);
   }
   // Run the mapping process, i.e. take the recorded material and map it
   // onto the mapping volume:

Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackParamsEstimationAlgorithm.hpp

@@ -63,37 +63,30 @@ class TrackParamsEstimationAlgorithm final : public IAlgorithm {
     /// Output prototrack collection - only tracks with successful parameter
     /// estimation are propagated (optional)
     std::string outputProtoTracks;
+
     /// Tracking geometry for surface lookup.
     std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry;
     /// Magnetic field variant.
     std::shared_ptr<const Acts::MagneticFieldProvider> magneticField;
+
     /// The minimum magnetic field to trigger the track parameters estimation
     double bFieldMin = 0.1 * Acts::UnitConstants::T;
-    /// Initial covariance matrix diagonal.
+
+    /// Initial sigmas for the track parameters.
     std::array<double, 6> initialSigmas = {
         1 * Acts::UnitConstants::mm,
         1 * Acts::UnitConstants::mm,
         1 * Acts::UnitConstants::degree,
         1 * Acts::UnitConstants::degree,
         0 * Acts::UnitConstants::e / Acts::UnitConstants::GeV,
         1 * Acts::UnitConstants::ns};
-    /// Initial q/p coefficient covariance matrix diagonal.
-    std::array<double, 6> initialSimgaQoverPCoefficients = {
-        0 * Acts::UnitConstants::mm /
-            (Acts::UnitConstants::e * Acts::UnitConstants::GeV),
-        0 * Acts::UnitConstants::mm /
-            (Acts::UnitConstants::e * Acts::UnitConstants::GeV),
-        0 * Acts::UnitConstants::degree /
-            (Acts::UnitConstants::e * Acts::UnitConstants::GeV),
-        0 * Acts::UnitConstants::degree /
-            (Acts::UnitConstants::e * Acts::UnitConstants::GeV),
-        0.1,
-        0 * Acts::UnitConstants::ns /
-            (Acts::UnitConstants::e * Acts::UnitConstants::GeV)};
+    /// Relative pt resolution used for the initial sigma of q/p.
+    double initialSigmaPtRel = 0.1;
     /// Inflate initial covariance.
     std::array<double, 6> initialVarInflation = {1., 1., 1., 1., 1., 1.};
     /// Inflate time covariance if no time measurement is available.
     double noTimeVarInflation = 100.;
+
     /// Particle hypothesis.
     Acts::ParticleHypothesis particleHypothesis =
         Acts::ParticleHypothesis::pion();

Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp

@@ -44,22 +44,32 @@ Acts::BoundSquareMatrix makeInitialCovariance(
 
   for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) {
     double sigma = config.initialSigmas[i];
+    double variance = sigma * sigma;
 
-    // Add momentum dependent uncertainties
-    sigma +=
-        config.initialSimgaQoverPCoefficients[i] * params[Acts::eBoundQOverP];
+    if (i == Acts::eBoundQOverP) {
+      // note that we rely on the fact that sigma theta is already computed
+      double varianceTheta = result(Acts::eBoundTheta, Acts::eBoundTheta);
 
-    double var = sigma * sigma;
+      // transverse momentum contribution
+      variance +=
+          std::pow(config.initialSigmaPtRel * params[Acts::eBoundQOverP], 2);
+
+      // theta contribution
+      variance +=
+          varianceTheta * std::pow(params[Acts::eBoundQOverP] *
+                                       std::tan(params[Acts::eBoundTheta]),
+                                   2);
+    }
 
     // Inflate the time uncertainty if no time measurement is available
     if (i == Acts::eBoundTime && !sp.t().has_value()) {
-      var *= config.noTimeVarInflation;
+      variance *= config.noTimeVarInflation;
     }
 
     // Inflate the initial covariance
-    var *= config.initialVarInflation[i];
+    variance *= config.initialVarInflation[i];
 
-    result(i, i) = var;
+    result(i, i) = variance;
   }
 
   return result;

Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -76,7 +76,7 @@ ActsExamples::ProcessCode ActsExamples::ParticleSmearing::execute(
       const auto theta = Acts::VectorHelpers::theta(particle.direction());
       const auto pt = particle.transverseMomentum();
       const auto p = particle.absoluteMomentum();
-      const auto q = particle.charge();
+      const auto qOverP = particle.qOverP();
       const auto particleHypothesis =
           m_cfg.particleHypothesis.value_or(particle.hypothesis());
 
@@ -87,16 +87,17 @@ ActsExamples::ProcessCode ActsExamples::ParticleSmearing::execute(
       const double sigmaZ0 =
           m_cfg.sigmaZ0 +
           m_cfg.sigmaZ0PtA * std::exp(-1.0 * std::abs(m_cfg.sigmaZ0PtB) * pt);
-      const double sigmaP = m_cfg.sigmaPRel * p;
-      // var(q/p) = (d(1/p)/dp)² * var(p) = (-1/p²)² * var(p)
-      const double sigmaQOverP = sigmaP / (p * p);
       // shortcuts for other resolutions
       const double sigmaT0 = m_cfg.sigmaT0;
       const double sigmaPhi = m_cfg.sigmaPhi;
       const double sigmaTheta = m_cfg.sigmaTheta;
+      const double sigmaQOverP =
+          std::sqrt(std::pow(m_cfg.sigmaPtRel * qOverP, 2) +
+                    std::pow(sigmaTheta * (qOverP * std::tan(theta)), 2));
 
       Acts::BoundVector params = Acts::BoundVector::Zero();
       // smear the position/time
+      // note that we smear d0 and z0 in the perigee frame
       params[Acts::eBoundLoc0] = sigmaD0 * stdNormal(rng);
       params[Acts::eBoundLoc1] = sigmaZ0 * stdNormal(rng);
       params[Acts::eBoundTime] = time + sigmaT0 * stdNormal(rng);
@@ -105,10 +106,8 @@ ActsExamples::ProcessCode ActsExamples::ParticleSmearing::execute(
           phi + sigmaPhi * stdNormal(rng), theta + sigmaTheta * stdNormal(rng));
       params[Acts::eBoundPhi] = newPhi;
       params[Acts::eBoundTheta] = newTheta;
-      // compute smeared absolute momentum vector
-      const double newP = std::max(0.0, p + sigmaP * stdNormal(rng));
-      const double qOverP = particleHypothesis.qOverP(newP, q);
-      params[Acts::eBoundQOverP] = qOverP;
+      // compute smeared q/p
+      params[Acts::eBoundQOverP] = qOverP + sigmaQOverP * stdNormal(rng);
 
       ACTS_VERBOSE("Smearing particle (pos, time, phi, theta, q/p):");
       ACTS_VERBOSE(" from: " << particle.position().transpose() << ", " << time
@@ -130,25 +129,45 @@ ActsExamples::ProcessCode ActsExamples::ParticleSmearing::execute(
       if (m_cfg.initialSigmas) {
         // use the initial sigmas if set
         for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) {
-          cov(i, i) = m_cfg.initialVarInflation[i] * (*m_cfg.initialSigmas)[i] *
-                      (*m_cfg.initialSigmas)[i];
+          double sigma = (*m_cfg.initialSigmas)[i];
+          double variance = sigma * sigma;
+
+          if (i == Acts::eBoundQOverP) {
+            // note that we rely on the fact that sigma theta is already
+            // computed
+            double varianceTheta = cov(Acts::eBoundTheta, Acts::eBoundTheta);
+
+            // transverse momentum contribution
+            variance += std::pow(
+                m_cfg.initialSigmaPtRel * params[Acts::eBoundQOverP], 2);
+
+            // theta contribution
+            variance += varianceTheta *
+                        std::pow(params[Acts::eBoundQOverP] *
+                                     std::tan(params[Acts::eBoundTheta]),
+                                 2);
+          }
+
+          // Inflate the initial covariance
+          variance *= m_cfg.initialVarInflation[i];
+
+          cov(i, i) = variance;
         }
       } else {
         // otherwise use the smearing sigmas
-        cov(Acts::eBoundLoc0, Acts::eBoundLoc0) =
-            m_cfg.initialVarInflation[Acts::eBoundLoc0] * sigmaD0 * sigmaD0;
-        cov(Acts::eBoundLoc1, Acts::eBoundLoc1) =
-            m_cfg.initialVarInflation[Acts::eBoundLoc1] * sigmaZ0 * sigmaZ0;
-        cov(Acts::eBoundTime, Acts::eBoundTime) =
-            m_cfg.initialVarInflation[Acts::eBoundTime] * sigmaT0 * sigmaT0;
-        cov(Acts::eBoundPhi, Acts::eBoundPhi) =
-            m_cfg.initialVarInflation[Acts::eBoundPhi] * sigmaPhi * sigmaPhi;
-        cov(Acts::eBoundTheta, Acts::eBoundTheta) =
-            m_cfg.initialVarInflation[Acts::eBoundTheta] * sigmaTheta *
-            sigmaTheta;
-        cov(Acts::eBoundQOverP, Acts::eBoundQOverP) =
-            m_cfg.initialVarInflation[Acts::eBoundQOverP] * sigmaQOverP *
-            sigmaQOverP;
+
+        Acts::BoundVector sigmas = Acts::BoundVector(
+            {sigmaD0, sigmaZ0, sigmaPhi, sigmaTheta, sigmaQOverP, sigmaT0});
+
+        for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) {
+          double sigma = sigmas[i];
+          double variance = sigma * sigma;
+
+          // Inflate the initial covariance
+          variance *= m_cfg.initialVarInflation[i];
+
+          cov(i, i) = variance;
+        }
       }
 
       parameters.emplace_back(perigee, params, cov, particleHypothesis);

Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.hpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2024 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -41,6 +41,11 @@ class ParticleSmearing final : public IAlgorithm {
     std::string inputParticles;
     /// Output smeared tracks parameters collection.
     std::string outputTrackParameters;
+
+    /// Random numbers service.
+    std::shared_ptr<const RandomNumbers> randomNumbers = nullptr;
+
+    // Smearing parameters
     /// Constant term of the d0 resolution.
     double sigmaD0 = 20 * Acts::UnitConstants::um;
     /// Pt-dependent d0 resolution of the form sigma_d0 = A*exp(-1.*abs(B)*pt).
@@ -57,18 +62,17 @@ class ParticleSmearing final : public IAlgorithm {
     double sigmaPhi = 1 * Acts::UnitConstants::degree;
     /// Theta angular resolution.
     double sigmaTheta = 1 * Acts::UnitConstants::degree;
-    /// Relative momentum resolution.
-    double sigmaPRel = 0.05;
-    /// Optional. Initial covariance matrix diagonal. Overwrites the default if
-    /// set.
-    std::optional<std::array<double, 6>> initialSigmas = std::array<double, 6>{
-        1 * Acts::UnitConstants::mm,     1 * Acts::UnitConstants::mm,
-        1 * Acts::UnitConstants::degree, 1 * Acts::UnitConstants::degree,
-        0.1 / Acts::UnitConstants::GeV,  1 * Acts::UnitConstants::ns};
+    /// Relative transverse momentum resolution.
+    double sigmaPtRel = 0.05;
+
+    /// Optional. Initial sigmas for the track parameters which overwrites the
+    /// smearing params if set.
+    std::optional<std::array<double, 6>> initialSigmas;
+    /// Relative pt resolution used for the initial sigma of q/p.
+    double initialSigmaPtRel = 0.1;
     /// Inflate the initial covariance matrix
-    std::array<double, 6> initialVarInflation = {1., 1., 1., 1., 1., 1.};
-    /// Random numbers service.
-    std::shared_ptr<const RandomNumbers> randomNumbers = nullptr;
+    std::array<double, 6> initialVarInflation = {1e4, 1e4, 1e4, 1e4, 1e4, 1e4};
+
     /// Optional particle hypothesis override.
     std::optional<Acts::ParticleHypothesis> particleHypothesis = std::nullopt;
   };

Examples/Io/Root/src/RootTrackSummaryWriter.cpp

@@ -425,7 +425,11 @@ ProcessCode RootTrackSummaryWriter::writeT(const AlgorithmContext& ctx,
              param[Acts::eBoundTime] - t_time};
 
       for (unsigned int i = 0; i < Acts::eBoundSize; ++i) {
-        pull[i] = res[i] / error[i];  // MARK: fpeMask(FLTINV, 1, #2284)
+        // MARK: fpeMaskBegin(FLTDIV, 1, #2348)
+        // MARK: fpeMaskBegin(FLTINV, 1, #2348)
+        pull[i] = res[i] / error[i];
+        // MARK: fpeMaskEnd(FLTINV)
+        // MARK: fpeMaskEnd(FLTDIV)
       }
     }