From ca4ed4f568320603f81b9a4f0c19a45fb6dd1d58 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 29 Nov 2018 15:54:44 +0000
Subject: [PATCH] fixes in particleConversion for sibyll, added decay table to
interface, added pid to sibstack
---
Documentation/Examples/cascade_example.cc | 50 +++++++++++++++--------
Framework/Cascade/SibStack.h | 8 +++-
Framework/Cascade/sibyll2.3c.h | 8 ++++
Processes/Sibyll/code_generator.py | 3 +-
Processes/Sibyll/sibyll_codes.dat | 3 +-
Processes/Sibyll/testSibyll.cc | 10 ++---
6 files changed, 57 insertions(+), 25 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 50a8f351c..2490cd344 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -42,19 +42,11 @@ public:
template <typename Particle>
double MinStepLength(Particle& p) const {
+
// beam particles for sibyll : 1, 2, 3 for p, pi, k
- corsika::particles::Code c_id = p.GetPID();
- corsika::process::sibyll::Code s_id = process::sibyll::ConvertToSibyll( p.GetPID() );
- corsika::particles::Code p_id_2 = process::sibyll::ConvertFromSibyll( corsika::process::sibyll::Code::Proton );
- cout << p_id_2 << endl;
- std::cout << "MinStepLength: particle input " << "corsika id: " << c_id << std::endl;
- auto test = static_cast<corsika::process::sibyll::SibyllCodeIntType>(s_id);
- std::cout << "MinStepLength: particle input " << "sibyll id: |" << (int)test << "|" <<std::endl;
- // std::cout << "MinStepLength: particle input " << "sibyll id: " << process::sibyll::ConvertToSibyll( p.GetPID() ) << std::endl;
- cout << p.GetPID() << " --> " << process::sibyll::ConvertToSibyllRaw( p.GetPID() ) << endl;
+ // read from cross section code table
+ int kBeam = 1;
- int kBeam = 1;
-
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
@@ -98,6 +90,9 @@ public:
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
+ cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() ) << endl;
+ if( process::sibyll::CanInteract( p.GetPID() ) ){
+
// get energy of particle from stack
/*
stack is in GeV in lab. frame
@@ -107,8 +102,9 @@ public:
*/
EnergyType E = p.GetEnergy();
EnergyType Ecm = sqrt( 2. * E * 0.93827_GeV );
- // FOR NOW: set beam to proton
- int kBeam = 13; //p.GetPID();
+
+ int kBeam = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
@@ -156,10 +152,15 @@ public:
//transform to lab. frame, primitve
const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();
// add to corsika stack
- s.NewParticle().SetEnergy( en_lab * 1_GeV );
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy( en_lab * 1_GeV );
+ pnew.SetPID( process::sibyll::ConvertFromSibyll( p.GetPID() ) );
+
}
}
+ }else
+ p.Delete();
}
@@ -183,14 +184,31 @@ public:
// test random number generator
std::cout << "ProcessSplit: " << " test sequence of random numbers." << std::endl;
int a = 0;
- for(int i=0; i<5; ++i)
+ for(int i=0; i<8; ++i)
std::cout << i << " " << s_rndm_(a) << std::endl;
//initialize Sibyll
sibyll_ini_();
// set particles stable / unstable
-
+ // use stack to loop over particles
+ setup::Stack ds;
+ ds.NewParticle().SetPID( Code::Proton );
+ ds.NewParticle().SetPID( Code::Neutron );
+ ds.NewParticle().SetPID( Code::PiPlus );
+ ds.NewParticle().SetPID( Code::PiMinus );
+ ds.NewParticle().SetPID( Code::KPlus );
+ ds.NewParticle().SetPID( Code::KMinus );
+ ds.NewParticle().SetPID( Code::K0Long );
+ ds.NewParticle().SetPID( Code::K0Short );
+
+ for( auto &p: ds){
+ int s_id = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+ // set particle stable by setting table value negative
+ cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")" << " stable in Sibyll .." << endl;
+ s_csydec_.idb[ s_id ] = -s_csydec_.idb[ s_id-1 ];
+ p.Delete();
+ }
}
int GetCount() { return fCount; }
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 053bcdeb1..d38bf0da9 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -6,6 +6,7 @@
#include <corsika/stack/Stack.h>
#include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
using namespace std;
using namespace corsika::stack;
@@ -17,8 +18,8 @@ class SibStackData {
void Clear() { s_plist_.np = 0; }
- int GetSize() const { return s_plist_.np-1; }
- int GetCapacity() const { return s_plist_.np-1; }
+ 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; }
@@ -45,6 +46,9 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
public:
void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
double 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())); }
+
};
diff --git a/Framework/Cascade/sibyll2.3c.h b/Framework/Cascade/sibyll2.3c.h
index 4a00010fc..0899d02ee 100644
--- a/Framework/Cascade/sibyll2.3c.h
+++ b/Framework/Cascade/sibyll2.3c.h
@@ -23,6 +23,14 @@ extern"C"{
}s_plist_;
+ extern struct {
+ double cbr[223+16+12+8];
+ int kdec[1338+6*(16+12+8)];
+ int lbarp[99];
+ int idb[99];
+ }s_csydec_;
+
+
// additional particle stack for the mother particles of unstable particles
// stable particles have entry zero
extern struct {
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
index 784b8e7d8..a87e644fa 100755
--- a/Processes/Sibyll/code_generator.py
+++ b/Processes/Sibyll/code_generator.py
@@ -28,7 +28,7 @@ def read_sibyll_codes(filename, pythia_db):
line = line.strip()
if line[0] == '#':
continue
- identifier, sib_code, canInteractFlag = line.split()
+ identifier, sib_code, canInteractFlag, xsctnId = line.split()
try:
pythia_db[identifier]["sibyll_code"] = int(sib_code)
pythia_db[identifier]["sibyll_canInteract"] = int(canInteractFlag)
@@ -148,3 +148,4 @@ if __name__ == "__main__":
print(generate_corsika2sibyll(pythia_db), file=f)
print(generate_known_particle(pythia_db), file=f)
print(generate_sibyll2corsika(pythia_db), file=f)
+ print(generate_interacting_particle(pythia_db), file=f)
diff --git a/Processes/Sibyll/sibyll_codes.dat b/Processes/Sibyll/sibyll_codes.dat
index 9b6479216..d04a8436b 100644
--- a/Processes/Sibyll/sibyll_codes.dat
+++ b/Processes/Sibyll/sibyll_codes.dat
@@ -12,7 +12,8 @@ NuTau 92 0 0
NuTauBar 93 0 0
Gamma 1 0 0
Pi0 6 1 2
-Rho0 27 1 2
+# rho0 could interact but sibyll has no cross section/interaction length. was used for gamma had int
+Rho0 27 0 0
K0Long 11 1 3
PiPlus 7 1 2
PiMinus 8 1 2
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 7e8049524..8302be009 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -23,20 +23,20 @@ TEST_CASE("Sibyll", "[processes]") {
SECTION("Sibyll -> Corsika") {
REQUIRE(corsika::particles::Electron::GetCode() ==
- process::sibyll::ConvertFromSibyll(process::sibyll::Code::Electron));
+ process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
}
SECTION("Corsika -> Sibyll") {
REQUIRE(process::sibyll::ConvertToSibyll(corsika::particles::Electron::GetCode()) ==
- process::sibyll::Code::Electron);
+ process::sibyll::SibyllCode::Electron);
REQUIRE(process::sibyll::ConvertToSibyllRaw(corsika::particles::Proton::GetCode()) ==
13 );
}
- SECTION("handledBySibyll") {
- REQUIRE(process::sibyll::HandledBySibyll(corsika::particles::Electron::GetCode()));
+ SECTION("KnownBySibyll") {
+ REQUIRE(process::sibyll::KnownBySibyll(corsika::particles::Electron::GetCode()));
REQUIRE_FALSE(
- process::sibyll::HandledBySibyll(corsika::particles::XiPrimeC0::GetCode()));
+ process::sibyll::KnownBySibyll(corsika::particles::XiPrimeC0::GetCode()));
}
}
--
GitLab