IAP GITLAB

Skip to content
Snippets Groups Projects
Commit a888bf32 authored by Ralf Ulrich's avatar Ralf Ulrich
Browse files

Merge branch '83-move-sibyll-stuff-from-framework-cascade-to-where-it-belongs' into 'master' 

Resolve "move SIBYLL stuff from Framework/Cascade to where it belongs"

Closes #83

See merge request AirShowerPhysics/corsika!30
parents 03310d22 f6bbc34f
No related branches found
No related tags found
1 merge request!30Resolve "move SIBYLL stuff from Framework/Cascade to where it belongs"
Pipeline #42 passed
Stage: build
Hide whitespace changes
Inline Side-by-side
Showing
with 416 additions and 414 deletions
CMakeLists.txt
+ 4
− 6
View file @ a888bf32
......@@ -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)
......
Documentation/Examples/CMakeLists.txt
+ 1
− 1
View file @ a888bf32
......@@ -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
......
Documentation/Examples/cascade_example.cc
+ 32
− 330
View file @ a888bf32
......@@ -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
......@@ -58,7 +53,7 @@ static int fEmCount;
static EnergyType fInvEnergy;
static int fInvCount;
class ProcessEMCut : public corsika::process::BaseProcess<ProcessEMCut> {
class ProcessEMCut : public corsika::process::ContinuousProcess<ProcessEMCut> {
public:
ProcessEMCut() {}
template <typename Particle>
......@@ -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;
......@@ -134,7 +129,24 @@ public:
}
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
const Code pid = p.GetPID();
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += p.GetEnergy();
fEmCount += 1;
p.Delete();
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += p.GetEnergy();
fInvCount += 1;
p.Delete();
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += p.GetEnergy();
p.Delete();
}
// cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
// cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
// cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
......@@ -156,28 +168,6 @@ public:
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack&) const {
cout << "ProcessCut: DoDiscrete: " << p.GetPID() << endl;
const Code pid = p.GetPID();
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += p.GetEnergy();
fEmCount += 1;
p.Delete();
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += p.GetEnergy();
fInvCount += 1;
p.Delete();
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += p.GetEnergy();
fCount += 1;
p.Delete();
}
}
void Init() {
fEmEnergy = 0. * 1_GeV;
fEmCount = 0;
......@@ -206,296 +196,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,10 +220,10 @@ int main() {
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
corsika::process::sibyll::ProcessDecay p2;
ProcessEMCut p3;
const auto sequence = /*p0 +*/ p1 + p2 + p3;
corsika::process::sibyll::Interaction sibyll;
corsika::process::sibyll::Decay decay;
ProcessEMCut cut;
const auto sequence = /*p0 +*/ sibyll + decay + cut;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
......@@ -543,9 +243,11 @@ 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;
p3.ShowResults();
//<< "GeV, particles below energy threshold =" << p1.GetCount()
<< endl;
cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
<< std::endl;
cut.ShowResults();
cout << "total energy (GeV): "
<< (p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy()) / 1_GeV << endl;
<< (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy()) / 1_GeV << endl;
}
Framework/Cascade/CMakeLists.txt
+ 1
− 24
View file @ a888bf32
......@@ -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
......
Framework/ProcessSequence/ContinuousProcess.h
+ 2
− 2
View file @ a888bf32
......@@ -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
......
Framework/ProcessSequence/InteractionProcess.h
+ 4
− 4
View file @ a888bf32
......@@ -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 {
......
Framework/ProcessSequence/ProcessSequence.h
+ 2
− 2
View file @ a888bf32
......@@ -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;
};
......
Processes/Sibyll/CMakeLists.txt
+ 8
− 1
View file @ a888bf32
......@@ -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 (
......
Processes/Sibyll/ProcessDecay.h → Processes/Sibyll/Decay.h
+ 51
− 30
View file @ a888bf32
#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
......
Processes/Sibyll/Interaction.h 0 → 100644
+ 283
− 0
View file @ a888bf32
#ifndef _corsika_process_sibyll_interaction_h_
#define _corsika_process_sibyll_interaction_h_
#include <corsika/process/InteractionProcess.h>
//#include <corsika/setup/SetupStack.h>
//#include <corsika/setup/SetupTrajectory.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>
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>> {
//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;
int a = 0;
for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl;
// initialize Sibyll
sibyll_ini_();
}
//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;
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);
/*
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());
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_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;
}
// 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;
}
}
return process::EProcessReturn::eOk;
}
};
}
#endif
Processes/Sibyll/ParticleConversion.h
+ 4
− 4
View file @ a888bf32
......@@ -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)];
}
......
Framework/Cascade/SibStack.h → Processes/Sibyll/SibStack.h
+ 12
− 9
View file @ a888bf32
#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
Processes/Sibyll/code_generator.py
+ 1
− 1
View file @ a888bf32
......@@ -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:
......
Framework/Cascade/gasdev.f → Processes/Sibyll/gasdev.f
+ 0
− 0
View file @ a888bf32
File moved
Framework/Cascade/rndm_dbl.f → Processes/Sibyll/rndm_dbl.f
+ 0
− 0
View file @ a888bf32
File moved
Processes/Sibyll/sibyll2.3c.cc 0 → 100644
+ 11
− 0
View file @ a888bf32
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/random/RNGManager.h>
double s_rndm_(int&) {
static corsika::random::RNG& rmng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
;
return rmng() / (double)rmng.max();
}
Framework/Cascade/sibyll2.3c.f → Processes/Sibyll/sibyll2.3c.f
+ 0
− 0
View file @ a888bf32
File moved
Framework/Cascade/sibyll2.3c.h → Processes/Sibyll/sibyll2.3c.h
+ 0
− 0
View file @ a888bf32
File moved
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum