Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:10:46

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Utilities/AngleHelpers.hpp"
0015 
0016 #include <cmath>
0017 #include <numbers>
0018 #include <random>
0019 #include <utility>
0020 
0021 namespace Acts::detail::Test {
0022 
0023 template <typename scalar_t, std::size_t kSize, typename generator_t>
0024 inline auto generateParameters(generator_t& rng)
0025     -> Eigen::Matrix<scalar_t, kSize, 1> {
0026   using Scalar = scalar_t;
0027   using ParametersVector = Eigen::Matrix<scalar_t, kSize, 1>;
0028 
0029   std::normal_distribution<Scalar> standardNormal(0, 1);
0030 
0031   ParametersVector params;
0032   for (auto i = 0u; i < kSize; ++i) {
0033     params[i] = standardNormal(rng);
0034   }
0035 
0036   return params;
0037 }
0038 
0039 template <typename scalar_t, std::size_t kSize, typename generator_t>
0040 inline auto generateCovariance(generator_t& rng)
0041     -> Eigen::Matrix<scalar_t, kSize, kSize> {
0042   using Scalar = scalar_t;
0043   using ParametersVector = Eigen::Matrix<scalar_t, kSize, 1>;
0044   using CovarianceMatrix = Eigen::Matrix<scalar_t, kSize, kSize>;
0045 
0046   std::normal_distribution<Scalar> standardNormal(0, 1);
0047   std::uniform_real_distribution<Scalar> distCorr(-1, 1);
0048 
0049   // generate standard deviations
0050   ParametersVector stddev;
0051   for (auto i = 0u; i < kSize; ++i) {
0052     stddev[i] = std::abs(standardNormal(rng));
0053   }
0054   // generate correlation matrix
0055   CovarianceMatrix corr;
0056   for (auto i = 0u; i < kSize; ++i) {
0057     corr(i, i) = 1;
0058     // only need generate the sub-diagonal elements
0059     for (auto j = 0u; j < i; ++j) {
0060       corr(i, j) = corr(j, i) = distCorr(rng);
0061     }
0062   }
0063   // construct the covariance matrix
0064   CovarianceMatrix cov = stddev.asDiagonal() * corr * stddev.asDiagonal();
0065 
0066   return cov;
0067 }
0068 
0069 /// Generate a random parameters vector and covariance matrix.
0070 ///
0071 /// @return std:::pair<ParametersVector, CovarianceMatrix>
0072 template <typename scalar_t, std::size_t kSize, typename generator_t>
0073 inline auto generateParametersCovariance(generator_t& rng)
0074     -> std::pair<Eigen::Matrix<scalar_t, kSize, 1>,
0075                  Eigen::Matrix<scalar_t, kSize, kSize>> {
0076   auto params = generateParameters<scalar_t, kSize, generator_t>(rng);
0077   auto cov = generateCovariance<scalar_t, kSize, generator_t>(rng);
0078   return {params, cov};
0079 }
0080 
0081 struct GenerateBoundDirectionOptions {
0082   /// Low, high (exclusive) for the transverse direction angle.
0083   double phiMin = -std::numbers::pi;
0084   double phiMax = std::numbers::pi;
0085 
0086   /// Low, high (inclusive) for  the longitudinal direction angle.
0087   ///
0088   /// This intentionally uses theta instead of eta so it can represent the
0089   /// full direction space with finite values.
0090   ///
0091   /// @note This is the standard generation, for detector performance
0092   /// classification, where a flat distribution in eta can be useful,
0093   /// this can be set by the etaUniform flag;
0094   ///
0095   double thetaMin = AngleHelpers::thetaFromEta(6.0);
0096   double thetaMax = AngleHelpers::thetaFromEta(-6.0);
0097 
0098   bool etaUniform = true;
0099 };
0100 
0101 template <typename generator_t>
0102 inline std::pair<double, double> generateBoundDirection(
0103     generator_t& rng, const GenerateBoundDirectionOptions& options) {
0104   using UniformReal = std::uniform_real_distribution<double>;
0105   assert(options.thetaMin >= 0.f);
0106   assert(options.thetaMax <= std::numbers::pi);
0107   assert(options.thetaMin <= options.thetaMax);
0108 
0109   // since we want to draw the direction uniform on the unit sphere, we must
0110   // draw from cos(theta) instead of theta. see e.g.
0111   // https://mathworld.wolfram.com/SpherePointPicking.html
0112   // Get cosThetaMin from thetaMax and vice versa, because cos is
0113   // monothonical decreasing between [0, pi]
0114   double cosThetaMin = std::cos(options.thetaMax);
0115   // ensure upper bound is included. see e.g.
0116   // https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution
0117   double cosThetaMax = std::nextafter(std::cos(options.thetaMin),
0118                                       std::numeric_limits<double>::max());
0119 
0120   // in case we force uniform eta generation
0121   double etaMin = Acts::AngleHelpers::etaFromTheta(options.thetaMax);
0122   double etaMax = Acts::AngleHelpers::etaFromTheta(options.thetaMin);
0123 
0124   UniformReal phiDist(options.phiMin, options.phiMax);
0125   UniformReal cosThetaDist(cosThetaMin, cosThetaMax);
0126   UniformReal etaDist(etaMin, etaMax);
0127 
0128   // draw parameters
0129   double phi = phiDist(rng);
0130 
0131   double theta = 0;
0132   if (!options.etaUniform) {
0133     const double cosTheta = cosThetaDist(rng);
0134     theta = std::acos(cosTheta);
0135   } else {
0136     const double eta = etaDist(rng);
0137     theta = AngleHelpers::thetaFromEta(eta);
0138   }
0139 
0140   return {phi, theta};
0141 }
0142 
0143 struct GenerateQoverPOptions {
0144   /// Low, high (exclusive) for absolute/transverse momentum.
0145   double pMin = 1 * UnitConstants::GeV;
0146   double pMax = 100 * UnitConstants::GeV;
0147 
0148   /// Indicate if the momentum referse to transverse momentum
0149   bool pTransverse = true;
0150 
0151   /// Indicate if the momentum should be uniformly distributed in log space.
0152   bool pLogUniform = false;
0153 
0154   /// Charge of the parameters.
0155   double charge = 1;
0156 
0157   /// Randomize the charge and flip the PDG particle number sign accordingly.
0158   bool randomizeCharge = true;
0159 };
0160 
0161 template <typename generator_t>
0162 inline double generateQoverP(generator_t& rng,
0163                              const GenerateQoverPOptions& options,
0164                              double theta) {
0165   using UniformIndex = std::uniform_int_distribution<std::uint8_t>;
0166   using UniformReal = std::uniform_real_distribution<double>;
0167 
0168   auto drawP = [&options](generator_t& rng_, double theta_) -> double {
0169     const double pTransverseScaling =
0170         options.pTransverse ? 1. / std::sin(theta_) : 1.;
0171 
0172     if (options.pLogUniform) {
0173       UniformReal pLogDist(std::log(options.pMin), std::log(options.pMax));
0174       return std::exp(pLogDist(rng_)) * pTransverseScaling;
0175     }
0176 
0177     UniformReal pDist(options.pMin, options.pMax);
0178     return pDist(rng_) * pTransverseScaling;
0179   };
0180 
0181   // choose between particle/anti-particle if requested
0182   // the upper limit of the distribution is inclusive
0183   UniformIndex particleTypeChoice(0u, options.randomizeCharge ? 1u : 0u);
0184   // (anti-)particle choice is one random draw but defines two properties
0185   const double qChoices[] = {
0186       options.charge,
0187       -options.charge,
0188   };
0189 
0190   // draw parameters
0191   const std::uint8_t type = particleTypeChoice(rng);
0192   const double q = qChoices[type];
0193 
0194   const double p = drawP(rng, theta);
0195   const double qOverP = (q != 0) ? q / p : 1 / p;
0196 
0197   return qOverP;
0198 }
0199 
0200 struct GenerateBoundParametersOptions {
0201   struct {
0202     double loc0Mean = 0 * UnitConstants::mm;
0203     double loc0Std = 1 * UnitConstants::mm;
0204 
0205     double loc1Mean = 0 * UnitConstants::mm;
0206     double loc1Std = 1 * UnitConstants::mm;
0207 
0208     double timeMean = 0 * UnitConstants::ns;
0209     double timeStd = 1 * UnitConstants::ns;
0210   } position;
0211 
0212   GenerateBoundDirectionOptions direction;
0213 
0214   GenerateQoverPOptions qOverP;
0215 };
0216 
0217 inline BoundVector someBoundParametersA() {
0218   return {0.0, 0.0, 0.0, std::numbers::pi / 2, 1.0, 0.0};
0219 }
0220 
0221 template <typename generator_t>
0222 inline BoundVector generateBoundParameters(
0223     generator_t& rng, const GenerateBoundParametersOptions& options) {
0224   std::normal_distribution<double> standardNormal(0, 1);
0225 
0226   const double loc0 = options.position.loc0Mean +
0227                       options.position.loc0Std * standardNormal(rng);
0228   const double loc1 = options.position.loc1Mean +
0229                       options.position.loc1Std * standardNormal(rng);
0230 
0231   auto [phi, theta] = generateBoundDirection(rng, options.direction);
0232   auto qOverP = generateQoverP(rng, options.qOverP, theta);
0233 
0234   const double time = options.position.timeMean +
0235                       options.position.timeStd * standardNormal(rng);
0236 
0237   return {loc0, loc1, phi, theta, qOverP, time};
0238 }
0239 
0240 template <typename generator_t>
0241 inline std::pair<BoundVector, BoundMatrix> generateBoundParametersCovariance(
0242     generator_t& rng, const GenerateBoundParametersOptions& options) {
0243   auto params = generateBoundParameters(rng, options);
0244   auto cov = generateCovariance<double, eBoundSize>(rng);
0245   return {params, cov};
0246 }
0247 
0248 struct GenerateFreeParametersOptions {
0249   struct {
0250     double xMean = 0 * UnitConstants::mm;
0251     double xStd = 1 * UnitConstants::mm;
0252 
0253     double yMean = 0 * UnitConstants::mm;
0254     double yStd = 1 * UnitConstants::mm;
0255 
0256     double zMean = 0 * UnitConstants::mm;
0257     double zStd = 1 * UnitConstants::mm;
0258 
0259     double timeMean = 0 * UnitConstants::ns;
0260     double timeStd = 1 * UnitConstants::ns;
0261   } position;
0262 
0263   GenerateBoundDirectionOptions direction;
0264 
0265   GenerateQoverPOptions qOverP;
0266 };
0267 
0268 inline FreeVector someFreeParametersA() {
0269   return {0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0};
0270 }
0271 
0272 template <typename generator_t>
0273 inline FreeVector generateFreeParameters(
0274     generator_t& rng, const GenerateFreeParametersOptions& options) {
0275   std::normal_distribution<double> standardNormal(0, 1);
0276 
0277   const double x =
0278       options.position.xMean + options.position.xStd * standardNormal(rng);
0279   const double y =
0280       options.position.yMean + options.position.yStd * standardNormal(rng);
0281   const double z =
0282       options.position.zMean + options.position.zStd * standardNormal(rng);
0283   const double time = options.position.timeMean +
0284                       options.position.timeStd * standardNormal(rng);
0285 
0286   auto [phi, theta] = generateBoundDirection(rng, options.direction);
0287 
0288   Vector3 direction = makeDirectionFromPhiTheta(phi, theta);
0289 
0290   auto qOverP = generateQoverP(rng, options.qOverP, theta);
0291 
0292   FreeVector freeParams;
0293   freeParams << x, y, z, time, direction, qOverP;
0294   return freeParams;
0295 }
0296 
0297 template <typename generator_t>
0298 inline std::pair<FreeVector, FreeMatrix> generateFreeParametersCovariance(
0299     generator_t& rng, const GenerateFreeParametersOptions& options) {
0300   auto params = generateFreeParameters(rng, options);
0301   auto cov = generateCovariance<double, eFreeSize>(rng);
0302   return {params, cov};
0303 }
0304 
0305 }  // namespace Acts::detail::Test