diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index c4f8274233ea273caf75e2410a74148534fea267..46809532afa29069d4c0d97f545ed59dddb7bbe1 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -15,350 +15,378 @@
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/stack/SuperStupidStack.hpp>
-#include <corsika/logging/Logging.h>
-
#include <algorithm>
#include <tuple>
#include <vector>
namespace corsika {
- /**
- * @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 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 corsika
diff --git a/corsika/stack/SuperStupidStack.hpp b/corsika/stack/SuperStupidStack.hpp
index d83dda4f3d8cbf3e2026d5d603c144e1d4eebe95..f35b41da77949341b9ee3bf9117f9fe45e819f9c 100644
--- a/corsika/stack/SuperStupidStack.hpp
+++ b/corsika/stack/SuperStupidStack.hpp
@@ -22,193 +22,259 @@
namespace corsika {
- typedef corsika::Vector<hepmomentum_d> MomentumVector;
-
- namespace super_stupid {
-
- /**
- * Example of a particle object on the stack.
- */
-
- template <typename StackIteratorInterface>
- class ParticleInterface : public ParticleBase<StackIteratorInterface> {
-
- protected:
- using corsika::ParticleBase<StackIteratorInterface>::GetStack;
- using corsika::ParticleBase<StackIteratorInterface>::GetStackData;
-
- public:
- using corsika::ParticleBase<StackIteratorInterface>::GetIndex;
-
- public:
- void SetParticleData(const std::tuple<corsika::Code, HEPEnergyType, MomentumVector,
- corsika::Point, TimeType>& v) {
- SetPID(std::get<0>(v));
- SetEnergy(std::get<1>(v));
- SetMomentum(std::get<2>(v));
- SetPosition(std::get<3>(v));
- SetTime(std::get<4>(v));
- }
- /*
- void SetParticleData(const corsika::Code vDataPID,
- const HEPEnergyType vDataE,
- const MomentumVector& vMomentum,
- const corsika::Point& vPosition,
- const TimeType vTime) {
- }*/
-
- void SetParticleData(ParticleInterface<StackIteratorInterface>&,
- const std::tuple<corsika::Code, HEPEnergyType, MomentumVector,
- corsika::Point, TimeType>& v) {
- SetPID(std::get<0>(v));
- SetEnergy(std::get<1>(v));
- SetMomentum(std::get<2>(v));
- SetPosition(std::get<3>(v));
- SetTime(std::get<4>(v));
- }
- /* void SetParticleData(ParticleInterface<StackIteratorInterface>&,
- const corsika::Code vDataPID,
- const HEPEnergyType vDataE,
- const MomentumVector& vMomentum,
- const corsika::Point& vPosition,
- const TimeType vTime) {
- SetPID(vDataPID);
- SetEnergy(vDataE);
- SetMomentum(vMomentum);
- SetPosition(vPosition);
- SetTime(vTime);
- }*/
-
- /// individual setters
- void SetPID(const corsika::Code id) { GetStackData().SetPID(GetIndex(), id); }
- void SetEnergy(const HEPEnergyType& e) { GetStackData().SetEnergy(GetIndex(), e); }
- void SetMomentum(const MomentumVector& v) {
- GetStackData().SetMomentum(GetIndex(), v);
- }
- void SetPosition(const corsika::Point& v) {
- GetStackData().SetPosition(GetIndex(), v);
- }
- void SetTime(const TimeType& v) { GetStackData().SetTime(GetIndex(), v); }
-
- /// individual getters
- corsika::Code GetPID() const { return GetStackData().GetPID(GetIndex()); }
- HEPEnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
- MomentumVector GetMomentum() const {
- return GetStackData().GetMomentum(GetIndex());
- }
- corsika::Point GetPosition() const {
- return GetStackData().GetPosition(GetIndex());
- }
- TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); }
- /**
- * @name derived quantities
- *
- * @{
- */
- corsika::Vector<dimensionless_d> GetDirection() const {
- return GetMomentum() / GetEnergy();
- }
- HEPMassType GetMass() const { return corsika::mass(GetPID()); }
- int16_t GetChargeNumber() const { return corsika::charge_number(GetPID()); }
- ///@}
- };
-
- /**
- * Memory implementation of the most simple (stupid) particle stack object.
- */
-
- class SuperStupidStackImpl {
-
- public:
- void Init() {}
- void Dump() const {}
-
- void Clear() {
- fDataPID.clear();
- fDataE.clear();
- fMomentum.clear();
- fPosition.clear();
- fTime.clear();
- }
-
- unsigned int GetSize() const { return fDataPID.size(); }
- unsigned int GetCapacity() const { return fDataPID.size(); }
-
- void SetPID(const unsigned int i, const corsika::Code id) { fDataPID[i] = id; }
- void SetEnergy(const unsigned int i, const HEPEnergyType e) { fDataE[i] = e; }
- void SetMomentum(const unsigned int i, const MomentumVector& v) {
- fMomentum[i] = v;
- }
- void SetPosition(const unsigned int i, const corsika::Point& v) {
- fPosition[i] = v;
- }
- void SetTime(const unsigned int i, const TimeType& v) { fTime[i] = v; }
-
- corsika::Code GetPID(const unsigned int i) const { return fDataPID[i]; }
- HEPEnergyType GetEnergy(const unsigned int i) const { return fDataE[i]; }
- MomentumVector GetMomentum(const unsigned int i) const { return fMomentum[i]; }
- corsika::Point GetPosition(const unsigned int i) const { return fPosition[i]; }
- TimeType GetTime(const unsigned int i) const { return fTime[i]; }
-
- /**
- * Function to copy particle at location i2 in stack to i1
- */
- void Copy(const unsigned int i1, const unsigned int i2) {
- fDataPID[i2] = fDataPID[i1];
- fDataE[i2] = fDataE[i1];
- fMomentum[i2] = fMomentum[i1];
- fPosition[i2] = fPosition[i1];
- fTime[i2] = fTime[i1];
- }
-
- /**
- * Function to copy particle at location i2 in stack to i1
- */
- void Swap(const unsigned int i1, const unsigned int i2) {
- std::swap(fDataPID[i2], fDataPID[i1]);
- std::swap(fDataE[i2], fDataE[i1]);
- std::swap(fMomentum[i2], fMomentum[i1]);
- std::swap(fPosition[i2], fPosition[i1]);
- std::swap(fTime[i2], fTime[i1]);
- }
-
- void IncrementSize() {
- using corsika::Code;
- using corsika::Point;
- fDataPID.push_back(Code::Unknown);
- fDataE.push_back(0 * electronvolt);
- CoordinateSystem& dummyCS =
- RootCoordinateSystem::getInstance().GetRootCoordinateSystem();
- fMomentum.push_back(MomentumVector(
- dummyCS, {0 * electronvolt, 0 * electronvolt, 0 * electronvolt}));
- fPosition.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
- fTime.push_back(0 * second);
- }
-
- void DecrementSize() {
- if (fDataE.size() > 0) {
- fDataPID.pop_back();
- fDataE.pop_back();
- fMomentum.pop_back();
- fPosition.pop_back();
- fTime.pop_back();
- }
- }
-
- private:
- /// the actual memory to store particle data
-
- std::vector<corsika::Code> fDataPID;
- std::vector<HEPEnergyType> fDataE;
- std::vector<MomentumVector> fMomentum;
- std::vector<corsika::Point> fPosition;
- std::vector<TimeType> fTime;
-
- }; // end class SuperStupidStackImpl
-
- typedef Stack<SuperStupidStackImpl, ParticleInterface> SuperStupidStack;
-
- } // namespace super_stupid
+
+/**
+ * Example of a particle object on the stack.
+ */
+
+template <typename StackIteratorInterface>
+struct ParticleInterface : public ParticleBase<StackIteratorInterface> {
+
+private:
+
+ typedef corsika::ParticleBase<StackIteratorInterface> super_type;
+
+public:
+
+ typedef corsika::Vector<corsika::units::si::hepmomentum_d> momentum_vector_type;
+
+ std::string as_string() const {
+ using namespace corsika::units::si;
+ return fmt::format("particle: i={}, PID={}, E={}GeV", super_type::GetIndex(),
+ particles::GetName(this->getPID()), this->getEnergy() / 1_GeV);
+ }
+
+ void setParticleData( std::tuple<corsika::Code, corsika::units::si::HEPEnergyType,
+ momentum_vector_type, corsika::Point, corsika::units::si::TimeType> const& v) {
+ this->setPID(std::get<0>(v));
+ this->setEnergy(std::get<1>(v));
+ this->setMomentum(std::get<2>(v));
+ this->setPosition(std::get<3>(v));
+ this->setTime(std::get<4>(v));
+ }
+
+ void setParticleData( ParticleInterface<StackIteratorInterface> const&,
+ std::tuple<corsika::Code, corsika::units::si::HEPEnergyType,
+ momentum_vector_type, corsika::Point, corsika::units::si::TimeType> const& v) {
+ this->setPID(std::get<0>(v));
+ this->setEnergy(std::get<1>(v));
+ this->setMomentum(std::get<2>(v));
+ this->setPosition(std::get<3>(v));
+ this->setTime(std::get<4>(v));
+ }
+
+ /// individual setters
+ void setPID(const corsika::Code id) {
+ super_type::GetStackData().setPID(super_type::GetIndex(), id);
+ }
+ void setEnergy(const corsika::units::si::HEPEnergyType& e) {
+ super_type::GetStackData().setEnergy(super_type::GetIndex(), e);
+ }
+ void setMomentum(const momentum_vector_type& v) {
+ super_type::GetStackData().setMomentum(super_type::GetIndex(), v);
+ }
+ void setPosition(const corsika::Point& v) {
+ super_type::GetStackData().setPosition(super_type::GetIndex(), v);
+ }
+ void setTime(const corsika::units::si::TimeType& v) {
+ super_type::GetStackData().setTime(super_type::GetIndex(), v);
+ }
+
+ /// individual getters
+ corsika::Code getPID() const {
+ return super_type::GetStackData().getPID(super_type::GetIndex());
+ }
+ corsika::units::si::HEPEnergyType getEnergy() const {
+ return super_type::GetStackData().getEnergy(super_type::GetIndex());
+ }
+ momentum_vector_type getMomentum() const {
+ return super_type::GetStackData().getMomentum(super_type::GetIndex());
+ }
+ corsika::Point getPosition() const {
+ return super_type::GetStackData().getPosition(super_type::GetIndex());
+ }
+ corsika::units::si::TimeType getTime() const {
+ return super_type::GetStackData().getTime(super_type::GetIndex());
+ }
+ /**
+ * @name derived quantities
+ *
+ * @{
+ */
+ corsika::Vector<corsika::units::si::dimensionless_d> getDirection() const {
+ return this->getMomentum() / this->getEnergy();
+ }
+
+ corsika::units::si::HEPMassType getMass() const {
+ return corsika::GetMass(this->getPID());
+ }
+
+ int16_t getChargeNumber() const {
+ return corsika::GetChargeNumber(this->getPID());
+ }
+ ///@}
+};
+
+/**
+ * Memory implementation of the most simple (stupid) particle stack object.
+ *
+ */
+
+class SuperStupidStackImpl {
+
+public:
+
+ typedef corsika::Vector<corsika::units::si::hepmomentum_d> momentum_type;
+ typedef std::vector<corsika::Code> code_vector_type;
+ typedef std::vector<corsika::units::si::HEPEnergyType> energy_vector_type;
+ typedef std::vector<corsika::Point> point_vector_type;
+ typedef std::vector<corsika::units::si::TimeType> time_vector_type;
+ typedef std::vector<momentum_type> momentum_vector_type;
+
+ SuperStupidStackImpl()=default;
+
+ SuperStupidStackImpl( SuperStupidStackImpl const& other)=default;
+
+ SuperStupidStackImpl( SuperStupidStackImpl && other)=default;
+
+
+ SuperStupidStackImpl& operator=( SuperStupidStackImpl const& other)=default;
+
+ SuperStupidStackImpl& operator=( SuperStupidStackImpl && other)=default;
+
+
+ void init() {}
+ void dump() const {}
+
+ void clear() {
+ dataPID_.clear();
+ dataE_.clear();
+ momentum_.clear();
+ position_.clear();
+ time_.clear();
+ }
+
+ unsigned int getSize() const { return dataPID_.size(); }
+ unsigned int getCapacity() const { return dataPID_.size(); }
+
+ void setPID(size_t i, const corsika::Code id) {
+ dataPID_[i] = id;
+ }
+ void setEnergy(size_t i, corsika::units::si::HEPEnergyType const& e) {
+ dataE_[i] = e;
+ }
+ void setMomentum(size_t i, momentum_type const& v) {
+ momentum_[i] = v;
+ }
+ void setPosition(size_t i, corsika::Point const& v) {
+ position_[i] = v;
+ }
+ void setTime(size_t i, corsika::units::si::TimeType const& v) {
+ time_[i] = v;
+ }
+
+ corsika::Code getPID(size_t i) const {
+ return dataPID_[i];
+ }
+
+ corsika::units::si::HEPEnergyType getEnergy(size_t i) const {
+ return dataE_[i];
+ }
+
+ momentum_type getMomentum(size_t i) const {
+ return momentum_[i];
+ }
+
+ corsika::Point getPosition(size_t i) const {
+ return position_[i];
+ }
+ corsika::units::si::TimeType getTime(size_t i) const {
+ return time_[i];
+ }
+
+ corsika::units::si::HEPEnergyType getDataE(size_t i) const {
+ return dataE_[i];
+ }
+
+ void setDataE(size_t i, corsika::units::si::HEPEnergyType const& dataE) {
+ dataE_[i] = dataE;
+ }
+
+ corsika::Code getDataPid(size_t i) const {
+ return dataPID_;
+ }
+
+ void setDataPid(size_t i, corsika::Code const& dataPid) {
+ dataPID_[i] = dataPid;
+ }
+ /**
+ * Function to copy particle at location i2 in stack to i1
+ */
+ void copy(size_t i1, size_t i2) {
+ dataPID_[i2] = dataPID_[i1];
+ dataE_[i2] = dataE_[i1];
+ momentum_[i2] = momentum_[i1];
+ position_[i2] = position_[i1];
+ time_[i2] = time_[i1];
+ }
+
+
+ /**
+ * FIXME: change to iterators.
+ * Function to copy particle at location i2 in stack to i1
+ */
+ void swap(size_t i1, size_t i2) {
+ std::swap(dataPID_[i2] , dataPID_[i1]);
+ std::swap(dataE_[i2] , dataE_[i1]);
+ std::swap(momentum_[i2], momentum_[i1]);
+ std::swap(position_[i2], position_[i1]);
+ std::swap(time_[i2] , time_[i1]);
+ }
+
+ void incrementSize() {
+ using corsika::Point;
+ using corsika::Code;
+
+ dataPID_.push_back(Code::Unknown);
+ dataE_.push_back(0 * corsika::units::si::electronvolt);
+
+ CoordinateSystem& dummyCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ momentum_.push_back(momentum_type( dummyCS,
+ {0 * corsika::units::si::electronvolt, 0 * corsika::units::si::electronvolt,
+ 0 * corsika::units::si::electronvolt}));
+
+ position_.push_back(
+ Point(dummyCS, {0 * corsika::units::si::meter, 0 * corsika::units::si::meter,
+ 0 * corsika::units::si::meter}));
+ time_.push_back(0 * corsika::units::si::second);
+ }
+
+ void decrementSize() {
+ if (dataE_.size() > 0) {
+ dataPID_.pop_back();
+ dataE_.pop_back();
+ momentum_.pop_back();
+ position_.pop_back();
+ time_.pop_back();
+ }
+ }
+
+
+private:
+
+ /// the actual memory to store particle data
+ code_vector_type dataPID_;
+ energy_vector_type dataE_;
+ std::vector<momentum_type> momentum_;
+ std::vector<corsika::Point> position_;
+ std::vector<corsika::units::si::TimeType> time_;
+
+}; // end class SuperStupidStackImpl
+
+
+typedef Stack<SuperStupidStackImpl, ParticleInterface> SuperStupidStack;
+
} // namespace corsika