Physics/NuclearInteraction/NuclearInteraction.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/NuclearInteraction/NuclearInteraction.hpp"
0010
0011 #include <algorithm>
0012 #include <cstddef>
0013 #include <cstdint>
0014 #include <iterator>
0015 #include <memory>
0016 #include <type_traits>
0017
0018 namespace ActsFatras {
0019
0020 const detail::NuclearInteractionParameters& NuclearInteraction::findParameters(
0021     double rnd,
0022     const detail::NuclearInteractionParametrisation& parametrisation,
0023     float particleMomentum) const {
0024   // Return lowest/highest if momentum outside the boundary
0025   if (particleMomentum <= parametrisation.front().first) {
0026     return parametrisation.front().second;
0027   }
0028   if (particleMomentum >= parametrisation.back().first) {
0029     return parametrisation.back().second;
0030   }
0031
0032   // Find the two neighbouring parametrisations
0033   const auto lowerBound = std::lower_bound(
0034       parametrisation.begin(), parametrisation.end(), particleMomentum,
0035       [](const std::pair<const float,
0036                          ActsFatras::detail::NuclearInteractionParameters>&
0037              params,
0038          const float mom) { return params.first < mom; });
0039   const float momentumUpperNeighbour = lowerBound->first;
0040   const float momentumLowerNeighbour = std::prev(lowerBound, 1)->first;
0041
0042   // Pick one randomly
0043   const float weight = (momentumUpperNeighbour - particleMomentum) /
0044                        (momentumUpperNeighbour - momentumLowerNeighbour);
0045   return (rnd < weight) ? std::prev(lowerBound, 1)->second : lowerBound->second;
0046 }  // namespace ActsFatras
0047
0048 unsigned int NuclearInteraction::sampleDiscreteValues(
0049     double rnd,
0050     const detail::NuclearInteractionParameters::CumulativeDistribution&
0051         distribution) const {
0052   // Fast exit
0053   if (distribution.second.empty()) {
0054     return 0;
0055   }
0056
0057   // Find the bin
0058   const uint32_t int_rnd = static_cast<uint32_t>(UINT32_MAX * rnd);
0059   const auto it = std::upper_bound(distribution.second.begin(),
0060                                    distribution.second.end(), int_rnd);
0061   std::size_t iBin =
0062       std::min((std::size_t)std::distance(distribution.second.begin(), it),
0063                distribution.second.size() - 1);
0064
0065   // Return the corresponding bin
0066   return static_cast<unsigned int>(distribution.first[iBin]);
0067 }
0068
0069 Particle::Scalar NuclearInteraction::sampleContinuousValues(
0070     double rnd,
0071     const detail::NuclearInteractionParameters::CumulativeDistribution&
0072         distribution,
0073     bool interpolate) const {
0074   // Fast exit
0075   if (distribution.second.empty()) {
0076     return std::numeric_limits<Scalar>::infinity();
0077   }
0078
0079   // Find the bin
0080   const uint32_t int_rnd = static_cast<uint32_t>(UINT32_MAX * rnd);
0081   // Fast exit for non-normalised CDFs like interaction probabiltiy
0082   if (int_rnd > distribution.second.back()) {
0083     return std::numeric_limits<Scalar>::infinity();
0084   }
0085   const auto it = std::upper_bound(distribution.second.begin(),
0086                                    distribution.second.end(), int_rnd);
0087   std::size_t iBin =
0088       std::min((std::size_t)std::distance(distribution.second.begin(), it),
0089                distribution.second.size() - 1);
0090
0091   if (interpolate) {
0092     // Interpolate between neighbouring bins and return a diced intermediate
0093     // value
0094     const uint32_t basecont = (iBin > 0 ? distribution.second[iBin - 1] : 0);
0095     const uint32_t dcont = distribution.second[iBin] - basecont;
0096     return distribution.first[iBin] +
0097            (distribution.first[iBin + 1] - distribution.first[iBin]) *
0098                (dcont > 0 ? (int_rnd - basecont) / dcont : 0.5);
0099   } else {
0100     return distribution.first[iBin];
0101   }
0102 }
0103
0104 unsigned int NuclearInteraction::finalStateMultiplicity(
0105     double rnd,
0106     const detail::NuclearInteractionParameters::CumulativeDistribution&
0107         distribution) const {
0108   return sampleDiscreteValues(rnd, distribution);
0109 }
0110
0111 std::pair<ActsFatras::Particle::Scalar, ActsFatras::Particle::Scalar>
0112 NuclearInteraction::globalAngle(ActsFatras::Particle::Scalar phi1,
0113                                 ActsFatras::Particle::Scalar theta1, float phi2,
0114                                 float theta2) const {
0115   // Rotation around the global y-axis
0116   Acts::SquareMatrix3 rotY = Acts::SquareMatrix3::Zero();
0117   rotY(0, 0) = std::cos(theta1);
0118   rotY(0, 2) = std::sin(theta1);
0119   rotY(1, 1) = 1.;
0120   rotY(2, 0) = -std::sin(theta1);
0121   rotY(2, 2) = std::cos(theta1);
0122
0123   // Rotation around the global z-axis
0124   Acts::SquareMatrix3 rotZ = Acts::SquareMatrix3::Zero();
0125   rotZ(0, 0) = std::cos(phi1);
0126   rotZ(0, 1) = -std::sin(phi1);
0127   rotZ(1, 0) = std::sin(phi1);
0128   rotZ(1, 1) = std::cos(phi1);
0129   rotZ(2, 2) = 1.;
0130
0131   // Rotate the direction vector of the second particle
0132   const Acts::Vector3 vector2(std::sin(theta2) * std::cos(phi2),
0133                               std::sin(theta2) * std::sin(phi2),
0134                               std::cos(theta2));
0135   const Acts::Vector3 vectorSum = rotZ * rotY * vector2;
0136
0137   // Calculate the global angles
0138   const float theta = std::acos(vectorSum.z() / vectorSum.norm());
0139   const float phi = std::atan2(vectorSum.y(), vectorSum.x());
0140
0141   return std::make_pair(phi, theta);
0142 }
0143
0144 bool NuclearInteraction::match(const Acts::ActsDynamicVector& momenta,
0145                                const Acts::ActsDynamicVector& invariantMasses,
0146                                float parametrizedMomentum) const {
0147   const unsigned int size = momenta.size();
0148   for (unsigned int i = 0; i < size; i++) {
0149     // Calculate the invariant masses
0150     const float momentum = momenta[i];
0151     const float invariantMass = invariantMasses[i];
0152     const float p1p2 = 2. * momentum * parametrizedMomentum;
0153     const float costheta = 1. - invariantMass * invariantMass / p1p2;
0154
0155     // Abort if an angle cannot be calculated
0156     if (std::abs(costheta) > 1) {
0157       return false;
0158     }
0159   }
0160   return true;
0161 }
0162 }  // namespace ActsFatras