From 450c7e0c823f81fbf6191edc579f952d8738c49c Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Mon, 17 Dec 2018 18:30:38 +0100
Subject: [PATCH] correct units in sibyll process interface
---
Documentation/Examples/cascade_example.cc | 10 +-
Framework/Cascade/Cascade.h | 8 +-
Framework/Geometry/Line.h | 2 -
Framework/ProcessSequence/ProcessSequence.h | 2 +-
Processes/Sibyll/Decay.h | 34 +-
Processes/Sibyll/Interaction.h | 426 ++++++++++----------
Processes/Sibyll/ParticleConversion.h | 2 +-
Processes/Sibyll/SibStack.h | 134 +++---
Processes/Sibyll/sibyll2.3c.cc | 1 -
9 files changed, 297 insertions(+), 322 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f8eb408b6..6d36ce8f0 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -200,7 +200,6 @@ public:
private:
};
-
int main() {
corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
@@ -224,13 +223,13 @@ int main() {
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
- corsika::process::sibyll::Interaction/*<setup::Stack>,setup::Trajectory>*/ p1;
+ corsika::process::sibyll::Interaction /*<setup::Stack>,setup::Trajectory>*/ p1;
corsika::process::sibyll::Decay p2;
ProcessEMCut p3;
const auto sequence = /*p0 +*/ p1 + p2 + p3;
setup::Stack stack;
- corsika::cascade::Cascade EAS(tracking, sequence, stack);
+ corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
stack.Clear();
auto particle = stack.NewParticle();
@@ -246,8 +245,9 @@ int main() {
particle.SetPosition(p);
EAS.Init();
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV
- //<< "GeV, particles below energy threshold =" << p1.GetCount()
+ cout << "Result: E0="
+ << E0 / 1_GeV
+ //<< "GeV, particles below energy threshold =" << p1.GetCount()
<< endl;
cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
<< std::endl;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 24f86db66..0edd50502 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -18,7 +18,6 @@
#include <corsika/random/UniformRealDistribution.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/environment/Environment.h>
#include <type_traits>
@@ -58,7 +57,7 @@ namespace corsika::cascade {
void Step(Particle& particle) {
using namespace corsika::units::si;
-
+
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
@@ -158,12 +157,7 @@ namespace corsika::cascade {
Stack& fStack;
corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
-<<<<<<< HEAD
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
-=======
- corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
-
->>>>>>> e9023467d9ae486f2436418f2004da612ebd02a7
};
} // namespace corsika::cascade
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index 5abf841ec..1b45f535a 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -31,8 +31,6 @@ namespace corsika::geometry {
, v0(pV0) {}
Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
-
- Point PositionFromArclength(corsika::units::si::LengthType l) const { return r0 + v0.normalized() * l; }
Point PositionFromArclength(corsika::units::si::LengthType l) const {
return r0 + v0.normalized() * l;
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 30a9c09af..d58c14b50 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -361,7 +361,7 @@ namespace corsika::process {
OPSEQ(DecayProcess, DecayProcess)
template <typename A, typename B>
- struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
+ struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
static const bool value = true;
};
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index cc13040b3..13524a1a0 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -1,9 +1,9 @@
#ifndef _include_corsika_process_sibyll_decay_h_
#define _include_corsika_process_sibyll_decay_h_
-#include <corsika/process/sibyll/SibStack.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/DecayProcess.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/SibStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/particles/ParticleProperties.h>
@@ -57,8 +57,8 @@ namespace corsika::process {
void setTrackedParticlesStable() {
/*
- Sibyll is hadronic generator
- only hadrons decay
+ Sibyll is hadronic generator
+ only hadrons decay
*/
// set particles unstable
setHadronsUnstable();
@@ -71,12 +71,12 @@ namespace corsika::process {
ds.NewParticle().SetPID(particles::Code::KMinus);
ds.NewParticle().SetPID(particles::Code::K0Long);
ds.NewParticle().SetPID(particles::Code::K0Short);
-
+
for (auto& p : ds) {
- int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- // set particle stable by setting table value negative
- s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
- p.Delete();
+ int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ // set particle stable by setting table value negative
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+ p.Delete();
}
}
@@ -84,7 +84,7 @@ namespace corsika::process {
public:
Decay() {}
void Init() {
- setHadronsUnstable();
+ setHadronsUnstable();
setTrackedParticlesStable();
}
@@ -110,11 +110,11 @@ namespace corsika::process {
friend void setHadronsUnstable();
template <typename Particle>
- double GetLifetime(Particle& p) const {
+ corsika::units::si::TimeType GetLifetime(Particle& p) const {
corsika::units::hep::EnergyType E = p.GetEnergy();
corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
- //const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
+ // const MassDensityType density = 1.25e-3 * kilogram / (1_cm * 1_cm * 1_cm);
const double gamma = E / m;
@@ -126,11 +126,11 @@ namespace corsika::process {
// return as column density
// const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
// cout << "Decay: MinStep: x0: " << x0 << endl;
- const double lifetime = gamma*t0 / 1_s;
- cout << "Decay: MinStep: tau: " << lifetime << endl;
- //int a = 1;
- //const double x = -x0 * log(s_rndm_(a));
- //cout << "Decay: next decay: " << x << endl;
+ corsika::units::si::TimeType const lifetime = gamma * t0;
+ cout << "Decay: MinStep: tau: " << lifetime << endl;
+ // int a = 1;
+ // const double x = -x0 * log(s_rndm_(a));
+ // cout << "Decay: next decay: " << x << endl;
return lifetime;
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index c0ceb87dd..87a2a0541 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -6,278 +6,262 @@
//#include <corsika/setup/SetupStack.h>
//#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
using namespace corsika;
using namespace corsika::process::sibyll;
using namespace corsika::units::si;
-
namespace corsika::process::sibyll {
// template <typename Stack, typename Track>
- //template <typename Stack>
- class Interaction : public corsika::process::InteractionProcess<Interaction> { // <Stack,Track>> {
+ // template <typename Stack>
+ class Interaction
+ : public corsika::process::InteractionProcess<Interaction> { // <Stack,Track>> {
- //typedef typename Stack::ParticleType Particle;
- //typedef typename corsika::setup::Stack::ParticleType Particle;
- //typedef corsika::setup::Trajectory Track;
+ // typedef typename Stack::ParticleType Particle;
+ // typedef typename corsika::setup::Stack::ParticleType Particle;
+ // typedef corsika::setup::Trajectory Track;
public:
-
Interaction() {}
~Interaction() {}
-
+
void Init() {
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("s_rndm");
-
+
// test random number generator
std::cout << "Interaction: "
- << " test sequence of random numbers." << std::endl;
+ << " test sequence of random numbers." << std::endl;
int a = 0;
for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl;
-
+
// initialize Sibyll
sibyll_ini_();
}
- //void setTrackedParticlesStable();
+ // void setTrackedParticlesStable();
template <typename Particle, typename Track>
- double GetInteractionLength(Particle& p, Track&) const {
-
- // coordinate system, get global frame of reference
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
- const particles::Code corsikaBeamId = p.GetPID();
-
- // beam particles for sibyll : 1, 2, 3 for p, pi, k
- // read from cross section code table
- int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
-
- bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
-
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- // target nuclei: A < 18
- // FOR NOW: assume target is oxygen
- int kTarget = 16;
-
- EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
- super_stupid::MomentumVector Ptot(
- rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
- // FOR NOW: assume target is at rest
- super_stupid::MomentumVector pTarget(
- rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
- Ptot += p.GetMomentum();
- Ptot += pTarget;
- // calculate cm. energy
- EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * constants::cSquared);
- double Ecm = sqs / 1_GeV;
-
- std::cout << "Interaction: "
- << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
- << " beam can interact:" << kBeam << endl
- << " beam XS code:" << kBeam << endl
- << " beam pid:" << p.GetPID() << endl
- << " target mass number:" << kTarget << std::endl;
-
- double int_length = 0;
- if (kInteraction) {
-
- double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
- double dumdif[3];
-
- if (kTarget == 1)
- sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
- else
- sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
-
- std::cout << "Interaction: "
- << "MinStep: sibyll return: " << prodCrossSection << std::endl;
- CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "Interaction: "
- << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+ corsika::units::si::GrammageType GetInteractionLength(Particle& p, Track&) const {
- const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
- std::cout << "Interaction: "
- << "nucleon mass " << nucleon_mass << std::endl;
- // calculate interaction length in medium
- int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
- // pick random step lenth
- std::cout << "Interaction: "
- << "interaction length (g/cm2): " << int_length << std::endl;
- } else
- int_length = std::numeric_limits<double>::infinity();
-
- /*
- what are the units of the output? slant depth or 3space length?
-
- */
- return int_length;
- //
- // int a = 0;
- // const double next_step = -int_length * log(s_rndm_(a));
- // std::cout << "Interaction: "
- // << "next step (g/cm2): " << next_step << std::endl;
- // return next_step;
- }
-
-
- template <typename Particle, typename Stack>
- corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const {
- cout << "ProcessSibyll: "
- << "DoInteraction: " << p.GetPID() << " interaction? "
- << process::sibyll::CanInteract(p.GetPID()) << endl;
- if (process::sibyll::CanInteract(p.GetPID())) {
- cout << "defining coordinates" << endl;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
- QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
- Point pOrig(rootCS, coordinates);
+ const particles::Code corsikaBeamId = p.GetPID();
+
+ // beam particles for sibyll : 1, 2, 3 for p, pi, k
+ // read from cross section code table
+ int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
+
+ bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
-
- here we need: GetTargetMassNumber() or GetTargetPID()??
- GetTargetMomentum() (zero in EAS)
- */
- // FOR NOW: set target to proton
- int kTarget = 1; // env.GetTargetParticle().GetPID();
-
- cout << "defining target momentum.." << endl;
- // FOR NOW: target is always at rest
- const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
- const auto pTarget = super_stupid::MomentumVector(
- rootCS, 0. * 1_GeV / constants::c, 0. * 1_GeV / constants::c,
- 0. * 1_GeV / constants::c);
- cout << "target momentum (GeV/c): "
- << pTarget.GetComponents() / 1_GeV * constants::c << endl;
- cout << "beam momentum (GeV/c): "
- << p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
-
- // get energy of particle from stack
- /*
- stack is in GeV in lab. frame
- convert to GeV in cm. frame
- (assuming proton at rest as target AND
- assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
- */
- // total energy: E_beam + E_target
- // in lab. frame: E_beam + m_target*c**2
- EnergyType E = p.GetEnergy();
- EnergyType Etot = E + Etarget;
- // total momentum
- super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
- // invariant mass, i.e. cm. energy
- EnergyType Ecm = sqrt(Etot * Etot -
- Ptot.squaredNorm() *
- constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
- /*
- 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
- */
- const double gamma = Etot / Ecm;
- const auto gambet = Ptot / (Ecm / constants::c);
+ */
+ // target nuclei: A < 18
+ // FOR NOW: assume target is oxygen
+ int kTarget = 16;
+
+ EnergyType Etot = p.GetEnergy() + kTarget * corsika::particles::Proton::GetMass();
+ super_stupid::MomentumVector Ptot(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ // FOR NOW: assume target is at rest
+ super_stupid::MomentumVector pTarget(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ Ptot += p.GetMomentum();
+ Ptot += pTarget;
+ // calculate cm. energy
+ EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm() * constants::cSquared);
+ double Ecm = sqs / 1_GeV;
std::cout << "Interaction: "
- << " DoDiscrete: gamma:" << gamma << endl;
- std::cout << "Interaction: "
- << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+ << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
+ << " beam can interact:" << kBeam << endl
+ << " beam XS code:" << kBeam << endl
+ << " beam pid:" << p.GetPID() << endl
+ << " target mass number:" << kTarget << std::endl;
- int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ if (kInteraction) {
- std::cout << "Interaction: "
- << " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
- << std::endl;
- if (E < 8.5_GeV || Ecm < 10_GeV) {
+ double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
+ double dumdif[3];
+
+ if (kTarget == 1)
+ sib_sigma_hp_(kBeam, Ecm, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
+ else
+ sib_sigma_hnuc_(kBeam, kTarget, Ecm, prodCrossSection, dummy);
+
+ std::cout << "Interaction: "
+ << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+ CrossSectionType sig = prodCrossSection * 1_mbarn;
+ std::cout << "Interaction: "
+ << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
+
+ const MassType nucleon_mass =
+ 0.93827_GeV / corsika::units::si::constants::cSquared;
+ std::cout << "Interaction: "
+ << "nucleon mass " << nucleon_mass << std::endl;
+ // calculate interaction length in medium
+ GrammageType int_length = kTarget * nucleon_mass / sig;
+ // pick random step lenth
std::cout << "Interaction: "
- << " DoInteraction: dropping particle.." << std::endl;
- p.Delete();
+ << "interaction length (g/cm2): " << int_length << std::endl;
+
+ return int_length;
} else {
- // Sibyll does not know about units..
- double sqs = Ecm / 1_GeV;
- // running sibyll, filling stack
- sibyll_(kBeam, kTarget, sqs);
- // running decays
- //setTrackedParticlesStable();
- 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;
-
- // SibStack does not know about momentum yet so we need counter to access momentum
- // array in Sibyll
- int i = -1;
- super_stupid::MomentumVector Ptot_final(
- rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
- for (auto& psib : ss) {
- ++i;
- // skip particles that have decayed in Sibyll
- if (abs(s_plist_.llist[i]) > 100) continue;
-
- // transform energy to lab. frame, primitve
- // compute beta_vec * p_vec
- // arbitrary Lorentz transformation based on sibyll routines
- const auto gammaBetaComponents = gambet.GetComponents();
- const auto pSibyllComponents = psib.GetMomentum().GetComponents();
- EnergyType en_lab = 0. * 1_GeV;
- MomentumType p_lab_components[3];
- en_lab = psib.GetEnergy() * gamma;
- EnergyType pnorm = 0. * 1_GeV;
- for (int j = 0; j < 3; ++j)
- pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * constants::c) /
- (gamma + 1.);
- pnorm += psib.GetEnergy();
-
- for (int j = 0; j < 3; ++j) {
- p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
- gammaBetaComponents[j] /
- constants::c;
- en_lab -=
- (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * constants::c;
- }
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+ }
- // add to corsika stack
- auto pnew = s.NewParticle();
- pnew.SetEnergy(en_lab);
- pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+ template <typename Particle, typename Stack>
+ corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) const {
+ cout << "ProcessSibyll: "
+ << "DoInteraction: " << p.GetPID() << " interaction? "
+ << process::sibyll::CanInteract(p.GetPID()) << endl;
+ if (process::sibyll::CanInteract(p.GetPID())) {
+ cout << "defining coordinates" << endl;
+ // coordinate system, get global frame of reference
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
+ QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+ Point pOrig(rootCS, coordinates);
+
+ /*
+ the target should be defined by the Environment,
+ ideally as full particle object so that the four momenta
+ and the boosts can be defined..
+
+ here we need: GetTargetMassNumber() or GetTargetPID()??
+ GetTargetMomentum() (zero in EAS)
+ */
+ // FOR NOW: set target to proton
+ int kTarget = 1; // env.GetTargetParticle().GetPID();
+
+ cout << "defining target momentum.." << endl;
+ // FOR NOW: target is always at rest
+ const EnergyType Etarget = 0. * 1_GeV + corsika::particles::Proton::GetMass();
+ const auto pTarget = super_stupid::MomentumVector(
+ rootCS, 0. * 1_GeV / constants::c, 0. * 1_GeV / constants::c,
+ 0. * 1_GeV / constants::c);
+ cout << "target momentum (GeV/c): "
+ << pTarget.GetComponents() / 1_GeV * constants::c << endl;
+ cout << "beam momentum (GeV/c): "
+ << p.GetMomentum().GetComponents() / 1_GeV * constants::c << endl;
+
+ // get energy of particle from stack
+ /*
+ stack is in GeV in lab. frame
+ convert to GeV in cm. frame
+ (assuming proton at rest as target AND
+ assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
+ */
+ // total energy: E_beam + E_target
+ // in lab. frame: E_beam + m_target*c**2
+ EnergyType E = p.GetEnergy();
+ EnergyType Etot = E + Etarget;
+ // total momentum
+ super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+ // invariant mass, i.e. cm. energy
+ EnergyType Ecm =
+ sqrt(Etot * Etot - Ptot.squaredNorm() *
+ constants::cSquared); // sqrt( 2. * E * 0.93827_GeV );
+ /*
+ 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
+ */
+ const double gamma = Etot / Ecm;
+ const auto gambet = Ptot / (Ecm / constants::c);
+
+ std::cout << "Interaction: "
+ << " DoDiscrete: gamma:" << gamma << endl;
+ std::cout << "Interaction: "
+ << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+
+ int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- corsika::geometry::QuantityVector<momentum_d> p_lab_c{
- p_lab_components[0], p_lab_components[1], p_lab_components[2]};
- pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c));
- Ptot_final += pnew.GetMomentum();
+ std::cout << "Interaction: "
+ << " DoInteraction: E(GeV):" << E / 1_GeV
+ << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
+ if (E < 8.5_GeV || Ecm < 10_GeV) {
+ std::cout << "Interaction: "
+ << " DoInteraction: dropping particle.." << std::endl;
+ p.Delete();
+ } else {
+ // Sibyll does not know about units..
+ double sqs = Ecm / 1_GeV;
+ // running sibyll, filling stack
+ sibyll_(kBeam, kTarget, sqs);
+ // running decays
+ // setTrackedParticlesStable();
+ 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;
+
+ // SibStack does not know about momentum yet so we need counter to access
+ // momentum array in Sibyll
+ int i = -1;
+ super_stupid::MomentumVector Ptot_final(rootCS, {0.0_Ns, 0.0_Ns, 0.0_Ns});
+ for (auto& psib : ss) {
+ ++i;
+ // skip particles that have decayed in Sibyll
+ if (abs(s_plist_.llist[i]) > 100) continue;
+
+ // transform energy to lab. frame, primitve
+ // compute beta_vec * p_vec
+ // arbitrary Lorentz transformation based on sibyll routines
+ const auto gammaBetaComponents = gambet.GetComponents();
+ const auto pSibyllComponents = psib.GetMomentum().GetComponents();
+ EnergyType en_lab = 0. * 1_GeV;
+ MomentumType p_lab_components[3];
+ en_lab = psib.GetEnergy() * gamma;
+ EnergyType pnorm = 0. * 1_GeV;
+ for (int j = 0; j < 3; ++j)
+ pnorm += (pSibyllComponents[j] * gammaBetaComponents[j] * constants::c) /
+ (gamma + 1.);
+ pnorm += psib.GetEnergy();
+
+ for (int j = 0; j < 3; ++j) {
+ p_lab_components[j] = pSibyllComponents[j] -
+ (-1) * pnorm * gammaBetaComponents[j] / constants::c;
+ en_lab -=
+ (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * constants::c;
+ }
+
+ // add to corsika stack
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy(en_lab);
+ pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+
+ corsika::geometry::QuantityVector<momentum_d> p_lab_c{
+ p_lab_components[0], p_lab_components[1], p_lab_components[2]};
+ pnew.SetMomentum(super_stupid::MomentumVector(rootCS, p_lab_c));
+ Ptot_final += pnew.GetMomentum();
+ }
+ // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV
+ // * constants::c << endl;
}
- // cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV *
- // constants::c << endl;
}
+ return process::EProcessReturn::eOk;
}
- return process::EProcessReturn::eOk;
- }
-
};
-
-}
+
+} // namespace corsika::process::sibyll
#endif
diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h
index e5e76c3b3..cd96f83d5 100644
--- a/Processes/Sibyll/ParticleConversion.h
+++ b/Processes/Sibyll/ParticleConversion.h
@@ -25,7 +25,7 @@ namespace corsika::process::sibyll {
#include <corsika/process/sibyll/Generated.inc>
- bool constexpr KnownBySibyll(corsika::particles::Code pCode) {
+ bool constexpr KnownBySibyll(corsika::particles::Code pCode) {
return isKnown[static_cast<corsika::particles::CodeIntType>(pCode)];
}
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index 4c0326a92..f79fe7da3 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -1,8 +1,8 @@
#ifndef _include_sibstack_h_
#define _include_sibstack_h_
-#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
@@ -14,75 +14,75 @@ using namespace corsika::units::si;
using namespace corsika::geometry;
namespace corsika::process::sibyll {
-
-class SibStackData {
-public:
- void Init();
+ class SibStackData {
+
+ public:
+ void Init();
- void Clear() { s_plist_.np = 0; }
+ void Clear() { s_plist_.np = 0; }
- int GetSize() const { return s_plist_.np; }
+ int GetSize() const { return s_plist_.np; }
#warning check actual capacity of sibyll stack
- int GetCapacity() const { return 8000; }
-
- void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
- void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
- void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
- auto tmp = v.GetComponents();
- for (int idx = 0; idx < 3; ++idx)
- s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
- }
-
- int GetId(const int i) const { return s_plist_.llist[i]; }
-
- EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
-
- super_stupid::MomentumVector GetMomentum(const int i) const {
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
- corsika::geometry::QuantityVector<momentum_d> components{
- s_plist_.p[0][i] * 1_GeV / constants::c,
- s_plist_.p[1][i] * 1_GeV / constants::c,
- s_plist_.p[2][i] * 1_GeV / constants::c};
- super_stupid::MomentumVector v1(rootCS, components);
- return v1;
- }
-
- 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 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);
- }
-};
-
-typedef Stack<SibStackData, ParticleInterface> SibStack;
-
-}
+ int GetCapacity() const { return 8000; }
+
+ void SetId(const int i, const int v) { s_plist_.llist[i] = v; }
+ void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV; }
+ void SetMomentum(const int i, const super_stupid::MomentumVector& v) {
+ auto tmp = v.GetComponents();
+ for (int idx = 0; idx < 3; ++idx)
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
+ }
+
+ int GetId(const int i) const { return s_plist_.llist[i]; }
+
+ EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+
+ super_stupid::MomentumVector GetMomentum(const int i) const {
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+ corsika::geometry::QuantityVector<momentum_d> components{
+ s_plist_.p[0][i] * 1_GeV / constants::c,
+ s_plist_.p[1][i] * 1_GeV / constants::c,
+ s_plist_.p[2][i] * 1_GeV / constants::c};
+ super_stupid::MomentumVector v1(rootCS, components);
+ return v1;
+ }
+
+ 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 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);
+ }
+ };
+
+ typedef Stack<SibStackData, ParticleInterface> SibStack;
+
+} // namespace corsika::process::sibyll
#endif
diff --git a/Processes/Sibyll/sibyll2.3c.cc b/Processes/Sibyll/sibyll2.3c.cc
index 336617fbf..77ad30527 100644
--- a/Processes/Sibyll/sibyll2.3c.cc
+++ b/Processes/Sibyll/sibyll2.3c.cc
@@ -2,7 +2,6 @@
#include <corsika/random/RNGManager.h>
-
double s_rndm_(int&) {
static corsika::random::RNG& rmng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
--
GitLab