complete interface change

Documentation/Examples/cascade_example.cc

@@ -256,8 +256,6 @@ int main() {
   double theta = 0.;
   double phi = 0.;
   {
-    auto particle = stack.NewParticle();
-    particle.SetPID(Code::Proton);
     HEPMomentumType P0 = sqrt(E0 * E0 - Proton::GetMass() * Proton::GetMass());
     auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
       return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
@@ -268,11 +266,8 @@ int main() {
     auto plab = stack::super_stupid::MomentumVector(rootCS, {px, py, pz});
     cout << "input angles: theta=" << theta << " phi=" << phi << endl;
     cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
-    particle.SetEnergy(E0);
-    particle.SetMomentum(plab);
-    particle.SetTime(0_ns);
-    Point p(rootCS, 0_m, 0_m, 0_m);
-    particle.SetPosition(p);
+    Point pos(rootCS, 0_m, 0_m, 0_m);
+    stack.AddParticle(Code::Proton, E0, plab, pos, 0_ns);
   }
 
   // define air shower object, run simulation

Documentation/Examples/stack_example.cc

@@ -11,6 +11,10 @@
 
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/stack/super_stupid/SuperStupidStack.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+
 #include <cassert>
 #include <iomanip>
 #include <iostream>
@@ -20,10 +24,12 @@ using namespace corsika::stack;
 using namespace std;
 
 void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
+  const geometry::CoordinateSystem& rootCS =
+      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
   for (int i = 0; i < 11; ++i) {
-    auto p = s.NewParticle();
-    p.SetPID(corsika::particles::Code::Electron);
-    p.SetEnergy(1.5_GeV * i);
+    s.AddParticle(corsika::particles::Code::Electron, 1.5_GeV * i,
+                  stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
+                  geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns);
   }
 }
 

Framework/Cascade/testCascade.cc

@@ -84,7 +84,7 @@ public:
   }
 
   template <typename Particle, typename T, typename Stack>
-  EProcessReturn DoContinuous(Particle& p, T&, Stack& s) {
+  EProcessReturn DoContinuous(Particle& p, T&, Stack&) {
     fCalls++;
     HEPEnergyType E = p.GetEnergy();
     if (E < fEcrit) {
@@ -92,13 +92,7 @@ public:
       fCount++;
     } else {
       p.SetEnergy(E / 2);
-      auto pnew = s.NewParticle();
-      // s.Copy(p, pnew); fix that .... todo
-      pnew.SetPID(p.GetPID());
-      pnew.SetTime(p.GetTime());
-      pnew.SetEnergy(E / 2);
-      pnew.SetPosition(p.GetPosition());
-      pnew.SetMomentum(p.GetMomentum());
+      p.AddSecondary(p.GetPID(), E / 2, p.GetMomentum(), p.GetPosition(), p.GetTime());
     }
     return EProcessReturn::eOk;
   }
@@ -133,14 +127,11 @@ TEST_CASE("Cascade", "[Cascade]") {
       RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
   stack.Clear();
-  auto particle = stack.NewParticle();
   HEPEnergyType E0 = 100_GeV;
-  particle.SetPID(particles::Code::Electron);
-  particle.SetEnergy(E0);
-  particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
-  particle.SetMomentum(
-      corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}));
-  particle.SetTime(0_ns);
+  stack.AddParticle(
+      particles::Code::Electron, E0,
+      corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
+      Point(rootCS, {0_m, 0_m, 10_km}), 0_ns);
   EAS.Init();
   EAS.Run();
 

Framework/StackInterface/ParticleBase.h

