From 456c4df36f048d75a3f6c9cf8a482b72e6fb4289 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 29 Nov 2018 09:30:29 +0000
Subject: [PATCH] added SibStack
---
Documentation/Examples/cascade_example.cc | 33 ++++++++++----
Framework/Cascade/SibStack.h | 53 +++++++++++++++++++++++
2 files changed, 78 insertions(+), 8 deletions(-)
create mode 100644 Framework/Cascade/SibStack.h
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index c4d7a7e8e..8b11f9ddf 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -19,6 +19,7 @@
#include <corsika/random/RNGManager.h>
#include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/cascade/SibStack.h>
//#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
@@ -103,6 +104,7 @@ public:
*/
// FOR NOW: set target to proton
int kTarget = 1; //p.GetPID();
+
std::cout << "ProcessSplit: " << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
if (E < 8.5_GeV || Ecm < 10_GeV ) {
std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
@@ -111,32 +113,45 @@ public:
} else {
// Sibyll does not know about units..
double sqs = Ecm / 1_GeV;
- // running sibyll
+ // running sibyll, filling stack
sibyll_( kBeam, kTarget, sqs);
// running decays
//decsib_();
// print final state
int print_unit = 6;
sib_list_( print_unit );
-
+
// delete current particle
p.Delete();
+ // add particles from sibyll to stack
+ // link to sibyll stack
+ SibStack ss;
+ /*
+ get transformation between Stack-frame and SibStack-frame
+ for EAS Stack-frame is lab. frame, could be different for CRMC-mode
+ the transformation should be derived from the input momenta
+ in general transformation is rotation + boost
+ */
const EnergyType proton_mass = 0.93827_GeV;
const double gamma = ( E + proton_mass ) / ( Ecm );
const double gambet = sqrt( E * E - proton_mass * proton_mass ) / Ecm;
-
- // add particles from sibyll to stack
- for(int i=0; i<s_plist_.np; ++i){
+
+ // SibStack does not know about momentum yet so we need counter to access momentum array in Sibyll
+ int i = -1;
+ for (auto &p: ss){
+ ++i;
//transform to lab. frame, primitve
- const double en_lab = gambet * s_plist_.p[2][i] + gamma * s_plist_.p[3][i];
+ const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();
// add to corsika stack
s.NewParticle().SetEnergy( en_lab * 1_GeV );
}
+
}
}
-
- void Init() //{ fCount = 0; }
+
+
+ void Init()
{
fCount = 0;
@@ -163,6 +178,7 @@ public:
sibyll_ini_();
// set particles stable / unstable
+
}
int GetCount() { return fCount; }
@@ -177,6 +193,7 @@ int main(){
const auto sequence = p0 + p1;
setup::Stack stack;
+
corsika::cascade::Cascade EAS(sequence, stack);
stack.Clear();
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
new file mode 100644
index 000000000..053bcdeb1
--- /dev/null
+++ b/Framework/Cascade/SibStack.h
@@ -0,0 +1,53 @@
+#ifndef _include_sibstack_h_
+#define _include_sibstack_h_
+
+#include <vector>
+#include <string>
+
+#include <corsika/stack/Stack.h>
+#include <corsika/cascade/sibyll2.3c.h>
+
+using namespace std;
+using namespace corsika::stack;
+
+class SibStackData {
+
+ public:
+ void Init();
+
+ void Clear() { s_plist_.np = 0; }
+
+ int GetSize() const { return s_plist_.np-1; }
+ int GetCapacity() const { return s_plist_.np-1; }
+
+
+ void SetId(const int i, const int v) { s_plist_.llist[i]=v; }
+ void SetEnergy(const int i, const double v) { s_plist_.p[3][i]=v; }
+
+ int GetId(const int i) const { return s_plist_.llist[i]; }
+ double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
+
+ void Copy(const int i1, const int i2) {
+ s_plist_.llist[i1] = s_plist_.llist[i2];
+ s_plist_.p[3][i1] = s_plist_.p[3][i2];
+ }
+
+ protected:
+ void IncrementSize() { s_plist_.np++; }
+ void DecrementSize() { if ( s_plist_.np>0) { s_plist_.np--; } }
+};
+
+
+template<typename StackIteratorInterface>
+class ParticleInterface : public ParticleBase<StackIteratorInterface> {
+ using ParticleBase<StackIteratorInterface>::GetStackData;
+ using ParticleBase<StackIteratorInterface>::GetIndex;
+ public:
+ void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
+ double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+};
+
+
+typedef Stack<SibStackData, ParticleInterface> SibStack;
+
+#endif
--
GitLab