Merge branch '69-extend-current-superstupidstack-to-work-for-first-paticle-cascade' into 'master'

Resolve "Extend current SuperStupidStack to work for first paticle cascade"

Closes #69

See merge request AirShowerPhysics/corsika!19

Stack/SuperStupidStack/CMakeLists.txt

@@ -12,6 +12,7 @@ target_link_libraries (
   CORSIKAstackinterface
   CORSIKAunits
   CORSIKAparticles
+  CORSIKAgeometry
   )
 
 target_include_directories (
@@ -27,3 +28,24 @@ install (
   DESTINATION
   include/${SuperStupidStack_NAMESPACE}
   )
+
+# ----------------
+# code unit testing
+add_executable (
+  testSuperStupidStack
+  testSuperStupidStack.cc
+  )
+
+target_link_libraries (
+  testSuperStupidStack
+#  CORSIKAutls
+  ProcessStackInspector
+  CORSIKAgeometry
+  CORSIKAunits
+  CORSIKAthirdparty # for catch2
+  )
+
+add_test (
+  NAME testSuperStupidStack
+  COMMAND testSuperStupidStack
+  )

Stack/SuperStupidStack/SuperStupidStack.h

@@ -16,7 +16,12 @@
 #include <corsika/stack/Stack.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/geometry/CoordinateSystem.h> // remove 
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+
 #include <vector>
+#include <algorithm>
 
 namespace corsika::stack {
 
@@ -24,24 +29,39 @@ namespace corsika::stack {
 
     using corsika::particles::Code;
     using corsika::units::si::EnergyType;
-    using corsika::units::si::operator""_GeV; // literals;
+    using corsika::units::si::TimeType;
+    using corsika::units::si::second;
+    using corsika::units::si::meter;
+    using corsika::units::si::joule;
+    using corsika::units::si::energy_d;
+    using corsika::geometry::Point;
+    using corsika::geometry::Vector;
+
+#warning replace this with a proper momentum vector:
+    typedef Vector<energy_d> MomentumVector; // should be momentum_d !!!
 
     /**
      * Example of a particle object on the stack.
      */
 
     template <typename StackIteratorInterface>
-    class ParticleInterface : public ParticleBase<StackIteratorInterface> {
-
+    class ParticleInterface : public ParticleBase<StackIteratorInterface> {    
+      
       using ParticleBase<StackIteratorInterface>::GetStackData;
       using ParticleBase<StackIteratorInterface>::GetIndex;
 
     public:
       void SetPID(const Code id) { GetStackData().SetPID(GetIndex(), id); }
       void SetEnergy(const EnergyType& e) { GetStackData().SetEnergy(GetIndex(), e); }
+      void SetMomentum(const MomentumVector& v) { GetStackData().SetMomentum(GetIndex(), v); }
+      void SetPosition(const Point& v) { GetStackData().SetPosition(GetIndex(), v); }
+      void SetTime(const TimeType& v) { GetStackData().SetTime(GetIndex(), v); }
 
       Code GetPID() const { return GetStackData().GetPID(GetIndex()); }
       EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+      MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
+      Point GetPosition() const { return GetStackData().GetPosition(GetIndex()); }
+      TimeType GetTime() const { return GetStackData().GetTime(GetIndex()); }
     };
 
     /**
@@ -61,42 +81,60 @@ namespace corsika::stack {
 
       int GetSize() const { return fDataPID.size(); }
       int GetCapacity() const { return fDataPID.size(); }
-
+      
       void SetPID(const int i, const Code id) { fDataPID[i] = id; }
       void SetEnergy(const int i, const EnergyType e) { fDataE[i] = e; }
-
+      void SetMomentum(const int i, const MomentumVector& v) { fMomentum[i] = v; }
+      void SetPosition(const int i, const Point& v) { fPosition[i] = v; }
+      void SetTime(const int i, const TimeType& v) { fTime[i] = v; }
+      
       Code GetPID(const int i) const { return fDataPID[i]; }
       EnergyType GetEnergy(const int i) const { return fDataE[i]; }
+      MomentumVector GetMomentum(const int i) const { return fMomentum[i]; }
+      Point GetPosition(const int i) const { return fPosition[i]; }
+      TimeType GetTime(const int i) const { return fTime[i]; }
 
       /**
        *   Function to copy particle at location i2 in stack to i1
        */
       void Copy(const int i1, const int i2) {
-        fDataE[i2] = fDataE[i1];
         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 int i1, const int i2) {
-        EnergyType tE = fDataE[i2];
-        Code tC = fDataPID[i2];
-        fDataE[i2] = fDataE[i1];
-        fDataPID[i2] = fDataPID[i1];
-        fDataE[i1] = tE;
-        fDataPID[i1] = tC;
+	std::swap(fDataPID[i2], fDataPID[i1]);
+	std::swap(fDataE[i2], fDataE[i1]);
+	std::swap(fMomentum[i2], fMomentum[i1]); // should be Momentum !!!!
+	std::swap(fPosition[i2], fPosition[i1]);
+	std::swap(fTime[i2], fTime[i1]);
       }
 
     protected:
       void IncrementSize() {
-        fDataE.push_back(0_GeV);
         fDataPID.push_back(Code::Unknown);
+        fDataE.push_back(0 * joule);
+#warning this here makes no sense: see issue #48
+	auto const dummyCS = corsika::geometry::CoordinateSystem::CreateRootCS();
+	fMomentum.push_back(MomentumVector(dummyCS,
+					   {0 * joule, 0 * joule, 0 * joule}));	
+	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();
-          fDataPID.pop_back();
+	  fMomentum.pop_back();
+	  fPosition.pop_back();
+	  fTime.pop_back();
         }
       }
 
@@ -105,6 +143,9 @@ namespace corsika::stack {
 
       std::vector<Code> fDataPID;
       std::vector<EnergyType> fDataE;
+      std::vector<Vector<corsika::units::si::energy_d>> fMomentum; // should be Momentum !!!!
+      std::vector<Point> fPosition;
+      std::vector<TimeType> fTime;
 
     }; // end class SuperStupidStackImpl
 

Stack/SuperStupidStack/testSuperStupidStack.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/stack/super_stupid/SuperStupidStack.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::super_stupid;
+
+
+#include <iostream>
+using namespace std;
+
+
+TEST_CASE("SuperStupidStack", "[stack]") {
+
+  SECTION("read+write") {
+
+    SuperStupidStack s;
+    auto p = s.NewParticle();
+    p.SetPID(corsika::particles::Code::Electron);
+    p.SetEnergy(1.5_GeV);
+    auto const dummyCS = corsika::geometry::CoordinateSystem::CreateRootCS();
+    p.SetMomentum(MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 * joule}));	
+    p.SetPosition(Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}));
+    p.SetTime(100_s);
+    
+    // read
+    REQUIRE(s.GetSize() == 1);
+    auto pout = s.GetNextParticle();
+    REQUIRE(pout.GetPID() == corsika::particles::Code::Electron);
+    REQUIRE(pout.GetEnergy() == 1.5_GeV);
+#warning Fix the next two lines:
+    //REQUIRE(pout.GetMomentum() == MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 * joule}));
+    //REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}));
+    REQUIRE(pout.GetTime() == 100_s);
+  }
+  
+  SECTION("write+delete") {
+
+    SuperStupidStack s;
+    for (int i=0; i<99; ++i)
+      s.NewParticle();
+
+    REQUIRE(s.GetSize() == 99);
+
+    for (int i=0; i<99; ++i)
+      s.GetNextParticle().Delete();
+
+
+    REQUIRE(s.GetSize() == 0);
+  }
+}