refactor: Rework random parameter and covariance generation (acts-pro…

…ject#3777)

- Reworks random bound and free track parameter generation
  - This should guarantee to obtain valid parameters acts-project#3756
- Split general random parameter and covariance generation
- Also, don't use covariance to draw parameters in general random generation

Core/include/Acts/EventData/detail/GenerateParameters.hpp

@@ -8,33 +8,48 @@
 
 #pragma once
 
+#include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
-#include "Acts/Utilities/detail/periodic.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Utilities/AngleHelpers.hpp"
 
 #include <cmath>
+#include <numbers>
 #include <random>
 #include <utility>
 
 namespace Acts::detail::Test {
 
-/// Generate a random parameters vector and covariance matrix.
-///
-/// @return std:::pair<ParametersVector, CovarianceMatrix>
 template <typename scalar_t, std::size_t kSize, typename generator_t>
-inline auto generateParametersCovariance(generator_t& rng)
-    -> std::pair<Eigen::Matrix<scalar_t, kSize, 1>,
-                 Eigen::Matrix<scalar_t, kSize, kSize>> {
+inline auto generateParameters(generator_t& rng)
+    -> Eigen::Matrix<scalar_t, kSize, 1> {
+  using Scalar = scalar_t;
+  using ParametersVector = Eigen::Matrix<scalar_t, kSize, 1>;
+
+  std::normal_distribution<Scalar> standardNormal(0, 1);
+
+  ParametersVector params;
+  for (auto i = 0u; i < kSize; ++i) {
+    params[i] = standardNormal(rng);
+  }
+
+  return params;
+}
+
+template <typename scalar_t, std::size_t kSize, typename generator_t>
+inline auto generateCovariance(generator_t& rng)
+    -> Eigen::Matrix<scalar_t, kSize, kSize> {
   using Scalar = scalar_t;
   using ParametersVector = Eigen::Matrix<scalar_t, kSize, 1>;
   using CovarianceMatrix = Eigen::Matrix<scalar_t, kSize, kSize>;
 
-  std::normal_distribution<Scalar> distNormal(0, 1);
+  std::normal_distribution<Scalar> standardNormal(0, 1);
   std::uniform_real_distribution<Scalar> distCorr(-1, 1);
 
   // generate standard deviations
   ParametersVector stddev;
   for (auto i = 0u; i < kSize; ++i) {
-    stddev[i] = std::abs(distNormal(rng));
+    stddev[i] = std::abs(standardNormal(rng));
   }
   // generate correlation matrix
   CovarianceMatrix corr;
@@ -48,32 +63,216 @@ inline auto generateParametersCovariance(generator_t& rng)
   // construct the covariance matrix
   CovarianceMatrix cov = stddev.asDiagonal() * corr * stddev.asDiagonal();
 
-  // generate random parameters
-  // this is ignoring the correlations; since this does not need to generate
-  // credible data, this should be fine.
-  ParametersVector params;
-  for (auto i = 0u; i < kSize; ++i) {
-    params[i] = stddev[i] * distNormal(rng);
+  return cov;
+}
+
+/// Generate a random parameters vector and covariance matrix.
+///
+/// @return std:::pair<ParametersVector, CovarianceMatrix>
+template <typename scalar_t, std::size_t kSize, typename generator_t>
+inline auto generateParametersCovariance(generator_t& rng)
+    -> std::pair<Eigen::Matrix<scalar_t, kSize, 1>,
+                 Eigen::Matrix<scalar_t, kSize, kSize>> {
+  auto params = generateParameters<scalar_t, kSize, generator_t>(rng);
+  auto cov = generateCovariance<scalar_t, kSize, generator_t>(rng);
+  return {params, cov};
+}
+
+struct GenerateBoundDirectionOptions {
+  /// Low, high (exclusive) for the transverse direction angle.
+  double phiMin = -std::numbers::pi;
+  double phiMax = std::numbers::pi;
+
+  /// Low, high (inclusive) for  the longitudinal direction angle.
+  ///
+  /// This intentionally uses theta instead of eta so it can represent the
+  /// full direction space with finite values.
+  ///
+  /// @note This is the standard generation, for detector performance
+  /// classification, where a flat distribution in eta can be useful,
+  /// this can be set by the etaUniform flag;
+  ///
+  double thetaMin = AngleHelpers::thetaFromEta(6.0);
+  double thetaMax = AngleHelpers::thetaFromEta(-6.0);
+
+  bool etaUniform = true;
+};
+
+template <typename generator_t>
+inline std::pair<double, double> generateBoundDirection(
+    generator_t& rng, const GenerateBoundDirectionOptions& options) {
+  using UniformReal = std::uniform_real_distribution<double>;
+
+  // since we want to draw the direction uniform on the unit sphere, we must
+  // draw from cos(theta) instead of theta. see e.g.
+  // https://mathworld.wolfram.com/SpherePointPicking.html
+  double cosThetaMin = std::cos(options.thetaMin);
+  // ensure upper bound is included. see e.g.
+  // https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution
+  double cosThetaMax = std::nextafter(std::cos(options.thetaMax),
+                                      std::numeric_limits<double>::max());
+
+  // in case we force uniform eta generation
+  double etaMin = Acts::AngleHelpers::etaFromTheta(options.thetaMin);
+  double etaMax = Acts::AngleHelpers::etaFromTheta(options.thetaMax);
+
+  UniformReal phiDist(options.phiMin, options.phiMax);
+  UniformReal cosThetaDist(cosThetaMin, cosThetaMax);
+  UniformReal etaDist(etaMin, etaMax);
+
+  // draw parameters
+  double phi = phiDist(rng);
+
+  double theta = 0;
+  if (!options.etaUniform) {
+    const double cosTheta = cosThetaDist(rng);
+    theta = std::acos(cosTheta);
+  } else {
+    const double eta = etaDist(rng);
+    theta = AngleHelpers::thetaFromEta(eta);
   }
 
+  return {phi, theta};
+}
+
+struct GenerateQoverPOptions {
+  /// Low, high (exclusive) for absolute/transverse momentum.
+  double pMin = 1 * UnitConstants::GeV;
+  double pMax = 100 * UnitConstants::GeV;
+
+  /// Indicate if the momentum referse to transverse momentum
+  bool pTransverse = true;
+
+  /// Charge of the parameters.
+  double charge = 1;
+
+  /// Randomize the charge and flip the PDG particle number sign accordingly.
+  bool randomizeCharge = true;
+};
+
+template <typename generator_t>
+inline double generateQoverP(generator_t& rng,
+                             const GenerateQoverPOptions& options,
+                             double theta) {
+  using UniformIndex = std::uniform_int_distribution<std::uint8_t>;
+  using UniformReal = std::uniform_real_distribution<double>;
+
+  // choose between particle/anti-particle if requested
+  // the upper limit of the distribution is inclusive
+  UniformIndex particleTypeChoice(0u, options.randomizeCharge ? 1u : 0u);
+  // (anti-)particle choice is one random draw but defines two properties
+  const double qChoices[] = {
+      options.charge,
+      -options.charge,
+  };
+  UniformReal pDist(options.pMin, options.pMax);
+
+  // draw parameters
+  const std::uint8_t type = particleTypeChoice(rng);
+  const double q = qChoices[type];
+
+  const double p =
+      pDist(rng) * (options.pTransverse ? 1. / std::sin(theta) : 1.);
+  const double qOverP = (q != 0) ? q / p : 1 / p;
+
+  return qOverP;
+}
+
+struct GenerateBoundParametersOptions {
+  struct {
+    double loc0Mean = 0 * UnitConstants::mm;
+    double loc0Std = 1 * UnitConstants::mm;
+
+    double loc1Mean = 0 * UnitConstants::mm;
+    double loc1Std = 1 * UnitConstants::mm;
+
+    double timeMean = 0 * UnitConstants::ns;
+    double timeStd = 1 * UnitConstants::ns;
+  } position;
+
+  GenerateBoundDirectionOptions direction;
+
+  GenerateQoverPOptions qOverP;
+};
+
+template <typename generator_t>
+inline BoundVector generateBoundParameters(
+    generator_t& rng, const GenerateBoundParametersOptions& options) {
+  std::normal_distribution<double> standardNormal(0, 1);
+
+  const double loc0 = options.position.loc0Mean +
+                      options.position.loc0Std * standardNormal(rng);
+  const double loc1 = options.position.loc1Mean +
+                      options.position.loc1Std * standardNormal(rng);
+
+  auto [phi, theta] = generateBoundDirection(rng, options.direction);
+  auto qOverP = generateQoverP(rng, options.qOverP, theta);
+
+  const double time = options.position.timeMean +
+                      options.position.timeStd * standardNormal(rng);
+
+  return {loc0, loc1, phi, theta, qOverP, time};
+}
+
+template <typename generator_t>
+inline std::pair<BoundVector, BoundMatrix> generateBoundParametersCovariance(
+    generator_t& rng, const GenerateBoundParametersOptions& options) {
+  auto params = generateBoundParameters(rng, options);
+  auto cov = generateCovariance<ActsScalar, eBoundSize>(rng);
   return {params, cov};
 }
 
-/// Generate a random bound parameters vector and covariance matrix.
+struct GenerateFreeParametersOptions {
+  struct {
+    double xMean = 0 * UnitConstants::mm;
+    double xStd = 1 * UnitConstants::mm;
+
+    double yMean = 0 * UnitConstants::mm;
+    double yStd = 1 * UnitConstants::mm;
+
+    double zMean = 0 * UnitConstants::mm;
+    double zStd = 1 * UnitConstants::mm;
+
+    double timeMean = 0 * UnitConstants::ns;
+    double timeStd = 1 * UnitConstants::ns;
+  } position;
+
+  GenerateBoundDirectionOptions direction;
+
+  GenerateQoverPOptions qOverP;
+};
+
 template <typename generator_t>
-inline auto generateBoundParametersCovariance(generator_t& rng) {
-  auto parCov = generateParametersCovariance<ActsScalar, eBoundSize>(rng);
-  auto [phi, theta] = detail::normalizePhiTheta(parCov.first[eBoundPhi],
-                                                parCov.first[eBoundTheta]);
-  parCov.first[eBoundPhi] = phi;
-  parCov.first[eBoundTheta] = theta;
-  return parCov;
+inline FreeVector generateFreeParameters(
+    generator_t& rng, const GenerateFreeParametersOptions& options) {
+  std::normal_distribution<double> standardNormal(0, 1);
+
+  const double x =
+      options.position.xMean + options.position.xStd * standardNormal(rng);
+  const double y =
+      options.position.yMean + options.position.yStd * standardNormal(rng);
+  const double z =
+      options.position.zMean + options.position.zStd * standardNormal(rng);
+  const double time = options.position.timeMean +
+                      options.position.timeStd * standardNormal(rng);
+
+  auto [phi, theta] = generateBoundDirection(rng, options.direction);
+
+  Vector3 direction = makeDirectionFromPhiTheta(phi, theta);
+
+  auto qOverP = generateQoverP(rng, options.qOverP, theta);
+
+  FreeVector freeParams;
+  freeParams << x, y, z, time, direction, qOverP;
+  return freeParams;
 }
 
-/// Generate a random free parameters vector and covariance matrix.
 template <typename generator_t>
-inline auto generateFreeParametersCovariance(generator_t& rng) {
-  return generateParametersCovariance<ActsScalar, eFreeSize>(rng);
+inline std::pair<FreeVector, FreeMatrix> generateFreeParametersCovariance(
+    generator_t& rng, const GenerateFreeParametersOptions& options) {
+  auto params = generateFreeParameters(rng, options);
+  auto cov = generateCovariance<ActsScalar, eFreeSize>(rng);
+  return {params, cov};
 }
 
 }  // namespace Acts::detail::Test

Core/include/Acts/EventData/detail/MultiTrajectoryTestsCommon.hpp

@@ -331,7 +331,7 @@ class MultiTrajectoryTestsCommon {
     auto tsb = traj.getTrackState(index);
     // then modify one and check that the other was modified as well
     {
-      auto [par, cov] = generateBoundParametersCovariance(rng);
+      auto [par, cov] = generateBoundParametersCovariance(rng, {});
       tsb.predicted() = par;
       tsb.predictedCovariance() = cov;
       BOOST_CHECK_EQUAL(tsa.predicted(), par);
@@ -340,7 +340,7 @@ class MultiTrajectoryTestsCommon {
       BOOST_CHECK_EQUAL(tsb.predictedCovariance(), cov);
     }
     {
-      auto [par, cov] = generateBoundParametersCovariance(rng);
+      auto [par, cov] = generateBoundParametersCovariance(rng, {});
       tsb.filtered() = par;
       tsb.filteredCovariance() = cov;
       BOOST_CHECK_EQUAL(tsa.filtered(), par);
@@ -349,7 +349,7 @@ class MultiTrajectoryTestsCommon {
       BOOST_CHECK_EQUAL(tsb.filteredCovariance(), cov);
     }
     {
-      auto [par, cov] = generateBoundParametersCovariance(rng);
+      auto [par, cov] = generateBoundParametersCovariance(rng, {});
       tsb.smoothed() = par;
       tsb.smoothedCovariance() = cov;
       BOOST_CHECK_EQUAL(tsa.smoothed(), par);
@@ -377,7 +377,7 @@ class MultiTrajectoryTestsCommon {
     }
     {
       // reset measurements w/ full parameters
-      auto [measPar, measCov] = generateBoundParametersCovariance(rng);
+      auto [measPar, measCov] = generateBoundParametersCovariance(rng, {});
       tsb.allocateCalibrated(eBoundSize);
       tsb.template calibrated<eBoundSize>() = measPar;
       tsb.template calibratedCovariance<eBoundSize>() = measCov;
@@ -390,7 +390,7 @@ class MultiTrajectoryTestsCommon {
     }
     {
       // reset only the effective measurements
-      auto [measPar, measCov] = generateBoundParametersCovariance(rng);
+      auto [measPar, measCov] = generateBoundParametersCovariance(rng, {});
       std::size_t nMeasurements = tsb.effectiveCalibrated().rows();
       auto effPar = measPar.head(nMeasurements);
       auto effCov = measCov.topLeftCorner(nMeasurements, nMeasurements);

Core/include/Acts/EventData/detail/TestTrackState.hpp

@@ -65,7 +65,7 @@ struct TestTrackState {
     }
 
     // create track parameters
-    auto [trkPar, trkCov] = generateBoundParametersCovariance(rng);
+    auto [trkPar, trkCov] = generateBoundParametersCovariance(rng, {});
     // trkPar[eBoundPhi] = 45_degree;
     // trkPar[eBoundTheta] = 90_degree;
     // trkPar[eBoundQOverP] = 5.;

Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp

@@ -89,7 +89,7 @@ ParametricParticleGenerator::operator()(RandomEngine& rng) {
       sinTheta = std::sqrt(1 - cosTheta * cosTheta);
     } else {
       const double eta = etaDist(rng);
-      const double theta = 2 * std::atan(std::exp(-eta));
+      const double theta = Acts::AngleHelpers::thetaFromEta(eta);
       sinTheta = std::sin(theta);
       cosTheta = std::cos(theta);
     }

Tests/UnitTests/Examples/EventData/MeasurementTests.cpp

@@ -85,7 +85,8 @@ BOOST_DATA_TEST_CASE(VariableBoundOneEmplace, bd::make(boundIndices), index) {
 BOOST_AUTO_TEST_CASE(VariableBoundAll) {
   MeasurementContainer container;
 
-  auto [params, cov] = generateBoundParametersCovariance(rng);
+  auto [params, cov] =
+      generateParametersCovariance<ActsScalar, eBoundSize>(rng);
 
   FixedBoundMeasurementProxy<eBoundSize> meas =
       container.makeMeasurement<eBoundSize>(geoId);
@@ -106,7 +107,8 @@ BOOST_AUTO_TEST_CASE(VariableBoundAll) {
 BOOST_AUTO_TEST_CASE(VariableBoundAllEmplace) {
   MeasurementContainer container;
 
-  auto [params, cov] = generateBoundParametersCovariance(rng);
+  auto [params, cov] =
+      generateParametersCovariance<ActsScalar, eBoundSize>(rng);
 
   FixedBoundMeasurementProxy<eBoundSize> meas =
       container.emplaceMeasurement<eBoundSize>(
@@ -144,7 +146,8 @@ BOOST_AUTO_TEST_CASE(VariableBoundReassign) {
   BOOST_CHECK(!meas.contains(eBoundQOverP));
 
   // reassign w/ all parameters
-  auto [paramsN, covN] = generateBoundParametersCovariance(rng);
+  auto [paramsN, covN] =
+      generateParametersCovariance<ActsScalar, eBoundSize>(rng);
 
   meas = container.makeMeasurement(eBoundSize, geoId);
   meas.setSubspaceIndices(std::array{eBoundLoc0, eBoundLoc1, eBoundTime,

Tests/UnitTests/Examples/Io/Json/JsonDigitizationConfigTests.cpp

@@ -57,7 +57,7 @@ struct Fixture {
 
     // generate random track parameters
     auto [par, cov] =
-        Acts::detail::Test::generateBoundParametersCovariance(rng);
+        Acts::detail::Test::generateBoundParametersCovariance(rng, {});
     boundParams = par;
 
     freeParams =