test/utils/landau_distribution.hpp

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 // Project include(s).
0012 #include "detray/definitions/algebra.hpp"
0013 #include "detray/definitions/containers.hpp"
0014 #include "detray/definitions/math.hpp"
0015
0016 // System include(s).
0017 #include <limits>
0018 #include <random>
0019
0020 namespace detray {
0021
0022 /// Landau distribution based on CERNLIB G110
0023 /// Computer Physics Communications 31(1984) 97—111)
0024 template <concepts::scalar scalar_t>
0025 class landau_distribution {
0026  public:
0027   using scalar_type = scalar_t;
0028
0029   /// Generate a random number following a Landau distribution
0030   /// with location parameter mu and scale parameter sigma:
0031   ///      Landau( (x-mu)/sigma )
0032   /// Note that mu is not the mpv(most probable value) of the Landa
0033   /// distribution and sigma is not the standard deviation of the distribution
0034   /// which is not defined. For mu  =0 and sigma=1, the mpv = -0.22278
0035   ///
0036   /// The Landau random number generation is implemented using the
0037   /// function landau_quantile(x,sigma), which provides
0038   /// the inverse of the landau cumulative distribution.
0039   /// landau_quantile has been converted from CERNLIB ranlan(G110).
0040   ///
0041   /// Reference[1]: ROOT TRandom.cxx
0042   /// Reference[2]: ACTS LandauDistribution.cxx
0043   template <typename generator_t>
0044   scalar_type operator()(generator_t &generator, const scalar_type location,
0045                          const scalar_type scale) const {
0046     const auto z = std::uniform_real_distribution<scalar_type>()(generator);
0047     // LANDAU quantile : algorithm from CERNLIB G110 ranlan
0048     // Converted by Rene Brun from CERNLIB routine ranlan(G110),
0049     // Moved and adapted to QuantFuncMathCore by B. List 29.4.2010
0050     return location + scale * quantile(z);
0051   }
0052
0053  private:
0054   scalar_type quantile(const scalar_type z) const {
0055     static const darray<double, 982> f{
0056         0.,        0.,        0.,        0.,        0.,        -2.244733,
0057         -2.204365, -2.168163, -2.135219, -2.104898, -2.076740, -2.050397,
0058         -2.025605, -2.002150, -1.979866, -1.958612, -1.938275, -1.918760,
0059         -1.899984, -1.881879, -1.864385, -1.847451, -1.831030, -1.815083,
0060         -1.799574, -1.784473, -1.769751, -1.755383, -1.741346, -1.727620,
0061         -1.714187, -1.701029, -1.688130, -1.675477, -1.663057, -1.650858,
0062         -1.638868, -1.627078, -1.615477, -1.604058, -1.592811, -1.581729,
0063         -1.570806, -1.560034, -1.549407, -1.538919, -1.528565, -1.518339,
0064         -1.508237, -1.498254, -1.488386, -1.478628, -1.468976, -1.459428,
0065         -1.449979, -1.440626, -1.431365, -1.422195, -1.413111, -1.404112,
0066         -1.395194, -1.386356, -1.377594, -1.368906, -1.360291, -1.351746,
0067         -1.343269, -1.334859, -1.326512, -1.318229, -1.310006, -1.301843,
0068         -1.293737, -1.285688, -1.277693, -1.269752, -1.261863, -1.254024,
0069         -1.246235, -1.238494, -1.230800, -1.223153, -1.215550, -1.207990,
0070         -1.200474, -1.192999, -1.185566, -1.178172, -1.170817, -1.163500,
0071         -1.156220, -1.148977, -1.141770, -1.134598, -1.127459, -1.120354,
0072         -1.113282, -1.106242, -1.099233, -1.092255, -1.085306, -1.078388,
0073         -1.071498, -1.064636, -1.057802, -1.050996, -1.044215, -1.037461,
0074         -1.030733, -1.024029, -1.017350, -1.010695, -1.004064, -.997456,
0075         -.990871,  -.984308,  -.977767,  -.971247,  -.964749,  -.958271,
0076         -.951813,  -.945375,  -.938957,  -.932558,  -.926178,  -.919816,
0077         -.913472,  -.907146,  -.900838,  -.894547,  -.888272,  -.882014,
0078         -.875773,  -.869547,  -.863337,  -.857142,  -.850963,  -.844798,
0079         -.838648,  -.832512,  -.826390,  -.820282,  -.814187,  -.808106,
0080         -.802038,  -.795982,  -.789940,  -.783909,  -.777891,  -.771884,
0081         -.765889,  -.759906,  -.753934,  -.747973,  -.742023,  -.736084,
0082         -.730155,  -.724237,  -.718328,  -.712429,  -.706541,  -.700661,
0083         -.694791,  -.688931,  -.683079,  -.677236,  -.671402,  -.665576,
0084         -.659759,  -.653950,  -.648149,  -.642356,  -.636570,  -.630793,
0085         -.625022,  -.619259,  -.613503,  -.607754,  -.602012,  -.596276,
0086         -.590548,  -.584825,  -.579109,  -.573399,  -.567695,  -.561997,
0087         -.556305,  -.550618,  -.544937,  -.539262,  -.533592,  -.527926,
0088         -.522266,  -.516611,  -.510961,  -.505315,  -.499674,  -.494037,
0089         -.488405,  -.482777,  -.477153,  -.471533,  -.465917,  -.460305,
0090         -.454697,  -.449092,  -.443491,  -.437893,  -.432299,  -.426707,
0091         -.421119,  -.415534,  -.409951,  -.404372,  -.398795,  -.393221,
0092         -.387649,  -.382080,  -.376513,  -.370949,  -.365387,  -.359826,
0093         -.354268,  -.348712,  -.343157,  -.337604,  -.332053,  -.326503,
0094         -.320955,  -.315408,  -.309863,  -.304318,  -.298775,  -.293233,
0095         -.287692,  -.282152,  -.276613,  -.271074,  -.265536,  -.259999,
0096         -.254462,  -.248926,  -.243389,  -.237854,  -.232318,  -.226783,
0097         -.221247,  -.215712,  -.210176,  -.204641,  -.199105,  -.193568,
0098         -.188032,  -.182495,  -.176957,  -.171419,  -.165880,  -.160341,
0099         -.154800,  -.149259,  -.143717,  -.138173,  -.132629,  -.127083,
0100         -.121537,  -.115989,  -.110439,  -.104889,  -.099336,  -.093782,
0101         -.088227,  -.082670,  -.077111,  -.071550,  -.065987,  -.060423,
0102         -.054856,  -.049288,  -.043717,  -.038144,  -.032569,  -.026991,
0103         -.021411,  -.015828,  -.010243,  -.004656,  .000934,   .006527,
0104         .012123,   .017722,   .023323,   .028928,   .034535,   .040146,
0105         .045759,   .051376,   .056997,   .062620,   .068247,   .073877,
0106         .079511,   .085149,   .090790,   .096435,   .102083,   .107736,
0107         .113392,   .119052,   .124716,   .130385,   .136057,   .141734,
0108         .147414,   .153100,   .158789,   .164483,   .170181,   .175884,
0109         .181592,   .187304,   .193021,   .198743,   .204469,   .210201,
0110         .215937,   .221678,   .227425,   .233177,   .238933,   .244696,
0111         .250463,   .256236,   .262014,   .267798,   .273587,   .279382,
0112         .285183,   .290989,   .296801,   .302619,   .308443,   .314273,
0113         .320109,   .325951,   .331799,   .337654,   .343515,   .349382,
0114         .355255,   .361135,   .367022,   .372915,   .378815,   .384721,
0115         .390634,   .396554,   .402481,   .408415,   .414356,   .420304,
0116         .426260,   .432222,   .438192,   .444169,   .450153,   .456145,
0117         .462144,   .468151,   .474166,   .480188,   .486218,   .492256,
0118         .498302,   .504356,   .510418,   .516488,   .522566,   .528653,
0119         .534747,   .540850,   .546962,   .553082,   .559210,   .565347,
0120         .571493,   .577648,   .583811,   .589983,   .596164,   .602355,
0121         .608554,   .614762,   .620980,   .627207,   .633444,   .639689,
0122         .645945,   .652210,   .658484,   .664768,   .671062,   .677366,
0123         .683680,   .690004,   .696338,   .702682,   .709036,   .715400,
0124         .721775,   .728160,   .734556,   .740963,   .747379,   .753807,
0125         .760246,   .766695,   .773155,   .779627,   .786109,   .792603,
0126         .799107,   .805624,   .812151,   .818690,   .825241,   .831803,
0127         .838377,   .844962,   .851560,   .858170,   .864791,   .871425,
0128         .878071,   .884729,   .891399,   .898082,   .904778,   .911486,
0129         .918206,   .924940,   .931686,   .938446,   .945218,   .952003,
0130         .958802,   .965614,   .972439,   .979278,   .986130,   .992996,
0131         .999875,   1.006769,  1.013676,  1.020597,  1.027533,  1.034482,
0132         1.041446,  1.048424,  1.055417,  1.062424,  1.069446,  1.076482,
0133         1.083534,  1.090600,  1.097681,  1.104778,  1.111889,  1.119016,
0134         1.126159,  1.133316,  1.140490,  1.147679,  1.154884,  1.162105,
0135         1.169342,  1.176595,  1.183864,  1.191149,  1.198451,  1.205770,
0136         1.213105,  1.220457,  1.227826,  1.235211,  1.242614,  1.250034,
0137         1.257471,  1.264926,  1.272398,  1.279888,  1.287395,  1.294921,
0138         1.302464,  1.310026,  1.317605,  1.325203,  1.332819,  1.340454,
0139         1.348108,  1.355780,  1.363472,  1.371182,  1.378912,  1.386660,
0140         1.394429,  1.402216,  1.410024,  1.417851,  1.425698,  1.433565,
0141         1.441453,  1.449360,  1.457288,  1.465237,  1.473206,  1.481196,
0142         1.489208,  1.497240,  1.505293,  1.513368,  1.521465,  1.529583,
0143         1.537723,  1.545885,  1.554068,  1.562275,  1.570503,  1.578754,
0144         1.587028,  1.595325,  1.603644,  1.611987,  1.620353,  1.628743,
0145         1.637156,  1.645593,  1.654053,  1.662538,  1.671047,  1.679581,
0146         1.688139,  1.696721,  1.705329,  1.713961,  1.722619,  1.731303,
0147         1.740011,  1.748746,  1.757506,  1.766293,  1.775106,  1.783945,
0148         1.792810,  1.801703,  1.810623,  1.819569,  1.828543,  1.837545,
0149         1.846574,  1.855631,  1.864717,  1.873830,  1.882972,  1.892143,
0150         1.901343,  1.910572,  1.919830,  1.929117,  1.938434,  1.947781,
0151         1.957158,  1.966566,  1.976004,  1.985473,  1.994972,  2.004503,
0152         2.014065,  2.023659,  2.033285,  2.042943,  2.052633,  2.062355,
0153         2.072110,  2.081899,  2.091720,  2.101575,  2.111464,  2.121386,
0154         2.131343,  2.141334,  2.151360,  2.161421,  2.171517,  2.181648,
0155         2.191815,  2.202018,  2.212257,  2.222533,  2.232845,  2.243195,
0156         2.253582,  2.264006,  2.274468,  2.284968,  2.295507,  2.306084,
0157         2.316701,  2.327356,  2.338051,  2.348786,  2.359562,  2.370377,
0158         2.381234,  2.392131,  2.403070,  2.414051,  2.425073,  2.436138,
0159         2.447246,  2.458397,  2.469591,  2.480828,  2.492110,  2.503436,
0160         2.514807,  2.526222,  2.537684,  2.549190,  2.560743,  2.572343,
0161         2.583989,  2.595682,  2.607423,  2.619212,  2.631050,  2.642936,
0162         2.654871,  2.666855,  2.678890,  2.690975,  2.703110,  2.715297,
0163         2.727535,  2.739825,  2.752168,  2.764563,  2.777012,  2.789514,
0164         2.802070,  2.814681,  2.827347,  2.840069,  2.852846,  2.865680,
0165         2.878570,  2.891518,  2.904524,  2.917588,  2.930712,  2.943894,
0166         2.957136,  2.970439,  2.983802,  2.997227,  3.010714,  3.024263,
0167         3.037875,  3.051551,  3.065290,  3.079095,  3.092965,  3.106900,
0168         3.120902,  3.134971,  3.149107,  3.163312,  3.177585,  3.191928,
0169         3.206340,  3.220824,  3.235378,  3.250005,  3.264704,  3.279477,
0170         3.294323,  3.309244,  3.324240,  3.339312,  3.354461,  3.369687,
0171         3.384992,  3.400375,  3.415838,  3.431381,  3.447005,  3.462711,
0172         3.478500,  3.494372,  3.510328,  3.526370,  3.542497,  3.558711,
0173         3.575012,  3.591402,  3.607881,  3.624450,  3.641111,  3.657863,
0174         3.674708,  3.691646,  3.708680,  3.725809,  3.743034,  3.760357,
0175         3.777779,  3.795300,  3.812921,  3.830645,  3.848470,  3.866400,
0176         3.884434,  3.902574,  3.920821,  3.939176,  3.957640,  3.976215,
0177         3.994901,  4.013699,  4.032612,  4.051639,  4.070783,  4.090045,
0178         4.109425,  4.128925,  4.148547,  4.168292,  4.188160,  4.208154,
0179         4.228275,  4.248524,  4.268903,  4.289413,  4.310056,  4.330832,
0180         4.351745,  4.372794,  4.393982,  4.415310,  4.436781,  4.458395,
0181         4.480154,  4.502060,  4.524114,  4.546319,  4.568676,  4.591187,
0182         4.613854,  4.636678,  4.659662,  4.682807,  4.706116,  4.729590,
0183         4.753231,  4.777041,  4.801024,  4.825179,  4.849511,  4.874020,
0184         4.898710,  4.923582,  4.948639,  4.973883,  4.999316,  5.024942,
0185         5.050761,  5.076778,  5.102993,  5.129411,  5.156034,  5.182864,
0186         5.209903,  5.237156,  5.264625,  5.292312,  5.320220,  5.348354,
0187         5.376714,  5.405306,  5.434131,  5.463193,  5.492496,  5.522042,
0188         5.551836,  5.581880,  5.612178,  5.642734,  5.673552,  5.704634,
0189         5.735986,  5.767610,  5.799512,  5.831694,  5.864161,  5.896918,
0190         5.929968,  5.963316,  5.996967,  6.030925,  6.065194,  6.099780,
0191         6.134687,  6.169921,  6.205486,  6.241387,  6.277630,  6.314220,
0192         6.351163,  6.388465,  6.426130,  6.464166,  6.502578,  6.541371,
0193         6.580553,  6.620130,  6.660109,  6.700495,  6.741297,  6.782520,
0194         6.824173,  6.866262,  6.908795,  6.951780,  6.995225,  7.039137,
0195         7.083525,  7.128398,  7.173764,  7.219632,  7.266011,  7.312910,
0196         7.360339,  7.408308,  7.456827,  7.505905,  7.555554,  7.605785,
0197         7.656608,  7.708035,  7.760077,  7.812747,  7.866057,  7.920019,
0198         7.974647,  8.029953,  8.085952,  8.142657,  8.200083,  8.258245,
0199         8.317158,  8.376837,  8.437300,  8.498562,  8.560641,  8.623554,
0200         8.687319,  8.751955,  8.817481,  8.883916,  8.951282,  9.019600,
0201         9.088889,  9.159174,  9.230477,  9.302822,  9.376233,  9.450735,
0202         9.526355,  9.603118,  9.681054,  9.760191,  9.840558,  9.922186,
0203         10.005107, 10.089353, 10.174959, 10.261958, 10.350389, 10.440287,
0204         10.531693, 10.624646, 10.719188, 10.815362, 10.913214, 11.012789,
0205         11.114137, 11.217307, 11.322352, 11.429325, 11.538283, 11.649285,
0206         11.762390, 11.877664, 11.995170, 12.114979, 12.237161, 12.361791,
0207         12.488946, 12.618708, 12.751161, 12.886394, 13.024498, 13.165570,
0208         13.309711, 13.457026, 13.607625, 13.761625, 13.919145, 14.080314,
0209         14.245263, 14.414134, 14.587072, 14.764233, 14.945778, 15.131877,
0210         15.322712, 15.518470, 15.719353, 15.925570, 16.137345, 16.354912,
0211         16.578520, 16.808433, 17.044929, 17.288305, 17.538873, 17.796967,
0212         18.062943, 18.337176, 18.620068, 18.912049, 19.213574, 19.525133,
0213         19.847249, 20.180480, 20.525429, 20.882738, 21.253102, 21.637266,
0214         22.036036, 22.450278, 22.880933, 23.329017, 23.795634, 24.281981,
0215         24.789364, 25.319207, 25.873062, 26.452634, 27.059789, 27.696581,
0216         28.365274, 29.068370, 29.808638, 30.589157, 31.413354, 32.285060,
0217         33.208568, 34.188705, 35.230920, 36.341388, 37.527131, 38.796172,
0218         40.157721, 41.622399, 43.202525, 44.912465, 46.769077, 48.792279,
0219         51.005773, 53.437996, 56.123356, 59.103894};
0220
0221     if (z <= 0) {
0222       // @TODO: Is using max OK?
0223       return -std::numeric_limits<scalar_type>::max();
0224     }
0225     if (z >= 1) {
0226       return std::numeric_limits<scalar_type>::max();
0227     }
0228
0229     double ranlan = 0.;
0230     double u = 0.;
0231     double v = 0.;
0232     u = 1000. * z;
0233     auto k = static_cast<int>(u);
0234     u -= static_cast<double>(k);
0235     if (k >= 70 && k < 800) {
0236       auto i{static_cast<std::size_t>(k)};
0237       ranlan = f[i - 1] + u * (f[i] - f[i - 1]);
0238     } else if (k >= 7 && k <= 980) {
0239       auto i{static_cast<std::size_t>(k)};
0240       ranlan = f[i - 1] +
0241                u * (f[i] - f[i - 1] -
0242                     0.25 * (1 - u) * (f[i + 1] - f[i] - f[i - 1] + f[i - 2]));
0243     } else if (k < 7) {
0244       v = math::log(z);
0245       u = 1 / v;
0246       ranlan = ((0.99858950 + (3.45213058E1 + 1.70854528E1 * u) * u) /
0247                 (1. + (3.41760202E1 + 4.01244582 * u) * u)) *
0248                (-math::log(-0.91893853 - v) - 1);
0249     } else {
0250       u = 1 - z;
0251       v = u * u;
0252       if (z <= 0.999) {
0253         ranlan = (1.00060006 + 2.63991156E2 * u + 4.37320068E3 * v) /
0254                  ((1. + 2.57368075E2 * u + 3.41448018E3 * v) * u);
0255       } else {
0256         ranlan = (1.00001538 + 6.07514119E3 * u + 7.34266409E5 * v) /
0257                  ((1. + 6.06511919E3 * u + 6.94021044E5 * v) * u);
0258       }
0259     }
0260     return static_cast<scalar_type>(ranlan);
0261   }
0262 };
0263
0264 }  // namespace detray