added nuclearstackextension

Stack/NuclearStackExtension/CMakeLists.txt

+
+set (NuclearStackExtension_HEADERS NuclearStackExtension.h)
+set (NuclearStackExtension_NAMESPACE corsika/stack/nuclear_extension)
+
+add_library (NuclearStackExtension INTERFACE)
+
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (NuclearStackExtension ${NuclearStackExtension_NAMESPACE} ${NuclearStackExtension_HEADERS})
+
+target_link_libraries (
+  NuclearStackExtension
+  INTERFACE
+  CORSIKAstackinterface
+  CORSIKAunits
+  CORSIKAparticles
+  CORSIKAgeometry
+  )
+
+target_include_directories (
+  NuclearStackExtension
+  INTERFACE
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include>
+  )
+
+install (
+  FILES
+  ${NuclearStackExtension_HEADERS}
+  DESTINATION
+  include/${NuclearStackExtension_NAMESPACE}
+  )
+
+# ----------------
+# code unit testing
+add_executable (
+  testNuclearStackExtension
+  testNuclearStackExtension.cc
+  )
+
+target_link_libraries (
+  testNuclearStackExtension
+  SuperStupidStack
+  NuclearStackExtension
+  #  CORSIKAutls
+  ProcessStackInspector
+  CORSIKAgeometry
+  CORSIKAunits
+  CORSIKAthirdparty # for catch2
+  )
+CORSIKA_ADD_TEST(testNuclearStackExtension)

Stack/NuclearStackExtension/NuclearStackExtension.h

