IAP GITLAB

Commit b5f4c7a0 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus 🖖
Browse files

workshop example

parent 45dc98c8
Pipeline #887 passed with stages
in 5 minutes and 58 seconds
......@@ -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
......
/*
* (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;
}
......@@ -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
......
Supports Markdown
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