File indexing completed on 2025-01-18 09:10:46
0001
0002
0003
0004
0005
0006
0007
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
0050 ParametersVector stddev;
0051 for (auto i = 0u; i < kSize; ++i) {
0052 stddev[i] = std::abs(standardNormal(rng));
0053 }
0054
0055 CovarianceMatrix corr;
0056 for (auto i = 0u; i < kSize; ++i) {
0057 corr(i, i) = 1;
0058
0059 for (auto j = 0u; j < i; ++j) {
0060 corr(i, j) = corr(j, i) = distCorr(rng);
0061 }
0062 }
0063
0064 CovarianceMatrix cov = stddev.asDiagonal() * corr * stddev.asDiagonal();
0065
0066 return cov;
0067 }
0068
0069
0070
0071
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
0083 double phiMin = -std::numbers::pi;
0084 double phiMax = std::numbers::pi;
0085
0086
0087
0088
0089
0090
0091
0092
0093
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
0110
0111
0112
0113
0114 double cosThetaMin = std::cos(options.thetaMax);
0115
0116
0117 double cosThetaMax = std::nextafter(std::cos(options.thetaMin),
0118 std::numeric_limits<double>::max());
0119
0120
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
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
0145 double pMin = 1 * UnitConstants::GeV;
0146 double pMax = 100 * UnitConstants::GeV;
0147
0148
0149 bool pTransverse = true;
0150
0151
0152 bool pLogUniform = false;
0153
0154
0155 double charge = 1;
0156
0157
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
0182
0183 UniformIndex particleTypeChoice(0u, options.randomizeCharge ? 1u : 0u);
0184
0185 const double qChoices[] = {
0186 options.charge,
0187 -options.charge,
0188 };
0189
0190
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 }