src/Physics/BetheHeitler.cpp

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 #include "ActsFatras/Physics/ElectroMagnetic/BetheHeitler.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/PdgParticle.hpp"
0013 #include "Acts/Utilities/UnitVectors.hpp"
0014 #include "ActsFatras/EventData/Barcode.hpp"
0015 #include "ActsFatras/EventData/GenerationProcess.hpp"
0016
0017 #include <cmath>
0018 #include <numbers>
0019
0020 namespace ActsFatras {
0021
0022 Particle BetheHeitler::bremPhoton(const Particle &particle, double gammaE,
0023                                   double rndPsi, double rndTheta1,
0024                                   double rndTheta2, double rndTheta3) const {
0025   // ------------------------------------------------------
0026   // simple approach
0027   // (a) simulate theta uniform within the opening angle of the relativistic
0028   // Hertz dipole
0029   //      theta_max = 1/gamma
0030   // (b)Following the Geant4 approximation from L. Urban -> encapsulate that
0031   // later
0032   //      the azimutal angle
0033
0034   const double psi = 2. * std::numbers::pi * rndPsi;
0035
0036   // the start of the equation
0037   double theta = 0.;
0038   if (uniformHertzDipoleAngle) {
0039     // the simplest simulation
0040     theta = particle.mass() / particle.energy() * rndTheta1;
0041   } else {
0042     // ----->
0043     theta = particle.mass() / particle.energy();
0044     // follow
0045     constexpr double a = 0.625;  // 5/8
0046     const double u = -std::log(rndTheta2 * rndTheta3) / a;
0047     theta *= (rndTheta1 < 0.25) ? u : u / 3.;  // 9./(9.+27) = 0.25
0048   }
0049
0050   const Acts::Vector3 particleDirection = particle.direction();
0051
0052   // construct the combined rotation to the scattered direction
0053   const Acts::RotationMatrix3 rotation(
0054       // rotation of the scattering deflector axis relative to the reference
0055       Acts::AngleAxis3(psi, particleDirection) *
0056       // rotation by the scattering angle around the deflector axis
0057       Acts::AngleAxis3(theta, Acts::createCurvilinearUnitU(particleDirection)));
0058   const Acts::Vector3 photonDirection = rotation * particleDirection;
0059
0060   Particle photon(particle.particleId().makeDescendant(0),
0061                   Acts::PdgParticle::eGamma);
0062   photon.setProcess(GenerationProcess::eBremsstrahlung)
0063       .setPosition4(particle.fourPosition())
0064       .setDirection(photonDirection)
0065       .setAbsoluteMomentum(gammaE)
0066       .setReferenceSurface(particle.referenceSurface());
0067   return photon;
0068 }
0069
0070 }  // namespace ActsFatras