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/ProcessType.hpp"
0016
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <numbers>
0020 #include <utility>
0021
0022 ActsFatras::Particle ActsFatras::BetheHeitler::bremPhoton(
0023     const Particle &particle, double gammaE, 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   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     double u = -log(rndTheta2 * rndTheta3) / a;
0047     theta *= (rndTheta1 < 0.25) ? u : u / 3.;  // 9./(9.+27) = 0.25
0048   }
0049
0050   Acts::Vector3 particleDirection = particle.direction();
0051   Acts::Vector3 photonDirection = particleDirection;
0052
0053   // construct the combined rotation to the scattered direction
0054   Acts::RotationMatrix3 rotation(
0055       // rotation of the scattering deflector axis relative to the reference
0056       Acts::AngleAxis3(psi, particleDirection) *
0057       // rotation by the scattering angle around the deflector axis
0058       Acts::AngleAxis3(theta, Acts::makeCurvilinearUnitU(particleDirection)));
0059   photonDirection.applyOnTheLeft(rotation);
0060
0061   Particle photon(particle.particleId().makeDescendant(0),
0062                   Acts::PdgParticle::eGamma);
0063   photon.setProcess(ActsFatras::ProcessType::eBremsstrahlung)
0064       .setPosition4(particle.fourPosition())
0065       .setDirection(photonDirection)
0066       .setAbsoluteMomentum(gammaE)
0067       .setReferenceSurface(particle.referenceSurface());
0068   return photon;
0069 }