diff --git a/corsika/detail/modules/sibyll/Decay.inl b/corsika/detail/modules/sibyll/Decay.inl
index e1e3978b29932ecbe2c42cae5d1513205e69d899..eff18c5408e16e3ffa629ad2591cf82f28a6ddd6 100644
--- a/corsika/detail/modules/sibyll/Decay.inl
+++ b/corsika/detail/modules/sibyll/Decay.inl
@@ -90,7 +90,7 @@ namespace corsika::sibyll {
const double gamma = E / m;
- const TimeType t0 = corsika::GetLifetime(vP.GetPID());
+ const TimeType t0 = corsika::lifetime(vP.GetPID());
auto const lifetime = gamma * t0;
const auto mkin =
@@ -123,7 +123,7 @@ namespace corsika::sibyll {
vP.GetMomentum(),
// setting particle mass with Corsika values, may be inconsistent
// with sibyll internal values
- corsika::GetMass(pCode));
+ corsika::mass(pCode));
// remember position
Point const decayPoint = vP.GetPosition();
TimeType const t0 = vP.GetTime();
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index 1aa26e4159e763b5d7aa81fcf5b69c7d5d68ad9f..cbd8b9f3a2826aed6f76b91b36bea8edc8e32073 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -87,13 +87,13 @@ namespace corsika::sibyll {
"Interaction: GetCrossSection: CoM energy outside range for Sibyll!");
}
const double dEcm = CoMenergy / 1_GeV;
- if (corsika::IsNucleus(TargetId)) {
- const int iTarget = corsika::GetNucleusA(TargetId);
+ if (corsika::is_nucleus(TargetId)) {
+ const int iTarget = corsika::nucleus_A(TargetId);
if (iTarget > maxTargetMassNumber_ || iTarget == 0)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
- } else if (TargetId == corsika::Proton::GetCode()) {
+ } else if (TargetId == corsika::Code::Proton) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
} else {
// no interaction in sibyll possible, return infinite cross section? or throw?
@@ -186,7 +186,7 @@ namespace corsika::sibyll {
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< corsika::sibyll::CanInteract(corsikaBeamId) << std::endl;
- if (corsika::IsNucleus(corsikaBeamId)) {
+ if (corsika::is_nucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
}
@@ -276,8 +276,8 @@ namespace corsika::sibyll {
allowed air in atmosphere also contains some Argon.
*/
int targetSibCode = -1;
- if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
- if (targetCode == corsika::Proton::GetCode()) targetSibCode = 1;
+ if (is_nucleus(targetCode)) targetSibCode = nucleus_A(targetCode);
+ if (targetCode == corsika::Code::Proton) targetSibCode = 1;
std::cout << "Interaction: sibyll code: " << targetSibCode << std::endl;
if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
throw std::runtime_error(
diff --git a/corsika/detail/modules/sibyll/NuclearInteraction.inl b/corsika/detail/modules/sibyll/NuclearInteraction.inl
index 0cd03c264c90531829eeb1d30597407346b12a57..d4399233c023f31f021f46361031d3b5f3aa1a36 100644
--- a/corsika/detail/modules/sibyll/NuclearInteraction.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteraction.inl
@@ -50,7 +50,7 @@ namespace corsika::sibyll {
for (unsigned int i = 0; i < GetNEnergyBins(); ++i) {
std::cout << " " << i << " ";
for (auto& n : pNuclei) {
- auto const j = GetNucleusA(n);
+ auto const j = corsika::nucleus_A(n);
std::cout << " " << std::setprecision(5) << std::setw(8)
<< cnucsignuc_.sigma[j - 1][k][i];
}
@@ -85,7 +85,7 @@ namespace corsika::sibyll {
for (auto& ptarg : allElementsInUniverse) {
++k;
std::cout << "NuclearInteraction: init target component: " << ptarg << std::endl;
- const int ib = GetNucleusA(ptarg);
+ const int ib = corsika::nucleus_A(ptarg);
if (!hadronicInteraction_.IsValidTarget(ptarg)) {
std::cout
<< "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id="
@@ -217,7 +217,7 @@ namespace corsika::sibyll {
if (corsikaBeamId != corsika::Code::Nucleus) {
// check if target-style nucleus (enum), these are not allowed as projectile
- if (corsika::IsNucleus(corsikaBeamId))
+ if (corsika::is_nucleus(corsikaBeamId))
throw std::runtime_error(
"NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
"projectiles should use NuclearStackExtension!");
@@ -326,9 +326,7 @@ namespace corsika::sibyll {
"NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
"should use NuclearStackExtension!");
- auto const ProjMass =
- vP.GetNuclearZ() * corsika::Proton::GetMass() +
- (vP.GetNuclearA() - vP.GetNuclearZ()) * corsika::Neutron::GetMass();
+ auto const ProjMass = corsika::nucleus_mass(vP.GetNuclearA(), vP.GetNuclearZ());
std::cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << std::endl;
count_++;
@@ -414,7 +412,7 @@ namespace corsika::sibyll {
// sample target nucleon number
//
// proton stand-in for nucleon
- const auto beamId = corsika::Proton::GetCode();
+ const auto beamId = corsika::Code::Proton;
auto const* const currentNode = vP.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
@@ -447,8 +445,8 @@ namespace corsika::sibyll {
allowed air in atmosphere also contains some Argon.
*/
int kATarget = -1;
- if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
- if (targetCode == corsika::Proton::GetCode()) kATarget = 1;
+ if (corsika::is_nucleus(targetCode)) kATarget = corsika::nucleus_A(targetCode);
+ else if (targetCode == corsika::Code::Proton) kATarget = 1;
std::cout << "NuclearInteraction: nuclib target code: " << kATarget << std::endl;
if (!hadronicInteraction_.IsValidTarget(targetCode))
throw std::runtime_error("target outside range. ");
@@ -459,7 +457,7 @@ namespace corsika::sibyll {
<< std::endl;
// get nucleon-nucleon cross section
// (needed to determine number of nucleon-nucleon scatterings)
- const auto protonId = corsika::Proton::GetCode();
+ const auto protonId = corsika::Code::Proton;
const auto [prodCrossSection, elaCrossSection] =
hadronicInteraction_.GetCrossSection(protonId, protonId, EcmNN);
const double sigProd = prodCrossSection / 1_mb;
@@ -519,10 +517,7 @@ namespace corsika::sibyll {
else
specCode = corsika::Code::Nucleus;
- // TODO: mass of nuclei?
- const HEPMassType mass =
- corsika::Proton::GetMass() * nuclZ +
- (nuclA - nuclZ) * corsika::Neutron::GetMass(); // this neglects binding energy
+ const HEPMassType mass = corsika::nucleus_mass(nuclA, nuclZ);
std::cout << "NuclearInteraction: adding fragment: " << specCode << std::endl;
std::cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << std::endl;
@@ -565,7 +560,7 @@ namespace corsika::sibyll {
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
- const double mass_ratio = corsika::GetMass(elaNucCode) / ProjMass;
+ const double mass_ratio = corsika::mass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
vP.AddSecondary(std::tuple<corsika::Code, si::HEPEnergyType,
@@ -578,7 +573,7 @@ namespace corsika::sibyll {
std::cout << "calculate inelastic nucleon-nucleon interactions.." << std::endl;
for (int j = 0; j < nInelNucleons; ++j) {
// TODO: sample neutron or proton
- auto pCode = corsika::Proton::GetCode();
+ auto pCode = corsika::Code::Proton;
// temporarily add to stack, will be removed after interaction in DoInteraction
std::cout << "inelastic interaction no. " << j << std::endl;
auto inelasticNucleon = vP.AddSecondary(
diff --git a/corsika/media/NuclearComposition.hpp b/corsika/media/NuclearComposition.hpp
index 94eef28088dad629805e71e4e4723e963a90d7c2..f7f4bd5a25391901ff59846afdfbe67798ce6763 100644
--- a/corsika/media/NuclearComposition.hpp
+++ b/corsika/media/NuclearComposition.hpp
@@ -72,10 +72,10 @@ namespace corsika {
, fAvgMassNumber(std::inner_product(
pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
std::plus<double>(), [](auto const compID, auto const fraction) -> double {
- if (IsNucleus(compID)) {
- return GetNucleusA(compID) * fraction;
+ if (is_nucleus(compID)) {
+ return nucleus_A(compID) * fraction;
} else {
- return GetMass(compID) / ConvertSIToHEP(constants::u) * fraction;
+ return mass(compID) / ConvertSIToHEP(constants::u) * fraction;
}
})) {
assert(pComponents.size() == pFractions.size());
diff --git a/corsika/modules/sibyll/Interaction.hpp b/corsika/modules/sibyll/Interaction.hpp
index 387d3e7f2666673b96cfa0d94b5f835545656da6..8f0869e51e47849a644af56248bb734a2243951d 100644
--- a/corsika/modules/sibyll/Interaction.hpp
+++ b/corsika/modules/sibyll/Interaction.hpp
@@ -46,8 +46,7 @@ namespace corsika::sibyll {
HEPEnergyType GetMinEnergyCoM() const { return minEnergyCoM_; }
HEPEnergyType GetMaxEnergyCoM() const { return maxEnergyCoM_; }
bool IsValidTarget(corsika::Code TargetId) const {
- return (corsika::GetNucleusA(TargetId) < maxTargetMassNumber_) &&
- corsika::IsNucleus(TargetId);
+ return corsika::is_nucleus(TargetId) && (corsika::nucleus_A(TargetId) < maxTargetMassNumber_);
}
std::tuple<CrossSectionType, CrossSectionType> GetCrossSection(
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index ac433e4de5724fe94aa83b4ea9284cb1c39be787..c4f8274233ea273caf75e2410a74148534fea267 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -23,372 +23,342 @@
namespace corsika {
-/**
- * @namespace nuclear_extension
- *
- * Add A and Z data to existing stack (currently SuperStupidStack) of particle
- * properties. This is done via inheritance, not via CombinedStack since the nuclear
- * data is stored ONLY when needed (for nuclei) and not for all particles. Thus, this is
- * a new, derived Stack object.
- *
- * Only for Code::Nucleus particles A and Z are stored, not for all
- * normal elementary particles.
- *
- * Thus in your code, make sure to always check <code>
- * particle.GetPID()==Code::Nucleus </code> before attempting to
- * read any nuclear information.
- *
- *
- */
-
-
-/**
- * @class NuclearParticleInterface
- *
- * Define ParticleInterface for NuclearStackExtension Stack derived from
- * ParticleInterface of Inner stack class
- */
-template < template <typename> class InnerParticleInterface, typename StackIteratorInterface>
-struct NuclearParticleInterface : public InnerParticleInterface<StackIteratorInterface> {
-
- typedef InnerParticleInterface<StackIteratorInterface> super_type;
-
-
-
-public:
-
-
- typedef std::tuple<
- corsika::Code, corsika::units::si::HEPEnergyType,
- momentum_type, corsika::Point,
- corsika::units::si::TimeType> particle_data_type;
-
- typedef std::tuple<
- corsika::Code, corsika::units::si::HEPEnergyType,
- momentum_type, corsika::Point,
- corsika::units::si::TimeType,
- unsigned short, unsigned short> altenative_particle_data_type;
-
-
- typedef corsika::Vector<corsika::units::si::hepmomentum_d> momentum_type;
-
- void setParticleData(particle_data_type const& v) {
-
- if (std::get<0>(v) == corsika::Code::Nucleus) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- super_type::setParticleData(v);
- setNucleusRef(-1); // this is not a nucleus
- }
-
- void setParticleData( altenative_particle_data_type const& v)
- {
- const unsigned short A = std::get<5>(v);
- const unsigned short Z = std::get<6>(v);
- if (std::get<0>(v) != corsika::Code::Nucleus || A == 0 || Z == 0) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
- setNucleusRef(super_type::GetStackData().getNucleusNextRef()); // store this nucleus data ref
- setNuclearA(A);
- setNuclearZ(Z);
- super_type::setParticleData(particle_data_type{std::get<0>(v), std::get<1>(v),
- std::get<2>(v), std::get<3>(v), std::get<4>(v)});
- }
-
- void setParticleData( super_type& p, particle_data_type const& v)
- {
- if (std::get<0>(v) == corsika::Code::Nucleus) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- super_type::setParticleData(p, particle_data_type{std::get<0>(v), std::get<1>(v),
- std::get<2>(v), std::get<3>(v), std::get<4>(v)});
-
- setNucleusRef(-1); // this is not a nucleus
- }
-
- void setParticleData( super_type& p, altenative_particle_data_type const& v) {
-
- const unsigned short A = std::get<5>(v);
- const unsigned short Z = std::get<6>(v);
-
- if (std::get<0>(v) != corsika::Code::Nucleus || A == 0 || Z == 0) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- setNucleusRef(super_type::GetStackData().getNucleusNextRef()); // store this nucleus data ref
- setNuclearA(A);
- setNuclearZ(Z);
- super_type::setParticleData(p, particle_data_type{std::get<0>(v), std::get<1>(v),
- std::get<2>(v), std::get<3>(v),
- std::get<4>(v)});
- }
-
- std::string as_string() const {
- return fmt::format(
- "{}, nuc({})", super_type::as_string(),
- (isNucleus() ? fmt::format("A={}, Z={}", getNuclearA(), getNuclearZ())
- : "n/a"));
- }
-
- /**
- * @name individual setters
- * @{
- */
- void setNuclearA(const unsigned short vA) {
- super_type::GetStackData().setNuclearA(super_type::GetIndex(), vA);
- }
- void setNuclearZ(const unsigned short vZ) {
- super_type::GetStackData().setNuclearZ(super_type::GetIndex(), vZ);
- }
- /// @}
-
- /**
- * @name individual getters
- * @{
- */
- int getNuclearA() const { return super_type::GetStackData().getNuclearA(super_type::GetIndex()); }
- int getNuclearZ() const { return super_type::GetStackData().getNuclearZ(super_type::GetIndex()); }
- /// @}
-
- /**
- * Overwrite normal GetParticleMass function with nuclear version
- */
- corsika::units::si::HEPMassType getMass() const {
- if (super_type::GetPID() ==
- corsika::Code::Nucleus)
- return corsika::GetNucleusMass(getNuclearA(), getNuclearZ());
- return super_type::getMass();
- }
- /**
- * Overwirte normal GetChargeNumber function with nuclear version
- **/
- int16_t getChargeNumber() const {
- if (super_type::GetPID() ==
- corsika::Code::Nucleus)
- return getNuclearZ();
- return super_type::getChargeNumber();
- }
-
- int getNucleusRef() const {
- return super_type::GetStackData().getNucleusRef(GetIndex());
- } // LCOV_EXCL_LINE
-
-protected:
-
- void setNucleusRef(const int vR) {
- super_type::GetStackData().setNucleusRef(super_type::GetIndex(), vR);
- }
-
- bool isNucleus() const {
- return super_type::GetStackData().isNucleus(super_type::GetIndex());
- }
-};
-
-/**
- * @class NuclearStackExtension
- *
- * Memory implementation of adding nuclear inforamtion to the
- * existing particle stack defined in class InnerStackImpl.
- *
- * Inside the NuclearStackExtension class there is a dedicated
- * fNucleusRef index, where fNucleusRef[i] is referring to the
- * correct A and Z for a specific particle index i. fNucleusRef[i]
- * == -1 means that this is not a nucleus, and a subsequent call to
- * GetNucleusA would produce an exception.
- */
-template <typename InnerStackImpl>
-class NuclearStackExtensionImpl : public InnerStackImpl {
-
- typedef InnerStackImpl super_type;
-
-public:
-
- typedef std::vector<int> nucleus_ref_type;
- typedef std::vector<unsigned short> nuclear_a_type;
- typedef std::vector<unsigned short> nuclear_z_type;
-
-
- NuclearStackExtensionImpl()= default;
-
- NuclearStackExtensionImpl( NuclearStackExtensionImpl<InnerStackImpl> const&)= default;
-
- NuclearStackExtensionImpl( NuclearStackExtensionImpl<InnerStackImpl> &&)= default;
-
- NuclearStackExtensionImpl<InnerStackImpl>&
- operator=( NuclearStackExtensionImpl<InnerStackImpl> const&)= default;
-
- NuclearStackExtensionImpl<InnerStackImpl>&
- operator=( NuclearStackExtensionImpl<InnerStackImpl> &&)= default;
-
- void init() {
- super_type::init();
- }
-
- void dump() {
- super_type::dump();
- }
-
- void clear() {
- super_type::clear();
- nucleusRef_.clear();
- nuclearA_.clear();
- nuclearZ_.clear();
- }
-
- unsigned int getSize() const {
- return nucleusRef_.size();
- }
-
- unsigned int getCapacity() const {
- return nucleusRef_.capacity();
- }
-
- void setNuclearA(const unsigned int i, const unsigned short vA) {
- nuclearA_[getNucleusRef(i)] = vA;
- }
-
- void setNuclearZ(const unsigned int i, const unsigned short vZ) {
- nuclearZ_[getNucleusRef(i)] = vZ;
- }
-
- void setNucleusRef(const unsigned int i, const int v) {
- nucleusRef_[i] = v;
- }
-
- int getNuclearA(const unsigned int i) const {
- return nuclearA_[getNucleusRef(i)];
- }
-
- int getNuclearZ(const unsigned int i) const {
- return nuclearZ_[getNucleusRef(i)];
- }
- // this function will create new storage for Nuclear Properties, and return the
- // reference to it
- int getNucleusNextRef() {
- nuclearA_.push_back(0);
- nuclearZ_.push_back(0);
- return nuclearA_.size() - 1;
- }
-
- int getNucleusRef(const unsigned int i) const {
- if (nucleusRef_[i] >= 0) return nucleusRef_[i];
- std::ostringstream err;
- err << "NuclearStackExtension: no nucleus at ref=" << i;
- throw std::runtime_error(err.str());
- }
-
- bool isNucleus(const unsigned int i) const { return nucleusRef_[i] >= 0; }
-
- /**
- * Function to copy particle at location i1 in stack to i2
- */
- void copy(const unsigned int i1, const unsigned int i2) {
- // index range check
- if (i1 >= getSize() || i2 >= getSize()) {
- std::ostringstream err;
- err << "NuclearStackExtension: trying to access data beyond size of stack!";
- throw std::runtime_error(err.str());
- }
- // copy internal particle data p[i2] = p[i1]
- super_type::copy(i1, i2);
- // check if any of p[i1] or p[i2] was a Code::Nucleus
- const int ref1 = nucleusRef_[i1];
- const int ref2 = nucleusRef_[i2];
- if (ref2 < 0) {
- if (ref1 >= 0) {
- // i1 is nucleus, i2 is not
- nucleusRef_[i2] = getNucleusNextRef();
- nuclearA_[nucleusRef_[i2]] = nuclearA_[ref1];
- nuclearZ_[nucleusRef_[i2]] = nuclearZ_[ref1];
- } else {
- // neither i1 nor i2 are nuclei
- }
- } else {
- if (ref1 >= 0) {
- // both are nuclei, i2 is overwritten with nucleus i1
- // fNucleusRef stays the same, but A and Z data is overwritten
- nuclearA_[ref2] = nuclearA_[ref1];
- nuclearZ_[ref2] = nuclearZ_[ref1];
- } else {
- // i2 is overwritten with non-nucleus i1
- nucleusRef_[i2] = -1; // flag as non-nucleus
- nuclearA_.erase(nuclearA_.cbegin() + ref2); // remove data for i2
- nuclearZ_.erase(nuclearZ_.cbegin() + ref2); // remove data for i2
- const int n = nucleusRef_.size(); // update fNucleusRef: indices above ref2
- // must be decremented by 1
- for (int i = 0; i < n; ++i) {
- if (nucleusRef_[i] > ref2) { nucleusRef_[i] -= 1; }
- }
- }
- }
- }
-
- /**
- * Function to copy particle at location i2 in stack to i1
- */
- void swap(const unsigned int i1, const unsigned int i2) {
- // index range check
- if (i1 >= getSize() || i2 >= getSize()) {
- std::ostringstream err;
- err << "NuclearStackExtension: trying to access data beyond size of stack!";
- throw std::runtime_error(err.str());
- }
- // swap original particle data
- super_type::swap(i1, i2);
- // swap corresponding nuclear reference data
- std::swap(nucleusRef_[i2], nucleusRef_[i1]);
- }
-
- void incrementSize() {
- super_type::incrementSize();
- nucleusRef_.push_back(-1);
- }
-
- void decrementSize() {
- super_type::decrementSize();
- if (nucleusRef_.size() > 0) {
- const int ref = nucleusRef_.back();
- nucleusRef_.pop_back();
- if (ref >= 0) {
- nuclearA_.erase(nuclearA_.begin() + ref);
- nuclearZ_.erase(nuclearZ_.begin() + ref);
- const int n = nucleusRef_.size();
- for (int i = 0; i < n; ++i) {
- if (nucleusRef_[i] >= ref) { nucleusRef_[i] -= 1; }
- }
- }
- }
- }
-
-private:
- /// the actual memory to store particle data
-
- nucleus_ref_type nucleusRef_;
- nuclear_a_type nuclearA_;
- nuclear_z_type nuclearZ_;
-
-}; // end class NuclearStackExtensionImpl
-
-template <typename InnerStack, template <typename> typename _PI>
-using NuclearStackExtension =
- Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, _PI> ;
-
-//
-template <typename StackIter>
-using ExtendedParticleInterfaceType = NuclearParticleInterface<SuperStupidStack, StackIter>;
-
-// the particle data stack with extra nuclear information:
-using ParticleDataStack = NuclearStackExtension<SuperStupidStack, ExtendedParticleInterfaceType>;
-
-
+ /**
+ * @namespace nuclear_extension
+ *
+ * Add A and Z data to existing stack of particle properties.
+ *
+ * Only for Code::Nucleus particles A and Z are stored, not for all
+ * normal elementary particles.
+ *
+ * Thus in your code, make sure to always check <code>
+ * particle.GetPID()==Code::Nucleus </code> before attempting to
+ * read any nuclear information.
+ *
+ *
+ */
+
+ typedef corsika::Vector<hepmomentum_d> MomentumVector;
+
+ namespace nuclear_extension {
+
+ /**
+ * @class NuclearParticleInterface
+ *
+ * Define ParticleInterface for NuclearStackExtension Stack derived from
+ * ParticleInterface of Inner stack class
+ */
+ template <template <typename> typename InnerParticleInterface,
+ typename StackIteratorInterface>
+ class NuclearParticleInterface
+ : public InnerParticleInterface<StackIteratorInterface> {
+
+ protected:
+ using InnerParticleInterface<StackIteratorInterface>::GetStackData;
+ using InnerParticleInterface<StackIteratorInterface>::GetIndex;
+
+ public:
+ void SetParticleData(
+ const std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, TimeType>& v) {
+ if (std::get<0>(v) == corsika::Code::Nucleus) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
+ }
+ InnerParticleInterface<StackIteratorInterface>::SetParticleData(v);
+ SetNucleusRef(-1); // this is not a nucleus
+ }
+
+ void SetParticleData(
+ const std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, TimeType, unsigned short, unsigned short>& v) {
+ const unsigned short A = std::get<5>(v);
+ const unsigned short Z = std::get<6>(v);
+ if (std::get<0>(v) != corsika::Code::Nucleus || A == 0 || Z == 0) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
+ }
+ SetNucleusRef(GetStackData().GetNucleusNextRef()); // store this nucleus data ref
+ SetNuclearA(A);
+ SetNuclearZ(Z);
+ InnerParticleInterface<StackIteratorInterface>::SetParticleData(
+ std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, TimeType>{std::get<0>(v), std::get<1>(v),
+ std::get<2>(v), std::get<3>(v),
+ std::get<4>(v)});
+ }
+
+ void SetParticleData(
+ InnerParticleInterface<StackIteratorInterface>& p,
+ const std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, TimeType>& v) {
+ if (std::get<0>(v) == corsika::Code::Nucleus) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
+ }
+ InnerParticleInterface<StackIteratorInterface>::SetParticleData(
+ p, std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, TimeType>{std::get<0>(v), std::get<1>(v),
+ std::get<2>(v), std::get<3>(v),
+ std::get<4>(v)});
+ SetNucleusRef(-1); // this is not a nucleus
+ }
+
+ void SetParticleData(
+ InnerParticleInterface<StackIteratorInterface>& p,
+ const std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, TimeType, unsigned short, unsigned short>& v) {
+ const unsigned short A = std::get<5>(v);
+ const unsigned short Z = std::get<6>(v);
+ if (std::get<0>(v) != corsika::Code::Nucleus || A == 0 || Z == 0) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
+ }
+ SetNucleusRef(GetStackData().GetNucleusNextRef()); // store this nucleus data ref
+ SetNuclearA(A);
+ SetNuclearZ(Z);
+ InnerParticleInterface<StackIteratorInterface>::SetParticleData(
+ p, std::tuple<corsika::Code, HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, TimeType>{std::get<0>(v), std::get<1>(v),
+ std::get<2>(v), std::get<3>(v),
+ std::get<4>(v)});
+ }
+
+ /**
+ * @name individual setters
+ * @{
+ */
+ void SetNuclearA(const unsigned short vA) {
+ GetStackData().SetNuclearA(GetIndex(), vA);
+ }
+ void SetNuclearZ(const unsigned short vZ) {
+ GetStackData().SetNuclearZ(GetIndex(), vZ);
+ }
+ /// @}
+
+ /**
+ * @name individual getters
+ * @{
+ */
+ int GetNuclearA() const { return GetStackData().GetNuclearA(GetIndex()); }
+ int GetNuclearZ() const { return GetStackData().GetNuclearZ(GetIndex()); }
+ /// @}
+
+ /**
+ * Overwrite normal GetParticleMass function with nuclear version
+ */
+ HEPMassType GetMass() const {
+ if (InnerParticleInterface<StackIteratorInterface>::GetPID() ==
+ corsika::Code::Nucleus)
+ return corsika::nucleus_mass(GetNuclearA(), GetNuclearZ());
+ return InnerParticleInterface<StackIteratorInterface>::GetMass();
+ }
+ /**
+ * Overwirte normal GetChargeNumber function with nuclear version
+ **/
+ int16_t GetChargeNumber() const {
+ if (InnerParticleInterface<StackIteratorInterface>::GetPID() ==
+ corsika::Code::Nucleus)
+ return GetNuclearZ();
+ return InnerParticleInterface<StackIteratorInterface>::GetChargeNumber();
+ }
+
+ int GetNucleusRef() const { return GetStackData().GetNucleusRef(GetIndex()); }
+
+ protected:
+ void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); }
+ };
+
+ /**
+ * @class NuclearStackExtension
+ *
+ * Memory implementation of adding nuclear inforamtion to the
+ * existing particle stack defined in class InnerStackImpl.
+ *
+ * Inside the NuclearStackExtension class there is a dedicated
+ * fNucleusRef index, where fNucleusRef[i] is referring to the
+ * correct A and Z for a specific particle index i. fNucleusRef[i]
+ * == -1 means that this is not a nucleus, and a subsequent call to
+ * GetNucleusA would produce an exception.
+ */
+ template <typename InnerStackImpl>
+ class NuclearStackExtensionImpl : public InnerStackImpl {
+
+ public:
+ void Init() { InnerStackImpl::Init(); }
+ void Dump() { InnerStackImpl::Dump(); }
+
+ void Clear() {
+ InnerStackImpl::Clear();
+ fNucleusRef.clear();
+ fNuclearA.clear();
+ fNuclearZ.clear();
+ }
+
+ unsigned int GetSize() const { return fNucleusRef.size(); }
+ unsigned int GetCapacity() const { return fNucleusRef.capacity(); }
+
+ void SetNuclearA(const unsigned int i, const unsigned short vA) {
+ fNuclearA[GetNucleusRef(i)] = vA;
+ }
+ void SetNuclearZ(const unsigned int i, const unsigned short vZ) {
+ fNuclearZ[GetNucleusRef(i)] = vZ;
+ }
+ void SetNucleusRef(const unsigned int i, const int v) { fNucleusRef[i] = v; }
+
+ int GetNuclearA(const unsigned int i) const { return fNuclearA[GetNucleusRef(i)]; }
+ int GetNuclearZ(const unsigned int i) const { return fNuclearZ[GetNucleusRef(i)]; }
+ // this function will create new storage for Nuclear Properties, and return the
+ // reference to it
+ int GetNucleusNextRef() {
+ fNuclearA.push_back(0);
+ fNuclearZ.push_back(0);
+ return fNuclearA.size() - 1;
+ }
+
+ int GetNucleusRef(const unsigned int i) const {
+ if (fNucleusRef[i] >= 0) return fNucleusRef[i];
+ std::ostringstream err;
+ err << "NuclearStackExtension: no nucleus at ref=" << i;
+ throw std::runtime_error(err.str());
+ }
+
+ /**
+ * Function to copy particle at location i1 in stack to i2
+ */
+ void Copy(const unsigned int i1, const unsigned int i2) {
+ // index range check
+ if (i1 >= GetSize() || i2 >= GetSize()) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: trying to access data beyond size of stack!";
+ throw std::runtime_error(err.str());
+ }
+ // copy internal particle data p[i2] = p[i1]
+ InnerStackImpl::Copy(i1, i2);
+ // check if any of p[i1] or p[i2] was a Code::Nucleus
+ const int ref1 = fNucleusRef[i1];
+ const int ref2 = fNucleusRef[i2];
+ if (ref2 < 0) {
+ if (ref1 >= 0) {
+ // i1 is nucleus, i2 is not
+ fNucleusRef[i2] = GetNucleusNextRef();
+ fNuclearA[fNucleusRef[i2]] = fNuclearA[ref1];
+ fNuclearZ[fNucleusRef[i2]] = fNuclearZ[ref1];
+ } else {
+ // neither i1 nor i2 are nuclei
+ }
+ } else {
+ if (ref1 >= 0) {
+ // both are nuclei, i2 is overwritten with nucleus i1
+ // fNucleusRef stays the same, but A and Z data is overwritten
+ fNuclearA[ref2] = fNuclearA[ref1];
+ fNuclearZ[ref2] = fNuclearZ[ref1];
+ } else {
+ // i2 is overwritten with non-nucleus i1
+ fNucleusRef[i2] = -1; // flag as non-nucleus
+ fNuclearA.erase(fNuclearA.cbegin() + ref2); // remove data for i2
+ fNuclearZ.erase(fNuclearZ.cbegin() + ref2); // remove data for i2
+ const int n = fNucleusRef.size(); // update fNucleusRef: indices above ref2
+ // must be decremented by 1
+ for (int i = 0; i < n; ++i) {
+ if (fNucleusRef[i] > ref2) { fNucleusRef[i] -= 1; }
+ }
+ }
+ }
+ }
+
+ /**
+ * Function to copy particle at location i2 in stack to i1
+ */
+ void Swap(const unsigned int i1, const unsigned int i2) {
+ // index range check
+ if (i1 >= GetSize() || i2 >= GetSize()) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: trying to access data beyond size of stack!";
+ throw std::runtime_error(err.str());
+ }
+ // swap original particle data
+ InnerStackImpl::Swap(i1, i2);
+ // swap corresponding nuclear reference data
+ std::swap(fNucleusRef[i2], fNucleusRef[i1]);
+ }
+
+ void IncrementSize() {
+ InnerStackImpl::IncrementSize();
+ fNucleusRef.push_back(-1);
+ }
+
+ void DecrementSize() {
+ InnerStackImpl::DecrementSize();
+ if (fNucleusRef.size() > 0) {
+ const int ref = fNucleusRef.back();
+ fNucleusRef.pop_back();
+ if (ref >= 0) {
+ fNuclearA.erase(fNuclearA.begin() + ref);
+ fNuclearZ.erase(fNuclearZ.begin() + ref);
+ const int n = fNucleusRef.size();
+ for (int i = 0; i < n; ++i) {
+ if (fNucleusRef[i] >= ref) { fNucleusRef[i] -= 1; }
+ }
+ }
+ }
+ }
+
+ private:
+ /// the actual memory to store particle data
+
+ std::vector<int> fNucleusRef;
+ std::vector<unsigned short> fNuclearA;
+ std::vector<unsigned short> fNuclearZ;
+
+ }; // end class NuclearStackExtensionImpl
+
+ // template<typename StackIteratorInterface>
+ // using NuclearParticleInterfaceType<StackIteratorInterface> =
+ // NuclearParticleInterface< ,StackIteratorInterface>
+
+ // works, but requires stupd _PI class
+ // template<typename SS> using TEST =
+ // NuclearParticleInterface<corsika::super_stupid::SuperStupidStack::PIType,
+ // SS>;
+ template <typename InnerStack, template <typename> typename _PI>
+ using NuclearStackExtension =
+ Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, _PI>;
+
+ // ----
+
+ // I'm dont't manage to do this properly.......
+ /*
+ template<typename TT, typename SS> using TESTi = typename
+ NuclearParticleInterface<TT::template PIType, SS>::ExtendedParticleInterface;
+ template<typename TT, typename SS> using TEST1 = TESTi<TT, SS>;
+ template<typename SS> using TEST2 = TEST1<typename
+ corsika::super_stupid::SuperStupidStack, SS>;
+
+ using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename
+ InnerStack::StackImpl>, TEST2>;
+ */
+ /*
+ // .... this should be fixed ....
+
+ template <typename InnerStack, typename SS=StackIteratorInterface>
+ //using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename
+ InnerStack::StackImpl>, NuclearParticleInterface<typename InnerStack::template PIType,
+ StackIteratorInterface>::ExtendedParticleInterface>; using NuclearStackExtension =
+ Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, TEST1<typename
+ corsika::super_stupid::SuperStupidStack, SS> >;
+
+ //template <typename InnerStack>
+ // using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename
+ InnerStack::StackImpl>, TEST<typename
+ corsika::super_stupid::SuperStupidStack::PIType>>;
+ //using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename
+ InnerStack::StackImpl>, TEST>;
+ */
+
+ } // namespace nuclear_extension
} // namespace corsika
diff --git a/corsika/stack/SuperStupidStack.hpp b/corsika/stack/SuperStupidStack.hpp
index ffa0e4cb20b33a2534624e3d6b0462ddd4931af7..d83dda4f3d8cbf3e2026d5d603c144e1d4eebe95 100644
--- a/corsika/stack/SuperStupidStack.hpp
+++ b/corsika/stack/SuperStupidStack.hpp
@@ -108,8 +108,8 @@ namespace corsika {
corsika::Vector<dimensionless_d> GetDirection() const {
return GetMomentum() / GetEnergy();
}
- HEPMassType GetMass() const { return corsika::GetMass(GetPID()); }
- int16_t GetChargeNumber() const { return corsika::GetChargeNumber(GetPID()); }
+ HEPMassType GetMass() const { return corsika::mass(GetPID()); }
+ int16_t GetChargeNumber() const { return corsika::charge_number(GetPID()); }
///@}
};
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index d203c652336b446160e8c834fe171997ad2f94df..c6666ab41844505182f05a0e18bd347f511bc180 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -23,26 +23,26 @@ using namespace corsika::sibyll;
TEST_CASE("Sibyll", "[processes]") {
SECTION("Sibyll -> Corsika") {
- REQUIRE(Electron::GetCode() ==
+ REQUIRE(Code::Electron ==
corsika::sibyll::ConvertFromSibyll(corsika::sibyll::SibyllCode::Electron));
}
SECTION("Corsika -> Sibyll") {
- REQUIRE(corsika::sibyll::ConvertToSibyll(Electron::GetCode()) ==
+ REQUIRE(corsika::sibyll::ConvertToSibyll(Code::Electron) ==
corsika::sibyll::SibyllCode::Electron);
- REQUIRE(corsika::sibyll::ConvertToSibyllRaw(Proton::GetCode()) == 13);
+ REQUIRE(corsika::sibyll::ConvertToSibyllRaw(Code::Proton) == 13);
}
SECTION("canInteractInSibyll") {
- REQUIRE(corsika::sibyll::CanInteract(Proton::GetCode()));
+ REQUIRE(corsika::sibyll::CanInteract(Code::Proton));
REQUIRE(corsika::sibyll::CanInteract(Code::XiCPlus));
- REQUIRE_FALSE(corsika::sibyll::CanInteract(Electron::GetCode()));
- REQUIRE_FALSE(corsika::sibyll::CanInteract(SigmaC0::GetCode()));
+ REQUIRE_FALSE(corsika::sibyll::CanInteract(Code::Electron));
+ REQUIRE_FALSE(corsika::sibyll::CanInteract(Code::SigmaC0));
- REQUIRE_FALSE(corsika::sibyll::CanInteract(Nucleus::GetCode()));
- REQUIRE_FALSE(corsika::sibyll::CanInteract(Helium::GetCode()));
+ REQUIRE_FALSE(corsika::sibyll::CanInteract(Code::Nucleus));
+ REQUIRE_FALSE(corsika::sibyll::CanInteract(Code::Helium));
}
SECTION("cross-section type") {
@@ -101,7 +101,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::setup::Stack stack;
const HEPEnergyType E0 = 100_GeV;
- HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass());
+ HEPMomentumType P0 = sqrt(E0 * E0 - Proton::mass() * Proton::mass());
auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
corsika::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
@@ -123,7 +123,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
setup::Stack stack;
const HEPEnergyType E0 = 400_GeV;
- HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass());
+ HEPMomentumType P0 = sqrt(E0 * E0 - Proton::mass() * Proton::mass());
auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
corsika::Point pos(cs, 0_m, 0_m, 0_m);
@@ -147,7 +147,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
setup::Stack stack;
const HEPEnergyType E0 = 10_GeV;
- HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass());
+ HEPMomentumType P0 = sqrt(E0 * E0 - Proton::mass() * Proton::mass());
auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
corsika::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(