cpu/material/bethe_equation.cpp

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 // Project include(s).
0010 #include "detray/definitions/pdg_particle.hpp"
0011 #include "detray/material/interaction.hpp"
0012 #include "detray/material/material.hpp"
0013 #include "detray/material/material_slab.hpp"
0014 #include "detray/material/predefined_materials.hpp"
0015
0016 // Detray test include(s)
0017 #include "detray/test/framework/types.hpp"
0018 #include "detray/test/utils/landau_distribution.hpp"
0019 #include "detray/test/utils/random_scatterer.hpp"
0020 #include "detray/test/utils/statistics.hpp"
0021
0022 // GTest include(s).
0023 #include <gtest/gtest.h>
0024
0025 // System include(s).
0026 #include <random>
0027
0028 using namespace detray;
0029
0030 using test_algebra = test::algebra;
0031 using scalar = test::scalar;
0032
0033 // Test class for energy loss with Bethe function
0034 // Input tuple: < material, particle type, momentum, expected output >
0035 // https://pdg.lbl.gov/2022/AtomicNuclearProperties for muon dEdX and range
0036 class EnergyLossBetheValidation
0037     : public ::testing::TestWithParam<
0038           std::tuple<material<scalar>, pdg_particle<scalar>, scalar, scalar>> {
0039 };
0040
0041 // This tests the material functionalities
0042 TEST_P(EnergyLossBetheValidation, bethe_energy_loss) {
0043   // Interaction object
0044   interaction<scalar> I;
0045
0046   // Zero incidence angle
0047   const scalar cos_inc_ang{1.f};
0048
0049   // Material slab with the 1 cm thickness
0050   material_slab<scalar> slab(std::get<0>(GetParam()), 1.f * unit<scalar>::cm);
0051
0052   // Particle
0053   pdg_particle<scalar> ptc = std::get<1>(GetParam());
0054
0055   // Path segment in the material
0056   const scalar path_segment{slab.path_segment(cos_inc_ang)};
0057
0058   // qOverP
0059   const scalar qop{ptc.charge() / std::get<2>(GetParam())};
0060
0061   // Bethe Stopping power in MeV * cm^2 / g
0062   const scalar dEdx{I.compute_energy_loss_bethe_bloch(
0063                         path_segment, slab.get_material(), ptc, {ptc, qop}) /
0064                     path_segment / slab.get_material().mass_density() /
0065                     (unit<scalar>::MeV * unit<scalar>::cm2 / unit<scalar>::g)};
0066
0067   const scalar expected_dEdx = std::get<3>(GetParam());
0068
0069   // Check if difference is within 5% error
0070   EXPECT_NEAR((expected_dEdx - dEdx) / dEdx, 0.f, 0.05f);
0071 }
0072
0073 INSTANTIATE_TEST_SUITE_P(
0074     detray_material_Bethe_0p1GeV_H2Liquid, EnergyLossBetheValidation,
0075     ::testing::Values(std::make_tuple(hydrogen_liquid<scalar>(), muon<scalar>(),
0076                                       0.1003f * unit<scalar>::GeV, 6.539f)));
0077
0078 INSTANTIATE_TEST_SUITE_P(
0079     detray_material_Bethe_1GeV_H2Liquid, EnergyLossBetheValidation,
0080     ::testing::Values(std::make_tuple(hydrogen_liquid<scalar>(), muon<scalar>(),
0081                                       1.101f * unit<scalar>::GeV, 4.182f)));
0082
0083 INSTANTIATE_TEST_SUITE_P(
0084     detray_material_Bethe_10GeV_H2Liquid, EnergyLossBetheValidation,
0085     ::testing::Values(std::make_tuple(hydrogen_liquid<scalar>(), muon<scalar>(),
0086                                       10.11f * unit<scalar>::GeV, 4.777f)));
0087
0088 INSTANTIATE_TEST_SUITE_P(
0089     detray_material_Bethe_100GeV_H2Liquid, EnergyLossBetheValidation,
0090     ::testing::Values(std::make_tuple(hydrogen_liquid<scalar>(), muon<scalar>(),
0091                                       100.1f * unit<scalar>::GeV, 5.305f)));
0092
0093 INSTANTIATE_TEST_SUITE_P(
0094     detray_material_Bethe_0p1GeV_HeGas, EnergyLossBetheValidation,
0095     ::testing::Values(std::make_tuple(helium_gas<scalar>(), muon<scalar>(),
0096                                       0.1003f * unit<scalar>::GeV, 3.082f)));
0097
0098 /*
0099 //@fixme: Test fails with He Gas and 1 GeV muons (18 % difference)
0100 INSTANTIATE_TEST_SUITE_P(
0101     detray_material_Bethe_1GeV_HeGas, EnergyLossBetheValidation,
0102     ::testing::Values(std::make_tuple(helium_gas<scalar>(), muon<scalar>(),
0103                                       1.101f * unit<scalar>::GeV, 2.133f)));
0104 */
0105
0106 INSTANTIATE_TEST_SUITE_P(
0107     detray_material_Bethe_10GeV_HeGas, EnergyLossBetheValidation,
0108     ::testing::Values(std::make_tuple(helium_gas<scalar>(), muon<scalar>(),
0109                                       10.11f * unit<scalar>::GeV, 2.768f)));
0110
0111 INSTANTIATE_TEST_SUITE_P(
0112     detray_material_Bethe_100GeV_HeGas, EnergyLossBetheValidation,
0113     ::testing::Values(std::make_tuple(helium_gas<scalar>(), muon<scalar>(),
0114                                       100.1f * unit<scalar>::GeV, 3.188f)));
0115
0116 INSTANTIATE_TEST_SUITE_P(
0117     detray_material_Bethe_0p1GeV_Al, EnergyLossBetheValidation,
0118     ::testing::Values(std::make_tuple(aluminium<scalar>(), muon<scalar>(),
0119                                       0.1003f * unit<scalar>::GeV, 2.533f)));
0120
0121 INSTANTIATE_TEST_SUITE_P(
0122     detray_material_Bethe_1GeV_Al, EnergyLossBetheValidation,
0123     ::testing::Values(std::make_tuple(aluminium<scalar>(), muon<scalar>(),
0124                                       1.101f * unit<scalar>::GeV, 1.744f)));
0125
0126 INSTANTIATE_TEST_SUITE_P(
0127     detray_material_Bethe_10GeV_Al, EnergyLossBetheValidation,
0128     ::testing::Values(std::make_tuple(aluminium<scalar>(), muon<scalar>(),
0129                                       10.11f * unit<scalar>::GeV, 2.097f)));
0130
0131 INSTANTIATE_TEST_SUITE_P(
0132     detray_material_Bethe_100GeV_Al, EnergyLossBetheValidation,
0133     ::testing::Values(std::make_tuple(aluminium<scalar>(), muon<scalar>(),
0134                                       100.1f * unit<scalar>::GeV, 2.360f)));
0135
0136 INSTANTIATE_TEST_SUITE_P(
0137     detray_material_Bethe_0p1GeV_Si, EnergyLossBetheValidation,
0138     ::testing::Values(std::make_tuple(silicon<scalar>(), muon<scalar>(),
0139                                       0.1003f * unit<scalar>::GeV, 2.608f)));
0140
0141 INSTANTIATE_TEST_SUITE_P(
0142     detray_material_Bethe_1GeV_Si, EnergyLossBetheValidation,
0143     ::testing::Values(std::make_tuple(silicon<scalar>(), muon<scalar>(),
0144                                       1.101f * unit<scalar>::GeV, 1.803f)));
0145
0146 INSTANTIATE_TEST_SUITE_P(
0147     detray_material_Bethe_10GeV_Si, EnergyLossBetheValidation,
0148     ::testing::Values(std::make_tuple(silicon<scalar>(), muon<scalar>(),
0149                                       10.11f * unit<scalar>::GeV, 2.177f)));
0150
0151 INSTANTIATE_TEST_SUITE_P(
0152     detray_material_Bethe_100GeV_Si, EnergyLossBetheValidation,
0153     ::testing::Values(std::make_tuple(silicon<scalar>(), muon<scalar>(),
0154                                       100.1f * unit<scalar>::GeV, 2.451f)));
0155
0156 INSTANTIATE_TEST_SUITE_P(detray_material_Bethe_0p1GeV_Si_with_DED,
0157                          EnergyLossBetheValidation,
0158                          ::testing::Values(std::make_tuple(
0159                              silicon_with_ded<scalar>(), muon<scalar>(),
0160                              0.1003f * unit<scalar>::GeV, 2.608f)));
0161
0162 INSTANTIATE_TEST_SUITE_P(detray_material_Bethe_1GeV_Si_with_DED,
0163                          EnergyLossBetheValidation,
0164                          ::testing::Values(std::make_tuple(
0165                              silicon_with_ded<scalar>(), muon<scalar>(),
0166                              1.101f * unit<scalar>::GeV, 1.803f)));
0167
0168 INSTANTIATE_TEST_SUITE_P(detray_material_Bethe_10GeV_Si_with_DED,
0169                          EnergyLossBetheValidation,
0170                          ::testing::Values(std::make_tuple(
0171                              silicon_with_ded<scalar>(), muon<scalar>(),
0172                              10.11f * unit<scalar>::GeV, 2.177f)));
0173
0174 INSTANTIATE_TEST_SUITE_P(detray_material_Bethe_100GeV_Si_with_DED,
0175                          EnergyLossBetheValidation,
0176                          ::testing::Values(std::make_tuple(
0177                              silicon_with_ded<scalar>(), muon<scalar>(),
0178                              100.1f * unit<scalar>::GeV, 2.451f)));
0179
0180 INSTANTIATE_TEST_SUITE_P(
0181     detray_material_Bethe_0p1GeV_ArLiquid, EnergyLossBetheValidation,
0182     ::testing::Values(std::make_tuple(argon_liquid<scalar>(), muon<scalar>(),
0183                                       0.1003f * unit<scalar>::GeV, 2.34f)));
0184
0185 // @fixme ~6% discrepancy
0186 /*INSTANTIATE_TEST_SUITE_P(
0187     detray_material_Bethe_1GeV_ArLiquid, EnergyLossBetheValidation,
0188     ::testing::Values(std::make_tuple(argon_liquid<scalar>(), muon<scalar>(),
0189                                       1.101f * unit<scalar>::GeV, 1.644f)));
0190 */
0191
0192 INSTANTIATE_TEST_SUITE_P(
0193     detray_material_Bethe_10GeV_ArLiquid, EnergyLossBetheValidation,
0194     ::testing::Values(std::make_tuple(argon_liquid<scalar>(), muon<scalar>(),
0195                                       10.11f * unit<scalar>::GeV, 2.003f)));
0196
0197 INSTANTIATE_TEST_SUITE_P(
0198     detray_material_Bethe_100GeV_ArLiquid, EnergyLossBetheValidation,
0199     ::testing::Values(std::make_tuple(argon_liquid<scalar>(), muon<scalar>(),
0200                                       100.1f * unit<scalar>::GeV, 2.258f)));
0201
0202 INSTANTIATE_TEST_SUITE_P(
0203     detray_material_Bethe_0p1GeV_Fe, EnergyLossBetheValidation,
0204     ::testing::Values(std::make_tuple(iron<scalar>(), muon<scalar>(),
0205                                       0.1003f * unit<scalar>::GeV, 2.274f)));
0206
0207 // @fixme ~6% discrepancy
0208 /*INSTANTIATE_TEST_SUITE_P(
0209     detray_material_Bethe_1GeV_Fe, EnergyLossBetheValidation,
0210     ::testing::Values(std::make_tuple(iron<scalar>(), muon<scalar>(),
0211                                       1.101f * unit<scalar>::GeV, 1.581f)));
0212 */
0213
0214 INSTANTIATE_TEST_SUITE_P(
0215     detray_material_Bethe_10GeV_Fe, EnergyLossBetheValidation,
0216     ::testing::Values(std::make_tuple(iron<scalar>(), muon<scalar>(),
0217                                       10.11f * unit<scalar>::GeV, 1.942f)));
0218
0219 INSTANTIATE_TEST_SUITE_P(
0220     detray_material_Bethe_100GeV_Fe, EnergyLossBetheValidation,
0221     ::testing::Values(std::make_tuple(iron<scalar>(), muon<scalar>(),
0222                                       100.1f * unit<scalar>::GeV, 2.207f)));
0223
0224 INSTANTIATE_TEST_SUITE_P(
0225     detray_material_Bethe_0p1GeV_Fe_with_DED, EnergyLossBetheValidation,
0226     ::testing::Values(std::make_tuple(iron_with_ded<scalar>(), muon<scalar>(),
0227                                       0.1003f * unit<scalar>::GeV, 2.274f)));
0228
0229 INSTANTIATE_TEST_SUITE_P(
0230     detray_material_Bethe_1GeV_Fe_with_DED, EnergyLossBetheValidation,
0231     ::testing::Values(std::make_tuple(iron_with_ded<scalar>(), muon<scalar>(),
0232                                       1.101f * unit<scalar>::GeV, 1.581f)));
0233
0234 INSTANTIATE_TEST_SUITE_P(
0235     detray_material_Bethe_10GeV_Fe_with_DED, EnergyLossBetheValidation,
0236     ::testing::Values(std::make_tuple(iron_with_ded<scalar>(), muon<scalar>(),
0237                                       10.11f * unit<scalar>::GeV, 1.942f)));
0238
0239 INSTANTIATE_TEST_SUITE_P(
0240     detray_material_Bethe_100GeV_Fe_with_DED, EnergyLossBetheValidation,
0241     ::testing::Values(std::make_tuple(iron_with_ded<scalar>(), muon<scalar>(),
0242                                       100.1f * unit<scalar>::GeV, 2.207f)));
0243
0244 INSTANTIATE_TEST_SUITE_P(
0245     detray_material_Bethe_0p1GeV_Cu, EnergyLossBetheValidation,
0246     ::testing::Values(std::make_tuple(copper<scalar>(), muon<scalar>(),
0247                                       0.1003f * unit<scalar>::GeV, 2.198f)));
0248
0249 // @fixme
0250 /*INSTANTIATE_TEST_SUITE_P(
0251     detray_material_Bethe_1GeV_Cu, EnergyLossBetheValidation,
0252     ::testing::Values(std::make_tuple(copper<scalar>(), muon<scalar>(),
0253                                       1.101f * unit<scalar>::GeV, 1.532f)));
0254 */
0255
0256 INSTANTIATE_TEST_SUITE_P(
0257     detray_material_Bethe_10GeV_Cu, EnergyLossBetheValidation,
0258     ::testing::Values(std::make_tuple(copper<scalar>(), muon<scalar>(),
0259                                       10.11f * unit<scalar>::GeV, 1.891f)));
0260
0261 INSTANTIATE_TEST_SUITE_P(
0262     detray_material_Bethe_100GeV_Cu, EnergyLossBetheValidation,
0263     ::testing::Values(std::make_tuple(copper<scalar>(), muon<scalar>(),
0264                                       100.1f * unit<scalar>::GeV, 2.155f)));
0265
0266 INSTANTIATE_TEST_SUITE_P(
0267     detray_material_Bethe_0p1GeV_Cu_with_DED, EnergyLossBetheValidation,
0268     ::testing::Values(std::make_tuple(copper_with_ded<scalar>(), muon<scalar>(),
0269                                       0.1003f * unit<scalar>::GeV, 2.198f)));
0270
0271 INSTANTIATE_TEST_SUITE_P(
0272     detray_material_Bethe_1GeV_Cu_with_DED, EnergyLossBetheValidation,
0273     ::testing::Values(std::make_tuple(copper_with_ded<scalar>(), muon<scalar>(),
0274                                       1.101f * unit<scalar>::GeV, 1.532f)));
0275
0276 INSTANTIATE_TEST_SUITE_P(
0277     detray_material_Bethe_10GeV_Cu_with_DED, EnergyLossBetheValidation,
0278     ::testing::Values(std::make_tuple(copper_with_ded<scalar>(), muon<scalar>(),
0279                                       10.11f * unit<scalar>::GeV, 1.891f)));
0280
0281 INSTANTIATE_TEST_SUITE_P(
0282     detray_material_Bethe_100GeV_Cu_with_DED, EnergyLossBetheValidation,
0283     ::testing::Values(std::make_tuple(copper_with_ded<scalar>(), muon<scalar>(),
0284                                       100.1f * unit<scalar>::GeV, 2.155f)));
0285
0286 // @fixme ~10% discrepancy
0287 /*INSTANTIATE_TEST_SUITE_P(
0288     detray_material_Bethe_0p1GeV_CsI, EnergyLossBetheValidation,
0289     ::testing::Values(std::make_tuple(cesium_iodide<scalar>(), muon<scalar>(),
0290                                       0.1003f * unit<scalar>::GeV, 1.869f)));
0291 */
0292
0293 // @fixme ~10% discrepancy
0294 /*INSTANTIATE_TEST_SUITE_P(
0295     detray_material_Bethe_1GeV_CsI, EnergyLossBetheValidation,
0296     ::testing::Values(std::make_tuple(cesium_iodide<scalar>(), muon<scalar>(),
0297                                       1.101f * unit<scalar>::GeV, 1.391f)));
0298 */
0299
0300 INSTANTIATE_TEST_SUITE_P(
0301     detray_material_Bethe_10GeV_CsI, EnergyLossBetheValidation,
0302     ::testing::Values(std::make_tuple(cesium_iodide<scalar>(), muon<scalar>(),
0303                                       10.11f * unit<scalar>::GeV, 1.755f)));
0304
0305 INSTANTIATE_TEST_SUITE_P(
0306     detray_material_Bethe_100GeV_CsI, EnergyLossBetheValidation,
0307     ::testing::Values(std::make_tuple(cesium_iodide<scalar>(), muon<scalar>(),
0308                                       100.1f * unit<scalar>::GeV, 2.012f)));
0309
0310 INSTANTIATE_TEST_SUITE_P(detray_material_Bethe_0p1GeV_CsI_with_DED,
0311                          EnergyLossBetheValidation,
0312                          ::testing::Values(std::make_tuple(
0313                              cesium_iodide_with_ded<scalar>(), muon<scalar>(),
0314                              0.1003f * unit<scalar>::GeV, 1.869f)));
0315
0316 INSTANTIATE_TEST_SUITE_P(detray_material_Bethe_1GeV_CsI_with_DED,
0317                          EnergyLossBetheValidation,
0318                          ::testing::Values(std::make_tuple(
0319                              cesium_iodide_with_ded<scalar>(), muon<scalar>(),
0320                              1.101f * unit<scalar>::GeV, 1.391f)));
0321
0322 INSTANTIATE_TEST_SUITE_P(detray_material_Bethe_10GeV_CsI_with_DED,
0323                          EnergyLossBetheValidation,
0324                          ::testing::Values(std::make_tuple(
0325                              cesium_iodide_with_ded<scalar>(), muon<scalar>(),
0326                              10.11f * unit<scalar>::GeV, 1.755f)));
0327
0328 INSTANTIATE_TEST_SUITE_P(detray_material_Bethe_100GeV_CsI_with_DED,
0329                          EnergyLossBetheValidation,
0330                          ::testing::Values(std::make_tuple(
0331                              cesium_iodide_with_ded<scalar>(), muon<scalar>(),
0332                              100.1f * unit<scalar>::GeV, 2.012f)));
0333
0334 // Test class for MUON energy loss with Landau function
0335 // Input tuple: < material / particle type / energy / expected energy loss  /
0336 // expected fwhm  >
0337 class EnergyLossLandauValidation
0338     : public ::testing::TestWithParam<std::tuple<
0339           material<scalar>, pdg_particle<scalar>, scalar, scalar, scalar>> {};
0340
0341 TEST_P(EnergyLossLandauValidation, landau_energy_loss) {
0342   // Interaction object
0343   interaction<scalar> I;
0344
0345   // Zero incidence angle
0346   const scalar cos_inc_ang{1.f};
0347
0348   // Material
0349   const auto mat = std::get<0>(GetParam());
0350
0351   // Thickness
0352   const scalar thickness = 0.17f * unit<scalar>::cm;
0353
0354   // Material slab with a unit thickness
0355   material_slab<scalar> slab(mat, thickness);
0356
0357   // Particle
0358   pdg_particle<scalar> ptc = std::get<1>(GetParam());
0359
0360   // Path segment in the material
0361   const scalar path_segment{slab.path_segment(cos_inc_ang)};
0362
0363   // Energy
0364   const scalar E = std::get<2>(GetParam());
0365
0366   // p
0367   const scalar p = std::sqrt(E * E - ptc.mass() * ptc.mass());
0368
0369   // q
0370   const scalar q = ptc.charge();
0371
0372   // qOverP
0373   const scalar qop{q / p};
0374
0375   // Landau Energy loss in MeV
0376   const scalar Landau_MeV{I.compute_energy_loss_landau(path_segment,
0377                                                        slab.get_material(), ptc,
0378                                                        {ptc, qop}) /
0379                           unit<scalar>::MeV};
0380
0381   // Check if difference is within 5% error
0382   EXPECT_TRUE(std::abs(std::get<3>(GetParam()) - Landau_MeV) / Landau_MeV <
0383               0.05f);
0384
0385   // Landau Energy loss Fluctuation
0386   const scalar fwhm_MeV{I.compute_energy_loss_landau_fwhm(path_segment,
0387                                                           slab.get_material(),
0388                                                           ptc, {ptc, qop}) /
0389                         unit<scalar>::MeV};
0390
0391   // Check if difference is within 10% error
0392   EXPECT_TRUE(std::abs(std::get<4>(GetParam()) - fwhm_MeV) / fwhm_MeV < 0.1f);
0393 }
0394
0395 // Expected output from Fig 33.7 in RPP2018
0396 INSTANTIATE_TEST_SUITE_P(
0397     detray_material_Landau_10GeV_Silicon, EnergyLossLandauValidation,
0398     ::testing::Values(std::make_tuple(silicon<scalar>(), muon<scalar>(),
0399                                       10.f * unit<scalar>::GeV, 0.525f,
0400                                       0.13f)));
0401
0402 // Input tuple: < energy >
0403 class LandauDistributionValidation
0404     : public ::testing::TestWithParam<
0405           std::tuple<pdg_particle<scalar>, scalar>> {};
0406
0407 TEST_P(LandauDistributionValidation, landau_distribution) {
0408   // Random generator
0409   std::random_device rd{};
0410   std::mt19937_64 generator{rd()};
0411   generator.seed(0u);
0412
0413   // Interaction object
0414   interaction<scalar> I;
0415
0416   // Zero incidence angle
0417   const scalar cos_inc_ang{1.f};
0418
0419   // Material
0420   const auto mat = silicon<scalar>();
0421
0422   // Thickness
0423   const auto thickness = 1.f * unit<scalar>::mm;
0424
0425   // Material slab with a unit thickness
0426   material_slab<scalar> slab(mat, thickness);
0427
0428   // Particle
0429   pdg_particle<scalar> ptc = std::get<0>(GetParam());
0430
0431   // Path segment in the material
0432   const scalar path_segment{slab.path_segment(cos_inc_ang)};
0433
0434   // Energy
0435   const scalar E = std::get<1>(GetParam());
0436
0437   // p
0438   const scalar p = std::sqrt(E * E - ptc.mass() * ptc.mass());
0439
0440   // q
0441   const scalar q = ptc.charge();
0442
0443   // qOverP
0444   const scalar qop{q / p};
0445
0446   // Bethe energy loss
0447   const scalar dE{I.compute_energy_loss_bethe_bloch(
0448       path_segment, slab.get_material(), ptc, {ptc, qop})};
0449
0450   // Landau Energy loss
0451   const scalar mpv{I.compute_energy_loss_landau(
0452       path_segment, slab.get_material(), ptc, {ptc, qop})};
0453
0454   // Landau Energy loss Sigma
0455   const scalar sigma{I.compute_energy_loss_landau_sigma(
0456       path_segment, slab.get_material(), ptc, {ptc, qop})};
0457
0458   // Landau Sampling
0459   landau_distribution<scalar> landau;
0460
0461   std::vector<scalar> samples;
0462   std::size_t n_samples{10000u};
0463   for (std::size_t i = 0u; i < n_samples; i++) {
0464     const scalar sa = landau(generator, mpv, sigma);
0465     samples.push_back(sa);
0466   }
0467
0468   // Mean energy loss from Landau samples
0469   const scalar sample_mean = statistics::mean(samples);
0470
0471   // Make sure that the difference is within 50% (...)
0472   // @note: We should not expect a good consistency between landau
0473   // distribution samples and bethe bloch function because the arguments
0474   // for landau_distribution are not the most probable value and sigma
0475   // (sigma is not even defined in Landau distribution). We might need to put
0476   // a smart scaling factor
0477   EXPECT_NEAR((dE - sample_mean) / dE, 0.f, 0.5f);
0478
0479   // Test attenuate function
0480   std::vector<scalar> energies;
0481
0482   const scalar mass = ptc.mass();
0483
0484   for (std::size_t i = 0u; i < n_samples; i++) {
0485     const scalar new_p = actor::random_scatterer<test_algebra>().attenuate(
0486         mpv, sigma, mass, p, generator);
0487     ASSERT_TRUE(new_p < p);
0488
0489     const scalar new_E = std::sqrt(new_p * new_p + mass * mass);
0490     energies.push_back(new_E);
0491   }
0492
0493   const scalar E_mean = statistics::mean(energies);
0494   const scalar E_expected = E - dE;
0495   EXPECT_TRUE(E_expected < E);
0496
0497   // Make sure that the difference is within 65%
0498   EXPECT_NEAR((E_mean - E_expected) / dE, 0.f, 0.65f);
0499 }
0500
0501 INSTANTIATE_TEST_SUITE_P(
0502     detray_material_Landau, LandauDistributionValidation,
0503     ::testing::Values(std::make_tuple(muon<scalar>(), 1.f * unit<scalar>::GeV),
0504                       std::make_tuple(muon<scalar>(), 10.f * unit<scalar>::GeV),
0505                       std::make_tuple(muon<scalar>(),
0506                                       100.f * unit<scalar>::GeV)));