From d10704f0e3716ddc6c569f12635345ed198ad88b Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 6 Dec 2018 10:59:14 +0000
Subject: [PATCH] added actual decay to ProcessDecay
---
Framework/Cascade/SibStack.h | 4 +++-
Processes/Sibyll/ProcessDecay.h | 34 ++++++++++++++++++++++++++++++++-
2 files changed, 36 insertions(+), 2 deletions(-)
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 57a3f9c16..2275cdf49 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -6,6 +6,7 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/units/PhysicalUnits.h>
#include <corsika/stack/Stack.h>
using namespace std;
@@ -64,11 +65,12 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
+ void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v ); }
EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::sibyll::SibyllCode GetPID() const { return static_cast<corsika::process::sibyll::SibyllCode> (GetStackData().GetId(GetIndex())); }
super_stupid::MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
+ void SetMomentum(const super_stupid::MomentumVector& v) { GetStackData().SetMomentum(GetIndex(), v); }
};
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index ecae303e7..ab11fc7e4 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -5,6 +5,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/cascade/SibStack.h>
//using namespace corsika::particles;
@@ -104,7 +105,38 @@ namespace corsika::process {
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
-
+ SibStack ss;
+ ss.Clear();
+ // copy particle to sibyll stack
+ auto pin = ss.NewParticle();
+ pin.SetPID( process::sibyll::ConvertToSibyllRaw( p.GetPID() ) );
+ pin.SetEnergy( p.GetEnergy() );
+ pin.SetMomentum( p.GetMomentum() );
+ // set all particles/hadrons unstable
+ setHadronsUnstable();
+ // call sibyll decay
+ std::cout << "calling Sibyll decay routine.." << std::endl;
+ decsib_();
+ // print output
+ int print_unit = 6;
+ sib_list_( print_unit );
+ // copy particles from sibyll stack to corsika
+ int i = -1;
+ for (auto &psib: ss){
+ ++i;
+ // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
+ if( abs(s_plist_.llist[ i ]) > 100 ) continue;
+ // add to corsika stack
+ cout << "decay product: " << process::sibyll::ConvertFromSibyll( psib.GetPID() ) << endl;
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy( psib.GetEnergy() );
+ pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
+ pnew.SetMomentum( psib.GetMomentum() );
+ }
+ // empty sibyll stack
+ ss.Clear();
+ // remove original particle from stack
+ p.Delete();
}
template <typename Particle, typename Stack>
--
GitLab