From cbcccd22812b75c5f72623006ca90f074641f9e1 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 26 Dec 2018 22:47:03 +0000
Subject: [PATCH] repaired interface to sibyll decay routine to properly pass
particle masses, added positrons as EM particle to EMCut
---
Documentation/Examples/cascade_example.cc | 3 +++
Processes/Sibyll/Decay.h | 18 ++++++++++++------
2 files changed, 15 insertions(+), 6 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f6d04cfaf..101aa15ee 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -75,6 +75,9 @@ public:
case Code::Electron:
is_em = true;
break;
+ case Code::Positron:
+ is_em = true;
+ break;
case Code::Gamma:
is_em = true;
break;
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 93516c141..19f270eb2 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -132,9 +132,15 @@ namespace corsika::process {
corsika::particles::GetLifetime(p.GetPID());
auto const lifetime = gamma * t0;
+ const auto mkin = (E * E - p.GetMomentum().squaredNorm());//delta_mass(p.GetMomentum(), E, m);
cout << "Decay: code: " << p.GetPID() << endl;
cout << "Decay: MinStep: t0: " << t0 << endl;
cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
+ cout << "Decay: momentum: " << p.GetMomentum().GetComponents() / 1_GeV << " GeV" << endl;
+ cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV << " " << m / 1_GeV*m / 1_GeV << endl;
+ //cout << "Decay: sib mass: " << s_mass1_.am2[ process::sibyll::ConvertToSibyllRaw(p.GetPID()) ] << endl;
+ auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ cout << "Decay: sib mass: " << get_sibyll_mass2( sib_id ) << endl;
cout << "Decay: MinStep: gamma: " << gamma << endl;
cout << "Decay: MinStep: tau: " << lifetime << endl;
@@ -154,10 +160,13 @@ namespace corsika::process {
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
+#warning setting particle mass with Corsika values, may be inconsistent with sibyll internal values
+ pin.SetMass( corsika::particles::GetMass( pCode ) );
// remember position
Point decayPoint = p.GetPosition();
TimeType t0 = p.GetTime();
- // remove original particle from corsika stack
+ // remove original particle from corsika stack
p.Delete();
// set all particles/hadrons unstable
// setHadronsUnstable();
@@ -170,15 +179,12 @@ namespace corsika::process {
// 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;
+ if ( psib.HasDecayed() ) 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()));
--
GitLab