From f551ac31d3a77cd3ca03b6037d57d2e9d3b49c69 Mon Sep 17 00:00:00 2001
From: Alan Coleman <alanco@umich.edu>
Date: Thu, 7 Dec 2023 11:25:29 +0000
Subject: [PATCH] Resolve "synchronize output format between examples"
---
CMakeLists.txt | 7 +
README.md | 30 +-
applications/CMakeLists.txt | 13 +
applications/README.md | 7 +
.../c8_air_shower.cpp | 68 ++---
.../process/SwitchProcessSequence.inl | 16 +-
.../modules/energy_loss/BetheBlochPDG.inl | 4 +-
.../modules/proposal/ContinuousProcess.inl | 4 +-
.../modules/proposal/InteractionModel.inl | 3 +-
.../modules/proposal/ProposalProcessBase.inl | 2 +-
corsika/detail/modules/radio/RadioProcess.inl | 5 +-
corsika/modules/Random.hpp | 2 +-
corsika/modules/writers/EnergyLossWriter.hpp | 2 +-
.../modules/writers/LongitudinalWriter.hpp | 4 +-
corsika/modules/writers/WriterOff.hpp | 2 +-
examples/CMakeLists.txt | 128 ++++----
examples/README.md | 9 +
examples/cascade_example.cpp | 158 ----------
examples/{ => cascade_examples}/em_shower.cpp | 96 +++---
examples/{ => cascade_examples}/mars.cpp | 68 +++--
.../mc_conex.cpp} | 150 ++++-----
.../radio_em_shower.cpp | 199 ++++++------
examples/{ => cascade_examples}/water.cpp | 142 +++++----
examples/cascade_proton_example.cpp | 147 ---------
.../boundary_crossing.cpp} | 10 +-
.../{ => framework_examples}/environment.cpp | 0
.../geometry.cpp} | 0
.../helix_trajectory.cpp} | 0
.../known_particles.cpp} | 11 +-
.../stack.cpp} | 0
.../static_sequence.cpp} | 2 +-
.../{ => physics_examples}/stopping_power.cpp | 9 +-
.../synchrotron_clover_leaf.cpp} | 22 +-
.../synchrotron_test_C8tracking.cpp | 0
.../synchrotron_test_manual_tracking.cpp | 0
examples/vertical_EAS.cpp | 286 ------------------
36 files changed, 566 insertions(+), 1040 deletions(-)
create mode 100644 applications/CMakeLists.txt
create mode 100644 applications/README.md
rename examples/corsika.cpp => applications/c8_air_shower.cpp (92%)
create mode 100644 examples/README.md
delete mode 100644 examples/cascade_example.cpp
rename examples/{ => cascade_examples}/em_shower.cpp (74%)
rename examples/{ => cascade_examples}/mars.cpp (88%)
rename examples/{hybrid_MC.cpp => cascade_examples/mc_conex.cpp} (81%)
rename examples/{ => cascade_examples}/radio_em_shower.cpp (58%)
rename examples/{ => cascade_examples}/water.cpp (72%)
delete mode 100644 examples/cascade_proton_example.cpp
rename examples/{boundary_example.cpp => framework_examples/boundary_crossing.cpp} (94%)
rename examples/{ => framework_examples}/environment.cpp (100%)
rename examples/{geometry_example.cpp => framework_examples/geometry.cpp} (100%)
rename examples/{helix_example.cpp => framework_examples/helix_trajectory.cpp} (100%)
rename examples/{particle_list_example.cpp => framework_examples/known_particles.cpp} (85%)
rename examples/{stack_example.cpp => framework_examples/stack.cpp} (100%)
rename examples/{staticsequence_example.cpp => framework_examples/static_sequence.cpp} (99%)
rename examples/{ => physics_examples}/stopping_power.cpp (92%)
rename examples/{clover_leaf.cpp => physics_examples/synchrotron_clover_leaf.cpp} (90%)
rename examples/{ => physics_examples}/synchrotron_test_C8tracking.cpp (100%)
rename examples/{ => physics_examples}/synchrotron_test_manual_tracking.cpp (100%)
delete mode 100644 examples/vertical_EAS.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 29ebd5356..65b91e2a7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -262,6 +262,13 @@ add_subdirectory (modules/conex)
add_subdirectory (modules/epos)
add_subdirectory (modules/fluka)
+#
+#+++++++++++++++++++++++++++++++
+# standard applications
+#
+
+add_subdirectory(applications)
+
#+++++++++++++++++++++++++++++++
# unit testing
#
diff --git a/README.md b/README.md
index 506ac4d93..105f090a6 100644
--- a/README.md
+++ b/README.md
@@ -102,8 +102,6 @@ make install
### FLUKA support
-Warning: may only work when the next version of FLUKA is released (as of 2023.06.15)
-
For legal reasons we do not distribute/bundle FLUKA together with CORSIKA 8.
If you want to use FLUKA as low-energy hadronic interaction model, you have to download
it separately from (http://www.fluka.org/), which requires registering there as FLUKA user.
@@ -142,22 +140,40 @@ make install
To run the Unit Tests, just type `ctest` in your build area.
-## Running examples
+## Running applications and examples
+
+### Standard applications
+
+Applications for standard use-cases are located in the `applications` directory.
+These are example scripts that can be used directly or slightly modified for your use case.
+See [applications/README.md] for more.
+The applications are compiled automatically after running `make` and will appear your `corsika-build/bin` directory.
+After running `make install` the binaries will also be copied into your `corsika-install/bin` directory as well.
+
+
+For example, from inside your `corsika-install/bin` directory, run
+```shell
+c8_air_shower --pdg 2212 -E 1e5 -f my_shower
+```
+This will run a vertical 100 TeV proton shower and will create and put the output into `./my_shower`.
+
### Building the examples
-From your top corsika build directory, (the one that includes `corsika-build` and `corsika-install`) type
+Unlike the applications, the examples must be compiled as a second step.
+From your top corsika directory, (the one that includes `corsika-build` and `corsika-install`) run
```shell
cmake -Dcorsika_DIR=$PWD/corsika-build -S ./corsika/examples -B ./corsika-build-examples
cd corsika-build-examples
make -j4 #The number should match the number of available cores on your machine
```
-From any directory, run the program `corsika-build-examples/bin/corsika`. As an example, you can run with the following flags:
+You can run the examples by using the binaries in `corsika-build-examples/bin/`.
+For example:
```shell
-corsika-build-examples/bin/corsika --pdg 2212 -E 1e5 -f my_shower
+corsika-build-examples/bin/known_particles
```
-This will run a vertical 100 TeV proton shower and will create and put the output into `./my_shower`.
+This will print out all of the particles that are known by CORSIKA.
### Generating doxygen documentation
diff --git a/applications/CMakeLists.txt b/applications/CMakeLists.txt
new file mode 100644
index 000000000..27dac14b5
--- /dev/null
+++ b/applications/CMakeLists.txt
@@ -0,0 +1,13 @@
+SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
+
+if(WITH_FLUKA)
+ message("compiling c8_air_shower.cpp with FLUKA")
+else()
+ message("compiling c8_air_shower.cpp with UrQMD")
+endif()
+add_executable (c8_air_shower c8_air_shower.cpp)
+target_link_libraries (c8_air_shower CORSIKA8)
+
+install (
+ TARGETS c8_air_shower DESTINATION bin
+)
diff --git a/applications/README.md b/applications/README.md
new file mode 100644
index 000000000..24b892df2
--- /dev/null
+++ b/applications/README.md
@@ -0,0 +1,7 @@
+# CORSIKA 8 Applications
+
+This directory contains standard applications which are typical for astroparticle physics solutions.
+They are "physics-complete" and are suitable for generating simulations that can be used in publications.
+
+For example, `c8_air_shower.cpp` would be a similar binary to what would be built by CORSIKA 7 and will simulate
+air showers in a curved atmosphere.
diff --git a/examples/corsika.cpp b/applications/c8_air_shower.cpp
similarity index 92%
rename from examples/corsika.cpp
rename to applications/c8_air_shower.cpp
index df5547107..fd7b02538 100644
--- a/examples/corsika.cpp
+++ b/applications/c8_air_shower.cpp
@@ -20,7 +20,6 @@
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/process/DynamicInteractionProcess.hpp>
-#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/random/RNGManager.hpp>
@@ -35,7 +34,6 @@
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/Environment.hpp>
-#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/GeomagneticModel.hpp>
#include <corsika/media/GladstoneDaleRefractiveIndex.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -50,7 +48,6 @@
#include <corsika/modules/Epos.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
-#include <corsika/modules/OnShellCheck.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
@@ -94,8 +91,17 @@ using namespace std;
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
using EnvType = Environment<EnvironmentInterface>;
-
-using Particle = setup::Stack<EnvType>::particle_type;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = setup::Tracking;
+using Particle = StackType::particle_type;
+
+//
+// This is the main example script which runs EAS with fairly standard settings
+// w.r.t. what was implemented in CORSIKA 7. Users may want to change some of the
+// specifics (observation altitude, magnetic field, energy cuts, etc.), but this
+// example is the most physics-complete one and should be used for full simulations
+// of particle cascades in air
+//
long registerRandomStreams(long seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
@@ -111,9 +117,9 @@ long registerRandomStreams(long seed) {
if (seed == 0) {
std::random_device rd;
seed = rd();
- std::cout << "random seed (auto) " << seed << std::endl;
+ CORSIKA_LOG_INFO("random seed (auto) {}", seed);
} else {
- std::cout << "random seed " << seed << std::endl;
+ CORSIKA_LOG_INFO("random seed {}", seed);
}
RNGManager<>::getInstance().setSeed(seed);
return seed;
@@ -346,6 +352,8 @@ int main(int argc, char** argv) {
// we make the axis much longer than the inj-core distance since the
// profile will go beyond the core, depending on zenith angle
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
+ auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
+ uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
/* === END: CONSTRUCT GEOMETRY === */
double const emthinfrac = app["--emthin"]->as<double>();
@@ -364,22 +372,22 @@ int main(int argc, char** argv) {
OutputManager output(app["--filename"]->as<std::string>(), seed, args.str(), outputDir);
// register energy losses as output
- EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
+ EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
output.add("energyloss", dEdX);
- DynamicInteractionProcess<setup::Stack<EnvType>> heModel;
+ DynamicInteractionProcess<StackType> heModel;
// have SIBYLL always for PROPOSAL photo-hadronic interactions
auto sibyll = std::make_shared<corsika::sibyll::Interaction>(env);
if (auto const modelStr = app["--hadronModel"]->as<std::string>();
modelStr == "SIBYLL-2.3d") {
- heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{sibyll};
+ heModel = DynamicInteractionProcess<StackType>{sibyll};
} else if (modelStr == "QGSJet-II.04") {
- heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{
+ heModel = DynamicInteractionProcess<StackType>{
std::make_shared<corsika::qgsjetII::Interaction>()};
} else if (modelStr == "EPOS-LHC") {
- heModel = DynamicInteractionProcess<setup::Stack<EnvType>>{
+ heModel = DynamicInteractionProcess<StackType>{
std::make_shared<corsika::epos::Interaction>()};
} else {
CORSIKA_LOG_CRITICAL("invalid choice \"{}\"; also check argument parser", modelStr);
@@ -426,7 +434,7 @@ int main(int argc, char** argv) {
auto emContinuous =
make_select(EMHadronSwitch(), emContinuousBethe, emContinuousProposal);
- LongitudinalWriter profile{showerAxis, 200, 10_g / square(1_cm)};
+ LongitudinalWriter profile{showerAxis, nAxisBins, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
@@ -438,7 +446,7 @@ int main(int argc, char** argv) {
#endif
InteractionCounter leIntCounted{leIntModel};
- StackInspector<setup::Stack<EnvType>> stackInspect(10000, false, E0);
+ StackInspector<StackType> stackInspect(10000, false, E0);
// assemble all processes into an ordered process list
struct EnergySwitch {
@@ -452,30 +460,20 @@ int main(int argc, char** argv) {
// observation plane
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+ ObservationPlane<TrackingType, ParticleWriterParquet> observationLevel{
obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
true, // plane should "absorb" particles
false}; // do not print z-coordinate
// register ground particle output
output.add("particles", observationLevel);
- PrimaryWriter<setup::Tracking, ParticleWriterParquet> primaryWriter(observationLevel);
+ PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
output.add("primary", primaryWriter);
int ring_number{app["--ring"]->as<int>()};
- std::cout << "Ring number : " << ring_number << std::endl;
auto const radius_{ring_number * 25_m};
- // std::cout << "Radius = " << radius_ << std::endl;
const int rr_ = static_cast<int>(radius_ / 1_m);
- // if (ring_number == 0) {
- // // assemble the final process sequence without radio
- // auto sequence = make_sequence(stackInspect, hadronSequence,
- // decaySequence, emCascade,
- // emContinuous, longprof, observationLevel,
- // thinning, cut);
- // } else {
-
// Radio antennas and relevant information
// the antenna time variables
const TimeType duration_{4e-7_s};
@@ -491,14 +489,11 @@ int main(int argc, char** argv) {
auto const injectionPosY_{injectionPos.getCoordinates().getY()};
auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
- std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
if (ring_number != 0) {
// setup CoREAS antennas - use the for loop for star shape pattern
for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
auto phiRad_1 = phi_1 / 180. * M_PI;
- // auto phi_1 = 0;
- // auto phiRad_1 = phi_1 / 180. * M_PI;
auto const point_1{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_1),
showerCoreY_ + radius_ * sin(phiRad_1),
constants::EarthRadius::Mean)};
@@ -514,8 +509,6 @@ int main(int argc, char** argv) {
// setup ZHS antennas - use the for loop for star shape pattern
for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
auto phiRad_ = phi_ / 180. * M_PI;
- // auto phi_ = 0; phi_;
- // auto phiRad_ = phi_ / 180. * M_PI;
auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
showerCoreY_ + radius_ * sin(phiRad_),
constants::EarthRadius::Mean)};
@@ -557,18 +550,19 @@ int main(int argc, char** argv) {
// create the cascade object using the default stack and tracking
// implementation
- setup::Tracking tracking;
- setup::Stack<EnvType> stack;
+ TrackingType tracking;
+ StackType stack;
Cascade EAS(env, tracking, sequence, output, stack);
// print our primary parameters all in one place
if (app["--pdg"]->count() > 0) {
- CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as<int>());
+ CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as<int>());
} else {
- CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
+ CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
}
- CORSIKA_LOG_INFO("Primary Energy: {}", E0);
- CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
+ CORSIKA_LOG_INFO("Primary Energy: {}", E0);
+ CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
+ CORSIKA_LOG_INFO("Primary Direction: {}", plab.getNorm());
CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates());
CORSIKA_LOG_INFO("Shower Axis Length: {}", (showerCore - injectionPos).getNorm() * 1.2);
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index 0e4584a76..364e307c4 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -77,10 +77,10 @@ namespace corsika {
template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TParticle>
- inline ProcessReturn SwitchProcessSequence<
- TCondition, TSequence, USequence, IndexStart, IndexProcess1,
- IndexProcess2>::doContinuous(Step<TParticle>& step,
- ContinuousProcessIndex const idLimit) {
+ inline ProcessReturn SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
+ IndexProcess1, IndexProcess2>::
+ doContinuous(Step<TParticle>& step,
+ [[maybe_unused]] ContinuousProcessIndex const idLimit) {
if (select_(step.getParticlePre())) {
if constexpr (process1_type::is_process_sequence) {
return A_.doContinuous(step, idLimit);
@@ -218,10 +218,10 @@ namespace corsika {
template <typename TCondition, typename TSequence, typename USequence, int IndexStart,
int IndexProcess1, int IndexProcess2>
template <typename TParticle>
- CrossSectionType SwitchProcessSequence<
- TCondition, TSequence, USequence, IndexStart, IndexProcess1,
- IndexProcess2>::getCrossSection(TParticle const& projectile, Code const targetId,
- FourMomentum const& targetP4) const {
+ CrossSectionType SwitchProcessSequence<TCondition, TSequence, USequence, IndexStart,
+ IndexProcess1, IndexProcess2>::
+ getCrossSection(TParticle const& projectile, [[maybe_unused]] Code const targetId,
+ [[maybe_unused]] FourMomentum const& targetP4) const {
if (select_(projectile)) {
if constexpr (is_interaction_process_v<process1_type>) {
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 6847d90a3..fed757cd4 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -162,10 +162,8 @@ namespace corsika {
HEPEnergyType const dE = getTotalEnergyLoss(step.getParticlePre(), dX);
// if (dE > HEPEnergyType::zero())
// dE = -dE;
- [[maybe_unused]] const auto Ekin = step.getEkinPre();
- auto EkinNew = Ekin + dE;
CORSIKA_LOG_TRACE("EnergyLoss dE={} MeV, Ekin={} GeV, EkinNew={} GeV", dE / 1_MeV,
- Ekin / 1_GeV, EkinNew / 1_GeV);
+ step.getEkinPre() / 1_GeV, (step.getEkinPre() + dE) / 1_GeV);
step.add_dEkin(dE);
// also send to output
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 36f277270..a49a846ba 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -48,8 +48,8 @@ namespace corsika::proposal {
template <typename TEnvironment, typename... TOutputArgs>
inline ContinuousProcess<TOutput>::ContinuousProcess(TEnvironment const& _env,
TOutputArgs&&... args)
- : TOutput(args...)
- , ProposalProcessBase(_env) {}
+ : ProposalProcessBase(_env)
+ , TOutput(args...) {}
template <typename TOutput>
template <typename TParticle>
diff --git a/corsika/detail/modules/proposal/InteractionModel.inl b/corsika/detail/modules/proposal/InteractionModel.inl
index 05c44ba2f..288b36c51 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -56,7 +56,8 @@ namespace corsika::proposal {
template <typename TStackView>
inline ProcessReturn
InteractionModel<THadronicLEModel, THadronicHEModel>::doInteraction(
- TStackView& view, Code const projectileId, FourMomentum const& projectileP4) {
+ TStackView& view, Code const projectileId,
+ [[maybe_unused]] FourMomentum const& projectileP4) {
auto const projectile = view.getProjectile();
diff --git a/corsika/detail/modules/proposal/ProposalProcessBase.inl b/corsika/detail/modules/proposal/ProposalProcessBase.inl
index 0ea18a4bf..2a72811ee 100644
--- a/corsika/detail/modules/proposal/ProposalProcessBase.inl
+++ b/corsika/detail/modules/proposal/ProposalProcessBase.inl
@@ -46,7 +46,7 @@ namespace corsika::proposal {
}
return lowest_table_value;
- };
+ }
template <typename TEnvironment>
inline ProposalProcessBase::ProposalProcessBase(TEnvironment const& _env) {
diff --git a/corsika/detail/modules/radio/RadioProcess.inl b/corsika/detail/modules/radio/RadioProcess.inl
index 9211af14b..bad68ac85 100644
--- a/corsika/detail/modules/radio/RadioProcess.inl
+++ b/corsika/detail/modules/radio/RadioProcess.inl
@@ -57,7 +57,8 @@ namespace corsika {
template <typename Particle, typename Track>
inline LengthType
RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::getMaxStepLength(
- const Particle& vParticle, const Track& vTrack) const {
+ [[maybe_unused]] const Particle& vParticle,
+ [[maybe_unused]] const Track& vTrack) const {
return meter * std::numeric_limits<double>::infinity();
}
@@ -176,4 +177,4 @@ namespace corsika {
return config;
}
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/modules/Random.hpp b/corsika/modules/Random.hpp
index a254aa229..d31bd35df 100644
--- a/corsika/modules/Random.hpp
+++ b/corsika/modules/Random.hpp
@@ -23,7 +23,7 @@ namespace corsika {
inline void rng_func(corsika::default_prng_type& rng, double* dest, std::size_t N) {
std::uniform_real_distribution<double> udist(0.0, 1.0);
std::generate(dest, std::next(dest, N), std::bind(udist, std::ref(rng)));
- };
+ }
} // namespace detail
inline void connect_random_stream(corsika::default_prng_type& rng,
diff --git a/corsika/modules/writers/EnergyLossWriter.hpp b/corsika/modules/writers/EnergyLossWriter.hpp
index b516ab626..c247b424e 100644
--- a/corsika/modules/writers/EnergyLossWriter.hpp
+++ b/corsika/modules/writers/EnergyLossWriter.hpp
@@ -156,7 +156,7 @@ namespace corsika {
/**
* Return the configuration of this output.
*/
- YAML::Node getConfig() const;
+ YAML::Node getConfig() const override;
private:
ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage
diff --git a/corsika/modules/writers/LongitudinalWriter.hpp b/corsika/modules/writers/LongitudinalWriter.hpp
index caa49f9af..0d6eb2822 100644
--- a/corsika/modules/writers/LongitudinalWriter.hpp
+++ b/corsika/modules/writers/LongitudinalWriter.hpp
@@ -122,12 +122,12 @@ namespace corsika {
/**
* Return a summary.
*/
- YAML::Node getSummary() const;
+ YAML::Node getSummary() const override;
/**
* Return the configuration of this output.
*/
- YAML::Node getConfig() const;
+ YAML::Node getConfig() const override;
number_profile::ProfileData const& getProfile(
number_profile::ProfileIndex index) const {
diff --git a/corsika/modules/writers/WriterOff.hpp b/corsika/modules/writers/WriterOff.hpp
index 3d752cc9a..da418d54d 100644
--- a/corsika/modules/writers/WriterOff.hpp
+++ b/corsika/modules/writers/WriterOff.hpp
@@ -43,7 +43,7 @@ namespace corsika {
template <typename... TArgs>
void write(TArgs&&...) {}
- virtual YAML::Node getConfig() const { return YAML::Node(); }
+ virtual YAML::Node getConfig() const override { return YAML::Node(); }
}; // class WriterOff
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index d3aa7ecc2..614bfb9bc 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -9,88 +9,88 @@ find_package (corsika CONFIG REQUIRED)
# this is only for CORSIKA_REGISTER_EXAMPLE
include ("${CMAKE_CURRENT_SOURCE_DIR}/corsikaExamples.cmake")
-add_executable (helix_example helix_example.cpp)
-target_link_libraries (helix_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (helix_example)
-add_executable (geometry_example geometry_example.cpp)
-target_link_libraries (geometry_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (geometry_example)
+###################
+## cascade_examples
+###################
-add_executable (stack_example stack_example.cpp)
-target_link_libraries (stack_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (stack_example)
+add_executable (em_shower cascade_examples/em_shower.cpp)
+target_link_libraries (em_shower CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS 100 8472)
-add_executable (cascade_example cascade_example.cpp)
-target_link_libraries (cascade_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (cascade_example)
+if (WITH_FLUKA)
+ add_executable (mars cascade_examples/mars.cpp)
+ target_link_libraries (mars CORSIKA8::CORSIKA8)
+ message("FLUKA found, will make 'mars' example")
+else()
+ message("FLUKA not found, will not make 'mars' example")
+endif()
-add_executable (boundary_example boundary_example.cpp)
-target_link_libraries (boundary_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (boundary_example)
+add_executable (mc_conex cascade_examples/mc_conex.cpp)
+target_link_libraries (mc_conex CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (mc_conex RUN_OPTIONS 4 2 10000.)
-add_executable (cascade_proton_example cascade_proton_example.cpp)
-target_link_libraries (cascade_proton_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (cascade_proton_example)
+add_executable (radio_em_shower cascade_examples/radio_em_shower.cpp)
+target_link_libraries (radio_em_shower CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (radio_em_shower RUN_OPTIONS 10 1121673489)
-add_executable (vertical_EAS vertical_EAS.cpp)
-target_link_libraries (vertical_EAS CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
+if (WITH_FLUKA)
+ add_executable (water cascade_examples/water.cpp)
+ target_link_libraries (water CORSIKA8::CORSIKA8)
+ message("FLUKA found, will make 'water' example")
+else()
+ message("FLUKA not found, will not make 'water' example")
+endif()
-add_executable (stopping_power stopping_power.cpp)
-target_link_libraries (stopping_power CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (stopping_power)
-add_executable (staticsequence_example staticsequence_example.cpp)
-target_link_libraries (staticsequence_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (staticsequence_example)
+#####################
+## framework_examples
+#####################
-add_executable (particle_list_example particle_list_example.cpp)
-target_link_libraries (particle_list_example CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (particle_list_example)
+add_executable (boundary_crossing framework_examples/boundary_crossing.cpp)
+target_link_libraries (boundary_crossing CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (boundary_crossing)
-add_executable (em_shower em_shower.cpp)
-target_link_libraries (em_shower CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS 100 8472)
+add_executable (environment framework_examples/environment.cpp)
+target_link_libraries (environment CORSIKA8::CORSIKA8)
-add_executable (hybrid_MC hybrid_MC.cpp)
-target_link_libraries (hybrid_MC CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
+add_executable (geometry framework_examples/geometry.cpp)
+target_link_libraries (geometry CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (geometry)
+add_executable (helix_trajectory framework_examples/helix_trajectory.cpp)
+target_link_libraries (helix_trajectory CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (helix_trajectory)
-add_executable (corsika corsika.cpp)
-# for ICRC 2023
-option(WITH_FLUKA "use FLUKA instead of UrQMD in corsika.cpp")
-if (WITH_FLUKA)
- target_compile_definitions(corsika PRIVATE WITH_FLUKA)
- message("compiling corsika.cpp with FLUKA")
-else()
- message("compiling corsika.cpp with UrQMD")
-endif()
-target_link_libraries (corsika CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (corsika RUN_OPTIONS --energy=10 --zenith=10 --nevent=2 --filename=corsika_test --pdg=2212)
+add_executable (known_particles framework_examples/known_particles.cpp)
+target_link_libraries (known_particles CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (known_particles)
-add_executable (mars mars.cpp)
-target_link_libraries (mars CORSIKA8::CORSIKA8)
+add_executable (stack framework_examples/stack.cpp)
+target_link_libraries (stack CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (stack)
-add_executable (water water.cpp)
-target_link_libraries (water CORSIKA8::CORSIKA8)
+add_executable (static_sequence framework_examples/static_sequence.cpp)
+target_link_libraries (static_sequence CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (static_sequence)
-add_executable (environment environment.cpp)
-target_link_libraries (environment CORSIKA8::CORSIKA8)
-add_executable (synchrotron_test_manual_tracking synchrotron_test_manual_tracking.cpp)
-target_link_libraries (synchrotron_test_manual_tracking CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (synchrotron_test_manual_tracking)
+###################
+## physics_examples
+###################
-add_executable (synchrotron_test_C8tracking synchrotron_test_C8tracking.cpp)
+add_executable (stopping_power physics_examples/stopping_power.cpp)
+target_link_libraries (stopping_power CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (stopping_power)
+
+add_executable (synchrotron_clover_leaf physics_examples/synchrotron_clover_leaf.cpp)
+target_link_libraries (synchrotron_clover_leaf CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (synchrotron_clover_leaf)
+
+add_executable (synchrotron_test_C8tracking physics_examples/synchrotron_test_C8tracking.cpp)
target_link_libraries (synchrotron_test_C8tracking CORSIKA8::CORSIKA8)
CORSIKA_REGISTER_EXAMPLE (synchrotron_test_C8tracking)
-add_executable (clover_leaf clover_leaf.cpp)
-target_link_libraries (clover_leaf CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (clover_leaf)
-
-add_executable (radio_em_shower radio_em_shower.cpp)
-target_link_libraries (radio_em_shower CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (radio_em_shower RUN_OPTIONS 10 1121673489)
+add_executable (synchrotron_test_manual_tracking physics_examples/synchrotron_test_manual_tracking.cpp)
+target_link_libraries (synchrotron_test_manual_tracking CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (synchrotron_test_manual_tracking)
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 000000000..ba9a661db
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,9 @@
+# CORSIKA 8 Examples
+
+This directory contains several example scripts that can help you get started to simulate your own showers and/or to begin developing your own interfaces within the C8 framework.
+The sub-directories contain the several types of examples:
+
+ * **cascade_examples**: scripts to set up running full cascade simulations for various use cases. This is a good starting point if you want to simulate showers for your own experiments.
+ * **framework_examples**: show how to use and interact with the internal framework of C8
+ * **physics_examples**: demonstrate various tests of the particle interactions and/or cascades
+
\ No newline at end of file
diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp
deleted file mode 100644
index 351809bf7..000000000
--- a/examples/cascade_example.cpp
+++ /dev/null
@@ -1,158 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/process/SwitchProcessSequence.hpp>
-#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/core/EnergyMomentumOperations.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-
-#include <corsika/output/OutputManager.hpp>
-#include <corsika/modules/writers/SubWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriter.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/HomogeneousMedium.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/media/UniformMagneticField.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <corsika/modules/BetheBlochPDG.hpp>
-#include <corsika/modules/StackInspector.hpp>
-#include <corsika/modules/Sibyll.hpp>
-#include <corsika/modules/ParticleCut.hpp>
-#include <corsika/modules/TrackWriter.hpp>
-
-#include <iostream>
-#include <limits>
-
-using namespace corsika;
-using namespace std;
-
-//
-// The example main program for a particle cascade
-//
-int main() {
-
- logging::set_level(logging::level::warn);
-
- CORSIKA_LOG_INFO("cascade_example");
-
- LengthType const height_atmosphere = 112.8_km;
-
- feenableexcept(FE_INVALID);
- // initialize random number sequence(s)
- RNGManager<>::getInstance().registerRandomStream("cascade");
-
- // setup environment, geometry
- using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
- using EnvType = Environment<EnvironmentInterface>;
- EnvType env;
- auto& universe = *(env.getUniverse());
-
- CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
-
- auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
-
- using MyHomogeneousModel =
- MediumPropertyModel<UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>;
-
- // fraction of oxygen
- double const fox = 0.20946;
- auto const props = world->setModelProperties<MyHomogeneousModel>(
- Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 0_T),
- 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition({Code::Nitrogen, Code::Oxygen}, {1. - fox, fox}));
-
- auto innerMedium = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
-
- innerMedium->setModelProperties(props);
- world->addChild(std::move(innerMedium));
- universe.addChild(std::move(world));
-
- // setup particle stack, and add primary particle
- setup::Stack<EnvType> stack;
- stack.clear();
- const int nuclA = 4;
- const int nuclZ = int(nuclA / 2.15 + 0.7);
- const Code beamCode = get_nucleus_code(nuclA, nuclZ);
- const HEPMassType mass = get_nucleus_mass(beamCode);
- const HEPEnergyType E0 = nuclA * 1_TeV;
- double theta = 0.;
- double phi = 0.;
-
- Point const injectionPos(
- rootCS, 0_m, 0_m,
- height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
-
- OutputManager output("cascade_outputs");
-
- theta *= M_PI / 180.;
- phi *= M_PI / 180.;
- DirectionVector const direction(
- rootCS, {sin(theta) * cos(phi), sin(theta) * sin(phi), -cos(theta)});
-
- ShowerAxis const showerAxis{injectionPos, direction * 100_km, env};
- EnergyLossWriter dEdX{showerAxis};
- output.add("energyloss", dEdX);
-
- {
- HEPMomentumType const P0 = calculate_momentum(E0, mass);
- auto plab = direction * P0;
- CORSIKA_LOG_INFO("input particle: {}", beamCode);
- CORSIKA_LOG_INFO("input angles: theta={} phi={}", theta, phi);
- CORSIKA_LOG_INFO("input momentum: {}", plab.getComponents() / 1_GeV);
- stack.addParticle(std::make_tuple(
- beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)), direction,
- injectionPos, 0_ns));
- }
-
- // setup processes, decays and interactions
- setup::Tracking tracking;
- StackInspector<setup::Stack<EnvType>> stackInspect(100, true, E0);
-
- RNGManager<>::getInstance().registerRandomStream("sibyll");
-
- corsika::sibyll::Interaction sibyll{env};
- corsika::sibyll::Decay decay;
-
- // cascade with only HE model ==> HE cut
- ParticleCut<SubWriter<decltype(dEdX)>> cut(80_GeV, true, dEdX);
- BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
-
- TrackWriter trackWriter;
- output.add("tracks", trackWriter); // register TrackWriter
-
- // assemble all processes into an ordered process list
- auto sequence = make_sequence(stackInspect, sibyll, decay, eLoss, trackWriter, cut);
-
- // define air shower object, run simulation
- Cascade EAS(env, tracking, sequence, output, stack);
-
- output.startOfLibrary();
- EAS.run();
- output.endOfLibrary();
-
- const HEPEnergyType Efinal = dEdX.getEnergyLost();
- CORSIKA_LOG_INFO(
- "\n"
- "total cut energy (GeV) : {}\n"
- "relative difference (%): {}\n"
- "total dEdX energy (GeV): {}\n"
- "relative difference (%): {}\n",
- Efinal / 1_GeV, (Efinal / E0 - 1) * 100, dEdX.getEnergyLost() / 1_GeV,
- dEdX.getEnergyLost() / E0 * 100);
-}
diff --git a/examples/em_shower.cpp b/examples/cascade_examples/em_shower.cpp
similarity index 74%
rename from examples/em_shower.cpp
rename to examples/cascade_examples/em_shower.cpp
index cd06faa12..a4d6c3471 100644
--- a/examples/em_shower.cpp
+++ b/examples/cascade_examples/em_shower.cpp
@@ -18,10 +18,10 @@
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
@@ -29,11 +29,9 @@
#include <corsika/media/Environment.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/UniformMagneticField.hpp>
-#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/PROPOSAL.hpp>
@@ -62,13 +60,19 @@ void registerRandomStreams(int seed) {
if (seed == 0) {
std::random_device rd;
seed = rd();
- cout << "new random seed (auto) " << seed << endl;
+ CORSIKA_LOG_INFO("random seed (auto) {} ", seed);
+ } else {
+ CORSIKA_LOG_INFO("random seed {} ", seed);
}
RNGManager<>::getInstance().setSeed(seed);
}
+using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using EnvType = Environment<EnvironmentInterface>;
template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = setup::Tracking;
int main(int argc, char** argv) {
@@ -87,29 +91,23 @@ int main(int argc, char** argv) {
registerRandomStreams(seed);
// setup environment, geometry
- using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
- using EnvType = Environment<EnvironmentInterface>;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
// build a Linsley US Standard atmosphere into `env`
+ MagneticFieldVector bField{rootCS, 50_uT, 0_T, 0_T};
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
- env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
- MagneticFieldVector{rootCS, 20.4_uT, 0_T, -43.23_uT});
+ env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, bField);
std::unordered_map<Code, HEPEnergyType> energy_resolution = {
- {Code::Electron, 2_MeV},
- {Code::Positron, 2_MeV},
- {Code::Photon, 2_MeV},
+ {Code::Electron, 5_MeV},
+ {Code::Positron, 5_MeV},
+ {Code::Photon, 5_MeV},
};
for (auto [pcode, energy] : energy_resolution)
set_energy_production_threshold(pcode, energy);
- // setup particle stack, and add primary particle
- setup::Stack<EnvType> stack;
- stack.clear();
-
const Code beamCode = Code::Electron;
auto const mass = get_mass(beamCode);
const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
@@ -123,10 +121,6 @@ int main(int argc, char** argv) {
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz});
- cout << "input particle: " << beamCode << endl;
- cout << "input angles: theta=" << theta << endl;
- cout << "input momentum: " << plab.getComponents() / 1_GeV
- << ", norm = " << plab.getNorm() << endl;
auto const observationHeight = 0.0_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
@@ -137,54 +131,67 @@ int main(int argc, char** argv) {
Point const injectionPos =
showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
- std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
-
- stack.addParticle(std::make_tuple(
- beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
- plab.normalized(), injectionPos, 0_ns));
-
- CORSIKA_LOG_INFO("shower axis length: {} ",
- (showerCore - injectionPos).getNorm() * 1.02);
-
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
-
- OutputManager output("em_shower_outputs");
-
- EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
- // register energy losses as output
- output.add("dEdX", dEdX);
+ auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
+ uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
+
+ CORSIKA_LOG_INFO("Primary particle: {}", beamCode);
+ CORSIKA_LOG_INFO("Zenith angle: {} (rad)", theta);
+ CORSIKA_LOG_INFO("Momentum: {} (GeV)", plab.getComponents() / 1_GeV);
+ CORSIKA_LOG_INFO("Propagation dir: {}", plab.getNorm());
+ CORSIKA_LOG_INFO("Injection point: {}", injectionPos.getCoordinates());
+ CORSIKA_LOG_INFO("shower axis length: {} ",
+ (showerCore - injectionPos).getNorm() * 1.02);
// setup processes, decays and interactions
+ EnergyLossWriter energyloss{showerAxis, dX, nAxisBins};
+ ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
+ energyloss);
- ParticleCut<SubWriter<decltype(dEdX)>> cut(2_MeV, 2_MeV, 100_GeV, 100_GeV, true, dEdX);
corsika::sibyll::Interaction sibyll{env};
corsika::sophia::InteractionModel sophia;
- HEPEnergyType heThresholdNN = 60_GeV;
+ HEPEnergyType heThresholdNN = 80_GeV;
corsika::proposal::Interaction emCascade(
env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
- corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
- // BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
+ corsika::proposal::ContinuousProcess<SubWriter<decltype(energyloss)>> emContinuous(
+ env, energyloss);
// NOT possible right now, due to interface differenc in PROPOSAL
// InteractionCounter emCascadeCounted(emCascade);
+ OutputManager output("em_shower_outputs");
+
+ output.add("energyloss", energyloss);
+
TrackWriter tracks;
output.add("tracks", tracks);
- // long. profile
- LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
+ LongitudinalWriter profile{showerAxis, nAxisBins, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+ ObservationPlane<TrackingType, ParticleWriterParquet> observationLevel{
obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
output.add("particles", observationLevel);
+ PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
+ output.add("primary", primaryWriter);
+
auto sequence = make_sequence(emCascade, emContinuous, longprof, observationLevel, cut);
// define air shower object, run simulation
- setup::Tracking tracking;
+ TrackingType tracking;
+
+ auto const primaryProperties = std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns);
+
+ // setup particle stack, and add primary particle
+ StackType stack;
+ stack.clear();
+ stack.addParticle(primaryProperties);
+ primaryWriter.recordPrimary(primaryProperties);
output.startOfLibrary();
Cascade EAS(env, tracking, sequence, output, stack);
@@ -194,7 +201,8 @@ int main(int argc, char** argv) {
EAS.run();
- HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+ HEPEnergyType const Efinal =
+ energyloss.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO(
"total energy budget (GeV): {}, "
@@ -202,4 +210,6 @@ int main(int argc, char** argv) {
Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
output.endOfLibrary();
+
+ return EXIT_SUCCESS;
}
diff --git a/examples/mars.cpp b/examples/cascade_examples/mars.cpp
similarity index 88%
rename from examples/mars.cpp
rename to examples/cascade_examples/mars.cpp
index 6a57b6ec3..55739516b 100644
--- a/examples/mars.cpp
+++ b/examples/cascade_examples/mars.cpp
@@ -30,6 +30,7 @@
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
@@ -55,7 +56,7 @@
#include <corsika/modules/Sophia.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
-#include <corsika/modules/UrQMD.hpp>
+#include <corsika/modules/FLUKA.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
@@ -74,8 +75,9 @@ using namespace std;
using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
-
-using Particle = setup::Stack<EnvType>::particle_type;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = setup::Tracking;
+using Particle = StackType::particle_type;
typedef decltype(1 * pascal) PressureType;
typedef decltype(1 * degree_celsius) TemperatureType;
@@ -111,7 +113,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("pythia");
- RNGManager<>::getInstance().registerRandomStream("urqmd");
+ RNGManager<>::getInstance().registerRandomStream("fluka");
RNGManager<>::getInstance().registerRandomStream("proposal");
if (seed == 0) {
std::random_device rd;
@@ -200,8 +202,6 @@ int main(int argc, char** argv) {
}
// check that we got either PDG or A/Z
- // this can be done with option_groups but the ordering
- // gets all messed up
if (app.count("--pdg") == 0) {
if ((app.count("-A") == 0) || (app.count("-Z") == 0)) {
CORSIKA_LOG_ERROR("If --pdg is not provided, then both -A and -Z are required.");
@@ -224,12 +224,10 @@ int main(int argc, char** argv) {
Medium::AirDry1Atm, // Mars, close enough
MagneticFieldVector{rootCS, 0_T, 0_uT, 0_T}); // Mars
- builder.setNuclearComposition( // Mars
- {{Code::Nitrogen, Code::Oxygen}, {1. / 3., 2. / 3.}}); // simplified
- //{{Code::Carbon, Code::Oxygen, // 95.97 CO2
- // Code::Nitrogen}, // 1.89 N2 + 1.93 Argon + 0.146 O2
- // {0.9597 / 3, 0.9597 * 2 / 3,
- // 1 - 0.9597}}); // values taken from AIRES manual, Ar removed for now
+ builder.setNuclearComposition( // Mars
+ {{Code::Carbon, Code::Oxygen, // 95.97 CO2
+ Code::Nitrogen}, // 1.89 N2 + 1.93 Argon + 0.146 O2
+ {0.9597 / 3, 0.9597 * 2 / 3, 1 - 0.9597}}); // values taken from AIRES manual
MarsAtmModel layer1(0.699e3 * pascal, 0.00009 / 1_m, 31.0 * degree_celsius,
0.000998 * 1 * degree_celsius / 1_m);
@@ -305,8 +303,10 @@ int main(int argc, char** argv) {
OutputManager output(app["--filename"]->as<std::string>());
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
+ auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
+ uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
- EnergyLossWriter dEdX{showerAxis};
+ EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
output.add("energyloss", dEdX);
HEPEnergyType const emcut = 1_GeV;
@@ -338,13 +338,13 @@ int main(int argc, char** argv) {
auto emContinuous =
make_select(EMHadronSwitch(), emContinuousBethe, emContinuousProposal);
- LongitudinalWriter longprof{showerAxis};
+ LongitudinalWriter longprof{showerAxis, nAxisBins, dX};
output.add("profile", longprof);
LongitudinalProfile<SubWriter<decltype(longprof)>> profile{longprof};
- corsika::urqmd::UrQMD urqmd;
- InteractionCounter urqmdCounted{urqmd};
- StackInspector<setup::Stack<EnvType>> stackInspect(5000, false, E0);
+ corsika::fluka::Interaction leIntModel{env};
+ InteractionCounter leIntCounted{leIntModel};
+ StackInspector<StackType> stackInspect(5000, false, E0);
// assemble all processes into an ordered process list
struct EnergySwitch {
@@ -354,7 +354,7 @@ int main(int argc, char** argv) {
bool operator()(Particle const& p) const { return (p.getKineticEnergy() < cutE_); }
};
auto hadronSequence =
- make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, sibyllCounted);
+ make_select(EnergySwitch(heHadronModelThreshold), leIntCounted, sibyllCounted);
// track writer
TrackWriter trackWriter;
@@ -362,11 +362,14 @@ int main(int argc, char** argv) {
// observation plane
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<setup::Tracking> observationLevel(
- obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
+ ObservationPlane<TrackingType> observationLevel(obsPlane,
+ DirectionVector(rootCS, {1., 0., 0.}));
// register the observation plane with the output
output.add("particles", observationLevel);
+ PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
+ output.add("primary", primaryWriter);
+
// assemble the final process sequence
auto sequence =
make_sequence(stackInspect, hadronSequence, decayPythia, emCascade, emContinuous,
@@ -375,18 +378,19 @@ int main(int argc, char** argv) {
// create the cascade object using the default stack and tracking
// implementation
- setup::Tracking tracking;
- setup::Stack<EnvType> stack;
+ TrackingType tracking;
+ StackType stack;
Cascade EAS(env, tracking, sequence, output, stack);
// print our primary parameters all in one place
if (app["--pdg"]->count() > 0) {
- CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as<int>());
+ CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as<int>());
} else {
- CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
+ CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
}
- CORSIKA_LOG_INFO("Primary Energy: {}", E0);
- CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
+ CORSIKA_LOG_INFO("Primary Energy: {}", E0);
+ CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
+ CORSIKA_LOG_INFO("Primary Direction: {}", plab.getNorm());
CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates());
CORSIKA_LOG_INFO("Shower Axis Length: {}", (showerCore - injectionPos).getNorm() * 1.2);
@@ -409,9 +413,12 @@ int main(int argc, char** argv) {
stack.clear();
// add the desired particle to the stack
- stack.addParticle(std::make_tuple(
+ auto const primaryProperties = std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
- plab.normalized(), injectionPos, 0_ns));
+ plab.normalized(), injectionPos, 0_ns);
+ stack.addParticle(primaryProperties);
+
+ primaryWriter.recordPrimary(primaryProperties);
// run the shower
EAS.run();
@@ -424,11 +431,6 @@ int main(int argc, char** argv) {
"relative difference (%): {}",
Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
- auto const hists = sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
-
- save_hist(hists.labHist(), labHist_file, true);
- save_hist(hists.CMSHist(), cMSHist_file, true);
-
// trigger the output manager to save this shower to disk
output.endOfShower();
}
diff --git a/examples/hybrid_MC.cpp b/examples/cascade_examples/mc_conex.cpp
similarity index 81%
rename from examples/hybrid_MC.cpp
rename to examples/cascade_examples/mc_conex.cpp
index cd819bdf1..cdee78eb4 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/cascade_examples/mc_conex.cpp
@@ -6,50 +6,43 @@
* the license.
*/
-/* clang-format off */
-// InteractionCounter used boost/histogram, which
-// fails if boost/type_traits have been included before. Thus, we have
-// to include it first...
-#include <corsika/framework/process/InteractionCounter.hpp>
-/* clang-format on */
-#include <corsika/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/process/SwitchProcessSequence.hpp>
-#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/Step.hpp>
-#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/process/InteractionCounter.hpp>
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+// #include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/output/OutputManager.hpp>
-#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/output/OutputManager.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/Environment.hpp>
-#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
-#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/CORSIKA7Atmospheres.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
-#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/UrQMD.hpp>
-#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/CONEX.hpp>
#include <corsika/setup/SetupStack.hpp>
@@ -63,6 +56,12 @@
using namespace corsika;
using namespace std;
+//
+// An example of running an EAS where the hadronic cascade is
+// handled by sibyll+URQMD and the EM cascade is treated with
+// CONEX + Bethe Bloch (as opposed to PROPOSAL).
+//
+
/**
* Random number stream initialization
*
@@ -70,17 +69,19 @@ using namespace std;
*/
void registerRandomStreams(uint64_t seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
- RNGManager<>::getInstance().registerRandomStream("qgsjet");
- RNGManager<>::getInstance().registerRandomStream("sibyll");
+ RNGManager<>::getInstance().registerRandomStream("conex");
RNGManager<>::getInstance().registerRandomStream("epos");
+ RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("pythia");
+ RNGManager<>::getInstance().registerRandomStream("qgsjet");
+ RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("urqmd");
- RNGManager<>::getInstance().registerRandomStream("proposal");
- RNGManager<>::getInstance().registerRandomStream("conex");
if (seed == 0) {
std::random_device rd;
seed = rd();
- CORSIKA_LOG_INFO("new random seed (auto) {}", seed);
+ CORSIKA_LOG_INFO("random seed (auto) {} ", seed);
+ } else {
+ CORSIKA_LOG_INFO("random seed {} ", seed);
}
RNGManager<>::getInstance().setSeed(seed);
}
@@ -141,8 +142,12 @@ private:
/**
* Selection of environment interface implementation:
*/
+using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using EnvType = Environment<EnvironmentInterface>;
template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = setup::Tracking;
int main(int argc, char** argv) {
@@ -157,7 +162,6 @@ int main(int argc, char** argv) {
" if no seed is given, a random seed is chosen");
return 1;
}
- feenableexcept(FE_INVALID);
uint64_t seed = 0;
if (argc > 4) seed = std::stol(std::string(argv[4]));
@@ -165,20 +169,15 @@ int main(int argc, char** argv) {
registerRandomStreams(seed);
// setup environment, geometry
- using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
- using EnvType = Environment<EnvironmentInterface>;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
// build a Linsley US Standard atmosphere into `env`
+ MagneticFieldVector bField{rootCS, 50_uT, 0_T, 0_T};
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
- env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
- MagneticFieldVector{rootCS, 0_T, 50_uT, 0_T});
+ env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, bField);
- // setup particle stack, and add primary particle
- setup::Stack<EnvType> stack;
- stack.clear();
unsigned short const A = std::stoi(std::string(argv[1]));
unsigned short const Z = std::stoi(std::string(argv[2]));
Code const beamCode = get_nucleus_code(A, Z);
@@ -194,12 +193,6 @@ int main(int argc, char** argv) {
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz});
- CORSIKA_LOG_INFO(
- "input particle: {}, "
- "input angles: theta={}, "
- "input momentum: {} GeV, "
- ", norm={}",
- beamCode, theta, plab.getComponents() / 1_GeV, plab.getNorm());
auto const observationHeight = 0_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
@@ -211,28 +204,25 @@ int main(int argc, char** argv) {
showerCore +
Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
- CORSIKA_LOG_INFO("point of injection: {} ", injectionPos.getCoordinates());
-
- stack.addParticle(std::make_tuple(
- Code::Proton, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
- plab.normalized(), injectionPos, 0_ns));
-
- CORSIKA_LOG_INFO("shower axis length: {} m",
- (showerCore - injectionPos).getNorm() * 1.02);
-
- OutputManager output("hybrid_MC_outputs");
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
+ auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
+ uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
+
+ CORSIKA_LOG_INFO("Primary particle: {}", beamCode);
+ CORSIKA_LOG_INFO("Zenith angle: {} (rad)", theta);
+ CORSIKA_LOG_INFO("Momentum: {} (GeV)", plab.getComponents() / 1_GeV);
+ CORSIKA_LOG_INFO("Propagation dir: {}", plab.getNorm());
+ CORSIKA_LOG_INFO("Injection point: {}", injectionPos.getCoordinates());
+ CORSIKA_LOG_INFO("shower axis length: {} ",
+ (showerCore - injectionPos).getNorm() * 1.02);
- // setup processes, decays and interactions
+ // SETUP WRITERS
- corsika::sibyll::Interaction sibyll{env};
- InteractionCounter sibyllCounted{sibyll};
-
- corsika::pythia8::Decay decayPythia;
+ OutputManager output("hybrid_MC_outputs");
// register energy losses as output
- EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
+ EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
output.add("energyloss", dEdX);
// create a track writer and register it with the output manager
@@ -242,19 +232,29 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(3_GeV, false, dEdX);
BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss(dEdX);
- LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
+ LongitudinalWriter profile{showerAxis, nAxisBins, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
+ Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+ ObservationPlane<TrackingType> observationLevel(obsPlane,
+ DirectionVector(rootCS, {1., 0., 0.}));
+ output.add("particles", observationLevel);
+
+ PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
+ output.add("primary", primaryWriter);
+
+ // SETUP PROCESSES, DECAYS, INTERACTIONS
+
+ corsika::sibyll::Interaction sibyll{env};
+ InteractionCounter sibyllCounted{sibyll};
+
+ corsika::pythia8::Decay decayPythia;
+
CONEXhybrid // SubWriter<decltype(dEdX>, SubWriter<decltype(profile)>>
conex_model(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton), dEdX,
profile);
- Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<setup::Tracking> observationLevel(
- obsPlane, DirectionVector(rootCS, {1., 0., 0.}), "particles.dat");
- output.add("obsplane", observationLevel);
-
corsika::urqmd::UrQMD urqmd_model;
InteractionCounter urqmdCounted{urqmd_model};
@@ -265,7 +265,7 @@ int main(int argc, char** argv) {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
- bool operator()(const setup::Stack<EnvType>::particle_type& p) const {
+ bool operator()(const StackType::particle_type& p) const {
return (p.getEnergy() < cutE_);
}
};
@@ -273,17 +273,27 @@ int main(int argc, char** argv) {
auto sequence = make_sequence(hadronSequence, decayPythia, eLoss, cut, conex_model,
longprof, observationLevel, trackCheck);
+ StackType stack;
+ stack.clear();
+
// define air shower object, run simulation
- setup::Tracking tracking;
+ TrackingType tracking;
Cascade EAS(env, tracking, sequence, output, stack);
+ output.startOfLibrary();
+
+ auto const primaryProperties = std::make_tuple(
+ Code::Proton, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns);
+
+ stack.addParticle(primaryProperties);
+ primaryWriter.recordPrimary(primaryProperties);
+
// to fix the point of first interaction, uncomment the following two lines:
// EAS.SetNodes();
// EAS.forceInteraction();
- output.startOfLibrary();
EAS.run();
- output.endOfLibrary();
const HEPEnergyType Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO(
@@ -296,5 +306,7 @@ int main(int argc, char** argv) {
save_hist(hists.labHist(), "inthist_lab_hybrid.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_hybrid.npz", true);
+ output.endOfLibrary();
+
CORSIKA_LOG_INFO("done");
}
diff --git a/examples/radio_em_shower.cpp b/examples/cascade_examples/radio_em_shower.cpp
similarity index 58%
rename from examples/radio_em_shower.cpp
rename to examples/cascade_examples/radio_em_shower.cpp
index c574755f7..65a0b1485 100644
--- a/examples/radio_em_shower.cpp
+++ b/examples/cascade_examples/radio_em_shower.cpp
@@ -21,10 +21,11 @@
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
-#include <corsika/output/OutputManager.hpp>
-#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
+#include <corsika/output/OutputManager.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
@@ -64,6 +65,13 @@
using namespace corsika;
using namespace std;
+//
+// An example of running an casacde which generates radio input for a
+// cascade that has hadronic interactions turned off. Currently, this
+// will produce output for only one antenna, but there are lines below,
+// which are commented out, to enable a full star-shape pattern of antennas
+//
+
void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("proposal");
@@ -72,14 +80,21 @@ void registerRandomStreams(int seed) {
if (seed == 0) {
std::random_device rd;
seed = rd();
- cout << "new random seed (auto) " << seed << endl;
+ CORSIKA_LOG_INFO("random seed (auto) {} ", seed);
+ } else {
+ CORSIKA_LOG_INFO("random seed {} ", seed);
}
RNGManager<>::getInstance().setSeed(seed);
}
+using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+using EnvType = Environment<EnvironmentInterface>;
template <typename TInterface>
using MyExtraEnv =
UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = setup::Tracking;
int main(int argc, char** argv) {
@@ -93,22 +108,21 @@ int main(int argc, char** argv) {
}
int seed{static_cast<int>(std::stof(std::string(argv[2])))};
- std::cout << "Seed: " << seed << std::endl;
+ CORSIKA_LOG_INFO("Seed {}", seed);
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams(seed);
// setup environment, geometry
- using EnvironmentInterface =
- IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
- using EnvType = Environment<EnvironmentInterface>;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
+ double const refractive_index = 1.000327;
+ MagneticFieldVector const bField{rootCS, 50_uT, 0_T, 0_T};
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
- env, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
- MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T});
+ env, AtmosphereId::LinsleyUSStd, center, refractive_index, Medium::AirDry1Atm,
+ bField);
std::unordered_map<Code, HEPEnergyType> energy_resolution = {
{Code::Electron, 5_MeV},
@@ -118,26 +132,19 @@ int main(int argc, char** argv) {
for (auto [pcode, energy] : energy_resolution)
set_energy_production_threshold(pcode, energy);
- // setup particle stack, and add primary particle
- setup::Stack<EnvType> stack;
- stack.clear();
- const Code beamCode = Code::Electron;
+ Code const beamCode = Code::Electron;
auto const mass = get_mass(beamCode);
- const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
- double theta = 0.;
+ HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[1]));
+ double const theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
- HEPMomentumType P0 = calculate_momentum(E0, mass);
+ HEPMomentumType const P0 = calculate_momentum(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
- auto plab = MomentumVector(rootCS, {px, py, pz});
- cout << "input particle: " << beamCode << endl;
- cout << "input angles: theta=" << theta << endl;
- cout << "input momentum: " << plab.getComponents() / 1_GeV
- << ", norm = " << plab.getNorm() << endl;
+ auto const plab = MomentumVector(rootCS, {px, py, pz});
auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
@@ -148,104 +155,103 @@ int main(int argc, char** argv) {
Point const injectionPos =
showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
- std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
-
- stack.addParticle(std::make_tuple(
- beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
- plab.normalized(), injectionPos, 0_ns));
-
- CORSIKA_LOG_INFO("shower axis length: {} ",
- (showerCore - injectionPos).getNorm() * 1.02);
-
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
+ auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
+ uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
+ // setup the radio antennas
TimeType const groundHitTime{(showerCore - injectionPos).getNorm() / constants::c};
- std::string outname_ = "radio_em_shower_outputs"; // + std::to_string(rr_);
- OutputManager output(outname_);
-
// Radio antennas and relevant information
// the antenna time variables
- const TimeType duration_{1e-6_s};
- const InverseTimeType sampleRate_{1e+9_Hz};
+ TimeType const duration{1e-6_s};
+ InverseTimeType const sampleRate{1e+9_Hz};
// the detector (aka antenna collection) for CoREAS and ZHS
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
AntennaCollection<TimeDomainAntenna> detectorZHS;
- auto const showerCoreX_{showerCore.getCoordinates().getX()};
- auto const showerCoreY_{showerCore.getCoordinates().getY()};
- auto const injectionPosX_{injectionPos.getCoordinates().getX()};
- auto const injectionPosY_{injectionPos.getCoordinates().getY()};
- auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
- auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
- std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
+ auto const showerCoreX{showerCore.getCoordinates().getX()};
+ auto const showerCoreY{showerCore.getCoordinates().getY()};
+ auto const injectionPosX{injectionPos.getCoordinates().getX()};
+ auto const injectionPosY{injectionPos.getCoordinates().getY()};
+ auto const injectionPosZ{injectionPos.getCoordinates().getZ()};
+ auto const triggerpoint{Point(rootCS, injectionPosX, injectionPosY, injectionPosZ)};
+
+ CORSIKA_LOG_INFO("Primary particle: {}", beamCode);
+ CORSIKA_LOG_INFO("Zenith angle: {} (rad)", theta);
+ CORSIKA_LOG_INFO("Momentum: {} (GeV)", plab.getComponents() / 1_GeV);
+ CORSIKA_LOG_INFO("Propagation dir: {}", plab.getNorm());
+ CORSIKA_LOG_INFO("Injection point: {}", injectionPos.getCoordinates());
+ CORSIKA_LOG_INFO("shower axis length: {} ",
+ (showerCore - injectionPos).getNorm() * 1.02);
+ CORSIKA_LOG_INFO("Trigger Point is: {}", triggerpoint);
// // setup CoREAS antennas - use the for loop for star shape pattern
- // for (auto radius_1 = 25_m; radius_1 <= 500_m; radius_1 += 25_m) {
- // for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
- auto radius_1 = 200_m;
- auto phi_1 = 45;
- auto phiRad_1 = phi_1 / 180. * M_PI;
- auto rr_1 = static_cast<int>(radius_1 / 1_m);
- auto const point_1{Point(rootCS, showerCoreX_ + radius_1 * cos(phiRad_1),
- showerCoreY_ + radius_1 * sin(phiRad_1),
- constants::EarthRadius::Mean)};
- std::cout << "Antenna point: " << point_1 << std::endl;
- auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
- std::string name_1 =
- "CoREAS_R=" + std::to_string(rr_1) + "_m--Phi=" + std::to_string(phi_1) + "degrees";
- TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_,
- sampleRate_, triggertime_1);
- detectorCoREAS.addAntenna(antenna_1);
+ // for (auto radius_coreas = 25_m; radius_coreas <= 500_m; radius_coreas += 25_m) {
+ // for (auto phi_coreas = 0; phi_coreas <= 315; phi_coreas += 45) {
+ auto radius_coreas = 200_m;
+ auto phi_coreas = 45;
+ auto phiRad_coreas = phi_coreas / 180. * M_PI;
+ auto rr_coreas = static_cast<int>(radius_coreas / 1_m);
+ auto const point_coreas{Point(rootCS, showerCoreX + radius_coreas * cos(phiRad_coreas),
+ showerCoreY + radius_coreas * sin(phiRad_coreas),
+ constants::EarthRadius::Mean)};
+ std::cout << "Antenna point: " << point_coreas << std::endl;
+ CORSIKA_LOG_INFO("Antenna point {}", injectionPos.getCoordinates());
+ auto triggertime_coreas{(triggerpoint - point_coreas).getNorm() / constants::c};
+ std::string name_coreas = "CoREAS_R=" + std::to_string(rr_coreas) +
+ "_m--Phi=" + std::to_string(phi_coreas) + "degrees";
+ TimeDomainAntenna antenna_coreas(name_coreas, point_coreas, rootCS, triggertime_coreas,
+ duration, sampleRate, triggertime_coreas);
+ detectorCoREAS.addAntenna(antenna_coreas);
// }
// }
// // setup ZHS antennas - use the for loop for star shape pattern
- // for (auto radius_ = 25_m; radius_ <= 500_m; radius_ += 25_m) {
- // for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
- auto radius_ = 200_m;
- auto phi_ = 45;
- auto phiRad_ = phi_ / 180. * M_PI;
- auto rr_ = static_cast<int>(radius_ / 1_m);
- auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
- showerCoreY_ + radius_ * sin(phiRad_),
- constants::EarthRadius::Mean)};
- auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
- std::string name_ =
- "ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
- TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_,
- triggertime_);
- detectorZHS.addAntenna(antenna_);
+ // for (auto radius_zhs = 25_m; radius_zhs <= 500_m; radius_zhs += 25_m) {
+ // for (auto phi_zhs = 0; phi_zhs <= 315; phi_zhs += 45) {
+ auto radius_zhs = 200_m;
+ auto phi_zhs = 45;
+ auto phiRad_zhs = phi_zhs / 180. * M_PI;
+ auto rr_zhs = static_cast<int>(radius_zhs / 1_m);
+ auto const point_zhs{Point(rootCS, showerCoreX + radius_zhs * cos(phiRad_zhs),
+ showerCoreY + radius_zhs * sin(phiRad_zhs),
+ constants::EarthRadius::Mean)};
+ auto triggertime_zhs{(triggerpoint - point_zhs).getNorm() / constants::c};
+ std::string name_zhs = "ZHS_R=" + std::to_string(rr_zhs) +
+ "_m--Phi=" + std::to_string(phi_zhs) + "degrees";
+ TimeDomainAntenna antenna_zhs(name_zhs, point_zhs, rootCS, triggertime_zhs, duration,
+ sampleRate, triggertime_zhs);
+ detectorZHS.addAntenna(antenna_zhs);
// }
// }
// setup processes, decays and interactions
-
- EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
- // register energy losses as output
- output.add("dEdX", dEdX);
-
- ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
-
- corsika::sophia::InteractionModel sophia;
+ EnergyLossWriter energyloss{showerAxis, dX, nAxisBins};
+ ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
+ energyloss);
corsika::sibyll::Interaction sibyll{env};
+ corsika::sophia::InteractionModel sophia;
HEPEnergyType heThresholdNN = 80_GeV;
corsika::proposal::Interaction emCascade(
env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
- corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
- // BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
+ corsika::proposal::ContinuousProcess<SubWriter<decltype(energyloss)>> emContinuous(
+ env, energyloss);
// NOT possible right now, due to interface differenc in PROPOSAL
// InteractionCounter emCascadeCounted(emCascade);
+ OutputManager output("radio_em_shower_outputs");
+
+ output.add("energyloss", energyloss);
+
TrackWriter tracks;
output.add("tracks", tracks);
- // long. profile
- LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
+ LongitudinalWriter profile{showerAxis, nAxisBins, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
@@ -269,15 +275,27 @@ int main(int argc, char** argv) {
output.add("ZHS", zhs);
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
+ ObservationPlane<TrackingType, ParticleWriterParquet> observationLevel{
obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
output.add("particles", observationLevel);
- // auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs);
+ PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
+ output.add("primary", primaryWriter);
+
auto sequence = make_sequence(emCascade, emContinuous, longprof, coreas, zhs, tracks,
observationLevel, cut);
// define air shower object, run simulation
- setup::Tracking tracking;
+ TrackingType tracking;
+
+ auto const primaryProperties = std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns);
+
+ // setup particle stack, and add primary particle
+ StackType stack;
+ stack.clear();
+ stack.addParticle(primaryProperties);
+ primaryWriter.recordPrimary(primaryProperties);
output.startOfLibrary();
Cascade EAS(env, tracking, sequence, output, stack);
@@ -287,7 +305,8 @@ int main(int argc, char** argv) {
EAS.run();
- HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+ HEPEnergyType const Efinal =
+ energyloss.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO(
"total energy budget (GeV): {}, "
@@ -295,4 +314,6 @@ int main(int argc, char** argv) {
Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
output.endOfLibrary();
+
+ return EXIT_SUCCESS;
}
diff --git a/examples/water.cpp b/examples/cascade_examples/water.cpp
similarity index 72%
rename from examples/water.cpp
rename to examples/cascade_examples/water.cpp
index 577ce4ad8..324e2e732 100644
--- a/examples/water.cpp
+++ b/examples/cascade_examples/water.cpp
@@ -24,15 +24,16 @@
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
-#include <corsika/modules/writers/SubWriter.hpp>
-#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
+#include <corsika/modules/writers/LongitudinalWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
+#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/Sophia.hpp>
-#include <corsika/modules/UrQMD.hpp>
+#include <corsika/modules/FLUKA.hpp>
#include <corsika/modules/tracking/TrackingStraight.hpp>
#include <corsika/output/OutputManager.hpp>
@@ -41,6 +42,7 @@
#include <corsika/stack/VectorStack.hpp>
#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
#include <CLI/App.hpp>
#include <CLI/Config.hpp>
@@ -50,9 +52,9 @@ using namespace corsika;
using IMediumType = IMediumPropertyModel<IMediumModel>;
using EnvType = Environment<IMediumType>;
-using StackActive = setup::Stack<EnvType>;
-using StackView = StackActive::stack_view_type;
-using Particle = StackActive::particle_type;
+using StackType = setup::Stack<EnvType>;
+using TrackingType = tracking_line::Tracking;
+using Particle = StackType::particle_type;
void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
@@ -61,7 +63,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("epos");
RNGManager<>::getInstance().registerRandomStream("pythia");
- RNGManager<>::getInstance().registerRandomStream("urqmd");
+ RNGManager<>::getInstance().registerRandomStream("fluka");
RNGManager<>::getInstance().registerRandomStream("proposal");
if (seed == 0) {
std::random_device rd;
@@ -75,12 +77,12 @@ void registerRandomStreams(int seed) {
int main(int argc, char** argv) {
// * process input
- Code primaryType;
- HEPEnergyType e0, eCut;
+ Code beamCode;
+ HEPEnergyType E0, eCut;
int A, Z, n_event;
int randomSeed;
std::string output_dir;
- CLI::App app{"Neutrino event generator"};
+ CLI::App app{"Cascade in water"};
// we start by definining a sub-group for the primary ID
auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary")
->check(CLI::Range(0, 26))
@@ -110,8 +112,25 @@ int main(int argc, char** argv) {
app.add_option("-s", randomSeed, "Seed for random number")
->check(CLI::NonNegativeNumber)
->default_val(0);
+
+ // parse the command line options into the variables
CLI11_PARSE(app, argc, argv);
+ std::string_view const loglevel = app["-v"]->as<std::string_view>();
+ if (loglevel == "warn") {
+ logging::set_level(logging::level::warn);
+ } else if (loglevel == "info") {
+ logging::set_level(logging::level::info);
+ } else if (loglevel == "debug") {
+ logging::set_level(logging::level::debug);
+ } else if (loglevel == "trace") {
+#ifndef _C8_DEBUG_
+ CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
+ return 1;
+#endif
+ logging::set_level(logging::level::trace);
+ }
+
// check that we got either PDG or A/Z
// this can be done with option_groups but the ordering
// gets all messed up
@@ -121,37 +140,24 @@ int main(int argc, char** argv) {
return 1;
}
}
+
+ // initialize random number sequence(s)
+ registerRandomStreams(randomSeed);
+
// check if we want to use a PDG code instead
if (app.count("--pdg") > 0) {
- primaryType = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
+ beamCode = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
} else {
// check manually for proton and neutrons
if ((A == 1) && (Z == 1))
- primaryType = Code::Proton;
+ beamCode = Code::Proton;
else if ((A == 1) && (Z == 0))
- primaryType = Code::Neutron;
+ beamCode = Code::Neutron;
else
- primaryType = get_nucleus_code(A, Z);
+ beamCode = get_nucleus_code(A, Z);
}
- e0 = app["-E"]->as<double>() * 1_GeV;
eCut = app["--eCut"]->as<double>() * 1_GeV;
- std::string_view const loglevel = app["-v"]->as<std::string_view>();
- if (loglevel == "warn") {
- logging::set_level(logging::level::warn);
- } else if (loglevel == "info") {
- logging::set_level(logging::level::info);
- } else if (loglevel == "debug") {
- logging::set_level(logging::level::debug);
- } else if (loglevel == "trace") {
-#ifndef _C8_DEBUG_
- CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
- return 1;
-#endif
- logging::set_level(logging::level::trace);
- }
-
- registerRandomStreams(randomSeed);
// * environment and universe
EnvType env;
@@ -163,27 +169,26 @@ int main(int argc, char** argv) {
Point const center{rootCS, 0_m, 0_m, 0_m};
auto sphere = std::make_unique<Sphere>(center, 100_m);
auto node = std::make_unique<VolumeTreeNode<IMediumType>>(std::move(sphere));
- // Hydrogen is not supported by UrQMD yet. See #456
- auto comp = NuclearComposition({{Code::Oxygen}, {1.0}});
+ NuclearComposition const nuclearComposition{{Code::Hydrogen, Code::Oxygen},
+ {2.0 / 3.0, 1.0 / 3.0}};
// density of sea water
auto density = 1.02_g / (1_cm * 1_cm * 1_cm);
auto water_medium =
std::make_shared<MediumPropertyModel<HomogeneousMedium<IMediumType>>>(
- Medium::WaterLiquid, density, comp);
+ Medium::WaterLiquid, density, nuclearComposition);
node->setModelProperties(water_medium);
universe->addChild(std::move(node));
}
- // * detector geometry
+ // * make downward-going shower axis and a observation plane in x-y-plane
auto injectorLength = 50_m;
- Point const injectorPos = Point(rootCS, {0_m, 0_m, injectorLength});
- auto const& injectCS = make_translation(rootCS, injectorPos.getCoordinates());
+ Point const injectionPos = Point(rootCS, {0_m, 0_m, injectorLength});
+ auto const& injectCS = make_translation(rootCS, injectionPos.getCoordinates());
DirectionVector upVec(rootCS, {0., 0., 1.});
DirectionVector leftVec(rootCS, {1., 0., 0.});
DirectionVector downVec(rootCS, {0., 0., -1.});
- // * observation plane
- std::vector<ObservationPlane<tracking_line::Tracking>> obsPlanes;
+ std::vector<ObservationPlane<TrackingType, ParticleWriterParquet>> obsPlanes;
const int nPlane = 5;
for (int i = 0; i < nPlane - 1; i++) {
Point planeCenter{injectCS, {0_m, 0_m, -(i + 1) * 3_m}};
@@ -193,12 +198,14 @@ int main(int argc, char** argv) {
Plane{Point{injectCS, {0_m, 0_m, -50_m}}, upVec}, leftVec, true);
// * longitutional profile
- ShowerAxis const showerAxis{injectorPos, 1.2 * injectorLength * downVec, env};
- LongitudinalWriter longiWriter{showerAxis, 5500, 1_g / square(1_cm)};
+ ShowerAxis const showerAxis{injectionPos, 1.2 * injectorLength * downVec, env};
+ auto const dX = 1_g / square(1_cm); // Binning of the writers along the shower axis
+ uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
+ LongitudinalWriter longiWriter{showerAxis, nAxisBins, dX};
LongitudinalProfile<SubWriter<decltype(longiWriter)>> longprof{longiWriter};
// * energy loss profile
- EnergyLossWriter dEdX{showerAxis, 1_g / square(1_cm), 5500};
+ EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
// * physical process list
// particle production threshold
@@ -207,19 +214,19 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(emCut, emCut, hadCut, hadCut, true, dEdX);
// hadronic interactions
- HEPEnergyType heHadronModelThreshold = 63.1_GeV;
+ HEPEnergyType heHadronModelThreshold = std::pow(10, 1.9) * 1_GeV;
corsika::sibyll::Interaction sibyll(env);
- corsika::urqmd::UrQMD urqmd;
- InteractionCounter urqmdCounted(urqmd);
+
+ corsika::fluka::Interaction leIntModel{env};
+ InteractionCounter leIntCounted{leIntModel};
struct EnergySwitch {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
};
- // auto lowModel = make_sequence(urqmd);
auto hadronSequence =
- make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, sibyll);
+ make_select(EnergySwitch(heHadronModelThreshold), leIntCounted, sibyll);
// decay process
corsika::pythia8::Decay decayPythia;
@@ -243,19 +250,41 @@ int main(int argc, char** argv) {
// hard coded
auto obsPlaneSequence =
make_sequence(obsPlanes[0], obsPlanes[1], obsPlanes[2], obsPlanes[3], obsPlanes[4]);
- output.add("longi_profile", longiWriter);
- output.add("energy_loss", dEdX);
+
+ PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(obsPlanes.back());
+ output.add("primary", primaryWriter);
+
+ output.add("profile", longiWriter);
+ output.add("energyloss", dEdX);
// * the final process sequence
auto sequence = make_sequence(physics_sequence, longprof, obsPlaneSequence, cut);
// * tracking and stack
- tracking_line::Tracking tracking;
- StackActive stack;
+ TrackingType tracking;
+ StackType stack;
// * cascade manager
Cascade EAS(env, tracking, sequence, output, stack);
+ E0 = app["-E"]->as<double>() * 1_GeV;
+ HEPEnergyType mass = get_mass(beamCode);
+ // convert Elab to Plab
+ HEPMomentumType P0 = calculate_momentum(E0, mass);
+ auto plab = MomentumVector(rootCS, P0 * downVec.getNorm());
+
+ // print our primary parameters all in one place
+ if (app["--pdg"]->count() > 0) {
+ CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as<int>());
+ } else {
+ CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
+ }
+ CORSIKA_LOG_INFO("Primary Energy: {}", E0);
+ CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
+ CORSIKA_LOG_INFO("Primary Direction: {}", plab.getNorm());
+ CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates());
+ CORSIKA_LOG_INFO("Shower Axis Length: {}", injectorLength);
+
// * main loop
output.startOfLibrary();
for (int i_shower = 0; i_shower < n_event; i_shower++) {
@@ -263,8 +292,11 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("Event: {} / {}", i_shower, n_event);
// * inject primary
- auto primary = stack.addParticle(std::make_tuple(
- primaryType, e0 - get_mass(primaryType), downVec, injectorPos, 0_ns));
+ auto const primaryProperties = std::make_tuple(
+ beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
+ plab.normalized(), injectionPos, 0_ns);
+ auto primary = stack.addParticle(primaryProperties);
+ stack.addParticle(primaryProperties);
EAS.run();
@@ -273,8 +305,8 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO(
"total energy budget (TeV): {:.2f} (dEdX={:.2f} ground={:.2f}), "
"relative difference (%): {:.3f}",
- e0 / 1_TeV, dEdX.getEnergyLost() / 1_TeV, obsPlaneFinal.getEnergyGround() / 1_TeV,
- (Efinal / e0 - 1.) * 100.);
+ E0 / 1_TeV, dEdX.getEnergyLost() / 1_TeV, obsPlaneFinal.getEnergyGround() / 1_TeV,
+ (Efinal / E0 - 1.) * 100.);
}
output.endOfLibrary();
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
deleted file mode 100644
index 68c7c72cd..000000000
--- a/examples/cascade_proton_example.cpp
+++ /dev/null
@@ -1,147 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/process/SwitchProcessSequence.hpp>
-#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/core/EnergyMomentumOperations.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-
-#include <corsika/output/OutputManager.hpp>
-#include <corsika/modules/writers/SubWriter.hpp>
-#include <corsika/modules/writers/EnergyLossWriter.hpp>
-
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/HomogeneousMedium.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/media/UniformMagneticField.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <corsika/modules/BetheBlochPDG.hpp>
-#include <corsika/modules/StackInspector.hpp>
-#include <corsika/modules/Sibyll.hpp>
-#include <corsika/modules/ParticleCut.hpp>
-#include <corsika/modules/TrackWriter.hpp>
-#include <corsika/modules/HadronicElasticModel.hpp>
-#include <corsika/modules/Pythia8.hpp>
-#include <corsika/modules/Sibyll.hpp>
-
-#include <iostream>
-#include <limits>
-#include <typeinfo>
-
-using namespace corsika;
-using namespace std;
-
-//
-// The example main program for a particle cascade
-//
-int main() {
-
- logging::set_level(logging::level::warn);
-
- std::cout << "cascade_proton_example" << std::endl;
-
- feenableexcept(FE_INVALID);
- // initialize random number sequence(s)
- RNGManager<>::getInstance().registerRandomStream("cascade");
-
- OutputManager output("cascade_proton_outputs");
-
- // setup environment, geometry
- using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
- using EnvType = Environment<EnvironmentInterface>;
- EnvType env;
- auto& universe = *(env.getUniverse());
- CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
-
- auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
-
- using MyHomogeneousModel =
- MediumPropertyModel<UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>;
-
- world->setModelProperties<MyHomogeneousModel>(
- Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 1_mT),
- 1_kg / (1_m * 1_m * 1_m), NuclearComposition({Code::Proton}, {1.}));
-
- universe.addChild(std::move(world));
-
- // setup particle stack, and add primary particle
- setup::Stack<EnvType> stack;
- stack.clear();
- const Code beamCode = Code::Proton;
- const HEPMassType mass = Proton::mass;
- const HEPEnergyType E0 = 200_GeV;
- double theta = 0.;
- double phi = 0.;
-
- Point injectionPos(rootCS, 0_m, 0_m, 0_m);
- {
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt(Elab * Elab - m * 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 = 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;
- stack.addParticle(std::make_tuple(
- beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
- plab.normalized(), injectionPos, 0_ns));
- }
-
- // setup processes, decays and interactions
- setup::Tracking tracking;
- StackInspector<setup::Stack<EnvType>> stackInspect(1000, true, E0);
-
- RNGManager<>::getInstance().registerRandomStream("pythia");
- RNGManager<>::getInstance().registerRandomStream("sibyll");
- corsika::sibyll::Interaction sibyll{env};
- // corsika::pythia8::Interaction pythia;
- corsika::pythia8::Decay decay;
-
- ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
- EnergyLossWriter dEdX{showerAxis};
- output.add("energyloss", dEdX);
-
- BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
- ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, true, dEdX);
- cut.printThresholds();
-
- // RNGManager::getInstance().registerRandomStream("HadronicElasticModel");
- // HadronicElasticModel::HadronicElasticInteraction
- // hadronicElastic(env);
-
- TrackWriter trackWriter;
- output.add("tracks", trackWriter); // register TrackWriter
-
- // assemble all processes into an ordered process list
- auto sequence = make_sequence(sibyll, decay, eLoss, trackWriter, stackInspect, cut);
-
- // define air shower object, run simulation
- Cascade EAS(env, tracking, sequence, output, stack);
- output.startOfLibrary();
- EAS.run();
- output.endOfLibrary();
-
- CORSIKA_LOG_INFO("Done");
-}
diff --git a/examples/boundary_example.cpp b/examples/framework_examples/boundary_crossing.cpp
similarity index 94%
rename from examples/boundary_example.cpp
rename to examples/framework_examples/boundary_crossing.cpp
index 895b5214e..d470a4c3c 100644
--- a/examples/boundary_example.cpp
+++ b/examples/framework_examples/boundary_crossing.cpp
@@ -37,6 +37,13 @@
using namespace corsika;
using namespace std;
+//
+// This example shows how to make a custom process which deletes particles that
+// cross a particular boundary (a sphere in this case)
+// For a plane boundary, this can be implemented by adding an ObservationPlane
+// or Observation Volume object instead (the standard "absorbing" geometry objects)
+//
+
template <bool deleteParticle>
struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
@@ -65,9 +72,6 @@ private:
std::ofstream file_;
};
-//
-// The example main program for a particle cascade
-//
int main() {
logging::set_level(logging::level::warn);
diff --git a/examples/environment.cpp b/examples/framework_examples/environment.cpp
similarity index 100%
rename from examples/environment.cpp
rename to examples/framework_examples/environment.cpp
diff --git a/examples/geometry_example.cpp b/examples/framework_examples/geometry.cpp
similarity index 100%
rename from examples/geometry_example.cpp
rename to examples/framework_examples/geometry.cpp
diff --git a/examples/helix_example.cpp b/examples/framework_examples/helix_trajectory.cpp
similarity index 100%
rename from examples/helix_example.cpp
rename to examples/framework_examples/helix_trajectory.cpp
diff --git a/examples/particle_list_example.cpp b/examples/framework_examples/known_particles.cpp
similarity index 85%
rename from examples/particle_list_example.cpp
rename to examples/framework_examples/known_particles.cpp
index 2c9016f52..4f6cf35ba 100644
--- a/examples/particle_list_example.cpp
+++ b/examples/framework_examples/known_particles.cpp
@@ -18,11 +18,14 @@ using namespace corsika;
using namespace std;
//
-// The example main program for a particle list
+// This example prints out all of the particles that are
+// known to CORSIKA 8 and the corresponding IDs in the
+// hadronic model codes.
//
+
int main() {
- logging::set_level(logging::level::warn);
+ logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
logging::info(
@@ -30,7 +33,7 @@ int main() {
"------------------------------------------\n"
" particles in CORSIKA\n"
"------------------------------------------\n");
- int const width = 20 + 10 + 10 + 10 + 15 + 15 + 17;
+ int const width = 20 + 10 + 10 + 10 + 16 + 16 + 17;
logging::info(
"Name | "
"PDG-id | "
@@ -46,7 +49,7 @@ int main() {
? to_string(corsika::sibyll::getSibyllMass(p) / 1_GeV)
: "");
auto const qgs_id = corsika::qgsjetII::convertToQgsjetII(p);
- logging::info("{:20} | {:10} | {:10} | {:10} | {:>15.5} | {:>15.5} |", p,
+ logging::info("{:20} | {:10} | {:10} | {:10} | {:>16.5} | {:>16.5} |", p,
static_cast<int>(get_PDG(p)),
(sib_id != corsika::sibyll::SibyllCode::Unknown
? to_string(static_cast<int>(sib_id))
diff --git a/examples/stack_example.cpp b/examples/framework_examples/stack.cpp
similarity index 100%
rename from examples/stack_example.cpp
rename to examples/framework_examples/stack.cpp
diff --git a/examples/staticsequence_example.cpp b/examples/framework_examples/static_sequence.cpp
similarity index 99%
rename from examples/staticsequence_example.cpp
rename to examples/framework_examples/static_sequence.cpp
index b9bbba59f..b3242dd5e 100644
--- a/examples/staticsequence_example.cpp
+++ b/examples/framework_examples/static_sequence.cpp
@@ -131,4 +131,4 @@ int main() {
modular();
return 0;
-}
\ No newline at end of file
+}
diff --git a/examples/stopping_power.cpp b/examples/physics_examples/stopping_power.cpp
similarity index 92%
rename from examples/stopping_power.cpp
rename to examples/physics_examples/stopping_power.cpp
index 8b39786e1..28a9ff7f0 100644
--- a/examples/stopping_power.cpp
+++ b/examples/physics_examples/stopping_power.cpp
@@ -9,7 +9,6 @@
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMediumModel.hpp>
-#include <corsika/media/ShowerAxis.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
@@ -30,7 +29,7 @@ using namespace std;
//
int main() {
- logging::set_level(logging::level::warn);
+ logging::set_level(logging::level::info);
CORSIKA_LOG_INFO("stopping_power");
@@ -52,7 +51,9 @@ int main() {
setup::Stack<EnvType> stack;
- std::ofstream file("dEdX.dat");
+ std::string fileName = "dEdX.dat";
+ CORSIKA_LOG_INFO("Writing to file {}", fileName);
+ std::ofstream file(fileName);
file << "# beta*gamma, dE/dX / MeV/(g/cm²)" << std::endl;
for (HEPEnergyType E0 = 200_MeV; E0 < 1_PeV; E0 *= 1.05) {
@@ -79,5 +80,5 @@ int main() {
HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm));
file << P0 / mass << "\t" << -dE / 1_MeV << std::endl;
}
- CORSIKA_LOG_INFO("finished writing dEdX.dat");
+ CORSIKA_LOG_INFO("finished writing {}", fileName);
}
diff --git a/examples/clover_leaf.cpp b/examples/physics_examples/synchrotron_clover_leaf.cpp
similarity index 90%
rename from examples/clover_leaf.cpp
rename to examples/physics_examples/synchrotron_clover_leaf.cpp
index f1044b9b3..a200ea101 100644
--- a/examples/clover_leaf.cpp
+++ b/examples/physics_examples/synchrotron_clover_leaf.cpp
@@ -24,7 +24,6 @@
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp>
-#include <corsika/media/ShowerAxis.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
@@ -35,14 +34,10 @@
#include <corsika/modules/radio/antennas/Antenna.hpp>
#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
-#include <corsika/modules/radio/propagators/NumericalIntegratingPropagator.hpp>
#include <corsika/modules/radio/propagators/DummyTestPropagator.hpp>
-#include <corsika/modules/TrackWriter.hpp>
+#include <corsika/modules/writers/PrimaryWriter.hpp>
-#include <corsika/modules/StackInspector.hpp>
-#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/TimeCut.hpp>
-// #include <corsika/modules/TrackWriter.hpp>
#include <iomanip>
#include <iostream>
@@ -57,6 +52,7 @@ using namespace std;
// A simple shower to get the electric field trace of an electron (& a positron)
// in order to reveal the "clover leaf" pattern of energy fluence to the ground.
//
+
int main() {
logging::set_level(logging::level::info);
@@ -70,7 +66,7 @@ int main() {
auto seed = rd();
RNGManager<>::getInstance().setSeed(seed);
- OutputManager output("1 electron - 1 positron");
+ OutputManager output("clover_leaf_outputs");
// create a suitable environment
using IModelInterface =
@@ -105,7 +101,6 @@ int main() {
auto const injectionPosY_{injectionPos.getCoordinates().getY()};
auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
- std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
// the antenna characteristics
const TimeType duration_{2e-6_s}; // 0.8e-4_s
@@ -113,7 +108,6 @@ int main() {
// the detectors
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
- // AntennaCollection<TimeDomainAntenna> detectorZHS;
std::string name_center = "CoREAS_R=0_m--Phi=0degrees";
auto triggertime_center{((triggerpoint_ - center).getNorm() / constants::c) - 500_ns};
@@ -127,7 +121,7 @@ int main() {
auto rr_1 = static_cast<int>(radius_1 / 1_m);
auto const point_1{Point(rootCS, centerX + radius_1 * cos(phiRad_1),
centerY + radius_1 * sin(phiRad_1), centerZ)};
- std::cout << "Antenna point: " << point_1 << std::endl;
+ CORSIKA_LOG_INFO("Antenna point: {}", point_1);
auto triggertime_1{((triggerpoint_ - point_1).getNorm() / constants::c) - 500_ns};
std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
"_m--Phi=" + std::to_string(phi_1) + "degrees";
@@ -186,16 +180,8 @@ int main() {
coreas(detectorCoREAS, SP);
output.add("CoREAS", coreas);
- // RadioProcess<decltype(detectorZHS), ZHS<decltype(detectorZHS),
- // decltype(SP)>, decltype(SP)>
- // zhs(detectorZHS, SP);
- // output.add("ZHS", zhs);
-
TimeCut cut(period / 4);
- // TrackWriter trackWriter;
- // output.add("tracks", trackWriter); // register TrackWriter
-
// assemble all processes into an ordered process list
auto sequence = make_sequence(coreas, cut);
diff --git a/examples/synchrotron_test_C8tracking.cpp b/examples/physics_examples/synchrotron_test_C8tracking.cpp
similarity index 100%
rename from examples/synchrotron_test_C8tracking.cpp
rename to examples/physics_examples/synchrotron_test_C8tracking.cpp
diff --git a/examples/synchrotron_test_manual_tracking.cpp b/examples/physics_examples/synchrotron_test_manual_tracking.cpp
similarity index 100%
rename from examples/synchrotron_test_manual_tracking.cpp
rename to examples/physics_examples/synchrotron_test_manual_tracking.cpp
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
deleted file mode 100644
index 4a64fc19b..000000000
--- a/examples/vertical_EAS.cpp
+++ /dev/null
@@ -1,286 +0,0 @@
-/*
- * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#define TRACE
-
-/* clang-format off */
-// InteractionCounter used boost/histogram, which
-// fails if boost/type_traits have been included before. Thus, we have
-// to include it first...
-#include <corsika/framework/process/InteractionCounter.hpp>
-/* clang-format on */
-#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/core/EnergyMomentumOperations.hpp>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/process/SwitchProcessSequence.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
-
-#include <corsika/modules/writers/EnergyLossWriter.hpp>
-#include <corsika/modules/writers/LongitudinalWriter.hpp>
-#include <corsika/modules/writers/SubWriter.hpp>
-#include <corsika/output/OutputManager.hpp>
-
-#include <corsika/media/CORSIKA7Atmospheres.hpp>
-#include <corsika/media/Environment.hpp>
-#include <corsika/media/FlatExponential.hpp>
-#include <corsika/media/GeomagneticModel.hpp>
-#include <corsika/media/HomogeneousMedium.hpp>
-#include <corsika/media/IMagneticFieldModel.hpp>
-#include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/media/NuclearComposition.hpp>
-#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/UniformMagneticField.hpp>
-
-#include <corsika/modules/BetheBlochPDG.hpp>
-#include <corsika/modules/Epos.hpp>
-#include <corsika/modules/LongitudinalProfile.hpp>
-#include <corsika/modules/ObservationPlane.hpp>
-#include <corsika/modules/PROPOSAL.hpp>
-#include <corsika/modules/ParticleCut.hpp>
-#include <corsika/modules/Pythia8.hpp>
-#include <corsika/modules/Sibyll.hpp>
-#include <corsika/modules/Sophia.hpp>
-#include <corsika/modules/StackInspector.hpp>
-#include <corsika/modules/TrackWriter.hpp>
-#include <corsika/modules/UrQMD.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <iomanip>
-#include <iostream>
-#include <limits>
-#include <string>
-
-using namespace corsika;
-using namespace std;
-
-using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
-using EnvType = Environment<EnvironmentInterface>;
-
-using Particle = setup::Stack<EnvType>::particle_type;
-
-void registerRandomStreams(int seed) {
- RNGManager<>::getInstance().registerRandomStream("cascade");
- RNGManager<>::getInstance().registerRandomStream("sibyll");
- // RNGManager<>::getInstance().registerRandomStream("sophia");
- RNGManager<>::getInstance().registerRandomStream("pythia");
- RNGManager<>::getInstance().registerRandomStream("urqmd");
- RNGManager<>::getInstance().registerRandomStream("proposal");
- if (seed == 0) {
- std::random_device rd;
- seed = rd();
- cout << "new random seed (auto) " << seed << endl;
- }
- RNGManager<>::getInstance().setSeed(seed);
-}
-
-template <typename T>
-using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
-
-// argv : 1.number of nucleons, 2.number of protons,
-// 3.total energy in GeV, 4.number of showers,
-// 5.seed (0 by default to generate random values for all)
-
-int main(int argc, char** argv) {
-
- logging::set_level(logging::level::warn);
-
- CORSIKA_LOG_INFO("vertical_EAS");
-
- if (argc < 4) {
- std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed] \n"
- " if A=0, Z is interpreted as PDG code \n"
- " if no seed is given, a random seed is chosen \n"
- << std::endl;
- return 1;
- }
- feenableexcept(FE_INVALID);
-
- int seed = 0;
-
- if (argc > 4) { seed = std::stoi(std::string(argv[4])); }
-
- // initialize random number sequence(s)
- registerRandomStreams(seed);
-
- // setup environment, geometry
- EnvType env;
- CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
- Point const center{rootCS, 0_m, 0_m, 0_m};
-
- // build a Linsley US Standard atmosphere into `env`
- create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
- env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
- MagneticFieldVector{rootCS, 20.4_uT, 0_T, -43.23_uT});
-
- // Uncomment if you want to use PROPOSAL
- // std::unordered_map<Code, HEPEnergyType> energy_resolution = {
- // {Code::Electron, 2_MeV},
- // {Code::Positron, 2_MeV},
- // {Code::Photon, 2_MeV},
- // };
- // for (auto [pcode, energy] : energy_resolution)
- // set_energy_production_threshold(pcode, energy);
-
- // pre-setup particle stack
- unsigned short const A = std::stoi(std::string(argv[1]));
- Code beamCode;
- unsigned short Z = 0;
- if (A > 0) {
- Z = std::stoi(std::string(argv[2]));
- if (A == 1 && Z == 0)
- beamCode = Code::Neutron;
- else if (A == 1 && Z == 1)
- beamCode = Code::Proton;
- else
- beamCode = get_nucleus_code(A, Z);
- } else {
- int pdg = std::stoi(std::string(argv[2]));
- beamCode = convert_from_PDG(PDGCode(pdg));
- }
- HEPEnergyType mass = get_mass(beamCode);
- HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
- double theta = 0.;
- double phi = 180.;
- auto const thetaRad = theta / 180. * constants::pi;
- auto const phiRad = phi / 180. * constants::pi;
-
- HEPMomentumType P0 = calculate_momentum(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(thetaRad, phiRad, P0);
- auto plab = 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
- << ", norm = " << plab.getNorm() << endl;
-
- auto const observationHeight = 0.0_km + constants::EarthRadius::Mean;
- auto const injectionHeight = 111.75_km + constants::EarthRadius::Mean;
- auto const t = -observationHeight * cos(thetaRad) +
- sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
- static_pow<2>(injectionHeight));
- Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
- Point const injectionPos =
- showerCore + DirectionVector{rootCS,
- {-sin(thetaRad) * cos(phiRad),
- -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
- t;
-
- std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
-
- // we make the axis much longer than the inj-core distance since the
- // profile will go beyond the core, depending on zenith angle
- std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
- << std::endl;
-
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env, false,
- 1000};
-
- // create the output manager that we then register outputs with
- OutputManager output("vertical_EAS_outputs");
-
- // EnergyLossWriterParquet dEdX_output{showerAxis, 10_g / square(1_cm), 200};
- // register energy losses as output
- // output.add("dEdX", dEdX_output);
- // register profile output
-
- // construct the overall energy loss writer and register it
- EnergyLossWriter dEdX{showerAxis};
- output.add("energyloss", dEdX);
-
- // construct the continuous energy loss model
- BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
-
- // construct a particle cut - cuts are set to values close to reality, put
- // higher values for faster runs
- ParticleCut<SubWriter<decltype(dEdX)>> cut{2_MeV, 2_MeV, 2_GeV, 300_MeV, true, dEdX};
-
- // setup longitudinal profile
- LongitudinalWriter longProf{showerAxis};
- output.add("profile", longProf);
-
- LongitudinalProfile<SubWriter<decltype(longProf)>> profile{longProf};
-
- // create a track writer and register it with the output manager
- TrackWriter<TrackWriterParquet> trackWriter;
- output.add("tracks", trackWriter);
-
- // setup processes, decays and interactions
-
- corsika::sibyll::Interaction sibyll{env};
- InteractionCounter sibyllCounted{sibyll};
-
- HEPEnergyType heThresholdNN = 60_GeV;
- // PROPOSAL is disabled for this example
- // corsika::sophia::InteractionModel sophia;
- // corsika::proposal::Interaction emCascade(env, sophia,
- // sibyll.getHadronInteractionModel(), heThresholdNN);
- // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
- // emContinuous(env, dEdX);
-
- corsika::pythia8::Decay decayPythia;
-
- corsika::urqmd::UrQMD urqmd;
- InteractionCounter urqmdCounted{urqmd};
- StackInspector<setup::Stack<EnvType>> stackInspect(50000, false, E0);
-
- // assemble all processes into an ordered process list
- struct EnergySwitch {
- HEPEnergyType cutE_;
- EnergySwitch(HEPEnergyType cutE)
- : cutE_(cutE) {}
- bool operator()(const Particle& p) const { return (p.getEnergy() < cutE_); }
- };
- auto hadronSequence = make_select(EnergySwitch(55_GeV), urqmdCounted, sibyllCounted);
-
- // directory for outputs
- string const labHist_file = "inthist_lab_verticalEAS.npz";
- string const cMSHist_file = "inthist_cms_verticalEAS.npz";
-
- // setup particle stack, and add primary particle
- setup::Stack<EnvType> stack;
- stack.clear();
-
- stack.addParticle(std::make_tuple(
- beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
- plab.normalized(), injectionPos, 0_ns));
-
- Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
- obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
- output.add("particles", observationLevel);
-
- auto sequence = make_sequence(stackInspect, hadronSequence, decayPythia, emContinuous,
- profile, observationLevel, cut);
-
- // define air shower object, run simulation
- setup::Tracking tracking;
- output.startOfLibrary();
- Cascade EAS(env, tracking, sequence, output, stack);
- EAS.run();
-
- auto const hists = sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
-
- save_hist(hists.labHist(), labHist_file, true);
- save_hist(hists.CMSHist(), cMSHist_file, true);
-
- output.endOfLibrary();
-}
--
GitLab