diff --git a/CMakeLists.txt b/CMakeLists.txt
index b6969bead512aacffc631f993a6a77d6e4f06bff..841f3a9d7b514c2cad0539cd85c2b9a55a07d2b7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -17,9 +17,6 @@ set (CMAKE_INSTALL_MESSAGE LAZY)
set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/CMakeModules)
include (CorsikaUtilities) # a few cmake function
-# enable warnings and disallow non-standard language
-add_definitions(-Wall -pedantic -Wextra -Wno-ignored-qualifiers)
-
set (CMAKE_CXX_STANDARD 17)
enable_testing ()
set (CTEST_OUTPUT_ON_FAILURE 1)
@@ -39,9 +36,10 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
-#set(CMAKE_CXX_FLAGS "-Wall -Wextra")
-set(CMAKE_CXX_FLAGS_DEBUG "-O0 -g")
-set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+# enable warnings and disallow non-standard language
+set(CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers")
+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS} -O0 -g")
+set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS} -O3")
# unit testing coverage, does not work yet
#include (CodeCoverage)
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 148df3493a54be8a98442a2f065a6c7a11858591..2459d4d9a7faf03b3a2d8fdb03b562db7848a4a7 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -25,7 +25,7 @@ add_executable (cascade_example cascade_example.cc)
target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging
CORSIKArandom
- CORSIKAsibyll
+ ProcessSibyll
CORSIKAcascade
ProcessStackInspector
CORSIKAprocesses
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index f9f344d4470a9a5638a00d8c041c87eb90a0b27b..f8eb408b62c4e199e8c8ca1cbf4a2a93677eb8ff 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -20,14 +20,10 @@
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
-#include <corsika/random/RNGManager.h>
-
-#include <corsika/cascade/SibStack.h>
-#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/geometry/Sphere.h>
-#include <corsika/process/sibyll/ParticleConversion.h>
-#include <corsika/process/sibyll/ProcessDecay.h>
+#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
#include <corsika/units/PhysicalUnits.h>
@@ -47,7 +43,6 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
-static int fCount = 0;
static EnergyType fEnergy = 0. * 1_GeV;
// FOR NOW: global static variables for ParticleCut process
@@ -121,7 +116,7 @@ public:
}
template <typename Particle>
- double MinStepLength(Particle& p, setup::Trajectory&) const {
+ double MaxStepLength(Particle& p, setup::Trajectory&) const {
const Code pid = p.GetPID();
if (isEmParticle(pid) || isInvisible(pid)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
@@ -173,7 +168,6 @@ public:
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += p.GetEnergy();
- fCount += 1;
p.Delete();
}
}
@@ -206,296 +200,6 @@ public:
private:
};
-class ProcessSplit : public corsika::process::InteractionProcess<ProcessSplit> {
-public:
- ProcessSplit() {}
-
- void setTrackedParticlesStable() const {
- /*
- Sibyll is hadronic generator
- only hadrons decay
- */
- // set particles unstable
- corsika::process::sibyll::setHadronsUnstable();
- // make tracked particles stable
- std::cout << "ProcessSplit: setting tracked hadrons stable.." << std::endl;
- setup::Stack ds;
- 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
- s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
- p.Delete();
- }
- }
-
- template <typename Particle, typename Track>
- double GetInteractionLength(Particle& p, Track&) const {
-
- // coordinate system, get global frame of reference
- CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
-
- const 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() * si::constants::cSquared);
- double Ecm = sqs / 1_GeV;
-
- std::cout << "ProcessSplit: "
- << "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 << "ProcessSplit: "
- << "MinStep: sibyll return: " << prodCrossSection << std::endl;
- CrossSectionType sig = prodCrossSection * 1_mbarn;
- std::cout << "ProcessSplit: "
- << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
-
- const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
- std::cout << "ProcessSplit: "
- << "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 << "ProcessSplit: "
- << "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 << "ProcessSplit: "
- // << "next step (g/cm2): " << next_step << std::endl;
- // return next_step;
- }
-
- template <typename Particle, typename Track, typename Stack>
- EProcessReturn DoContinuous(Particle&, Track&, Stack&) const {
- return EProcessReturn::eOk;
- }
-
- template <typename Particle, typename Stack>
- void 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 / si::constants::c, 0. * 1_GeV / si::constants::c,
- 0. * 1_GeV / si::constants::c);
- cout << "target momentum (GeV/c): "
- << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
- cout << "beam momentum (GeV/c): "
- << p.GetMomentum().GetComponents() / 1_GeV * si::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() *
- si::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 / si::constants::c);
-
- std::cout << "ProcessSplit: "
- << " DoDiscrete: gamma:" << gamma << endl;
- std::cout << "ProcessSplit: "
- << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
-
- int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
-
- std::cout << "ProcessSplit: "
- << " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
- << std::endl;
- if (E < 8.5_GeV || Ecm < 10_GeV) {
- std::cout << "ProcessSplit: "
- << " DoInteraction: dropping particle.." << std::endl;
- p.Delete();
- fCount++;
- } 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] * si::constants::c) /
- (gamma + 1.);
- pnorm += psib.GetEnergy();
-
- for (int j = 0; j < 3; ++j) {
- p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm *
- gammaBetaComponents[j] /
- si::constants::c;
- en_lab -=
- (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::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 *
- // si::constants::c << endl;
- }
- }
- }
-
- void Init() {
- fCount = 0;
-
- corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
-
- rmng.RegisterRandomStream("s_rndm");
-
- // test random number generator
- std::cout << "ProcessSplit: "
- << " 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_();
-
- setTrackedParticlesStable();
- }
-
- int GetCount() { return fCount; }
- EnergyType GetEnergy() { return fEnergy; }
-
-private:
-};
-
-double s_rndm_(int&) {
- static corsika::random::RNG& rmng =
- corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
- ;
- return rmng() / (double)rmng.max();
-}
int main() {
corsika::environment::Environment env; // dummy environment
@@ -520,8 +224,8 @@ int main() {
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
- ProcessSplit p1;
- corsika::process::sibyll::ProcessDecay p2;
+ 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;
@@ -543,8 +247,10 @@ int main() {
EAS.Init();
EAS.Run();
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;
+ //<< "GeV, particles below energy threshold =" << p1.GetCount()
+ << endl;
+ cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
+ << std::endl;
p3.ShowResults();
cout << "total energy (GeV): "
<< (p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy()) / 1_GeV << endl;
diff --git a/Framework/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt
index ad48313530b2df7dfcbda79cf8dcaf2791a0bf18..8da2afe065f5fce33f2fa3a9f9a695a930471a09 100644
--- a/Framework/Cascade/CMakeLists.txt
+++ b/Framework/Cascade/CMakeLists.txt
@@ -9,33 +9,10 @@ set (
set (
CORSIKAcascade_HEADERS
Cascade.h
- sibyll2.3c.h
- SibStack.h
)
-#set (
-# CORSIKAcascade_SOURCES
-# TrackingStep.cc
-# Cascade.cc
-# )
-
-set (
- CORSIKAsibyll_NAMESPACE
- corsika/cascade
- )
-
-set (
- CORSIKAsibyll_SOURCES
- sibyll2.3c.f
- gasdev.f
- )
-add_library (CORSIKAsibyll STATIC ${CORSIKAsibyll_SOURCES})
-
-#add_library (CORSIKAcascade STATIC ${CORSIKAcascade_SOURCES})
add_library (CORSIKAcascade INTERFACE)
-
-
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
#target_link_libraries (
@@ -70,7 +47,7 @@ target_link_libraries (
testCascade
# CORSIKAutls
CORSIKArandom
- CORSIKAsibyll
+ ProcessSibyll
CORSIKAcascade
ProcessStackInspector
CORSIKAstackinterface
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index a9bc21e5db1c6dcbb2ff8e13110ec957e6242e39..c578f235244e62bb54729f3bf5d6014706c01fa9 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/Framework/ProcessSequence/ContinuousProcess.h
@@ -33,8 +33,8 @@ namespace corsika::process {
// here starts the interface part
// -> enforce derived to implement DoContinuous...
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D&, T&, S&) const;
+ template <typename P, typename T, typename S>
+ inline EProcessReturn DoContinuous(P&, T&, S&) const;
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
index 40264e606a7de8ecf607fa8d5173190e9e4ca0e4..b0df1eebb5cc88cdb6a6b207906e0a1f2699bed1 100644
--- a/Framework/ProcessSequence/InteractionProcess.h
+++ b/Framework/ProcessSequence/InteractionProcess.h
@@ -34,11 +34,11 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce derived to implement DoInteraction...
- template <typename Particle, typename Stack>
- inline EProcessReturn DoInteraction(Particle&, Stack&) const;
+ template <typename P, typename S>
+ inline EProcessReturn DoInteraction(P&, S&) const;
- template <typename Particle, typename Track>
- inline double GetInteractionLength(Particle& p, Track& t) const;
+ template <typename P, typename T>
+ inline double GetInteractionLength(P&, T&) const;
template <typename Particle, typename Track>
inline double GetInverseInteractionLength(Particle& p, Track& t) const {
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index 8517d4345202f8fb667f830eab3c093937b2ad05..014630e1a46c076f5f8347686f1d047771c0b29f 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -348,8 +348,8 @@ namespace corsika::process {
OPSEQ(DecayProcess, ContinuousProcess)
OPSEQ(DecayProcess, DecayProcess)
- template <template <typename, typename> class T, typename A, typename B>
- struct is_process_sequence<T<A, B> > {
+ template <typename A, typename B>
+ struct is_process_sequence<corsika::process::ProcessSequence<A, B> > {
static const bool value = true;
};
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index a0985a1ad2304e7e5bd4f1024cf43b5cefaf8bc4..31d204714dba4aa72301e97ab3f5402e1abdeb73 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -15,12 +15,18 @@ add_custom_command (
set (
MODEL_SOURCES
ParticleConversion.cc
+ sibyll2.3c.f
+ sibyll2.3c.cc
+ gasdev.f
)
set (
MODEL_HEADERS
ParticleConversion.h
- ProcessDecay.h
+ sibyll2.3c.h
+ SibStack.h
+ Decay.h
+ Interaction.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
@@ -61,6 +67,7 @@ target_link_libraries (
CORSIKAparticles
CORSIKAunits
CORSIKAthirdparty
+ CORSIKAgeometry
)
target_include_directories (
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 701094e733ad463f673c8d73cdee0bdf34bdddab..cc13040b39101ccc0a6df13a93e1ae84fa696fbd 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -1,13 +1,12 @@
-#ifndef _include_ProcessDecay_h_
-#define _include_ProcessDecay_h_
+#ifndef _include_corsika_process_sibyll_decay_h_
+#define _include_corsika_process_sibyll_decay_h_
-#include <corsika/process/ContinuousProcess.h>
-
-#include <corsika/cascade/SibStack.h>
+#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/DecayProcess.h>
#include <corsika/setup/SetupTrajectory.h>
-// using namespace corsika::particles;
+#include <corsika/particles/ParticleProperties.h>
namespace corsika::process {
@@ -56,11 +55,37 @@ namespace corsika::process {
}
}
- class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+ void setTrackedParticlesStable() {
+ /*
+ Sibyll is hadronic generator
+ only hadrons decay
+ */
+ // set particles unstable
+ setHadronsUnstable();
+ // make tracked particles stable
+ std::cout << "Interaction: setting tracked hadrons stable.." << std::endl;
+ setup::Stack ds;
+ ds.NewParticle().SetPID(particles::Code::PiPlus);
+ ds.NewParticle().SetPID(particles::Code::PiMinus);
+ ds.NewParticle().SetPID(particles::Code::KPlus);
+ 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();
+ }
+ }
+
+ class Decay : public corsika::process::DecayProcess<Decay> {
public:
- ProcessDecay() {}
+ Decay() {}
void Init() {
- // setHadronsUnstable();
+ setHadronsUnstable();
+ setTrackedParticlesStable();
}
void setAllStable() {
@@ -85,31 +110,32 @@ namespace corsika::process {
friend void setHadronsUnstable();
template <typename Particle>
- double MinStepLength(Particle& p, setup::Trajectory&) const {
+ double GetLifetime(Particle& p) const {
corsika::units::hep::EnergyType E = p.GetEnergy();
corsika::units::hep::MassType m = corsika::particles::GetMass(p.GetPID());
- // env.GetDensity();
- 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;
- TimeType t0 = GetLifetime(p.GetPID());
- cout << "ProcessDecay: code: " << (p.GetPID()) << endl;
- cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
- cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
- cout << "ProcessDecay: MinStep: density: " << density << endl;
+ const TimeType t0 = particles::GetLifetime(p.GetPID());
+ cout << "Decay: code: " << (p.GetPID()) << endl;
+ cout << "Decay: MinStep: t0: " << t0 << endl;
+ cout << "Decay: MinStep: gamma: " << gamma << endl;
+ // cout << "Decay: MinStep: density: " << density << endl;
// return as column density
- const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
- cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
- int a = 1;
- const double x = -x0 * log(s_rndm_(a));
- cout << "ProcessDecay: next decay: " << x << endl;
- return x;
+ // 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;
+ return lifetime;
}
template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const {
+ void DoDecay(Particle& p, Stack& s) const {
SibStack ss;
ss.Clear();
// copy particle to sibyll stack
@@ -122,7 +148,7 @@ namespace corsika::process {
// set all particles/hadrons unstable
setHadronsUnstable();
// call sibyll decay
- std::cout << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
+ std::cout << "Decay: calling Sibyll decay routine.." << std::endl;
decsib_();
// print output
int print_unit = 6;
@@ -144,11 +170,6 @@ namespace corsika::process {
// empty sibyll stack
ss.Clear();
}
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
- }
};
} // namespace sibyll
} // namespace corsika::process
diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h
index 35bfbc617df904ba0e859fbf451ab1356092f207..e5e76c3b398e2d93a7af8b8fbe652dacf37ea8cb 100644
--- a/Processes/Sibyll/ParticleConversion.h
+++ b/Processes/Sibyll/ParticleConversion.h
@@ -25,11 +25,11 @@ namespace corsika::process::sibyll {
#include <corsika/process/sibyll/Generated.inc>
- bool KnownBySibyll(corsika::particles::Code pCode) {
+ bool constexpr KnownBySibyll(corsika::particles::Code pCode) {
return isKnown[static_cast<corsika::particles::CodeIntType>(pCode)];
}
- bool CanInteract(corsika::particles::Code pCode) {
+ bool constexpr CanInteract(corsika::particles::Code pCode) {
return canInteract[static_cast<corsika::particles::CodeIntType>(pCode)];
}
@@ -43,12 +43,12 @@ namespace corsika::process::sibyll {
return sibyll2corsika[static_cast<SibyllCodeIntType>(pCode) - minSibyll];
}
- int ConvertToSibyllRaw(corsika::particles::Code pCode) {
+ int constexpr ConvertToSibyllRaw(corsika::particles::Code pCode) {
return (int)static_cast<corsika::process::sibyll::SibyllCodeIntType>(
corsika::process::sibyll::ConvertToSibyll(pCode));
}
- int GetSibyllXSCode(corsika::particles::Code pCode) {
+ int constexpr GetSibyllXSCode(corsika::particles::Code pCode) {
return corsika2sibyllXStype[static_cast<corsika::particles::CodeIntType>(pCode)];
}
diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h
index a25316553da2182f6da9b2668f7f67336ede9a84..4c0326a9268ce61abc821e94bd4fd3861a985d5f 100644
--- a/Processes/Sibyll/SibStack.h
+++ b/Processes/Sibyll/SibStack.h
@@ -1,19 +1,20 @@
#ifndef _include_sibstack_h_
#define _include_sibstack_h_
-#include <string>
-#include <vector>
-
-#include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/stack/super_stupid/SuperStupidStack.h>
+
using namespace std;
using namespace corsika::stack;
-using namespace corsika::units;
+using namespace corsika::units::si;
using namespace corsika::geometry;
+namespace corsika::process::sibyll {
+
class SibStackData {
public:
@@ -30,7 +31,7 @@ public:
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 * si::constants::c;
+ s_plist_.p[idx][i] = tmp[idx] / 1_GeV * constants::c;
}
int GetId(const int i) const { return s_plist_.llist[i]; }
@@ -40,9 +41,9 @@ public:
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 / si::constants::c,
- s_plist_.p[1][i] * 1_GeV / si::constants::c,
- s_plist_.p[2][i] * 1_GeV / si::constants::c};
+ 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;
}
@@ -82,4 +83,6 @@ public:
typedef Stack<SibStackData, ParticleInterface> SibStack;
+}
+
#endif
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
index a593a1f2ce5a5c3a281c758e2404c72347286ec0..388668e1de0d5b04c8b5fa920308d20ab9982d39 100755
--- a/Processes/Sibyll/code_generator.py
+++ b/Processes/Sibyll/code_generator.py
@@ -89,7 +89,7 @@ def generate_sibyll2corsika(pythia_db) :
pDict[sib_code] = identifier
nPart = max(pDict.keys()) - min(pDict.keys()) + 1
- string += "std::array<corsika::particles::Code, {:d}> sibyll2corsika = {{\n".format(nPart)
+ string += "std::array<corsika::particles::Code, {:d}> constexpr sibyll2corsika = {{\n".format(nPart)
for iPart in range(nPart) :
if iPart in pDict: