/* * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * * See file AUTHORS for a list of contributors. * * This software is distributed under the terms of the GNU General Public * Licence version 3 (GPL Version 3). See file LICENSE for a full version of * the license. */ #ifndef _include_stack_nuclearstackextension_h_ #define _include_stack_nuclearstackextension_h_ #include <corsika/stack/super_stupid/SuperStupidStack.h> #include <corsika/particles/ParticleProperties.h> #include <corsika/stack/Stack.h> #include <corsika/units/PhysicalUnits.h> #include <corsika/geometry/Point.h> #include <corsika/geometry/Vector.h> #include <algorithm> #include <vector> namespace corsika::stack { namespace nuclear_extension { typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector; /** * Define ParticleInterface for NuclearStackExtension Stack derived from * ParticleInterface of Inner stack class */ template <template <typename> typename InnerParticleInterface, typename StackIteratorInterface> class NuclearParticleInterface : public InnerParticleInterface<StackIteratorInterface> { public: using ExtendedParticleInterface = NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>; protected: using InnerParticleInterface<StackIteratorInterface>::GetStackData; using InnerParticleInterface<StackIteratorInterface>::GetIndex; public: void SetParticleData(const corsika::particles::Code vDataPID, const corsika::units::si::HEPEnergyType vDataE, const MomentumVector& vMomentum, const corsika::geometry::Point& vPosition, const corsika::units::si::TimeType vTime, const int vA = 0, const int vZ = 0) { if (vDataPID == corsika::particles::Code::Nucleus) { if (vA == 0 || vZ == 0) { std::ostringstream err; err << "NuclearStackExtension: no A and Z specified for new Nucleus!"; throw std::runtime_error(err.str()); } SetNuclearRef( GetStackData().GetNucleusNextRef()); // store this nucleus data ref SetNuclearA(vA); SetNuclearZ(vZ); } else { SetNuclearRef(-1); // this is not a nucleus } InnerParticleInterface<StackIteratorInterface>:: SetParticleData(vDataPID, vDataE, vMomentum, vPosition, vTime); } void SetParticleData(InnerParticleInterface<StackIteratorInterface>&, const corsika::particles::Code vDataPID, const corsika::units::si::HEPEnergyType vDataE, const MomentumVector& vMomentum, const corsika::geometry::Point& vPosition, const corsika::units::si::TimeType vTime, const int vA = 0, const int vZ = 0) { SetParticleData(vDataPID, vDataE, vMomentum, vPosition, vTime, vA, vZ); } /** * @name individual setters * @{ */ void SetNuclearA(const int vA) { GetStackData().SetNuclearA(GetIndex(), vA); } void SetNuclearZ(const int vZ) { GetStackData().SetNuclearZ(GetIndex(), vZ); } /// @} /** * @name individual getters * @{ */ int GetNuclearA() const { return GetStackData().GetNuclearA(GetIndex()); } int GetNuclearZ() const { return GetStackData().GetNuclearZ(GetIndex()); } /// @} protected: void SetNuclearRef(const int vR) { GetStackData().SetNuclearRef(GetIndex(), vR); } }; /** * Memory implementation of the most simple (stupid) particle stack object. */ template <typename InnerStackImpl> class NuclearStackExtensionImpl : public InnerStackImpl { public: void Init() { InnerStackImpl::Init(); } void Clear() { InnerStackImpl::Clear(); fNuclearRef.clear(); fNuclearA.clear(); fNuclearZ.clear(); } unsigned int GetSize() const { return fNuclearRef.size(); } unsigned int GetCapacity() const { return fNuclearRef.size(); } void SetNuclearA(const unsigned int i, const int vA) { fNuclearA[GetNucleusRef(i)] = vA; } void SetNuclearZ(const unsigned int i, const int vZ) { fNuclearZ[GetNucleusRef(i)] = vZ; } void SetNuclearRef(const unsigned int i, const int v) { fNuclearRef[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 (fNuclearRef[i] >= 0) return fNuclearRef[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) { if (i1 >= GetSize() || i2 >= GetSize()) { std::ostringstream err; err << "NuclearStackExtension: trying to access data beyond size of stack!"; throw std::runtime_error(err.str()); } InnerStackImpl::Copy(i1, i2); const int ref1 = fNuclearRef[i1]; const int ref2 = fNuclearRef[i2]; if (ref2 < 0) { if (ref1 >= 0) { // i1 is nucleus, i2 is not fNuclearRef[i2] = GetNucleusNextRef(); fNuclearA[fNuclearRef[i2]] = fNuclearA[ref1]; fNuclearZ[fNuclearRef[i2]] = fNuclearZ[ref1]; } else { // neither i1 nor i2 are nuclei } } else { if (ref1 >= 0) { // both are nuclei, i2 is overwritten with nucleus i1 fNuclearA[ref2] = fNuclearA[ref1]; fNuclearZ[ref2] = fNuclearZ[ref1]; } else { // i2 is overwritten with non-nucleus i1 fNuclearA.erase(fNuclearA.cbegin() + ref2); fNuclearZ.erase(fNuclearZ.cbegin() + ref2); const int n = fNuclearRef.size(); for (int i = 0; i < n; ++i) { if (fNuclearRef[i] >= ref2) { fNuclearRef[i] -= 1; } } } } } /** * Function to copy particle at location i2 in stack to i1 */ void Swap(const unsigned int i1, const unsigned int i2) { if (i1 >= GetSize() || i2 >= GetSize()) { std::ostringstream err; err << "NuclearStackExtension: trying to access data beyond size of stack!"; throw std::runtime_error(err.str()); } InnerStackImpl::Swap(i1, i2); const int ref1 = fNuclearRef[i1]; const int ref2 = fNuclearRef[i2]; if (ref2 < 0) { if (ref1 >= 0) { // i1 is nucleus, i2 is not std::swap(fNuclearRef[i2], fNuclearRef[i1]); } else { // neither i1 nor i2 are nuclei } } else { if (ref1 >= 0) { // both are nuclei, i2 is overwritten with nucleus i1 std::swap(fNuclearRef[i2], fNuclearRef[i1]); } else { // i2 is overwritten with non-nucleus i1 std::swap(fNuclearRef[i2], fNuclearRef[i1]); } } } protected: void IncrementSize() { InnerStackImpl::IncrementSize(); fNuclearRef.push_back(-1); } void DecrementSize() { InnerStackImpl::DecrementSize(); if (fNuclearRef.size() > 0) { const int ref = fNuclearRef.back(); fNuclearRef.pop_back(); if (ref >= 0) { fNuclearA.erase(fNuclearA.begin() + ref); fNuclearZ.erase(fNuclearZ.begin() + ref); const int n = fNuclearRef.size(); for (int i = 0; i < n; ++i) { if (fNuclearRef[i] >= ref) { fNuclearRef[i] -= 1; } } } } } private: /// the actual memory to store particle data std::vector<int> fNuclearRef; std::vector<int> fNuclearA; std::vector<int> 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::stack::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::stack::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::stack::super_stupid::SuperStupidStack, SS> >; //template <typename InnerStack> // using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, TEST<typename corsika::stack::super_stupid::SuperStupidStack::PIType>>; //using NuclearStackExtension = Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, TEST>; */ } // namespace nuclear_extension } // namespace corsika::stack #endif