Physics/NuclearInteraction/NuclearInteraction.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020-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 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "Acts/Definitions/PdgParticle.hpp"
0014 #include "Acts/Material/MaterialSlab.hpp"
0015 #include "Acts/Utilities/UnitVectors.hpp"
0016 #include "Acts/Utilities/VectorHelpers.hpp"
0017 #include "ActsFatras/EventData/Barcode.hpp"
0018 #include "ActsFatras/EventData/Particle.hpp"
0019 #include "ActsFatras/EventData/ProcessType.hpp"
0020 #include "ActsFatras/Physics/NuclearInteraction/NuclearInteractionParameters.hpp"
0021
0022 #include <any>
0023 #include <cmath>
0024 #include <iterator>
0025 #include <limits>
0026 #include <optional>
0027 #include <random>
0028 #include <utility>
0029 #include <vector>
0030
0031 namespace ActsFatras {
0032
0033 /// This class provides a parametrised nuclear interaction. The thereby
0034 /// required parametrisation needs to be set and is not provided by default.
0035 ///
0036 /// @note This class differs between two different processes labelled as nuclear
0037 /// interaction. Either the initial particle survives (soft) or it gets
0038 /// destroyed (hard) by this process.
0039 struct NuclearInteraction {
0040   using Scalar = Particle::Scalar;
0041   /// The storage of the parameterisation
0042   detail::MultiParticleNuclearInteractionParametrisation
0043       multiParticleParameterisation;
0044   /// The number of trials to match momenta and inveriant masses
0045   //~ unsigned int nMatchingTrials = std::numeric_limits<unsigned int>::max();
0046   unsigned int nMatchingTrials = 100;
0047   unsigned int nMatchingTrialsTotal = 1000;
0048
0049   /// This method evaluates the nuclear interaction length L0.
0050   ///
0051   /// @tparam generator_t The random number generator type
0052   /// @param [in, out] generator The random number generator
0053   /// @param [in] particle The ingoing particle
0054   ///
0055   /// @return valid X0 limit and no limit on L0
0056   template <typename generator_t>
0057   std::pair<Scalar, Scalar> generatePathLimits(generator_t& generator,
0058                                                const Particle& particle) const {
0059     // Fast exit: No paramtrization provided
0060     if (multiParticleParameterisation.empty()) {
0061       return std::make_pair(std::numeric_limits<Scalar>::infinity(),
0062                             std::numeric_limits<Scalar>::infinity());
0063     }
0064     // Find the parametrisation that corresponds to the particle type
0065     for (const auto& particleParametrisation : multiParticleParameterisation) {
0066       if (particleParametrisation.first == particle.pdg()) {
0067         std::uniform_real_distribution<double> uniformDistribution{0., 1.};
0068
0069         // Get the parameters
0070         const detail::NuclearInteractionParametrisation&
0071             singleParticleParametrisation = particleParametrisation.second;
0072         const detail::NuclearInteractionParameters& parametrisation =
0073             findParameters(uniformDistribution(generator),
0074                            singleParticleParametrisation,
0075                            particle.absoluteMomentum());
0076
0077         // Set the L0 limit if not done already
0078         const auto& distribution =
0079             parametrisation.nuclearInteractionProbability;
0080         auto limits =
0081             std::make_pair(std::numeric_limits<Scalar>::infinity(),
0082                            sampleContinuousValues(
0083                                uniformDistribution(generator), distribution));
0084         return limits;
0085       }
0086     }
0087     return std::make_pair(std::numeric_limits<Scalar>::infinity(),
0088                           std::numeric_limits<Scalar>::infinity());
0089   }
0090
0091   /// This method performs a nuclear interaction.
0092   ///
0093   /// @tparam generator_t The random number generator type
0094   /// @param [in, out] generator The random number generator
0095   /// @param [in, out] particle The ingoing particle
0096   /// @param [out] generated Additional generated particles
0097   ///
0098   /// @return True if the particle was killed, false otherwise
0099   template <typename generator_t>
0100   bool run(generator_t& generator, Particle& particle,
0101            std::vector<Particle>& generated) const {
0102     // Fast exit: No paramtrization provided
0103     if (multiParticleParameterisation.empty()) {
0104       return false;
0105     }
0106
0107     // Find the parametrisation that corresponds to the particle type
0108     for (const auto& particleParametrisation : multiParticleParameterisation) {
0109       if (particleParametrisation.first == particle.pdg()) {
0110         std::uniform_real_distribution<double> uniformDistribution{0., 1.};
0111
0112         // Get the parameters
0113         const detail::NuclearInteractionParametrisation&
0114             singleParticleParametrisation = particleParametrisation.second;
0115         const detail::NuclearInteractionParameters& parametrisation =
0116             findParameters(uniformDistribution(generator),
0117                            singleParticleParametrisation,
0118                            particle.absoluteMomentum());
0119
0120         std::normal_distribution<double> normalDistribution{0., 1.};
0121         // Dice the interaction type
0122         const bool interactSoft =
0123             softInteraction(normalDistribution(generator),
0124                             parametrisation.softInteractionProbability);
0125
0126         // Get the final state multiplicity
0127         const unsigned int multiplicity = finalStateMultiplicity(
0128             uniformDistribution(generator),
0129             interactSoft ? parametrisation.softMultiplicity
0130                          : parametrisation.hardMultiplicity);
0131
0132         // Get the parameters for the kinematics
0133         const std::vector<detail::NuclearInteractionParameters::
0134                               ParametersWithFixedMultiplicity>&
0135             parametrisationOfType =
0136                 interactSoft ? parametrisation.softKinematicParameters
0137                              : parametrisation.hardKinematicParameters;
0138         const detail::NuclearInteractionParameters::
0139             ParametersWithFixedMultiplicity& parametrisationOfMultiplicity =
0140                 parametrisationOfType[multiplicity];
0141         if (!parametrisationOfMultiplicity.validParametrisation) {
0142           return false;
0143         }
0144
0145         // Get the kinematics
0146         const auto kinematics = sampleKinematics(
0147             generator, parametrisationOfMultiplicity, parametrisation.momentum);
0148         if (!kinematics.has_value()) {
0149           return run(generator, particle, generated);
0150         }
0151
0152         // Get the particle types
0153         const std::vector<int> pdgIds =
0154             samplePdgIds(generator, parametrisation.pdgMap, multiplicity,
0155                          particle.pdg(), interactSoft);
0156
0157         // Construct the particles
0158         const auto particles = convertParametersToParticles(
0159             generator, pdgIds, kinematics->first, kinematics->second, particle,
0160             parametrisation.momentum, interactSoft);
0161
0162         // Kill the particle in a hard process
0163         if (!interactSoft) {
0164           particle.setAbsoluteMomentum(0);
0165         }
0166
0167         generated.insert(generated.end(), particles.begin(), particles.end());
0168         return !interactSoft;
0169       }
0170     }
0171     // Fast exit if no parametrisation for the particle was provided
0172     return false;
0173   }
0174
0175  private:
0176   /// Retrieves the parametrisation for the particle
0177   ///
0178   /// @param [in] rnd A random number
0179   /// @param [in] parametrisation The storage of parametrisations
0180   /// @param [in] particleMomentum The particles momentum
0181   ///
0182   /// @return The parametrisation
0183   const detail::NuclearInteractionParameters& findParameters(
0184       double rnd,
0185       const detail::NuclearInteractionParametrisation& parametrisation,
0186       float particleMomentum) const;
0187
0188   /// Estimates the interaction type
0189   ///
0190   /// @param [in] rnd Random number
0191   /// @param [in] probability The probability for a soft interaction
0192   ///
0193   /// @return True if a soft interaction occurs
0194   inline bool softInteraction(double rnd, float probability) const {
0195     return rnd <= probability;
0196   }
0197
0198   /// Evaluates the multiplicity of the final state
0199   ///
0200   /// @param [in] rnd Random number
0201   /// @param [in] distribution The multiplicity distribution
0202   ///
0203   /// @return The final state multiplicity
0204   unsigned int finalStateMultiplicity(
0205       double rnd,
0206       const detail::NuclearInteractionParameters::CumulativeDistribution&
0207           distribution) const;
0208
0209   /// Evaluates the final state PDG IDs
0210   ///
0211   /// @tparam generator_t The random number generator type
0212   /// @param [in, out] generator The random number generator
0213   /// @param [in] pdgMap The branching probability map
0214   /// @param [in] multiplicity The final state multiplicity
0215   /// @param [in] particlePdg The PDG ID of the initial particle
0216   /// @param [in] soft Treat it as soft or hard nuclear interaction
0217   ///
0218   /// @return Vector containing the PDG IDs
0219   template <typename generator_t>
0220   std::vector<int> samplePdgIds(
0221       generator_t& generator,
0222       const detail::NuclearInteractionParameters::PdgMap& pdgMap,
0223       unsigned int multiplicity, int particlePdg, bool soft) const;
0224
0225   /// Evaluates the final state invariant masses
0226   ///
0227   /// @tparam generator_t The random number generator type
0228   /// @param [in, out] generator The random number generator
0229   /// @param [in] parametrisation Parametrisation of kinematic properties
0230   ///
0231   /// @return Vector containing the invariant masses
0232   template <typename generator_t>
0233   Acts::ActsDynamicVector sampleInvariantMasses(
0234       generator_t& generator,
0235       const detail::NuclearInteractionParameters::
0236           ParametersWithFixedMultiplicity& parametrisation) const;
0237
0238   /// Evaluates the final state momenta
0239   ///
0240   /// @tparam generator_t The random number generator type
0241   /// @param [in, out] generator The random number generator
0242   /// @param [in] parametrisation Parametrisation of kinematic properties
0243   /// @param [in] initialMomentum The initial momentum
0244   ///
0245   /// @return Vector containing the momenta
0246   template <typename generator_t>
0247   Acts::ActsDynamicVector sampleMomenta(
0248       generator_t& generator,
0249       const detail::NuclearInteractionParameters::
0250           ParametersWithFixedMultiplicity& parametrisation,
0251       float initialMomentum) const;
0252
0253   /// Tests whether the final state momenta and invariant masses are
0254   /// matching to each other to allow the evaluation of particle directions.
0255   ///
0256   /// @param [in] momenta The final state momenta
0257   /// @param [in] invariantMasses The final state invariant masses
0258   /// @param [in] parametrizedMomentum The momentum of the parametrized particle
0259   ///
0260   /// @return Decision whether the parameters can be matched to each other or
0261   /// not.
0262   bool match(const Acts::ActsDynamicVector& momenta,
0263              const Acts::ActsDynamicVector& invariantMasses,
0264              float parametrizedMomentum) const;
0265
0266   /// This method samples the kinematics of the final state particles
0267   ///
0268   /// @tparam generator_t The random number generator type
0269   /// @param [in, out] generator The random number generator
0270   /// @param [in] parameters The parametrisation
0271   /// @param [in] momentum The momentum of the parametrisation
0272   ///
0273   /// @return The final state momenta and invariant masses
0274   template <typename generator_t>
0275   std::optional<std::pair<Acts::ActsDynamicVector, Acts::ActsDynamicVector>>
0276   sampleKinematics(generator_t& generator,
0277                    const detail::NuclearInteractionParameters::
0278                        ParametersWithFixedMultiplicity& parameters,
0279                    float momentum) const;
0280
0281   /// Converts relative angles to absolute angles wrt the global
0282   /// coordinate system.
0283   /// @note It is assumed that the angles of the first particle are provided in
0284   /// the context of the global coordinate system whereas the angles of the
0285   /// second particle are provided relatively to the first particle.
0286   ///
0287   /// @param [in] phi1 The azimuthal angle of the first particle
0288   /// @param [in] theta1 The polar angle of the first particle
0289   /// @param [in] phi2 The azimuthal angle of the second particle
0290   /// @param [in] theta2 The polar angle of the second particle
0291   ///
0292   /// @return Azimuthal and polar angle of the second particle in the global
0293   /// coordinate system
0294   std::pair<ActsFatras::Particle::Scalar, ActsFatras::Particle::Scalar>
0295   globalAngle(ActsFatras::Particle::Scalar phi1,
0296               ActsFatras::Particle::Scalar theta1, float phi2,
0297               float theta2) const;
0298
0299   /// Converter from sampled numbers to a vector of particles
0300   ///
0301   /// @tparam generator_t The random number generator type
0302   /// @param [in, out] generator The random number generator
0303   /// @param [in] pdgId The PDG IDs
0304   /// @param [in] momenta The momenta
0305   /// @param [in] invariantMasses The invariant masses
0306   /// @param [in] initialParticle The initial particle
0307   /// @param [in] parametrizedMomentum Momentum of the parametrisation
0308   /// @param [in] soft Treat it as soft or hard nuclear interaction
0309   ///
0310   /// @return Vector containing the final state particles
0311   template <typename generator_t>
0312   std::vector<Particle> convertParametersToParticles(
0313       generator_t& generator, const std::vector<int>& pdgId,
0314       const Acts::ActsDynamicVector& momenta,
0315       const Acts::ActsDynamicVector& invariantMasses, Particle& initialParticle,
0316       float parametrizedMomentum, bool soft) const;
0317
0318   /// This function performs an inverse sampling to provide a discrete
0319   /// value from a distribution.
0320   ///
0321   /// @param [in] rnd A random number in [0,1]
0322   /// @param [in] distribution The distribution to sample from
0323   ///
0324   /// @return The sampled value
0325   unsigned int sampleDiscreteValues(
0326       double rnd,
0327       const detail::NuclearInteractionParameters::CumulativeDistribution&
0328           distribution) const;
0329
0330   /// This function performs an inverse sampling to provide a continuous
0331   /// value from a distribition.
0332   ///
0333   /// @param [in] rnd A random number in [0,1]
0334   /// @param [in] distribution The distribution to sample from
0335   /// @param [in] interpolate Flag to steer whether an interpolation between
0336   /// neighbouring bins should be performed instead of a bin lookup
0337   ///
0338   /// @return The sampled value
0339   Scalar sampleContinuousValues(
0340       double rnd,
0341       const detail::NuclearInteractionParameters::CumulativeDistribution&
0342           distribution,
0343       bool interpolate = false) const;
0344 };
0345
0346 template <typename generator_t>
0347 std::vector<int> NuclearInteraction::samplePdgIds(
0348     generator_t& generator,
0349     const detail::NuclearInteractionParameters::PdgMap& pdgMap,
0350     unsigned int multiplicity, int particlePdg, bool soft) const {
0351   // Fast exit in case of no final state particles
0352   if (multiplicity == 0) {
0353     return {};
0354   }
0355
0356   // The final state PDG IDs
0357   std::vector<int> pdgIds;
0358   pdgIds.reserve(multiplicity);
0359
0360   std::uniform_real_distribution<float> uniformDistribution{0., 1.};
0361
0362   // Find the producers probability distribution
0363   auto citProducer = pdgMap.cbegin();
0364   while (citProducer->first != particlePdg && citProducer != pdgMap.end()) {
0365     citProducer++;
0366   }
0367
0368   const std::vector<std::pair<int, float>>& mapInitial = citProducer->second;
0369   // Set the first particle depending on the interaction type
0370   if (soft) {
0371     // Store the initial particle if the interaction is soft
0372     pdgIds.push_back(particlePdg);
0373   } else {
0374     // Otherwise dice the particle
0375     const float rndInitial = uniformDistribution(generator);
0376
0377     pdgIds.push_back(
0378         std::lower_bound(mapInitial.begin(), mapInitial.end(), rndInitial,
0379                          [](const std::pair<int, float>& element,
0380                             float random) { return element.second < random; })
0381             ->first);
0382   }
0383
0384   // Set the remaining particles
0385   for (unsigned int i = 1; i < multiplicity; i++) {
0386     // Find the producers probability distribution from the last produced
0387     // particle
0388     citProducer = pdgMap.cbegin();
0389     while (citProducer->first != pdgIds[i - 1] && citProducer != pdgMap.end()) {
0390       citProducer++;
0391     }
0392
0393     // Set the next particle
0394     const std::vector<std::pair<int, float>>& map = citProducer->second;
0395     const float rnd = uniformDistribution(generator);
0396     pdgIds.push_back(
0397         std::lower_bound(map.begin(), map.end(), rnd,
0398                          [](const std::pair<int, float>& element,
0399                             float random) { return element.second < random; })
0400             ->first);
0401   }
0402   return pdgIds;
0403 }
0404
0405 template <typename generator_t>
0406 Acts::ActsDynamicVector NuclearInteraction::sampleInvariantMasses(
0407     generator_t& generator,
0408     const detail::NuclearInteractionParameters::ParametersWithFixedMultiplicity&
0409         parametrisation) const {
0410   // The resulting vector
0411   Acts::ActsDynamicVector parameters;
0412   const unsigned int size = parametrisation.eigenvaluesInvariantMass.size();
0413   parameters.resize(size);
0414
0415   // Sample in the eigenspace
0416   for (unsigned int i = 0; i < size; i++) {
0417     float variance = parametrisation.eigenvaluesInvariantMass[i];
0418     std::normal_distribution<Acts::ActsScalar> dist{
0419         parametrisation.meanInvariantMass[i], sqrtf(variance)};
0420     parameters[i] = dist(generator);
0421   }
0422   // Transform to multivariate normal distribution
0423   parameters = parametrisation.eigenvectorsInvariantMass * parameters;
0424
0425   // Perform the inverse sampling from the distributions
0426   for (unsigned int i = 0; i < size; i++) {
0427     const double cdf = (std::erff(parameters[i]) + 1) * 0.5;
0428     parameters[i] = sampleContinuousValues(
0429         cdf, parametrisation.invariantMassDistributions[i]);
0430   }
0431   return parameters;
0432 }
0433
0434 template <typename generator_t>
0435 Acts::ActsDynamicVector NuclearInteraction::sampleMomenta(
0436     generator_t& generator,
0437     const detail::NuclearInteractionParameters::ParametersWithFixedMultiplicity&
0438         parametrisation,
0439     float initialMomentum) const {
0440   // The resulting vector
0441   Acts::ActsDynamicVector parameters;
0442   const unsigned int size = parametrisation.eigenvaluesMomentum.size();
0443   parameters.resize(size);
0444
0445   // Sample in the eigenspace
0446   for (unsigned int i = 0; i < size; i++) {
0447     float variance = parametrisation.eigenvaluesMomentum[i];
0448     std::normal_distribution<Acts::ActsScalar> dist{
0449         parametrisation.meanMomentum[i], sqrtf(variance)};
0450     parameters[i] = dist(generator);
0451   }
0452
0453   // Transform to multivariate normal distribution
0454   parameters = parametrisation.eigenvectorsMomentum * parameters;
0455
0456   // Perform the inverse sampling from the distributions
0457   for (unsigned int i = 0; i < size; i++) {
0458     const float cdf = (std::erff(parameters[i]) + 1) * 0.5;
0459     parameters[i] =
0460         sampleContinuousValues(cdf, parametrisation.momentumDistributions[i]);
0461   }
0462
0463   // Scale the momenta
0464   Acts::ActsDynamicVector momenta = parameters.head(size - 1);
0465   const float sum = momenta.sum();
0466   const float scale = parameters.template tail<1>()(0, 0) / sum;
0467   momenta *= scale * initialMomentum;
0468   return momenta;
0469 }
0470
0471 template <typename generator_t>
0472 std::optional<std::pair<Acts::ActsDynamicVector, Acts::ActsDynamicVector>>
0473 NuclearInteraction::sampleKinematics(
0474     generator_t& generator,
0475     const detail::NuclearInteractionParameters::ParametersWithFixedMultiplicity&
0476         parameters,
0477     float momentum) const {
0478   unsigned int trials = 0;
0479   Acts::ActsDynamicVector invariantMasses =
0480       sampleInvariantMasses(generator, parameters);
0481   Acts::ActsDynamicVector momenta =
0482       sampleMomenta(generator, parameters, momentum);
0483   // Repeat momentum evaluation until the parameters match
0484   while (!match(momenta, invariantMasses, momentum)) {
0485     if (trials == nMatchingTrialsTotal) {
0486       return std::nullopt;
0487     }
0488     // Re-sampole invariant masses if no fitting momenta were found
0489     if (trials++ % nMatchingTrials == 0) {
0490       invariantMasses = sampleInvariantMasses(generator, parameters);
0491     } else {
0492       momenta = sampleMomenta(generator, parameters, momentum);
0493     }
0494   }
0495   return std::make_pair(momenta, invariantMasses);
0496 }
0497
0498 template <typename generator_t>
0499 std::vector<Particle> NuclearInteraction::convertParametersToParticles(
0500     generator_t& generator, const std::vector<int>& pdgId,
0501     const Acts::ActsDynamicVector& momenta,
0502     const Acts::ActsDynamicVector& invariantMasses, Particle& initialParticle,
0503     float parametrizedMomentum, bool soft) const {
0504   std::uniform_real_distribution<double> uniformDistribution{0., 1.};
0505   const auto& initialDirection = initialParticle.direction();
0506   const double phi = Acts::VectorHelpers::phi(initialDirection);
0507   const double theta = Acts::VectorHelpers::theta(initialDirection);
0508   const unsigned int size = momenta.size();
0509
0510   std::vector<Particle> result;
0511   result.reserve(size);
0512
0513   // Build the particles
0514   for (unsigned int i = 0; i < size; i++) {
0515     const float momentum = momenta[i];
0516     const float invariantMass = invariantMasses[i];
0517     const float p1p2 = 2. * momentum * parametrizedMomentum;
0518     const float costheta = 1. - invariantMass * invariantMass / p1p2;
0519
0520     const auto phiTheta =
0521         globalAngle(phi, theta, uniformDistribution(generator) * 2. * M_PI,
0522                     std::acos(costheta));
0523     const auto direction =
0524         Acts::makeDirectionFromPhiTheta(phiTheta.first, phiTheta.second);
0525
0526     Particle p = Particle(initialParticle.particleId().makeDescendant(i),
0527                           static_cast<Acts::PdgParticle>(pdgId[i]));
0528     p.setProcess(ProcessType::eNuclearInteraction)
0529         .setPosition4(initialParticle.fourPosition())
0530         .setAbsoluteMomentum(momentum)
0531         .setDirection(direction)
0532         .setReferenceSurface(initialParticle.referenceSurface());
0533
0534     // Store the particle
0535     if (i == 0 && soft) {
0536       initialParticle = p;
0537     } else {
0538       result.push_back(std::move(p));
0539     }
0540   }
0541
0542   return result;
0543 }
0544 }  // namespace ActsFatras