src/Physics/BetheHeitler.cpp

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 #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 <utility>
0020
0021 ActsFatras::Particle ActsFatras::BetheHeitler::bremPhoton(
0022     const Particle &particle, Scalar gammaE, Scalar rndPsi, Scalar rndTheta1,
0023     Scalar rndTheta2, Scalar rndTheta3) const {
0024   // ------------------------------------------------------
0025   // simple approach
0026   // (a) simulate theta uniform within the opening angle of the relativistic
0027   // Hertz dipole
0028   //      theta_max = 1/gamma
0029   // (b)Following the Geant4 approximation from L. Urban -> encapsulate that
0030   // later
0031   //      the azimutal angle
0032
0033   Scalar psi = 2. * M_PI * rndPsi;
0034
0035   // the start of the equation
0036   Scalar theta = 0.;
0037   if (uniformHertzDipoleAngle) {
0038     // the simplest simulation
0039     theta = particle.mass() / particle.energy() * rndTheta1;
0040   } else {
0041     // ----->
0042     theta = particle.mass() / particle.energy();
0043     // follow
0044     constexpr Scalar a = 0.625;  // 5/8
0045     Scalar u = -log(rndTheta2 * rndTheta3) / a;
0046     theta *= (rndTheta1 < 0.25) ? u : u / 3.;  // 9./(9.+27) = 0.25
0047   }
0048
0049   Vector3 particleDirection = particle.direction();
0050   Vector3 photonDirection = particleDirection;
0051
0052   // construct the combined rotation to the scattered direction
0053   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::makeCurvilinearUnitU(particleDirection)));
0058   photonDirection.applyOnTheLeft(rotation);
0059
0060   Particle photon(particle.particleId().makeDescendant(0),
0061                   Acts::PdgParticle::eGamma);
0062   photon.setProcess(ActsFatras::ProcessType::eBremsstrahlung)
0063       .setPosition4(particle.fourPosition())
0064       .setDirection(photonDirection)
0065       .setAbsoluteMomentum(gammaE)
0066       .setReferenceSurface(particle.referenceSurface());
0067   return photon;
0068 }