Physics/ElectroMagnetic/PhotonConversion.hpp

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