From b5f4c7a0811d9c0ce57022162cf95336e10d4a31 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Thu, 20 Jun 2019 02:30:17 +0200
Subject: [PATCH] workshop example
---
Documentation/Examples/CMakeLists.txt | 20 +++
Documentation/Examples/workshop_example.cc | 185 +++++++++++++++++++++
Tools/plot_tracks.sh | 9 +-
3 files changed, 213 insertions(+), 1 deletion(-)
create mode 100644 Documentation/Examples/workshop_example.cc
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 97acd456..2824dfc7 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -37,6 +37,26 @@ target_link_libraries (cascade_example
)
install (TARGETS cascade_example DESTINATION share/examples)
+add_executable(workshop_example workshop_example.cc)
+target_link_libraries(workshop_example
+SuperStupidStack
+ CORSIKAunits
+ CORSIKAlogging
+ CORSIKArandom
+ ProcessSibyll
+ CORSIKAcascade
+ ProcessTrackWriter
+ ProcessParticleCut
+ ProcessTrackingLine
+ CORSIKAprocesses
+ CORSIKAparticles
+ CORSIKAgeometry
+ CORSIKAenvironment
+ CORSIKAprocesssequence
+ )
+
+install (TARGETS workshop_example DESTINATION share/examples)
+
CORSIKA_ADD_TEST (boundary_example)
target_link_libraries (boundary_example
SuperStupidStack
diff --git a/Documentation/Examples/workshop_example.cc b/Documentation/Examples/workshop_example.cc
new file mode 100644
index 00000000..34baba8e
--- /dev/null
+++ b/Documentation/Examples/workshop_example.cc
@@ -0,0 +1,185 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/cascade/Cascade.h>
+#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/FlatExponential.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
+
+#include <corsika/geometry/Sphere.h>
+
+#include <corsika/process/sibyll/Decay.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+
+#include <corsika/process/track_writer/TrackWriter.h>
+
+#include <corsika/process/particle_cut/ParticleCut.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/utl/CorsikaFenv.h>
+
+#include <iostream>
+#include <limits>
+#include <typeinfo>
+
+using namespace corsika;
+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;
+
+using namespace std;
+using namespace corsika::units::si;
+
+// example boundary crossing process
+template <bool deleteParticle>
+struct MyBoundaryCrossingProcess
+ : public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
+
+ MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
+
+ template <typename Particle>
+ EProcessReturn DoBoundaryCrossing(Particle& p,
+ typename Particle::BaseNodeType const& from,
+ typename Particle::BaseNodeType const& to) {
+ std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
+
+ auto const& name = particles::GetName(p.GetPID());
+ auto const pos = p.GetPosition().GetCoordinates();
+
+ fFile << name << " " << pos[0] / 1_m << ' ' << pos[1] / 1_m << ' ' << pos[2] / 1_m
+ << '\n';
+
+ if constexpr (deleteParticle) { p.Delete(); }
+
+ return EProcessReturn::eOk;
+ }
+
+ void Init() {}
+
+private:
+ std::ofstream fFile;
+};
+
+//
+// The example main program for a particle cascade
+//
+int main() {
+ feenableexcept(FE_INVALID);
+
+ // initialize random number sequence(s)
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+
+ using EnvType = Environment<setup::IEnvironmentModel>;
+
+ // inheritance from IMediumModel to implement the pure virtual methods
+ using FlatExp = environment::FlatExponential<setup::IEnvironmentModel>;
+ using Homogeneous = environment::HomogeneousMedium<setup::IEnvironmentModel>;
+
+ EnvType env;
+ auto& universe = *(env.GetUniverse());
+
+ // obtain the root coordinate system provided by environment object
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ auto outerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_km}, 500_km);
+
+ outerMedium->SetModelProperties<FlatExp>(
+ Point{rootCS, 0_m, 0_m, 5_km}, Vector<dimensionless_d>{rootCS, {0., 0., 1.}},
+ 1_kg / cube(1_m), 8_km,
+ environment::NuclearComposition{
+ {particles::Code::Nitrogen, particles::Code::Oxygen},
+ {0.7847f, 1.f - 0.7847f}});
+
+ auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
+
+ auto const props = innerMedium->SetModelProperties<Homogeneous>(
+ 1_g / cube(1_m), environment::NuclearComposition{{particles::Code::Proton}, {1.f}});
+
+ universe.SetModelProperties(props);
+
+ outerMedium->AddChild(std::move(innerMedium));
+
+ universe.AddChild(std::move(outerMedium));
+
+ // setup processes, decays and interactions
+ tracking_line::TrackingLine tracking;
+
+ process::sibyll::Interaction sibyll;
+ process::sibyll::Decay decay;
+
+ process::particle_cut::ParticleCut cut(55_GeV);
+
+ process::track_writer::TrackWriter trackWriter("tracks.dat");
+ MyBoundaryCrossingProcess<false> boundaryCrossing("crossings.dat");
+
+ // assemble all processes into an ordered process list
+ auto sequence = sibyll << decay << cut << boundaryCrossing << trackWriter;
+
+ // setup particle stack, and add primary particles
+ setup::Stack stack;
+ stack.Clear();
+ const Code beamCode = Code::Proton;
+ const HEPMassType mass = particles::GetMass(Code::Proton);
+ const HEPEnergyType E0 = 500_TeV;
+
+ for (int N = 2, i = 0; i < N; ++i) {
+ auto const phi = 0;
+ auto const theta = i * 360 / N;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ auto const [px, py, pz] =
+ momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << " phi=" << phi << endl;
+ cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
+ Point pos(rootCS, 0_m, 0_m, 0_m);
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ beamCode, E0, plab, pos, 0_ns});
+ }
+
+ // define air shower object, run simulation
+ cascade::Cascade EAS(env, tracking, sequence, stack);
+ EAS.Init();
+ EAS.Run();
+
+ cout << "Result: E0=" << E0 / 1_GeV << endl;
+ cut.ShowResults();
+ const HEPEnergyType Efinal =
+ cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
+ cout << "total energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+}
diff --git a/Tools/plot_tracks.sh b/Tools/plot_tracks.sh
index 135b410d..888fe311 100755
--- a/Tools/plot_tracks.sh
+++ b/Tools/plot_tracks.sh
@@ -28,10 +28,17 @@ set output "$output"
set xlabel "x / m"
set ylabel "y / m"
set zlabel "z / m"
+
+d = 200
+
+set xrange [-d:d]
+set yrange [-d:d]
+set zrange [-15e3:15e3]
+
set title "CORSIKA 8 preliminary"
do for [t=0:360:1] {
- set view 60, t
+ set view 85, t
splot "$track_dat" u 3:4:5:6:7:8 w vectors nohead t ""
}
EOF
--
GitLab