Physics/ElectroMagnetic/PhotonConversion.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 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "Acts/Definitions/PdgParticle.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/Utilities/UnitVectors.hpp"
0016 #include "ActsFatras/EventData/Barcode.hpp"
0017 #include "ActsFatras/EventData/Particle.hpp"
0018 #include "ActsFatras/EventData/ProcessType.hpp"
0019
0020 #include <algorithm>
0021 #include <array>
0022 #include <cmath>
0023 #include <limits>
0024 #include <numbers>
0025 #include <random>
0026 #include <utility>
0027 #include <vector>
0028
0029 namespace ActsFatras {
0030
0031 /// This class handles the photon conversion. It evaluates the distance
0032 /// after which the interaction will occur and the final state due the
0033 /// interaction itself.
0034 class PhotonConversion {
0035  public:
0036   /// Scaling factor of children energy
0037   double childEnergyScaleFactor = 2.;
0038   /// Scaling factor for photon conversion probability
0039   double conversionProbScaleFactor = 0.98;
0040
0041   /// Method for evaluating the distance after which the photon
0042   /// conversion will occur.
0043   ///
0044   /// @tparam generator_t Type of the random number generator
0045   /// @param [in, out] generator The random number generator
0046   /// @param [in] particle The particle
0047   ///
0048   /// @return valid X0 limit and no limit on L0
0049   template <typename generator_t>
0050   std::pair<double, double> generatePathLimits(generator_t& generator,
0051                                                const Particle& particle) const;
0052
0053   /// This method evaluates the final state due to the photon conversion.
0054   ///
0055   /// @tparam generator_t Type of the random number generator
0056   /// @param [in, out] generator The random number generator
0057   /// @param [in, out] particle The interacting photon
0058   /// @param [out] generated List of generated particles
0059   ///
0060   /// @return True if the conversion occurred, else false
0061   template <typename generator_t>
0062   bool run(generator_t& generator, Particle& particle,
0063            std::vector<Particle>& generated) const;
0064
0065  private:
0066   /// This method constructs and returns the child particles.
0067   ///
0068   /// @param [in] photon The interacting photon
0069   /// @param [in] childEnergy The energy of one child particle
0070   /// @param [in] childDirection The direction of the child particle
0071   ///
0072   /// @return Array containing the produced leptons
0073   std::array<Particle, 2> generateChildren(
0074       const Particle& photon, double childEnergy,
0075       const Acts::Vector3& childDirection) const;
0076
0077   /// Generate the energy fraction of the first child particle.
0078   ///
0079   /// @tparam generator_t Type of the random number generator
0080   /// @param [in, out] generator The random number generator
0081   /// @param [in] gammaMom The momentum of the photon
0082   ///
0083   /// @return The energy of the child particle
0084   template <typename generator_t>
0085   double generateFirstChildEnergyFraction(generator_t& generator,
0086                                           double gammaMom) const;
0087
0088   /// Generate the direction of the child particles.
0089   ///
0090   /// @tparam generator_t Type of the random number generator
0091   /// @param [in, out] generator The random number generator
0092   /// @param [in] particle The photon
0093   ///
0094   /// @return The direction vector of the child particle
0095   template <typename generator_t>
0096   Acts::Vector3 generateChildDirection(generator_t& generator,
0097                                        const Particle& particle) const;
0098
0099   /// Helper methods for momentum evaluation
0100   /// @note These methods are taken from the Geant4 class
0101   /// G4PairProductionRelModel
0102   double screenFunction1(double delta) const;
0103   double screenFunction2(double delta) const;
0104
0105   /// Helper method for the electron mass. This is used to avoid multiple
0106   /// lookups in the internal data tables.
0107   static float electronMass() {
0108     return Acts::findMass(Acts::PdgParticle::eElectron).value();
0109   }
0110 };
0111
0112 inline double PhotonConversion::screenFunction1(double delta) const {
0113   // Compute the value of the screening function 3*PHI1(delta) - PHI2(delta)
0114   return (delta > 1.4) ? 42.038 - 8.29 * std::log(delta + 0.958)
0115                        : 42.184 - delta * (7.444 - 1.623 * delta);
0116 }
0117
0118 inline double PhotonConversion::screenFunction2(double delta) const {
0119   // Compute the value of the screening function 1.5*PHI1(delta)
0120   // +0.5*PHI2(delta)
0121   return (delta > 1.4) ? 42.038 - 8.29 * std::log(delta + 0.958)
0122                        : 41.326 - delta * (5.848 - 0.902 * delta);
0123 }
0124
0125 template <typename generator_t>
0126 std::pair<double, double> PhotonConversion::generatePathLimits(
0127     generator_t& generator, const Particle& particle) const {
0128   /// This method is based upon the Athena class PhotonConversionTool
0129
0130   // Fast exit if not a photon or the energy is too low
0131   if (particle.pdg() != Acts::PdgParticle::eGamma ||
0132       particle.absoluteMomentum() < (2 * electronMass())) {
0133     return {std::numeric_limits<double>::infinity(),
0134             std::numeric_limits<double>::infinity()};
0135   }
0136
0137   // Use for the moment only Al data - Yung Tsai - Rev.Mod.Particle Physics Vol.
0138   // 46, No.4, October 1974 optainef from a fit given in the momentum range 100
0139   // 10 6 2 1 0.6 0.4 0.2 0.1 GeV
0140
0141   //// Quadratic background function
0142   //  Double_t fitFunction(Double_t *x, Double_t *par) {
0143   //  return par[0] + par[1]*pow(x[0],par[2]);
0144   // }
0145   // EXT PARAMETER                                   STEP         FIRST
0146   // NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
0147   //  1  p0          -7.01612e-03   8.43478e-01   1.62766e-04   1.11914e-05
0148   //  2  p1           7.69040e-02   1.00059e+00   8.90718e-05  -8.41167e-07
0149   //  3  p2          -6.07682e-01   5.13256e+00   6.07228e-04  -9.44448e-07
0150   constexpr double p0 = -7.01612e-03;
0151   constexpr double p1 = 7.69040e-02;
0152   constexpr double p2 = -6.07682e-01;
0153
0154   // Calculate xi
0155   const double xi = p0 + p1 * std::pow(particle.absoluteMomentum(), p2);
0156
0157   std::uniform_real_distribution<double> uniformDistribution{0., 1.};
0158   // This is a transformation of eq. 3.75
0159   return {-9. / 7. *
0160               std::log(conversionProbScaleFactor *
0161                        (1 - uniformDistribution(generator))) /
0162               (1. - xi),
0163           std::numeric_limits<double>::infinity()};
0164 }
0165
0166 template <typename generator_t>
0167 double PhotonConversion::generateFirstChildEnergyFraction(
0168     generator_t& generator, double gammaMom) const {
0169   /// This method is based upon the Geant4 class G4PairProductionRelModel
0170
0171   /// @note This method is from the Geant4 class G4Element
0172   //
0173   //  Compute Coulomb correction factor (Phys Rev. D50 3-1 (1994) page 1254)
0174   constexpr double k1 = 0.0083;
0175   constexpr double k2 = 0.20206;
0176   constexpr double k3 = 0.0020;  // This term is missing in Athena
0177   constexpr double k4 = 0.0369;
0178   constexpr double alphaEM = 1. / 137.;
0179   constexpr double m_Z = 13.;  // Aluminium
0180   constexpr double az2 = (alphaEM * m_Z) * (alphaEM * m_Z);
0181   constexpr double az4 = az2 * az2;
0182   constexpr double coulombFactor =
0183       (k1 * az4 + k2 + 1. / (1. + az2)) * az2 - (k3 * az4 + k4) * az4;
0184
0185   const double logZ13 = std::log(m_Z) * 1. / 3.;
0186   const double FZ = 8. * (logZ13 + coulombFactor);
0187   const double deltaMax = std::exp((42.038 - FZ) * 0.1206) - 0.958;
0188
0189   const double deltaPreFactor = 136. / std::pow(m_Z, 1. / 3.);
0190   const double eps0 = electronMass() / gammaMom;
0191   const double deltaFactor = deltaPreFactor * eps0;
0192   const double deltaMin = 4. * deltaFactor;
0193
0194   // Compute the limits of eps
0195   const double epsMin =
0196       std::max(eps0, 0.5 - 0.5 * std::sqrt(1. - deltaMin / deltaMax));
0197   const double epsRange = 0.5 - epsMin;
0198
0199   // Sample the energy rate (eps) of the created electron (or positron)
0200   const double F10 = screenFunction1(deltaMin) - FZ;
0201   const double F20 = screenFunction2(deltaMin) - FZ;
0202   const double NormF1 = F10 * epsRange * epsRange;
0203   const double NormF2 = 1.5 * F20;
0204
0205   // We will need 3 uniform random number for each trial of sampling
0206   double greject = 0.;
0207   double eps = 0.;
0208   std::uniform_real_distribution<double> rndmEngine;
0209   do {
0210     if (NormF1 > rndmEngine(generator) * (NormF1 + NormF2)) {
0211       eps = 0.5 - epsRange * std::pow(rndmEngine(generator), 1. / 3.);
0212       const double delta = deltaFactor / (eps * (1. - eps));
0213       greject = (screenFunction1(delta) - FZ) / F10;
0214     } else {
0215       eps = epsMin + epsRange * rndmEngine(generator);
0216       const double delta = deltaFactor / (eps * (1. - eps));
0217       greject = (screenFunction2(delta) - FZ) / F20;
0218     }
0219   } while (greject < rndmEngine(generator));
0220   //  End of eps sampling
0221   return eps * childEnergyScaleFactor;
0222 }
0223
0224 template <typename generator_t>
0225 Acts::Vector3 PhotonConversion::generateChildDirection(
0226     generator_t& generator, const Particle& particle) const {
0227   /// This method is based upon the Athena class PhotonConversionTool
0228
0229   // Following the Geant4 approximation from L. Urban
0230   // the azimutal angle
0231   double theta = electronMass() / particle.energy();
0232
0233   std::uniform_real_distribution<double> uniformDistribution{0., 1.};
0234   const double u = -std::log(uniformDistribution(generator) *
0235                              uniformDistribution(generator)) *
0236                    1.6;
0237
0238   theta *= (uniformDistribution(generator) < 0.25)
0239                ? u
0240                : u * 1. / 3.;  // 9./(9.+27) = 0.25
0241
0242   // draw the random orientation angle
0243   const auto psi = std::uniform_real_distribution<double>(
0244       -std::numbers::pi, std::numbers::pi)(generator);
0245
0246   Acts::Vector3 direction = particle.direction();
0247   // construct the combined rotation to the scattered direction
0248   Acts::RotationMatrix3 rotation(
0249       // rotation of the scattering deflector axis relative to the reference
0250       Acts::AngleAxis3(psi, direction) *
0251       // rotation by the scattering angle around the deflector axis
0252       Acts::AngleAxis3(theta, Acts::createCurvilinearUnitU(direction)));
0253   direction.applyOnTheLeft(rotation);
0254   return direction;
0255 }
0256
0257 inline std::array<Particle, 2> PhotonConversion::generateChildren(
0258     const Particle& photon, double childEnergy,
0259     const Acts::Vector3& childDirection) const {
0260   using namespace Acts::UnitLiterals;
0261
0262   // Calculate the child momentum
0263   const double massChild = electronMass();
0264   const double momentum1 =
0265       std::sqrt(childEnergy * childEnergy - massChild * massChild);
0266
0267   // Use energy-momentum conservation for the other child
0268   const Acts::Vector3 vtmp =
0269       photon.fourMomentum().template segment<3>(Acts::eMom0) -
0270       momentum1 * childDirection;
0271   const double momentum2 = vtmp.norm();
0272
0273   // The daughter particles are created with the explicit electron mass used in
0274   // the calculations for consistency. Using the full Particle constructor with
0275   // charge and mass also avoids an additional lookup in the internal data
0276   // tables.
0277   std::array<Particle, 2> children = {
0278       Particle(photon.particleId().makeDescendant(0), Acts::eElectron, -1_e,
0279                electronMass())
0280           .setPosition4(photon.fourPosition())
0281           .setDirection(childDirection)
0282           .setAbsoluteMomentum(momentum1)
0283           .setProcess(ProcessType::ePhotonConversion)
0284           .setReferenceSurface(photon.referenceSurface()),
0285       Particle(photon.particleId().makeDescendant(1), Acts::ePositron, 1_e,
0286                electronMass())
0287           .setPosition4(photon.fourPosition())
0288           .setDirection(childDirection)
0289           .setAbsoluteMomentum(momentum2)
0290           .setProcess(ProcessType::ePhotonConversion)
0291           .setReferenceSurface(photon.referenceSurface()),
0292   };
0293   return children;
0294 }
0295
0296 template <typename generator_t>
0297 bool PhotonConversion::run(generator_t& generator, Particle& particle,
0298                            std::vector<Particle>& generated) const {
0299   // Fast exit if particle is not a photon
0300   if (particle.pdg() != Acts::PdgParticle::eGamma) {
0301     return false;
0302   }
0303
0304   // Fast exit if momentum is too low
0305   const double p = particle.absoluteMomentum();
0306   if (p < (2 * electronMass())) {
0307     return false;
0308   }
0309
0310   // Get one child energy
0311   const double childEnergy = p * generateFirstChildEnergyFraction(generator, p);
0312
0313   // Now get the deflection
0314   const Acts::Vector3 childDir = generateChildDirection(generator, particle);
0315
0316   // Produce the final state
0317   const std::array<Particle, 2> finalState =
0318       generateChildren(particle, childEnergy, childDir);
0319   generated.insert(generated.end(), finalState.begin(), finalState.end());
0320
0321   return true;
0322 }
0323
0324 }  // namespace ActsFatras