Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /acts/Core/src/Definitions/ParticleData.cpp was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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