Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-03-31 07:45:40

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 "Acts/Definitions/ParticleData.hpp"
0010 
0011 #include "Acts/Definitions/PdgParticle.hpp"
0012 
0013 #include <cassert>
0014 #include <cmath>
0015 #include <cstdint>
0016 #include <optional>
0017 #include <ostream>
0018 #include <utility>
0019 
0020 #include "ParticleDataTable.hpp"
0021 
0022 namespace {
0023 
0024 enum class Type {
0025   Charge,
0026   Mass,
0027   Name,
0028 };
0029 
0030 // Template function depending on the PDG and *type* of data
0031 // so we can use statically cached values
0032 template <typename T, Acts::PdgParticle pdg, Type type>
0033 std::optional<T> findCachedImpl(const std::map<std::int32_t, T>& map) {
0034   const static std::optional<T> value = [&map]() -> std::optional<T> {
0035     const auto it = map.find(static_cast<std::int32_t>(pdg));
0036     if (it == map.end()) {
0037       return std::nullopt;
0038     }
0039     return it->second;
0040   }();
0041 
0042   return value;
0043 }
0044 
0045 // Cache lookup for particle data
0046 // Uses a switch statement to map the PDG code to the correct cached value
0047 template <typename T, Type type>
0048 std::optional<T> findCached(Acts::PdgParticle pdg,
0049                             const std::map<std::int32_t, T>& map) {
0050   using enum Acts::PdgParticle;
0051   switch (pdg) {
0052     case eElectron:
0053       return findCachedImpl<T, eElectron, type>(map);
0054     case ePositron:
0055       return findCachedImpl<T, ePositron, type>(map);
0056     case eMuon:
0057       return findCachedImpl<T, eMuon, type>(map);
0058     case eAntiMuon:
0059       return findCachedImpl<T, eAntiMuon, type>(map);
0060     case eTau:
0061       return findCachedImpl<T, eTau, type>(map);
0062     case eAntiTau:
0063       return findCachedImpl<T, eAntiTau, type>(map);
0064     case eGamma:
0065       return findCachedImpl<T, eGamma, type>(map);
0066     case ePionZero:
0067       return findCachedImpl<T, ePionZero, type>(map);
0068     case ePionPlus:
0069       return findCachedImpl<T, ePionPlus, type>(map);
0070     case ePionMinus:
0071       return findCachedImpl<T, ePionMinus, type>(map);
0072     case eKaonPlus:
0073       return findCachedImpl<T, eKaonPlus, type>(map);
0074     case eKaonMinus:
0075       return findCachedImpl<T, eKaonMinus, type>(map);
0076     case eNeutron:
0077       return findCachedImpl<T, eNeutron, type>(map);
0078     case eAntiNeutron:
0079       return findCachedImpl<T, eAntiNeutron, type>(map);
0080     case eProton:
0081       return findCachedImpl<T, eProton, type>(map);
0082     case eAntiProton:
0083       return findCachedImpl<T, eAntiProton, type>(map);
0084     case eLead:
0085       return findCachedImpl<T, eLead, type>(map);
0086     case eKaon0Short:
0087       return findCachedImpl<T, eKaon0Short, type>(map);
0088     case eLambda0:
0089       return findCachedImpl<T, eLambda0, type>(map);
0090     default:
0091       return std::nullopt;
0092   }
0093 }
0094 
0095 }  // namespace
0096 
0097 std::optional<float> Acts::findCharge(Acts::PdgParticle pdg) {
0098   if (auto cached = findCached<float, Type::Charge>(pdg, particlesMapCharge());
0099       cached) {
0100     return cached;
0101   }
0102 
0103   const auto it = particlesMapCharge().find(pdg);
0104   if (it == particlesMapCharge().end()) {
0105     if (isNucleus(pdg)) {
0106       return Acts::findChargeOfNucleus(pdg);
0107     } else {
0108       return std::nullopt;
0109     }
0110   }
0111   return it->second;
0112 }
0113 
0114 float Acts::findChargeOfNucleus(Acts::PdgParticle pdg) {
0115   if (!isNucleus(pdg)) {
0116     throw std::invalid_argument("PDG must represent a nucleus");
0117   }
0118   auto pdgGround = makeNucleusGroundState(pdg);
0119   const auto it = particlesMapCharge().find(pdgGround);
0120   if (it == particlesMapCharge().end()) {
0121     // extract charge from PDG number
0122     auto ZandA = extractNucleusZandA(pdg);
0123     return static_cast<float>(ZandA.first);
0124   }
0125   return it->second;
0126 }
0127 
0128 std::optional<float> Acts::findMass(Acts::PdgParticle pdg) {
0129   if (auto cached = findCached<float, Type::Mass>(pdg, particlesMapMass());
0130       cached) {
0131     return cached;
0132   }
0133 
0134   const auto it = particlesMapMass().find(pdg);
0135   if (it == particlesMapMass().end()) {
0136     if (isNucleus(pdg)) {
0137       return Acts::findMassOfNucleus(pdg);
0138     } else {
0139       return std::nullopt;
0140     }
0141   }
0142   return it->second;
0143 }
0144 
0145 float Acts::findMassOfNucleus(Acts::PdgParticle pdg) {
0146   if (!isNucleus(pdg)) {
0147     throw std::invalid_argument("PDG must represent a nucleus");
0148   }
0149   auto pdgGround = makeNucleusGroundState(pdg);
0150   const auto it = particlesMapMass().find(pdgGround);
0151   if (it == particlesMapMass().end()) {
0152     // calculate mass using Bethe-Weizsacker formula
0153     return calculateNucleusMass(pdg);
0154   }
0155   return it->second;
0156 }
0157 
0158 float Acts::calculateNucleusMass(Acts::PdgParticle pdg) {
0159   if (!isNucleus(pdg)) {
0160     throw std::invalid_argument("PDG must represent a nucleus");
0161   }
0162 
0163   auto ZandA = extractNucleusZandA(pdg);
0164   int Z = std::abs(ZandA.first);
0165   int A = ZandA.second;
0166 
0167   float a_Vol = 15.260f * Acts::UnitConstants::MeV;
0168   float a_Surf = 16.267f * Acts::UnitConstants::MeV;
0169   float a_Col = 0.689f * Acts::UnitConstants::MeV;
0170   float a_Sym = 22.209f * Acts::UnitConstants::MeV;
0171   float a_Pair =
0172       (1 - 2 * (Z % 2)) * (1 - A % 2) * 10.076f * Acts::UnitConstants::MeV;
0173 
0174   float massP = 0.938272f * Acts::UnitConstants::GeV;
0175   float massN = 0.939565f * Acts::UnitConstants::GeV;
0176 
0177   float bindEnergy = 0.f;
0178   bindEnergy += a_Vol * A;
0179   bindEnergy -= a_Surf * std::pow(A, 2. / 3.);
0180   bindEnergy -= a_Col * Z * (Z - 1) * std::pow(A, -1. / 3.);
0181   bindEnergy -= a_Sym * std::pow(A - 2 * Z, 2.) / A;
0182   bindEnergy -= a_Pair * std::pow(A, -1. / 2.);
0183 
0184   return massP * Z + massN * (A - Z) - bindEnergy;
0185 }
0186 
0187 std::optional<std::string_view> Acts::findName(Acts::PdgParticle pdg) {
0188   if (auto cached =
0189           findCached<const char* const, Type::Name>(pdg, particlesMapName());
0190       cached) {
0191     return cached;
0192   }
0193 
0194   const auto it = particlesMapName().find(pdg);
0195   if (it == particlesMapName().end()) {
0196     if (isNucleus(pdg)) {
0197       return Acts::findNameOfNucleus(pdg);
0198     } else {
0199       return std::nullopt;
0200     }
0201   }
0202   return it->second;
0203 }
0204 
0205 std::optional<std::string_view> Acts::findNameOfNucleus(Acts::PdgParticle pdg) {
0206   if (!isNucleus(pdg)) {
0207     throw std::invalid_argument("PDG must represent a nucleus");
0208   }
0209 
0210   auto pdgGround = makeNucleusGroundState(pdg);
0211   const auto it = particlesMapName().find(pdgGround);
0212   if (it == particlesMapName().end()) {
0213     return std::nullopt;
0214   }
0215   return it->second;
0216 }
0217 
0218 std::optional<Acts::ParticleData> Acts::findParticleData(PdgParticle pdg) {
0219   auto charge = findCharge(pdg);
0220   auto mass = findMass(pdg);
0221   auto name = findName(pdg);
0222 
0223   if (charge.has_value() && mass.has_value() && name.has_value()) {
0224     return ParticleData{charge.value(), mass.value(), name.value()};
0225   }
0226 
0227   return std::nullopt;
0228 }
0229 
0230 std::ostream& Acts::operator<<(std::ostream& os, Acts::PdgParticle pdg) {
0231   const auto name = Acts::findName(pdg);
0232   if (name) {
0233     os << *name;
0234   } else {
0235     os << static_cast<std::int32_t>(pdg);
0236   }
0237   return os;
0238 }
0239 
0240 std::optional<std::string_view> Acts::pdgToShortAbsString(PdgParticle pdg) {
0241   pdg = makeAbsolutePdgParticle(pdg);
0242   if (pdg == eElectron) {
0243     return "e";
0244   }
0245   if (pdg == eMuon) {
0246     return "mu";
0247   }
0248   if (pdg == eTau) {
0249     return "t";
0250   }
0251   if (pdg == eGamma) {
0252     return "g";
0253   }
0254   if (pdg == ePionZero) {
0255     return "pi0";
0256   }
0257   if (pdg == ePionPlus) {
0258     return "pi";
0259   }
0260   if (pdg == eKaonPlus) {
0261     return "K";
0262   }
0263   if (pdg == eNeutron) {
0264     return "n";
0265   }
0266   if (pdg == eProton) {
0267     return "p";
0268   }
0269   if (pdg == eLead) {
0270     return "lead";
0271   }
0272   if (pdg == eKaon0Short) {
0273     return "Kaon0Short";
0274   }
0275   if (pdg == eLambda0) {
0276     return "Lambda0";
0277   }
0278   return std::nullopt;
0279 }