File indexing completed on 2026-03-31 07:45:40
0001
0002
0003
0004
0005
0006
0007
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
0031
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
0046
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 }
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
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
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 }