+
+/*
+ * (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 <vector>
+
+namespace corsika::stack {
+
+  namespace nuclear_extension {
+
+    using corsika::stack::super_stupid::MomentumVector;
+    
+    template <typename StackIteratorInterface>
+    class ParticleInterface
+        : public corsika::stack::super_stupid::ParticleInterface<StackIteratorInterface> {
+
+    protected:
+      //      using corsika::stack::ParticleBase<StackIteratorInterface>::GetStack;
+      //using corsika::stack::super_stupid::ParticleInterface<StackIteratorInterface>::GetStackData;
+      using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+      using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+
+    public:
+      void SetParticleData(const corsika::particles::Code vDataPID,
+                           const corsika::units::si::HEPEnergyType vDataE,
+                           const corsika::stack::super_stupid::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(corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData().GetNucleusNextRef()); // store this nucleus data ref
+	  SetNuclearA(vA);
+	  SetNuclearZ(vZ);
+	} else {
+	  SetNuclearRef(-1); // this is not a nucleus
+	}
+        corsika::stack::super_stupid::ParticleInterface<StackIteratorInterface>::SetParticleData(vDataPID,
+												 vDataE,
+												 vMomentum,
+												 vPosition,
+												 vTime);
+      }
+
+      void SetParticleData(ParticleInterface<StackIteratorInterface>& parent,
+                           const corsika::particles::Code vDataPID,
+                           const corsika::units::si::HEPEnergyType vDataE,
+                           const corsika::stack::super_stupid::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.
+     */
+
+    class NuclearStackExtensionImpl
+        : public corsika::stack::super_stupid::SuperStupidStackImpl {
+
+    public:
+      void Init() { corsika::stack::super_stupid::SuperStupidStackImpl::Init(); }
+
+      void Clear() {
+        corsika::stack::super_stupid::SuperStupidStackImpl::Clear();
+        fNuclearRef.clear();
+        fNuclearA.clear();
+        fNuclearZ.clear();
+      }
+
+      int GetSize() const { return fNuclearRef.size(); }
+      int GetCapacity() const { return fNuclearRef.size(); }
+
+      void SetNuclearA(const int i, const int vA) { fNuclearA[GetNucleusRef(i)] = vA; }
+      void SetNuclearZ(const int i, const int vZ) { fNuclearZ[GetNucleusRef(i)] = vZ; }
+      void SetNuclearRef(const int i, const int v) { fNuclearRef[i] = v; }
+      
+      int GetNuclearA(const int i) const { return fNuclearA[GetNucleusRef(i)]; }
+      int GetNuclearZ(const 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 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 int i1, const int i2) {
+        corsika::stack::super_stupid::SuperStupidStackImpl::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.begin() + ref2);
+	    fNuclearZ.erase(fNuclearZ.begin() + 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 int i1, const int i2) {
+        corsika::stack::super_stupid::SuperStupidStackImpl::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() {
+        corsika::stack::super_stupid::SuperStupidStackImpl::IncrementSize();
+        fNuclearRef.push_back(-1);
+      }
+
+      void DecrementSize() {
+        corsika::stack::super_stupid::SuperStupidStackImpl::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 SuperStupidStackImpl
+
+    typedef Stack<NuclearStackExtensionImpl, ParticleInterface> NuclearStackExtension;
+
+  } // namespace nuclear_extension
+} // namespace corsika::stack
+
+#endif

Stack/NuclearStackExtension/testNuclearStackExtension.cc

+
+/*
+ * (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.
+ */
+
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
+#include <corsika/units/PhysicalUnits.h>
+
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+                          // cpp file
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::stack::nuclear_extension;
+
+#include <iostream>
+using namespace std;
+
+TEST_CASE("NuclearStackExtension", "[stack]") {
+
+  geometry::CoordinateSystem& dummyCS =
+      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+  SECTION("write non nucleus") {
+    NuclearStackExtension s;
+    s.AddParticle(particles::Code::Electron, 1.5_GeV,
+                  MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+                  Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+    REQUIRE(s.GetSize() == 1);
+  }
+  
+  SECTION("write nucleus") {
+    NuclearStackExtension s;
+    s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
+                  MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+                  Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 10);
+    REQUIRE(s.GetSize() == 1);
+  }
+
+  SECTION("write invalid nucleus") {
+    NuclearStackExtension s;
+    REQUIRE_THROWS(s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
+				 MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+				 Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 0, 0));
+  }
+
+  SECTION("read non nucleus") {
+    NuclearStackExtension s;
+    s.AddParticle(particles::Code::Electron, 1.5_GeV,
+                  MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+                  Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+    const auto pout = s.GetNextParticle();
+    REQUIRE(pout.GetPID() == particles::Code::Electron);
+    REQUIRE(pout.GetEnergy() == 1.5_GeV);
+    REQUIRE(pout.GetTime() == 100_s);
+  }
+  
+  SECTION("read nucleus") {
+    NuclearStackExtension s;
+    s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
+                  MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+                  Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9);
+    const auto pout = s.GetNextParticle();
+    REQUIRE(pout.GetPID() == particles::Code::Nucleus);
+    REQUIRE(pout.GetEnergy() == 1.5_GeV);
+    REQUIRE(pout.GetTime() == 100_s);
+    REQUIRE(pout.GetNuclearA() == 10);
+    REQUIRE(pout.GetNuclearZ() == 9);
+  }
+
+  SECTION("read invalid nucleus") {
+    NuclearStackExtension s;
+    s.AddParticle(particles::Code::Electron, 1.5_GeV,
+                  MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+                  Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+    const auto pout = s.GetNextParticle();
+    REQUIRE_THROWS(pout.GetNuclearA());
+    REQUIRE_THROWS(pout.GetNuclearZ());
+  }
+  
+
+  SECTION("stack fill and cleanup") {
+
+    NuclearStackExtension s;
+    // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
+    for (int i = 0; i < 99; ++i) {
+      if ((i+1)%10 == 0) {
+	s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
+		      MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+		      Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i/2);
+      } else {
+	s.AddParticle(particles::Code::Electron, 1.5_GeV,
+		      MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+		      Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+      }
+    }
+
+    REQUIRE(s.GetSize() == 99);
+    
+    for (int i = 0; i < 99; ++i) s.GetNextParticle().Delete();
+
+    REQUIRE(s.GetSize() == 0);
+  }
+
+
+  SECTION("stack operations") {
+
+    NuclearStackExtension s;
+    // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
+    for (int i = 0; i < 99; ++i) {
+      if ((i+1)%10 == 0) {
+	s.AddParticle(particles::Code::Nucleus, i*15_GeV,
+		      MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+		      Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i/2);
+      } else {
+	s.AddParticle(particles::Code::Electron, i*1.5_GeV,
+		      MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+		      Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+      }
+    }
+
+    // copy
+    {
+      s.Copy(s.begin()+9, s.begin()+10);
+      const auto& p9 = s.cbegin() + 9;
+      const auto& p10 = s.cbegin() + 10;
+      
+      REQUIRE(p9.GetPID() == particles::Code::Nucleus);
+      REQUIRE(p9.GetEnergy() == 9*15_GeV);
+      REQUIRE(p9.GetTime() == 100_s);
+      REQUIRE(p9.GetNuclearA() == 9);
+      REQUIRE(p9.GetNuclearZ() == 9/2);
+      
+      REQUIRE(p10.GetPID() == particles::Code::Nucleus);
+      REQUIRE(p10.GetEnergy() == 9*15_GeV);
+      REQUIRE(p10.GetTime() == 100_s);
+      REQUIRE(p10.GetNuclearA() == 9);
+      REQUIRE(p10.GetNuclearZ() == 9/2);
+    }
+    
+    // swap
+    {
+      s.Swap(s.begin()+11, s.begin()+10);
+      const auto& p11 = s.cbegin() + 11; // now: nucleus
+      const auto& p10 = s.cbegin() + 10; // now: electron
+      
+      REQUIRE(p11.GetPID() == particles::Code::Nucleus);
+      REQUIRE(p11.GetEnergy() == 9*15_GeV);
+      REQUIRE(p11.GetTime() == 100_s);
+      REQUIRE(p11.GetNuclearA() == 9);
+      REQUIRE(p11.GetNuclearZ() == 9/2);
+      
+      REQUIRE(p10.GetPID() == particles::Code::Electron);
+      REQUIRE(p10.GetEnergy() == 11*1.5_GeV);
+      REQUIRE(p10.GetTime() == 100_s);
+    }
+
+    // swap two nuclei
+    {
+      s.Copy(s.begin()+29, s.begin()+59);
+      const auto& p29 = s.cbegin() + 29;
+      const auto& p59 = s.cbegin() + 59;
+      
+      REQUIRE(p29.GetPID() == particles::Code::Nucleus);
+      REQUIRE(p29.GetEnergy() == 59*15_GeV);
+      REQUIRE(p29.GetTime() == 100_s);
+      REQUIRE(p29.GetNuclearA() == 59);
+      REQUIRE(p29.GetNuclearZ() == 59/2);
+
+      REQUIRE(p59.GetPID() == particles::Code::Nucleus);
+      REQUIRE(p59.GetEnergy() == 29*15_GeV);
+      REQUIRE(p59.GetTime() == 100_s);
+      REQUIRE(p59.GetNuclearA() == 29);
+      REQUIRE(p59.GetNuclearZ() == 29/2);
+    }
+    
+    for (int i = 0; i < 99; ++i) s.DeleteLast();
+    REQUIRE(s.GetSize() == 0);
+
+  }
+
+}