diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 60e436a4e7574f233d642bd6673a635a6db25823..c5889f56b9b8e1c58424e36ec3f5c44b7079e14e 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -55,14 +55,17 @@ static int fEmCount;
static EnergyType fInvEnergy;
static int fInvCount;
-class ProcessEMCut : public corsika::process::ContinuousProcess<ProcessEMCut> {
+class ProcessCut : public corsika::process::ContinuousProcess<ProcessCut> {
+ EnergyType fECut;
+
public:
- ProcessEMCut() {}
+ ProcessCut(const EnergyType v)
+ : fECut(v) {}
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)
+ if (p.GetEnergy() < fECut || Ecm < 10_GeV)
return true;
else
return false;
@@ -134,24 +137,25 @@ public:
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) const {
- cout << "ProcessCut: DoContinuous: " << p.GetPID() << endl;
const Code pid = p.GetPID();
+ EnergyType energy = p.GetEnergy();
+ cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy << endl;
EProcessReturn ret = EProcessReturn::eOk;
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
- fEmEnergy += p.GetEnergy();
+ fEmEnergy += energy;
fEmCount += 1;
p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
- fInvEnergy += p.GetEnergy();
+ fInvEnergy += energy;
fInvCount += 1;
p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
- fEnergy += p.GetEnergy();
+ fEnergy += energy;
p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
}
@@ -178,20 +182,21 @@ public:
<< " ******************************" << endl;
}
- EnergyType GetInvEnergy() { return fInvEnergy; }
-
- EnergyType GetCutEnergy() { return fEnergy; }
-
- EnergyType GetEmEnergy() { return fEmEnergy; }
-
-private:
+ EnergyType GetInvEnergy() const { return fInvEnergy; }
+ EnergyType GetCutEnergy() const { return fEnergy; }
+ EnergyType GetEmEnergy() const { return fEmEnergy; }
};
+//
+// The example main program for a particle cascade
+//
int main() {
+ // initialize random number sequence(s)
corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
- corsika::environment::Environment env; // dummy environment
+ // setup environment, geometry
+ corsika::environment::Environment env;
auto& universe = *(env.GetUniverse());
auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
@@ -201,47 +206,47 @@ int main() {
using MyHomogeneousModel =
corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
- 1_g / (1_m * 1_m * 1_m),
+ 1_kg / (1_m * 1_m * 1_m),
corsika::environment::NuclearComposition(
- std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
+ std::vector<corsika::particles::Code>{corsika::particles::Code::Oxygen},
std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
- CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ // setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
-
corsika::process::sibyll::Interaction sibyll;
corsika::process::sibyll::Decay decay;
- ProcessEMCut cut;
+ ProcessCut cut(8_GeV);
+
+ // assemble all processes into an ordered process list
const auto sequence = p0 << sibyll << decay << cut;
+ // setup particle stack, and add primary particle
setup::Stack stack;
+ stack.Clear();
+ const EnergyType E0 = 1_TeV;
+ {
+ auto particle = stack.NewParticle();
+ particle.SetPID(Code::Proton);
+ hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
+ auto plab = stack::super_stupid::MomentumVector(rootCS, 0. * 1_GeV, 0. * 1_GeV, -P0);
+ particle.SetEnergy(E0);
+ particle.SetMomentum(plab);
+ particle.SetTime(0_ns);
+ Point p(rootCS, 0_m, 0_m, 10_km);
+ particle.SetPosition(p);
+ }
+ // define air shower object, run simulation
corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
-
- stack.Clear();
- auto particle = stack.NewParticle();
- EnergyType E0 = 100_GeV;
- hep::MomentumType P0 = sqrt(E0 * E0 - 0.93827_GeV * 0.93827_GeV);
- auto plab = stack::super_stupid::MomentumVector(rootCS, 0. * 1_GeV, 0. * 1_GeV, 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, particles below energy threshold =" << p1.GetCount()
- << endl;
- cout << "total energy below threshold (GeV): " //<< p1.GetEnergy() / 1_GeV
- << std::endl;
+
+ cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
cout << "total energy (GeV): "
<< (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy()) / 1_GeV << endl;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 31c16e2eab8164db2598c9149409cd5dc7900f76..6693f764c5a7ddacce02ca7a20f720bcf1a0ddaf 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -71,7 +71,7 @@ namespace corsika::cascade {
fProcessSequence.GetTotalInverseInteractionLength(particle, step);
// sample random exponential step length in grammage
- std::exponential_distribution expDist((1_m * 1_m / 1_g) / total_inv_lambda);
+ std::exponential_distribution expDist(total_inv_lambda * (1_g / (1_m * 1_m)));
GrammageType const next_interact = (1_g / (1_m * 1_m)) * expDist(fRNG);
std::cout << "total_inv_lambda=" << total_inv_lambda
@@ -99,7 +99,6 @@ namespace corsika::cascade {
<< ", next_decay=" << next_decay << std::endl;
// convert next_decay from time to length [m]
- // Environment::GetDistance(step, next_decay);
LengthType const distance_decay = next_decay * particle.GetMomentum().norm() /
particle.GetEnergy() *
corsika::units::constants::c;
@@ -108,10 +107,11 @@ namespace corsika::cascade {
auto const min_distance =
std::min({distance_interact, distance_decay, distance_max});
+ std::cout << " move particle by : " << min_distance << std::endl;
+
// here the particle is actually moved along the trajectory to new position:
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.PositionFromArclength(min_distance));
-
// .... also update time, momentum, direction, ...
// apply all continuous processes on particle + track
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 24fd807d2723343be6a827f0783fa8fd5fdd0383..4d3eca70f0df8eadde09469ad70ea326cdf7ece7 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -15,8 +15,11 @@ namespace corsika::process {
namespace sibyll {
class Decay : public corsika::process::DecayProcess<Decay> {
+ mutable int fCount = 0;
+
public:
Decay() {}
+ ~Decay() { std::cout << "Sibyll::Decay n=" << fCount << std::endl; }
void Init() {
setHadronsUnstable();
setTrackedParticlesStable();
@@ -128,19 +131,13 @@ namespace corsika::process {
const corsika::units::si::TimeType t0 =
corsika::particles::GetLifetime(p.GetPID());
- cout << "Decay: code: " << p.GetPID() << endl;
- cout << "Decay: MinStep: t0: " << t0 << endl;
- cout << "Decay: MinStep: energy: " << E / 1_GeV << 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 << "Decay: MinStep: x0: " << x0 << endl;
+ cout << "Decay: GetLifetime: \n"
+ << " code: " << p.GetPID() << endl;
+ cout << " t0: " << t0 << endl;
+ cout << " energy: " << E / 1_GeV << endl;
+ cout << " gamma: " << gamma << endl;
corsika::units::si::TimeType const lifetime = gamma * t0;
- cout << "Decay: MinStep: tau: " << lifetime << endl;
- // int a = 1;
- // const double x = -x0 * log(s_rndm_(a));
- // cout << "Decay: next decay: " << x << endl;
+ cout << " -> tau: " << lifetime << endl;
return lifetime;
}
@@ -148,6 +145,7 @@ namespace corsika::process {
void DoDecay(Particle& p, Stack& s) const {
using corsika::geometry::Point;
using namespace corsika::units::si;
+ fCount++;
SibStack ss;
ss.Clear();
// copy particle to sibyll stack
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fbaf7c5f519391af653da149d7c9aa3d173d6b0e..6897113ca6bf7b39fd2dcb724ff3bc99ee080a17 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -15,9 +15,11 @@ namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
+ mutable int fCount = 0;
+
public:
Interaction() {}
- ~Interaction() {}
+ ~Interaction() { std::cout << "Sibyll::Interaction n=" << fCount << std::endl; }
void Init() {
@@ -73,8 +75,8 @@ namespace corsika::process::sibyll {
const hep::EnergyType sqs = sqrt(Etot * Etot - Ptot.squaredNorm());
const double Ecm = sqs / 1_GeV;
- std::cout << "Interaction: "
- << "MinStep: input en: " << p.GetEnergy() / 1_GeV << endl
+ std::cout << "Interaction: LambdaInt: \n"
+ << " input energy: " << p.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kBeam << endl
<< " beam XS code:" << kBeam << endl
<< " beam pid:" << p.GetPID() << endl
@@ -137,9 +139,7 @@ namespace corsika::process::sibyll {
<< "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 =
+ const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point pOrig = p.GetPosition();
@@ -157,7 +157,6 @@ namespace corsika::process::sibyll {
const int kTarget = corsika::particles::Oxygen::
GetNucleusA(); // 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 = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
@@ -207,8 +206,9 @@ namespace corsika::process::sibyll {
<< " DoInteraction: should have dropped particle.." << std::endl;
// p.Delete(); delete later... different process
} else {
+ fCount++;
// Sibyll does not know about units..
- double sqs = Ecm / 1_GeV;
+ const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, kTarget, sqs);
// running decays
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index f92ee8eee67b1b242b5e248ba07b815e4ee0bbb3..6fdd4f6624edc1c522352e954ed70a9c52b4b8d4 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -12,6 +12,7 @@
#define _include_corsika_processes_TrackingLine_h_
#include <corsika/geometry/Point.h>
+#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
@@ -75,12 +76,19 @@ namespace corsika::process {
void Init() {}
auto GetTrack(Particle const& p) {
+ using std::cout;
+ using std::endl;
using namespace corsika::units::si;
+ using namespace corsika::geometry;
geometry::Vector<SpeedType::dimension_type> const velocity =
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
auto const currentPosition = p.GetPosition();
+ std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
+ << std::endl;
+ std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+
// to do: include effect of magnetic field
geometry::Line line(currentPosition, velocity);
@@ -125,6 +133,8 @@ namespace corsika::process {
min = *minIter;
}
+ std::cout << " t-intersect: " << min << std::endl;
+
return geometry::Trajectory<corsika::geometry::Line>(line, min);
}
};
diff --git a/Stack/SuperStupidStack/testSuperStupidStack.cc b/Stack/SuperStupidStack/testSuperStupidStack.cc
index 092baf5f9259ad658392a78b9a753c37ad74d3e4..6047fc8513ddaeef769ad6fc401b63f6be9339d0 100644
--- a/Stack/SuperStupidStack/testSuperStupidStack.cc
+++ b/Stack/SuperStupidStack/testSuperStupidStack.cc
@@ -45,10 +45,9 @@ TEST_CASE("SuperStupidStack", "[stack]") {
auto pout = s.GetNextParticle();
REQUIRE(pout.GetPID() == particles::Code::Electron);
REQUIRE(pout.GetEnergy() == 1.5_GeV);
-#warning Fix the next two lines:
- // REQUIRE(pout.GetMomentum() == MomentumVector(dummyCS, {1 * joule, 1 * joule, 1 *
- // joule})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 * meter, 1 *
- // meter}));
+ // REQUIRE(pout.GetMomentum() == stack::super_stupid::MomentumVector(dummyCS, {1_GeV,
+ // 1_GeV, 1_GeV})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 * meter,
+ // 1 * meter}));
REQUIRE(pout.GetTime() == 100_s);
}