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