Physics/EnergyLoss/BetheHeitler.hpp

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 #pragma once
0010
0011 #include "Fatras/Kernel/detail/RandomNumberDistributions.hpp"
0012
0013 namespace Fatras {
0014
0015 const double log_2 = std::log(2.);
0016
0017 /// The struct for the EnergyLoss physics list
0018 ///
0019 /// Bethe-Heitler for electron brem description as described here:
0020 /// "A Gaussian-mixture approximation of the Bethe–Heitler model of electron
0021 /// energy loss by bremsstrahlung" R. Frühwirth
0022 ///
0023 struct BetheHeitler {
0024
0025   /// The flag to include BetheHeitler process or not
0026   bool betheHeitler = true;
0027
0028   /// A scaling factor to
0029   double scaleFactor = 1.;
0030
0031   /// @brief Call operator for the Bethe-Heitler energy loss
0032   ///
0033   /// @tparam generator_t is a random number generator type
0034   /// @tparam detector_t is the detector information type
0035   /// @tparam particle_t is the particle information type
0036   ///
0037   /// @param[in] generator is the random number generator
0038   /// @param[in] detector the detector information
0039   /// @param[in] particle the particle which is being scattered
0040   ///
0041   /// @return eventually produced photons
0042   template <typename generator_t, typename detector_t, typename particle_t>
0043   std::vector<particle_t> operator()(generator_t &generator,
0044                                      const detector_t &detector,
0045                                      particle_t &particle) const {
0046
0047     // Do nothing if the flag is set to false
0048     if (not betheHeitler) {
0049       return {};
0050     }
0051
0052     double tInX0 = detector.thickness() / detector.material().X0();
0053
0054     // Take a random gamma-distributed value - depending on t/X0
0055     GammaDist gDist = GammaDist(tInX0 / log_2, 1.);
0056
0057     double u = gDist(generator);
0058     double z = std::exp(-1. * u);
0059     double sampledEnergyLoss = std::abs(scaleFactor * particle.E() * (z - 1.));
0060
0061     // apply the energy loss
0062     particle.energyLoss(sampledEnergyLoss);
0063
0064     // @TODO return photons, needs particle_creator_t
0065     return {};
0066   }
0067 };
0068
0069 } // namespace Fatras