Tests/Physics/EnergyLossTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018 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 http://mozilla.org/MPL/2.0/.
0008
0009 ///  Boost include(s)
0010 #define BOOST_TEST_MODULE AbortList Tests
0011
0012 #include <boost/test/included/unit_test.hpp>
0013 // leave blank line
0014
0015 #include <boost/test/data/test_case.hpp>
0016 // leave blank line
0017
0018 #include <boost/test/output_test_stream.hpp>
0019 // leave blank line
0020
0021 #include "Acts/Material/Material.hpp"
0022 #include "Acts/Material/MaterialProperties.hpp"
0023 #include "Fatras/Kernel/PhysicsList.hpp"
0024 #include "Fatras/Kernel/Process.hpp"
0025 #include "Fatras/Physics/EnergyLoss/BetheBloch.hpp"
0026 #include "Fatras/Physics/EnergyLoss/BetheHeitler.hpp"
0027 #include "Particle.hpp"
0028 #include <fstream>
0029 #include <random>
0030
0031 namespace bdata = boost::unit_test::data;
0032 namespace tt = boost::test_tools;
0033 namespace au = Acts::units;
0034
0035 namespace Fatras {
0036
0037 namespace Test {
0038
0039 // the generator
0040 typedef std::mt19937 Generator;
0041
0042 // standard generator
0043 Generator generator;
0044
0045 // some material
0046 Acts::Material berilium = Acts::Material(352.8, 407., 9.012, 4.,
0047                                          1.848 / (au::_cm * au::_cm * au::_cm));
0048
0049 /// The selector
0050 struct Selector {
0051
0052   /// call operator
0053   template <typename detector_t, typename particle_t>
0054   bool operator()(const detector_t, const particle_t &) const {
0055     return true;
0056   }
0057 };
0058
0059 /// write CSV flag
0060 bool write_csv = true;
0061
0062 std::ofstream os("EnergyLoss.csv", std::ofstream::out | std::ofstream::trunc);
0063
0064 /// Test the scattering implementation
0065 BOOST_DATA_TEST_CASE(
0066     EnergyLoss_test_,
0067     bdata::random(
0068         (bdata::seed = 20,
0069          bdata::distribution = std::uniform_real_distribution<>(0., 1.))) ^
0070         bdata::random(
0071             (bdata::seed = 21,
0072              bdata::distribution = std::uniform_real_distribution<>(0., 1.))) ^
0073         bdata::random(
0074             (bdata::seed = 22,
0075              bdata::distribution = std::uniform_real_distribution<>(0., 1.))) ^
0076         bdata::random((bdata::seed = 23,
0077                        bdata::distribution =
0078                            std::uniform_real_distribution<>(1.5, 10.5))) ^
0079         bdata::xrange(10000),
0080     x, y, z, p, index) {
0081
0082   Acts::MaterialProperties detector(berilium, 10. * Acts::units::_mm);
0083
0084   // create the particle and set the momentum
0085   /// position at 0.,0.,0
0086   Acts::Vector3D position{0., 0., 0.};
0087   // p of 1 GeV
0088   Acts::Vector3D momentum =
0089       p * Acts::units::_GeV * Acts::Vector3D(x, y, z).normalized();
0090   // positively charged
0091   double q = -1.;
0092   double m = 105.658367 * Acts::units::_MeV;        // muon mass
0093   const double me = 0.51099891 * Acts::units::_MeV; // electron mass
0094
0095   // create the particle
0096   Particle particle(position, momentum, m, q, 13, 1);
0097   BOOST_CHECK_EQUAL(particle.m(), m);
0098   BOOST_CHECK_EQUAL(particle.pdg(), 13);
0099
0100   // make the highland scatterer
0101   BetheBloch bbloch;
0102   BetheHeitler bheitler;
0103   double E = particle.E();
0104
0105   auto bbr = bbloch(generator, detector, particle);
0106   double eloss_io = E - particle.E();
0107
0108   // Check if the particle actually lost energy
0109   BOOST_CHECK(E >= particle.E());
0110   BOOST_CHECK(bbr.size() == 0);
0111
0112   // recreate the particle as an electron
0113   particle = Particle(position, momentum, me, q, 11, 1);
0114   BOOST_CHECK_EQUAL(particle.m(), me);
0115   BOOST_CHECK_EQUAL(particle.pdg(), 11);
0116
0117   E = particle.E();
0118
0119   auto bhr = bheitler(generator, detector, particle);
0120   double eloss_rad = E - particle.E();
0121   BOOST_CHECK(E >= particle.E());
0122   BOOST_CHECK(bhr.size() == 0);
0123
0124   // write out a csv file
0125   if (write_csv) {
0126     if (!index)
0127       os << "p,bethe_bloch,bethe_heitler" << '\n';
0128     os << particle.p() << "," << eloss_io << "," << eloss_rad << '\n';
0129   }
0130
0131   // Accept everything
0132   typedef Selector All;
0133   // Define the processes with selectors
0134   typedef Process<BetheBloch, All, All, All> BetheBlochProcess;
0135   typedef Process<BetheHeitler, All, All, All> BetheHeitlerProcess;
0136
0137   // now check the EnergyLoss as a PhysicsList
0138   typedef PhysicsList<BetheBlochProcess, BetheHeitlerProcess> EnergyLoss;
0139   EnergyLoss eLossPhysicsList;
0140
0141   std::vector<Particle> outgoing;
0142   BOOST_CHECK(!eLossPhysicsList(generator, detector, particle, outgoing));
0143   BOOST_CHECK(!outgoing.size());
0144 }
0145
0146 } // namespace Test
0147 } // namespace Fatras