fix: Fix Fatras photon conversion and cleanup interaction code

Author: andiwand

Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheBloch.hpp

@@ -56,7 +56,7 @@ struct BetheBloch {
     //      probable value and the Gaussian-equivalent sigma
     LandauDistribution lossDistribution(scaleFactorMPV * energyLoss,
                                         scaleFactorSigma * energyLossSigma);
-    const auto loss = lossDistribution(generator);
+    const double loss = lossDistribution(generator);
 
     // Apply the energy loss
     particle.correctEnergy(-loss);

Fatras/include/ActsFatras/Physics/ElectroMagnetic/BetheHeitler.hpp

@@ -61,9 +61,9 @@ struct BetheHeitler {
     detail::FpeSafeGammaDistribution gDist(
         slab.thicknessInX0() / std::numbers::ln2, 1.);
 
-    const auto u = gDist(generator);
-    const auto z = std::exp(-u);
-    const auto sampledEnergyLoss =
+    const double u = gDist(generator);
+    const double z = std::exp(-u);
+    const double sampledEnergyLoss =
         std::abs(scaleFactor * particle.energy() * (z - 1.));
 
     std::uniform_real_distribution<double> uDist(0., 1.);

Fatras/include/ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp

@@ -67,13 +67,13 @@ class PhotonConversion {
   /// This method constructs and returns the child particles.
   ///
   /// @param [in] photon The interacting photon
-  /// @param [in] childEnergy The energy of one child particle
-  /// @param [in] childDirection The direction of the child particle
+  /// @param [in] child1Energy The energy of the first child particle
+  /// @param [in] child1Direction The direction of the first child particle
   ///
   /// @return Array containing the produced leptons
   std::array<Particle, 2> generateChildren(
-      const Particle& photon, double childEnergy,
-      const Acts::Vector3& childDirection) const;
+      const Particle& photon, double child1Energy,
+      const Acts::Vector3& child1Direction) const;
 
   /// Generate the energy fraction of the first child particle.
   ///
@@ -225,50 +225,47 @@ double PhotonConversion::generateFirstChildEnergyFraction(
 template <typename generator_t>
 Acts::Vector3 PhotonConversion::generateChildDirection(
     generator_t& generator, const Particle& particle) const {
-  /// This method is based upon the Athena class PhotonConversionTool
-
-  // Following the Geant4 approximation from L. Urban
-  // the azimutal angle
-  double theta = electronMass() / particle.energy();
+  // This method is based upon the Athena class PhotonConversionTool
 
+  // Following the Geant4 approximation from L. Urban the azimutal angle
   std::uniform_real_distribution<double> uniformDistribution{0., 1.};
   const double u = -std::log(uniformDistribution(generator) *
                              uniformDistribution(generator)) *
                    1.6;
-
-  theta *= (uniformDistribution(generator) < 0.25)
-               ? u
-               : u * 1. / 3.;  // 9./(9.+27) = 0.25
+  const double theta = (electronMass() / particle.energy()) *
+                       ((uniformDistribution(generator) < 0.25)
+                            ? u
+                            : u * 1. / 3.);  // 9./(9.+27) = 0.25
 
   // draw the random orientation angle
   const auto psi = std::uniform_real_distribution<double>(
       -std::numbers::pi, std::numbers::pi)(generator);
 
-  Acts::Vector3 direction = particle.direction();
   // construct the combined rotation to the scattered direction
-  Acts::RotationMatrix3 rotation(
+  const Acts::RotationMatrix3 rotation(
       // rotation of the scattering deflector axis relative to the reference
-      Acts::AngleAxis3(psi, direction) *
+      Acts::AngleAxis3(psi, particle.direction()) *
       // rotation by the scattering angle around the deflector axis
-      Acts::AngleAxis3(theta, Acts::createCurvilinearUnitU(direction)));
-  direction.applyOnTheLeft(rotation);
-  return direction;
+      Acts::AngleAxis3(theta,
+                       Acts::createCurvilinearUnitU(particle.direction())));
+  return rotation * particle.direction();
 }
 
 inline std::array<Particle, 2> PhotonConversion::generateChildren(
-    const Particle& photon, double childEnergy,
-    const Acts::Vector3& childDirection) const {
+    const Particle& photon, double child1Energy,
+    const Acts::Vector3& child1Direction) const {
   using namespace Acts::UnitLiterals;
 
   // Calculate the child momentum
   const double massChild = electronMass();
-  const double momentum1 = Acts::fastCathetus(childEnergy, massChild);
+  const double absoluteMomentum1 = Acts::fastCathetus(child1Energy, massChild);
 
   // Use energy-momentum conservation for the other child
-  const Acts::Vector3 vtmp =
+  const Acts::Vector3 momentum2 =
       photon.fourMomentum().template segment<3>(Acts::eMom0) -
-      momentum1 * childDirection;
-  const double momentum2 = vtmp.norm();
+      absoluteMomentum1 * child1Direction;
+  const Acts::Vector3 child2Direction = momentum2.normalized();
+  const double absoluteMomentum2 = momentum2.norm();
 
   // The daughter particles are created with the explicit electron mass used in
   // the calculations for consistency. Using the full Particle constructor with
@@ -278,15 +275,15 @@ inline std::array<Particle, 2> PhotonConversion::generateChildren(
       Particle(photon.particleId().makeDescendant(0), Acts::eElectron, -1_e,
                electronMass())
           .setPosition4(photon.fourPosition())
-          .setDirection(childDirection)
-          .setAbsoluteMomentum(momentum1)
+          .setDirection(child1Direction)
+          .setAbsoluteMomentum(absoluteMomentum1)
           .setProcess(ProcessType::ePhotonConversion)
           .setReferenceSurface(photon.referenceSurface()),
       Particle(photon.particleId().makeDescendant(1), Acts::ePositron, 1_e,
                electronMass())
           .setPosition4(photon.fourPosition())
-          .setDirection(childDirection)
-          .setAbsoluteMomentum(momentum2)
+          .setDirection(child2Direction)
+          .setAbsoluteMomentum(absoluteMomentum2)
           .setProcess(ProcessType::ePhotonConversion)
           .setReferenceSurface(photon.referenceSurface()),
   };
@@ -308,14 +305,15 @@ bool PhotonConversion::run(generator_t& generator, Particle& particle,
   }
 
   // Get one child energy
-  const double childEnergy = p * generateFirstChildEnergyFraction(generator, p);
+  const double child1Energy =
+      p * generateFirstChildEnergyFraction(generator, p);
 
   // Now get the deflection
-  const Acts::Vector3 childDir = generateChildDirection(generator, particle);
+  const Acts::Vector3 child1Dir = generateChildDirection(generator, particle);
 
   // Produce the final state
   const std::array<Particle, 2> finalState =
-      generateChildren(particle, childEnergy, childDir);
+      generateChildren(particle, child1Energy, child1Dir);
   generated.insert(generated.end(), finalState.begin(), finalState.end());
 
   return true;

Fatras/include/ActsFatras/Physics/ElectroMagnetic/Scattering.hpp

@@ -53,20 +53,19 @@ struct GenericScattering {
     // drawn from the specific scattering model distribution.
 
     // draw the random orientation angle
-    const auto psi = std::uniform_real_distribution<double>(
+    const double psi = std::uniform_real_distribution<double>(
         -std::numbers::pi, std::numbers::pi)(generator);
     // draw the scattering angle
-    const auto theta = angle(generator, slab, particle);
+    const double theta = angle(generator, slab, particle);
 
-    Acts::Vector3 direction = particle.direction();
     // construct the combined rotation to the scattered direction
-    Acts::RotationMatrix3 rotation(
+    const Acts::RotationMatrix3 rotation(
         // rotation of the scattering deflector axis relative to the reference
-        Acts::AngleAxis3(psi, direction) *
+        Acts::AngleAxis3(psi, particle.direction()) *
         // rotation by the scattering angle around the deflector axis
-        Acts::AngleAxis3(theta, Acts::createCurvilinearUnitU(direction)));
-    direction.applyOnTheLeft(rotation);
-    particle.setDirection(direction);
+        Acts::AngleAxis3(theta,
+                         Acts::createCurvilinearUnitU(particle.direction())));
+    particle.setDirection(rotation * particle.direction());
 
     // scattering is non-destructive and produces no secondaries
     return {};

Fatras/include/ActsFatras/Physics/ElectroMagnetic/detail/GaussianMixture.hpp

@@ -42,10 +42,9 @@ struct GaussianMixture {
   double operator()(generator_t &generator, const Acts::MaterialSlab &slab,
                     Particle &particle) const {
     /// Calculate the highland formula first
-    double sigma = Acts::computeMultipleScatteringTheta0(
+    const double sigma = Acts::computeMultipleScatteringTheta0(
         slab, particle.absolutePdg(), particle.mass(), particle.qOverP(),
         particle.absoluteCharge());
-    double sigma2 = sigma * sigma;
 
     // Gauss distribution, will be sampled with generator
     std::normal_distribution<double> gaussDist(0., 1.);
@@ -56,30 +55,31 @@ struct GaussianMixture {
     // d_0'
     // beta² = (p/E)² = p²/(p² + m²) = 1/(1 + (m/p)²)
     // 1/beta² = 1 + (m/p)²
-    double mOverP = particle.mass() / particle.absoluteMomentum();
-    double beta2inv = 1 + mOverP * mOverP;
-    double dprime = slab.thicknessInX0() * beta2inv;
-    double log_dprime = std::log(dprime);
+    const double mOverP = particle.mass() / particle.absoluteMomentum();
+    const double beta2inv = 1 + mOverP * mOverP;
+    const double dprime = slab.thicknessInX0() * beta2inv;
+    const double log_dprime = std::log(dprime);
     // d_0''
-    double log_dprimeprime =
+    const double log_dprimeprime =
         std::log(std::pow(slab.material().Z(), 2.0 / 3.0) * dprime);
 
     // get epsilon
-    double epsilon =
+    const double epsilon =
         log_dprimeprime < 0.5
             ? gausMixEpsilon_a0 + gausMixEpsilon_a1 * log_dprimeprime +
                   gausMixEpsilon_a2 * log_dprimeprime * log_dprimeprime
             : gausMixEpsilon_b0 + gausMixEpsilon_b1 * log_dprimeprime +
                   gausMixEpsilon_b2 * log_dprimeprime * log_dprimeprime;
 
     // the standard sigma
-    double sigma1square = gausMixSigma1_a0 + gausMixSigma1_a1 * log_dprime +
-                          gausMixSigma1_a2 * log_dprime * log_dprime;
+    const double sigma1square = gausMixSigma1_a0 +
+                                gausMixSigma1_a1 * log_dprime +
+                                gausMixSigma1_a2 * log_dprime * log_dprime;
 
+    double sigma2 = Acts::square(sigma);
     // G4 optimised / native double Gaussian model
     if (optGaussianMixtureG4) {
-      sigma2 = 225. * dprime /
-               (particle.absoluteMomentum() * particle.absoluteMomentum());
+      sigma2 = 225. * dprime / Acts::square(particle.absoluteMomentum());
     }
     // throw the random number core/tail
     if (uniformDist(generator) < epsilon) {