@@ -46,6 +46,11 @@ namespace corsika::stack {
     /// will be invalidated by this operation
     void Delete() { GetIterator().GetStack().Delete(GetIterator()); }
 
+    template <typename... Args>
+    StackIterator AddSecondary(const Args... v) {
+      return GetStack().AddSecondary(GetIterator(), v...);
+    }
+
     //  protected: // todo should be proteced, but don't now how to 'friend Stack'
     /// Function to provide CRTP access to inheriting class (type)
     StackIterator& GetIterator() { return static_cast<StackIterator&>(*this); }

Framework/StackInterface/Stack.h

@@ -71,12 +71,14 @@ namespace corsika::stack {
 
     /// increase stack size, create new particle at end of stack
     template <typename... Args>
-    StackIterator AddParticle(Args... v) {
+    StackIterator AddParticle(const Args... v) {
       IncrementSize();
       return StackIterator(*this, GetSize() - 1, v...);
-      // auto p = StackIterator(*this, GetSize() - 1);
-      // p.SetParticleData(v...);
-      // return p;
+    }
+    template <typename... Args>
+    StackIterator AddSecondary(StackIterator& parent, const Args... v) {
+      IncrementSize();
+      return StackIterator(*this, GetSize() - 1, parent, v...);
     }
     void Copy(StackIterator& a, StackIterator& b) { Copy(a.GetIndex(), b.GetIndex()); }
     /// delete this particle

Framework/StackInterface/StackIterator.h

@@ -95,6 +95,14 @@ namespace corsika::stack {
       (**this).SetParticleData(args...);
     }
 
+    template <typename... Args>
+    StackIteratorInterface(StackType& data, const int index,
+                           StackIteratorInterface& parent, const Args... args)
+        : fIndex(index)
+        , fData(&data) {
+      (**this).SetParticleData(*parent, args...);
+    }
+
   public:
     StackIteratorInterface& operator++() {
       ++fIndex;
@@ -116,7 +124,7 @@ namespace corsika::stack {
     }
 
   protected:
-    int GetIndex() const { return fIndex; }
+    inline int GetIndex() const { return fIndex; }
     StackType& GetStack() { return *fData; }
     const StackType& GetStack() const { return *fData; }
     StackData& /*typename std::decay<StackData>::type&*/ GetStackData() {

Framework/StackInterface/testStackInterface.cc

@@ -61,10 +61,14 @@ class TestParticleInterface : public ParticleBase<StackIteratorInterface> {
   using ParticleBase<StackIteratorInterface>::GetStack;
   using ParticleBase<StackIteratorInterface>::GetStackData;
   using ParticleBase<StackIteratorInterface>::GetIndex;
+  using ParticleBase<StackIteratorInterface>::GetIterator;
 
 public:
   // one version
-  void AddSecondary(const double v) { GetStack().AddParticle(v); }
+  StackIteratorInterface& AddSecondary(const double v) {
+    GetStack().AddParticle(v);
+    return GetIterator();
+  }
   // another version
   void AddSecondary(const double v, const double p) { GetStack().AddParticle(v + p); }
 
@@ -82,7 +86,7 @@ TEST_CASE("Stack", "[Stack]") {
   // helper function for sum over stack data
   auto sum = [](const StackTest& stack) {
     double v = 0;
-    for (const auto&& p : stack) v += p.GetData();
+    for (const auto& p : stack) v += p.GetData();
     return v;
   };
 

Processes/NullModel/testNullModel.cc

@@ -21,6 +21,8 @@
 
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/stack/super_stupid/SuperStupidStack.h>
+
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 
@@ -30,16 +32,19 @@ using namespace corsika;
 
 TEST_CASE("NullModel", "[processes]") {
 
-  auto const& cs =
+  auto const& dummyCS =
       geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-  geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+  geometry::Point const origin(dummyCS, {0_m, 0_m, 0_m});
   geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(
-      cs, 0_m / second, 0_m / second, 1_m / second);
+      dummyCS, 0_m / second, 0_m / second, 1_m / second);
   geometry::Line line(origin, v);
   geometry::Trajectory<geometry::Line> track(line, 10_s);
 
   setup::Stack stack;
-  auto particle = stack.NewParticle();
+  auto particle = stack.AddParticle(
+      particles::Code::Electron, 1.5_GeV,
+      stack::super_stupid::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+      geometry::Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
 
   SECTION("interface") {
 

Processes/Sibyll/Decay.h

@@ -163,24 +163,22 @@ namespace corsika::process {
       }
 
       template <typename Particle, typename Stack>
-      void DoDecay(Particle& p, Stack& s) {
+      void DoDecay(Particle& p, Stack&) {
         using corsika::geometry::Point;
         using namespace corsika::units::si;
 
         fCount++;
         SibStack ss;
         ss.Clear();
-        // copy particle to sibyll stack
-        auto pin = ss.NewParticle();
         const corsika::particles::Code pCode = p.GetPID();
-        pin.SetPID(process::sibyll::ConvertToSibyllRaw(pCode));
-        pin.SetEnergy(p.GetEnergy());
-        pin.SetMomentum(p.GetMomentum());
-        // setting particle mass with Corsika values, may be inconsistent with sibyll
-        // internal values
-        // TODO: #warning setting particle mass with Corsika values, may be inconsistent
-        // with sibyll internal values
-        pin.SetMass(corsika::particles::GetMass(pCode));
+        // copy particle to sibyll stack
+        ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), p.GetEnergy(),
+                       p.GetMomentum(),
+                       // setting particle mass with Corsika values, may be inconsistent
+                       // with sibyll internal values
+                       // TODO: #warning setting particle mass with Corsika values, may be
+                       // inconsistent with sibyll internal values
+                       corsika::particles::GetMass(pCode));
         // remember position
         Point const decayPoint = p.GetPosition();
         TimeType const t0 = p.GetTime();
@@ -203,12 +201,8 @@ namespace corsika::process {
           // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
           if (psib.HasDecayed()) continue;
           // add to corsika stack
-          auto pnew = s.NewParticle();
-          pnew.SetEnergy(psib.GetEnergy());
-          pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
-          pnew.SetMomentum(psib.GetMomentum());
-          pnew.SetPosition(decayPoint);
-          pnew.SetTime(t0);
+          p.AddSecondary(process::sibyll::ConvertFromSibyll(psib.GetPID()),
+                         psib.GetEnergy(), psib.GetMomentum(), decayPoint, t0);
         }
         // empty sibyll stack
         ss.Clear();

Processes/Sibyll/Interaction.h

@@ -172,7 +172,7 @@ namespace corsika::process::sibyll {
      */
 
     template <typename Particle, typename Stack>
-    corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
+    corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&) {
 
       using namespace corsika::units;
       using namespace corsika::utl;
@@ -328,12 +328,9 @@ namespace corsika::process::sibyll {
             auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
 
             // add to corsika stack
-            auto pnew = s.NewParticle();
-            pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
-            pnew.SetEnergy(Plab.GetTimeLikeComponent());
-            pnew.SetMomentum(Plab.GetSpaceLikeComponents());
-            pnew.SetPosition(pOrig);
-            pnew.SetTime(tOrig);
+            auto pnew = p.AddSecondary(process::sibyll::ConvertFromSibyll(psib.GetPID()),
+                                       Plab.GetTimeLikeComponent(),
+                                       Plab.GetSpaceLikeComponents(), pOrig, tOrig);
 
             Plab_final += pnew.GetMomentum();
             Elab_final += pnew.GetEnergy();

Processes/Sibyll/SibStack.h

@@ -29,9 +29,7 @@ namespace corsika::process::sibyll {
     void Init();
 
     void Clear() { s_plist_.np = 0; }
-
     int GetSize() const { return s_plist_.np; }
-
     int GetCapacity() const { return 8000; }
 
     void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
@@ -43,7 +41,6 @@ namespace corsika::process::sibyll {
       using namespace corsika::units::si;
       s_plist_.p[4][i] = v / 1_GeV;
     }
-
     void SetMomentum(const int i, const MomentumVector& v) {
       using namespace corsika::units::si;
       auto tmp = v.GetComponents();
@@ -51,7 +48,6 @@ namespace corsika::process::sibyll {
     }
 
     int GetId(const int i) const { return s_plist_.llist[i]; }
-
     corsika::units::si::HEPEnergyType GetEnergy(const int i) const {
       using namespace corsika::units::si;
       return s_plist_.p[3][i] * 1_GeV;
@@ -60,7 +56,6 @@ namespace corsika::process::sibyll {
       using namespace corsika::units::si;
       return s_plist_.p[4][i] * 1_GeV;
     }
-
     MomentumVector GetMomentum(const int i) const {
       using corsika::geometry::CoordinateSystem;
       using corsika::geometry::QuantityVector;
@@ -93,6 +88,27 @@ namespace corsika::process::sibyll {
     using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
 
   public:
+    void SetParticleData(const int vID, // corsika::process::sibyll::SibyllCode vID,
+                         const corsika::units::si::HEPEnergyType vE,
+                         const MomentumVector& vP,
+                         const corsika::units::si::HEPMassType vM) {
+      SetPID(vID);
+      SetEnergy(vE);
+      SetMomentum(vP);
+      SetMass(vM);
+    }
+
+    void SetParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
+                         const int vID, //  corsika::process::sibyll::SibyllCode vID,
+                         const corsika::units::si::HEPEnergyType vE,
+                         const MomentumVector& vP,
+                         const corsika::units::si::HEPMassType vM) {
+      SetPID(vID);
+      SetEnergy(vE);
+      SetMomentum(vP);
+      SetMass(vM);
+    }
+
     void SetEnergy(const corsika::units::si::HEPEnergyType v) {
       GetStackData().SetEnergy(GetIndex(), v);
     }

Processes/Sibyll/testSibyll.cc

@@ -116,7 +116,12 @@ TEST_CASE("SibyllInterface", "[processes]") {
   SECTION("InteractionInterface") {
 
     setup::Stack stack;
-    auto particle = stack.NewParticle();
+    const HEPEnergyType E0 = 10_GeV;
+    HEPMomentumType P0 =
+        sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+    auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+    geometry::Point pos(cs, 0_m, 0_m, 0_m);
+    auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
 
     Interaction model(env);
 
@@ -130,19 +135,12 @@ TEST_CASE("SibyllInterface", "[processes]") {
   SECTION("DecayInterface") {
 
     setup::Stack stack;
-    auto particle = stack.NewParticle();
-    {
-      const HEPEnergyType E0 = 10_GeV;
-      particle.SetPID(particles::Code::Proton);
-      HEPMomentumType P0 =
-          sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
-      auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
-      particle.SetEnergy(E0);
-      particle.SetMomentum(plab);
-      particle.SetTime(0_ns);
-      geometry::Point p(cs, 0_m, 0_m, 0_m);
-      particle.SetPosition(p);
-    }
+    const HEPEnergyType E0 = 10_GeV;
+    HEPMomentumType P0 =
+        sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
+    auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+    geometry::Point pos(cs, 0_m, 0_m, 0_m);
+    auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
 
     Decay model;
 

Processes/StackInspector/testStackInspector.cc

@@ -30,16 +30,19 @@ using namespace corsika;
 
 TEST_CASE("StackInspector", "[processes]") {
 
-  auto const& cs =
+  auto const& rootCS =
       geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-  geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+  geometry::Point const origin(rootCS, {0_m, 0_m, 0_m});
   geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(
-      cs, 0_m / second, 0_m / second, 1_m / second);
+      rootCS, 0_m / second, 0_m / second, 1_m / second);
   geometry::Line line(origin, v);
   geometry::Trajectory<geometry::Line> track(line, 10_s);
 
   setup::Stack stack;
-  auto particle = stack.NewParticle();
+  auto particle = stack.AddParticle(
+      particles::Code::Electron, 10_GeV,
+      corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
+      geometry::Point(rootCS, {0_m, 0_m, 10_km}), 0_ns);
 
   SECTION("interface") {
 

Stack/SuperStupidStack/SuperStupidStack.h

@@ -43,23 +43,53 @@ namespace corsika::stack {
       using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
 
     public:
+      /// the factory function, this is how to create a new particle:
+      /*void AddSecondary(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) {
+        GetStack().AddParticle(vDataPID, vDataE, vMomentum, vPosition, vTime);
+        }*/
+      //
+      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) {
+        SetPID(vDataPID);
+        SetEnergy(vDataE);
+        SetMomentum(vMomentum);
+        SetPosition(vPosition);
+        SetTime(vTime);
+      }
+
+      void SetParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
+                           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) {
+        SetPID(vDataPID);
+        SetEnergy(vDataE);
+        SetMomentum(vMomentum);
+        SetPosition(vPosition);
+        SetTime(vTime);
+      }
+
       /// individual setters
       void SetPID(const corsika::particles::Code id) {
         GetStackData().SetPID(GetIndex(), id);
       }
-
       void SetEnergy(const corsika::units::si::HEPEnergyType& e) {
         GetStackData().SetEnergy(GetIndex(), e);
       }
-
       void SetMomentum(const MomentumVector& v) {
         GetStackData().SetMomentum(GetIndex(), v);
       }
-
       void SetPosition(const corsika::geometry::Point& v) {
         GetStackData().SetPosition(GetIndex(), v);
       }
-
       void SetTime(const corsika::units::si::TimeType& v) {
         GetStackData().SetTime(GetIndex(), v);
       }
@@ -68,23 +98,18 @@ namespace corsika::stack {
       corsika::particles::Code GetPID() const {
         return GetStackData().GetPID(GetIndex());
       }
-
       corsika::units::si::HEPEnergyType GetEnergy() const {
         return GetStackData().GetEnergy(GetIndex());
       }
-
       MomentumVector GetMomentum() const {
         return GetStackData().GetMomentum(GetIndex());
       }
-
       corsika::geometry::Point GetPosition() const {
         return GetStackData().GetPosition(GetIndex());
       }
-
       corsika::units::si::TimeType GetTime() const {
         return GetStackData().GetTime(GetIndex());
       }
-
       corsika::geometry::Vector<corsika::units::si::dimensionless_d> GetDirection()
           const {
         return GetMomentum() / GetEnergy();
@@ -110,30 +135,22 @@ namespace corsika::stack {
       }
 
       int GetSize() const { return fDataPID.size(); }
-
       int GetCapacity() const { return fDataPID.size(); }
 
       void SetPID(const int i, const corsika::particles::Code id) { fDataPID[i] = id; }
-
       void SetEnergy(const int i, const corsika::units::si::HEPEnergyType e) {
         fDataE[i] = e;
       }
       void SetMomentum(const int i, const MomentumVector& v) { fMomentum[i] = v; }
-
       void SetPosition(const int i, const corsika::geometry::Point& v) {
         fPosition[i] = v;
       }
-
       void SetTime(const int i, const corsika::units::si::TimeType& v) { fTime[i] = v; }
 
       corsika::particles::Code GetPID(const int i) const { return fDataPID[i]; }
-
       corsika::units::si::HEPEnergyType GetEnergy(const int i) const { return fDataE[i]; }
-
       MomentumVector GetMomentum(const int i) const { return fMomentum[i]; }
-
       corsika::geometry::Point GetPosition(const int i) const { return fPosition[i]; }
-
       corsika::units::si::TimeType GetTime(const int i) const { return fTime[i]; }
 
       /**

Stack/SuperStupidStack/testSuperStupidStack.cc

@@ -28,17 +28,15 @@ using namespace std;
 
 TEST_CASE("SuperStupidStack", "[stack]") {
 
+  geometry::CoordinateSystem& dummyCS =
+      geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
   SECTION("read+write") {
 
     SuperStupidStack s;
-    auto p = s.NewParticle();
-    p.SetPID(particles::Code::Electron);
-    p.SetEnergy(1.5_GeV);
-    geometry::CoordinateSystem& dummyCS =
-        geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-    p.SetMomentum(MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}));
-    p.SetPosition(Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}));
-    p.SetTime(100_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);
 
     // read
     REQUIRE(s.GetSize() == 1);
@@ -54,7 +52,10 @@ TEST_CASE("SuperStupidStack", "[stack]") {
   SECTION("write+delete") {
 
     SuperStupidStack s;
-    for (int i = 0; i < 99; ++i) s.NewParticle();
+    for (int i = 0; i < 99; ++i)
+      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);