diff --git a/CMakeLists.txt b/CMakeLists.txt
index 62d488fa2713a44964460d939aa6a5742b1840d8..b6969bead512aacffc631f993a6a77d6e4f06bff 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,6 +7,7 @@ project (
LANGUAGES CXX
)
+# as long as there still are modules using it:
enable_language (Fortran)
# ignore many irrelevant Up-to-date messages during install
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 49ab47047d476934a45de114f015cbfeea2c260b..f9f344d4470a9a5638a00d8c041c87eb90a0b27b 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -27,6 +27,8 @@
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/ProcessDecay.h>
+
#include <corsika/units/PhysicalUnits.h>
#include <iostream>
@@ -38,6 +40,7 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
+using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
@@ -45,16 +48,206 @@ 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
+// this is just wrong...
+static EnergyType fEmEnergy;
+static int fEmCount;
+
+static EnergyType fInvEnergy;
+static int fInvCount;
+
+class ProcessEMCut : public corsika::process::BaseProcess<ProcessEMCut> {
+public:
+ ProcessEMCut() {}
+ template <typename Particle>
+ bool isBelowEnergyCut(Particle& p) const {
+ // FOR NOW: center-of-mass energy hard coded
+ const EnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
+ if (p.GetEnergy() < 50_GeV || Ecm < 10_GeV)
+ return true;
+ else
+ return false;
+ }
+
+ bool isEmParticle(Code pCode) const {
+ bool is_em = false;
+ // FOR NOW: switch
+ switch (pCode) {
+ case Code::Electron:
+ is_em = true;
+ break;
+ case Code::Gamma:
+ is_em = true;
+ break;
+ default:
+ break;
+ }
+ return is_em;
+ }
+
+ void defineEmParticles() const {
+ // create bool array identifying em particles
+ }
+
+ bool isInvisible(Code pCode) const {
+ bool is_inv = false;
+ // FOR NOW: switch
+ switch (pCode) {
+ case Code::NuE:
+ is_inv = true;
+ break;
+ case Code::NuEBar:
+ is_inv = true;
+ break;
+ case Code::NuMu:
+ is_inv = true;
+ break;
+ case Code::NuMuBar:
+ is_inv = true;
+ break;
+ case Code::MuPlus:
+ is_inv = true;
+ break;
+ case Code::MuMinus:
+ is_inv = true;
+ break;
+
+ default:
+ break;
+ }
+ return is_inv;
+ }
+
+ template <typename Particle>
+ double MinStepLength(Particle& p, setup::Trajectory&) const {
+ const Code pid = p.GetPID();
+ if (isEmParticle(pid) || isInvisible(pid)) {
+ cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
+ return 0.;
+ } else {
+ double next_step = std::numeric_limits<double>::infinity();
+ cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
+ return next_step;
+ }
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+ // cout << "ProcessCut: DoContinous: " << p.GetPID() << endl;
+ // cout << " is em: " << isEmParticle( p.GetPID() ) << endl;
+ // cout << " is inv: " << isInvisible( p.GetPID() ) << endl;
+ // const Code pid = p.GetPID();
+ // if( isEmParticle( pid ) ){
+ // cout << "removing em. particle..." << endl;
+ // fEmEnergy += p.GetEnergy();
+ // fEmCount += 1;
+ // p.Delete();
+ // return EProcessReturn::eParticleAbsorbed;
+ // }
+ // if ( isInvisible( pid ) ){
+ // cout << "removing inv. particle..." << endl;
+ // fInvEnergy += p.GetEnergy();
+ // fInvCount += 1;
+ // p.Delete();
+ // return EProcessReturn::eParticleAbsorbed;
+ // }
+ 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;
+ fInvEnergy = 0. * 1_GeV;
+ fInvCount = 0;
+ fEnergy = 0. * 1_GeV;
+ // defineEmParticles();
+ }
+
+ void ShowResults() {
+ cout << " ******************************" << endl
+ << " ParticleCut: " << endl
+ << " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
+ << " no. of em. particles injected: " << fEmCount << endl
+ << " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
+ << " no. of inv. particles injected: " << fInvCount << endl
+ << " ******************************" << endl;
+ }
+
+ EnergyType GetInvEnergy() { return fInvEnergy; }
+
+ EnergyType GetCutEnergy() { return fEnergy; }
+
+ EnergyType GetEmEnergy() { return fEmEnergy; }
+
+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 = 1;
+ int kBeam = process::sibyll::GetSibyllXSCode(corsikaBeamId);
+
+ bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
/*
the target should be defined by the Environment,
@@ -63,33 +256,55 @@ public:
*/
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
- int kTarget = 1;
- double beamEnergy = p.GetEnergy() / 1_GeV;
+ 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: en: " << beamEnergy << " pid:" << kBeam << std::endl;
- double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
- double dumdif[3];
+ << "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;
- if (kTarget == 1)
- sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
- else
- sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy);
+ double int_length = 0;
+ if (kInteraction) {
- 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;
+ 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();
- 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
- double 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;
- // add exponential sampling
/*
what are the units of the output? slant depth or 3space length?
@@ -110,9 +325,38 @@ public:
template <typename Particle, typename Stack>
void DoInteraction(Particle& p, Stack& s) const {
- cout << "DoInteraction: " << p.GetPID() << " interaction? "
+ 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
/*
@@ -121,18 +365,30 @@ public:
(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 Ecm = sqrt(2. * E * 0.93827_GeV);
+ 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);
- int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ std::cout << "ProcessSplit: "
+ << " DoDiscrete: gamma:" << gamma << endl;
+ std::cout << "ProcessSplit: "
+ << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
- /*
- the target should be defined by the Environment,
- ideally as full particle object so that the four momenta
- and the boosts can be defined..
- */
- // FOR NOW: set target to proton
- int kTarget = 1; // p.GetPID();
+ int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
std::cout << "ProcessSplit: "
<< " DoInteraction: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
@@ -148,7 +404,8 @@ public:
// running sibyll, filling stack
sibyll_(kBeam, kTarget, sqs);
// running decays
- // decsib_();
+ setTrackedParticlesStable();
+ decsib_();
// print final state
int print_unit = 6;
sib_list_(print_unit);
@@ -159,49 +416,62 @@ public:
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
- /*
- 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
- in general transformation is rotation + boost
- */
- const EnergyType proton_mass = 0.93827_GeV;
- const double gamma = (E + proton_mass) / (Ecm);
- const double gambet = sqrt(E * E - proton_mass * proton_mass) / Ecm;
// SibStack does not know about momentum yet so we need counter to access momentum
// array in Sibyll
int i = -1;
- for (auto& p : ss) {
+ super_stupid::MomentumVector Ptot_final(
+ rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+ for (auto& psib : ss) {
++i;
- // transform to lab. frame, primitve
- const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();
+ // 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 * 1_GeV);
- pnew.SetPID(process::sibyll::ConvertFromSibyll(p.GetPID()));
+ 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;
}
- } else
- p.Delete();
+ }
}
void Init() {
fCount = 0;
- // define reference frame? --> defines boosts between corsika stack and model stack.
-
- // initialize random numbers for sibyll
- // FOR NOW USE SIBYLL INTERNAL !!!
- // rnd_ini_();
-
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("s_rndm");
- // // corsika::random::RNG srng;
- // auto srng = rmng.GetRandomStream("s_rndm");
-
// test random number generator
std::cout << "ProcessSplit: "
<< " test sequence of random numbers." << std::endl;
@@ -211,29 +481,11 @@ public:
// initialize Sibyll
sibyll_ini_();
- // set particles stable / unstable
- // use stack to loop over particles
- setup::Stack ds;
- ds.NewParticle().SetPID(Code::Proton);
- ds.NewParticle().SetPID(Code::Neutron);
- ds.NewParticle().SetPID(Code::PiPlus);
- ds.NewParticle().SetPID(Code::PiMinus);
- ds.NewParticle().SetPID(Code::KPlus);
- ds.NewParticle().SetPID(Code::KMinus);
- ds.NewParticle().SetPID(Code::K0Long);
- ds.NewParticle().SetPID(Code::K0Short);
-
- for (auto& p : ds) {
- int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- // set particle stable by setting table value negative
- cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")"
- << " stable in Sibyll .." << endl;
- s_csydec_.idb[s_id] = -s_csydec_.idb[s_id - 1];
- p.Delete();
- }
+ setTrackedParticlesStable();
}
int GetCount() { return fCount; }
+ EnergyType GetEnergy() { return fEnergy; }
private:
};
@@ -245,9 +497,6 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
-class A {};
-class B : public A {};
-
int main() {
corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
@@ -266,20 +515,37 @@ int main() {
universe.AddChild(std::move(theMedium));
+ CoordinateSystem& rootCS = RootCoordinateSystem::GetInstance().GetRootCS();
+
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
+
ProcessSplit p1;
- const auto sequence = p0 + p1;
+ corsika::process::sibyll::ProcessDecay p2;
+ ProcessEMCut p3;
+ const auto sequence = /*p0 +*/ p1 + p2 + p3;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 100_TeV;
+ EnergyType E0 = 100_GeV;
+ MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV) / si::constants::c;
+ auto plab = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c,
+ 0. * 1_GeV / si::constants::c, P0);
particle.SetEnergy(E0);
+ particle.SetMomentum(plab);
particle.SetPID(Code::Proton);
+ particle.SetTime(0_ns);
+ Point p(rootCS, 0_m, 0_m, 0_m);
+ particle.SetPosition(p);
EAS.Init();
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
+ 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();
+ cout << "total energy (GeV): "
+ << (p3.GetCutEnergy() + p3.GetInvEnergy() + p3.GetEmEnergy()) / 1_GeV << endl;
}
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 7e9bd3bac3bc743bea8fcb0c2b0d1d401e16efae..a25316553da2182f6da9b2668f7f67336ede9a84 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -7,9 +7,12 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
using namespace std;
using namespace corsika::stack;
+using namespace corsika::units;
+using namespace corsika::geometry;
class SibStackData {
@@ -19,13 +22,30 @@ public:
void Clear() { s_plist_.np = 0; }
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 double v) { s_plist_.p[3][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 * si::constants::c;
+ }
int GetId(const int i) const { return s_plist_.llist[i]; }
- double GetEnergy(const int i) const { return s_plist_.p[3][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 / si::constants::c,
+ s_plist_.p[1][i] * 1_GeV / si::constants::c,
+ s_plist_.p[2][i] * 1_GeV / si::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];
@@ -45,13 +65,19 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
- double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+ 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;
diff --git a/Framework/Geometry/BaseTrajectory.h b/Framework/Geometry/BaseTrajectory.h
index feedace208c89eefc0ecd02cd066e0dc5eadf1a4..5540fc89d62e240217b79346d86aa32b04004baf 100644
--- a/Framework/Geometry/BaseTrajectory.h
+++ b/Framework/Geometry/BaseTrajectory.h
@@ -46,7 +46,7 @@ namespace corsika::geometry {
corsika::units::si::LengthType) const = 0;
virtual LengthType ArcLength(corsika::units::si::TimeType t1,
- corsika::units::si::TimeType t2) const = 0;
+ corsika::units::si::TimeType t2) const = 0;
virtual corsika::units::si::TimeType GetDuration(
corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index 0fa07e64277f534a815e7e067c08e88e0901d783..2c8f6b8792cb204e652d1222b38f7849b70e0792 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -58,8 +58,8 @@ namespace corsika::geometry {
auto GetRadius() const { return radius; }
- corsika::units::si::LengthType ArcLength(
- corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
+ corsika::units::si::LengthType ArcLength(corsika::units::si::TimeType t1,
+ corsika::units::si::TimeType t2) const {
return (vPar + vPerp).norm() * (t2 - t1);
}
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
index f3cabbac7886225846fd8cae1ca7079237812749..5bfcea0216ee6fd496275003b307cccc3f804ee8 100644
--- a/Framework/Geometry/Trajectory.h
+++ b/Framework/Geometry/Trajectory.h
@@ -12,8 +12,8 @@
#ifndef _include_TRAJECTORY_H
#define _include_TRAJECTORY_H
-#include <corsika/units/PhysicalUnits.h>
#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
using corsika::units::si::LengthType;
using corsika::units::si::TimeType;
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index f2e4a9fee1aabbbb5196a597db36de06e6385594..a0985a1ad2304e7e5bd4f1024cf43b5cefaf8bc4 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -20,6 +20,7 @@ set (
set (
MODEL_HEADERS
ParticleConversion.h
+ ProcessDecay.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
new file mode 100644
index 0000000000000000000000000000000000000000..701094e733ad463f673c8d73cdee0bdf34bdddab
--- /dev/null
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -0,0 +1,156 @@
+#ifndef _include_ProcessDecay_h_
+#define _include_ProcessDecay_h_
+
+#include <corsika/process/ContinuousProcess.h>
+
+#include <corsika/cascade/SibStack.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+// using namespace corsika::particles;
+
+namespace corsika::process {
+
+ namespace sibyll {
+
+ void setHadronsUnstable() {
+ // name? also makes EM particles stable
+
+ // loop over all particles in sibyll
+ // should be changed to loop over human readable list
+ // i.e. corsika::particles::ListOfParticles()
+ std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
+ // make ALL particles unstable, then set EM stable
+ for (auto& p : corsika2sibyll) {
+ // std::cout << (int)p << std::endl;
+ const int sibCode = (int)p;
+ // skip unknown and antiparticles
+ if (sibCode < 1) continue;
+ // std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll(
+ // static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
+ s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
+ // std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] <<
+ // std::endl;
+ }
+ // set Leptons and Proton and Neutron stable
+ // use stack to loop over particles
+ setup::Stack ds;
+ ds.NewParticle().SetPID(corsika::particles::Code::Proton);
+ ds.NewParticle().SetPID(corsika::particles::Code::Neutron);
+ ds.NewParticle().SetPID(corsika::particles::Code::Electron);
+ ds.NewParticle().SetPID(corsika::particles::Code::Positron);
+ ds.NewParticle().SetPID(corsika::particles::Code::NuE);
+ ds.NewParticle().SetPID(corsika::particles::Code::NuEBar);
+ ds.NewParticle().SetPID(corsika::particles::Code::MuMinus);
+ ds.NewParticle().SetPID(corsika::particles::Code::MuPlus);
+ ds.NewParticle().SetPID(corsika::particles::Code::NuMu);
+ ds.NewParticle().SetPID(corsika::particles::Code::NuMuBar);
+
+ for (auto& p : ds) {
+ int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ // set particle stable by setting table value negative
+ // cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
+ // << " stable in Sibyll .." << endl;
+ s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+ p.Delete();
+ }
+ }
+
+ class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+ public:
+ ProcessDecay() {}
+ void Init() {
+ // setHadronsUnstable();
+ }
+
+ void setAllStable() {
+ // name? also makes EM particles stable
+
+ // loop over all particles in sibyll
+ // should be changed to loop over human readable list
+ // i.e. corsika::particles::ListOfParticles()
+ for (auto& p : corsika2sibyll) {
+ // std::cout << (int)p << std::endl;
+ const int sibCode = (int)p;
+ // skip unknown and antiparticles
+ if (sibCode < 1) continue;
+ std::cout << "Sibyll: Decay: setting "
+ << ConvertFromSibyll(static_cast<SibyllCode>(sibCode)) << " stable"
+ << std::endl;
+ s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
+ std::cout << "decay table value: " << s_csydec_.idb[sibCode - 1] << std::endl;
+ }
+ }
+
+ friend void setHadronsUnstable();
+
+ template <typename Particle>
+ double MinStepLength(Particle& p, setup::Trajectory&) 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 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;
+ // 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;
+ }
+
+ template <typename Particle, typename Stack>
+ void DoDiscrete(Particle& p, Stack& s) const {
+ SibStack ss;
+ ss.Clear();
+ // copy particle to sibyll stack
+ auto pin = ss.NewParticle();
+ pin.SetPID(process::sibyll::ConvertToSibyllRaw(p.GetPID()));
+ pin.SetEnergy(p.GetEnergy());
+ pin.SetMomentum(p.GetMomentum());
+ // remove original particle from corsika stack
+ p.Delete();
+ // set all particles/hadrons unstable
+ setHadronsUnstable();
+ // call sibyll decay
+ std::cout << "ProcessDecay: calling Sibyll decay routine.." << std::endl;
+ decsib_();
+ // 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;
+ // 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()));
+ pnew.SetMomentum(psib.GetMomentum());
+ }
+ // 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
+
+#endif
diff --git a/Processes/Sibyll/code_generator.py b/Processes/Sibyll/code_generator.py
index 77c1f5dd8216ad48ebf5c4aa476ef57695bb7fea..a593a1f2ce5a5c3a281c758e2404c72347286ec0 100755
--- a/Processes/Sibyll/code_generator.py
+++ b/Processes/Sibyll/code_generator.py
@@ -159,8 +159,6 @@ if __name__ == "__main__":
pythia_db = load_pythiadb(sys.argv[1])
read_sibyll_codes(sys.argv[2], pythia_db)
- print (str(pythia_db))
-
with open("Generated.inc", "w") as f:
print("// this file is automatically generated\n// edit at your own risk!\n", file=f)
print(generate_sibyll_enum(pythia_db), file=f)
diff --git a/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index 266005cde54295a790a8748e91614bf6f9d58b37..ad90442848976c4cc779977a2c2f1890cedb1091 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -10,9 +10,9 @@
*/
#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/units/PhysicalUnits.h>
-#include <corsika/particles/ParticleProperties.h>
#include <corsika/logging/Logger.h>
@@ -28,7 +28,8 @@ using namespace corsika::process::stack_inspector;
template <typename Stack>
StackInspector<Stack>::StackInspector(const bool aReport)
- : fReport(aReport) {}
+ : fReport(aReport)
+ , fCountStep(0) {}
template <typename Stack>
StackInspector<Stack>::~StackInspector() {}
@@ -36,7 +37,6 @@ StackInspector<Stack>::~StackInspector() {}
template <typename Stack>
process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
Stack& s) const {
- static int countStep = 0;
if (!fReport) return EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
EnergyType Etot = 0_GeV;
@@ -51,8 +51,8 @@ process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Tr
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
<< " pos=" << pos << endl;
}
- countStep++;
- cout << "StackInspector: nStep=" << countStep << " stackSize=" << s.GetSize()
+ fCountStep++;
+ cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << s.GetSize()
<< " Estack=" << Etot / 1_GeV << " GeV" << endl;
return EProcessReturn::eOk;
}
@@ -63,7 +63,9 @@ double StackInspector<Stack>::MaxStepLength(Particle&, setup::Trajectory&) const
}
template <typename Stack>
-void StackInspector<Stack>::Init() {}
+void StackInspector<Stack>::Init() {
+ fCountStep = 0;
+}
#include <corsika/setup/SetupStack.h>
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 58b9ec5ebd7652b127f719beafcedacbbf79d503..6032beac618b0941b01ffd76962c3bc13ad2300b 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -40,6 +40,7 @@ namespace corsika::process {
private:
bool fReport;
+ mutable int fCountStep = 0;
};
} // namespace stack_inspector
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index ef44f2390a5071f645e383bff23373d8ec240993..a3df326480c9417b736f6ac058dd70f6e8cce253 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -13,8 +13,8 @@
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Vector.h>
#include <corsika/setup/SetupTrajectory.h>
using corsika::setup::Trajectory;
@@ -33,66 +33,80 @@ using namespace std;
using namespace corsika::units::si;
struct DummyParticle {
- EnergyType fEnergy;
- Vector<momentum_d> fMomentum;
- Point fPosition;
-
- DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition) : fEnergy(pEnergy), fMomentum(pMomentum), fPosition(pPosition) {}
-
- auto GetEnergy() const { return fEnergy; }
- auto GetMomentum() const { return fMomentum; }
- auto GetPosition() const { return fPosition; }
+ EnergyType fEnergy;
+ Vector<momentum_d> fMomentum;
+ Point fPosition;
+
+ DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition)
+ : fEnergy(pEnergy)
+ , fMomentum(pMomentum)
+ , fPosition(pPosition) {}
+
+ auto GetEnergy() const { return fEnergy; }
+ auto GetMomentum() const { return fMomentum; }
+ auto GetPosition() const { return fPosition; }
};
struct DummyStack {
- using ParticleType = DummyParticle;
+ using ParticleType = DummyParticle;
};
TEST_CASE("TrackingLine") {
corsika::environment::Environment env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
-
+
tracking_line::TrackingLine<DummyStack> tracking(env);
SECTION("intersection with sphere") {
Point const origin(cs, {0_m, 0_m, 0_m});
Point const center(cs, {0_m, 0_m, 10_m});
Sphere const sphere(center, 1_m);
- Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
+ Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+ 0_m / second, 1_m / second);
Line line(origin, v);
-
+
geometry::Trajectory<Line> traj(line, 12345_s);
-
- auto const opt = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
+
+ auto const opt =
+ tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
REQUIRE(opt.has_value());
-
+
auto [t1, t2] = opt.value();
REQUIRE(t1 / 9_s == Approx(1));
REQUIRE(t2 / 11_s == Approx(1));
-
- auto const optNoIntersection = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
+
+ auto const optNoIntersection =
+ tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
REQUIRE_FALSE(optNoIntersection.has_value());
}
-
+
SECTION("maximally possible propagation") {
- auto& universe = *(env.GetUniverse());
-
- //~ std::cout << env.GetUniverse().get() << std::endl;
-
- DummyParticle p(1_J, Vector<momentum_d>(cs, 0*kilogram*meter/second, 0*kilogram*meter/second, 1*kilogram*meter/second), Point(cs, 0_m, 0_m,0_m));
-
- auto const radius = 20_m;
-
- auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
- universe.AddChild(std::move(theMedium));
-
- Point const origin(cs, {0_m, 0_m, 0_m});
- Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
- Line line(origin, v);
-
- auto const traj = tracking.GetTrack(p);
-
- REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius)).GetComponents(cs).norm().magnitude() == Approx(0).margin(1e-4));
+ auto& universe = *(env.GetUniverse());
+
+ //~ std::cout << env.GetUniverse().get() << std::endl;
+
+ DummyParticle p(
+ 1_J,
+ Vector<momentum_d>(cs, 0 * kilogram * meter / second,
+ 0 * kilogram * meter / second, 1 * kilogram * meter / second),
+ Point(cs, 0_m, 0_m, 0_m));
+
+ auto const radius = 20_m;
+
+ auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
+ universe.AddChild(std::move(theMedium));
+
+ Point const origin(cs, {0_m, 0_m, 0_m});
+ Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
+ 0_m / second, 1_m / second);
+ Line line(origin, v);
+
+ auto const traj = tracking.GetTrack(p);
+
+ REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
+ .GetComponents(cs)
+ .norm()
+ .magnitude() == Approx(0).margin(1e-4));
}
}
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index e822b64ec0dd4b393fd7285360c96879ce2ccc0f..c5de9d53c9e26d6f5afe88d9762ec88df3043a64 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -15,7 +15,7 @@
#include <corsika/environment/IMediumModel.h>
namespace corsika::setup {
- using IEnvironmentModel = corsika::environment::IMediumModel;
+ using IEnvironmentModel = corsika::environment::IMediumModel;
}
#endif
diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h
index 8c79934422420791c1d33395810e5b90afa4038e..a889d32a652753955ba2553012350e521b1d99dc 100644
--- a/Stack/SuperStupidStack/SuperStupidStack.h
+++ b/Stack/SuperStupidStack/SuperStupidStack.h
@@ -137,7 +137,7 @@ namespace corsika::stack {
void IncrementSize() {
fDataPID.push_back(Code::Unknown);
fDataE.push_back(0 * joule);
-#warning this here makes no sense: see issue #48
+ //#TODO this here makes no sense: see issue #48
geometry::CoordinateSystem& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCS();
fMomentum.push_back(MomentumVector(