IAP GITLAB

Skip to content
Snippets Groups Projects
Commit f5268d41 authored by ralfulrich's avatar ralfulrich
Browse files

make cascade_example a bit more exciting -> 1TeV, and added more printout

parent 76cf8e50
No related branches found
No related tags found
1 merge request!39Resolve "cascade example debugging"
Hide whitespace changes
Inline Side-by-side
Documentation/Examples/cascade_example.cc
+ 44
− 39
View file @ f5268d41
......@@ -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;
......
Framework/Cascade/Cascade.h
+ 3
− 3
View file @ f5268d41
......@@ -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
......
Processes/Sibyll/Decay.h
+ 10
− 12
View file @ f5268d41
......@@ -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
......
Processes/Sibyll/Interaction.h
+ 8
− 8
View file @ f5268d41
......@@ -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
......
Processes/TrackingLine/TrackingLine.h
+ 10
− 0
View file @ f5268d41
......@@ -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);
}
};
......
Stack/SuperStupidStack/testSuperStupidStack.cc
+ 3
− 4
View file @ f5268d41
......@@ -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);
}
......
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