From 234f612365b433601eab73ad2d79a4c68e91878e Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sun, 9 Oct 2022 22:55:54 +0200
Subject: [PATCH] use LE photon model in examples
---
corsika/modules/Random.hpp | 1 +
examples/corsika.cpp | 140 +++++++++++++++++++++----------------
examples/em_shower.cpp | 74 +++++++++++---------
3 files changed, 123 insertions(+), 92 deletions(-)
diff --git a/corsika/modules/Random.hpp b/corsika/modules/Random.hpp
index 0ab888ab8..a81bc1ad2 100644
--- a/corsika/modules/Random.hpp
+++ b/corsika/modules/Random.hpp
@@ -17,6 +17,7 @@
link this togehter, it will fail.
*/
#include <corsika/modules/sibyll/Random.hpp>
+#include <corsika/modules/sophia/Random.hpp>
#include <corsika/modules/epos/Random.hpp>
#include <corsika/modules/urqmd/Random.hpp>
#include <corsika/modules/qgsjetII/Random.hpp>
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index ca65c0c52..5d8eddec7 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -12,57 +12,58 @@
// 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/core/Logging.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
-#include <corsika/framework/utility/CorsikaFenv.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/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/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/LayeredSphericalAtmosphereBuilder.hpp>
-#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/CORSIKA7Atmospheres.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/OnShellCheck.hpp>
-#include <corsika/modules/StackInspector.hpp>
-#include <corsika/modules/TrackWriter.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
+#include <corsika/modules/QGSJetII.hpp>
#include <corsika/modules/Sibyll.hpp>
-#include <corsika/modules/Epos.hpp>
+#include <corsika/modules/Sophia.hpp>
+#include <corsika/modules/StackInspector.hpp>
+#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/UrQMD.hpp>
-#include <corsika/modules/PROPOSAL.hpp>
-#include <corsika/modules/QGSJetII.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <CLI/App.hpp>
-#include <CLI/Formatter.hpp>
#include <CLI/Config.hpp>
+#include <CLI/Formatter.hpp>
#include <iomanip>
#include <limits>
@@ -81,7 +82,8 @@
using namespace corsika;
using namespace std;
-using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using EnvironmentInterface =
+ IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
using Particle = setup::Stack<EnvType>::particle_type;
@@ -90,6 +92,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("qgsjet");
RNGManager<>::getInstance().registerRandomStream("sibyll");
+ RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("epos");
RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("urqmd");
@@ -107,7 +110,7 @@ void registerRandomStreams(int seed) {
template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
-int main(int argc, char** argv) {
+int main(int argc, char **argv) {
// the main command line description
CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."};
@@ -168,7 +171,8 @@ int main(int argc, char** argv) {
CLI11_PARSE(app, argc, argv);
if (app.count("--verbosity")) {
- std::string_view const loglevel = app["--verbosity"]->as<std::string_view>();
+ std::string_view const loglevel =
+ app["--verbosity"]->as<std::string_view>();
if (loglevel == "warn") {
logging::set_level(logging::level::warn);
} else if (loglevel == "info") {
@@ -189,7 +193,8 @@ int main(int argc, char** argv) {
// 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.");
+ CORSIKA_LOG_ERROR(
+ "If --pdg is not provided, then both -A and -Z are required.");
return 1;
}
}
@@ -199,7 +204,7 @@ int main(int argc, char** argv) {
/* === START: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
EnvType env;
- CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ CoordinateSystemPtr const &rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
GeomagneticModel wmm(center, corsika_data("GeoMag/WMM.COF"));
@@ -213,9 +218,10 @@ int main(int argc, char** argv) {
ofstream atmout("earth.dat");
for (LengthType h = 0_m; h < 110_km; h += 100_m) {
Point const ptest{rootCS, 0_m, 0_m, constants::EarthRadius::Mean + h};
- auto rho =
- env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity(
- ptest);
+ auto rho = env.getUniverse()
+ ->getContainingNode(ptest)
+ ->getModelProperties()
+ .getMassDensity(ptest);
atmout << h / 1_m << " " << rho / 1_kg * cube(1_m) << "\n";
}
atmout.close();
@@ -251,8 +257,8 @@ int main(int argc, char** argv) {
// convert the momentum to the zenith and azimuth angle of the primary
auto const [px, py, pz] =
- std::make_tuple(P0 * sin(thetaRad) * cos(phiRad), P0 * sin(thetaRad) * sin(phiRad),
- -P0 * cos(thetaRad));
+ std::make_tuple(P0 * sin(thetaRad) * cos(phiRad),
+ P0 * sin(thetaRad) * sin(phiRad), -P0 * cos(thetaRad));
auto plab = MomentumVector(rootCS, {px, py, pz});
/* === END: CONSTRUCT PRIMARY PARTICLE === */
@@ -264,14 +270,16 @@ int main(int argc, char** argv) {
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;
+ showerCore +
+ DirectionVector{rootCS,
+ {-sin(thetaRad) * cos(phiRad),
+ -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
+ t;
// 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};
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2,
+ env};
/* === END: CONSTRUCT GEOMETRY === */
// create the output manager that we then register outputs with
@@ -313,18 +321,23 @@ int main(int argc, char** argv) {
// decaySibyll.printDecayConfig();
+ // hadronic photon interactions in resonance region
+ corsika::sophia::InteractionModel sophia;
+
HEPEnergyType const emcut = 50_GeV;
HEPEnergyType const hadcut = 50_GeV;
- ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, dEdX);
+ ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true,
+ dEdX);
- // energy threshold for high energy hadronic model. Affects LE/HE switch for hadron
- // interactions and the hadronic photon model in proposal
+ // energy threshold for high energy hadronic model. Affects LE/HE switch for
+ // hadron interactions and the hadronic photon model in proposal
HEPEnergyType heHadronModelThreshold = 63.1_GeV;
- corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
- heHadronModelThreshold);
+ corsika::proposal::Interaction emCascade(
+ env, sophia, sibyll.getHadronInteractionModel(), heHadronModelThreshold);
// NOT available for PROPOSAL due to interface trouble:
// InteractionCounter emCascadeCounted(emCascade);
- // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env);
+ // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
+ // emContinuous(env);
BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
LongitudinalWriter profile{showerAxis, 200, 10_g / square(1_cm)};
@@ -338,12 +351,13 @@ int main(int argc, char** argv) {
// 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.getKineticEnergy() < cutE_); }
+ EnergySwitch(HEPEnergyType cutE) : cutE_(cutE) {}
+ bool operator()(const Particle &p) const {
+ return (p.getKineticEnergy() < cutE_);
+ }
};
- auto hadronSequence =
- make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, sibyllCounted);
+ auto hadronSequence = make_select(EnergySwitch(heHadronModelThreshold),
+ urqmdCounted, sibyllCounted);
auto decaySequence = make_sequence(decayPythia, decaySibyll);
TrackWriter trackWriter{tracks};
@@ -356,12 +370,13 @@ int main(int argc, char** argv) {
output.add("particles", observationLevel);
// assemble the final process sequence
- auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, cut,
- emCascade, emContinuous, // trackWriter,
+ auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence,
+ cut, emCascade, emContinuous, // trackWriter,
observationLevel, longprof);
/* === END: SETUP PROCESS LIST === */
- // create the cascade object using the default stack and tracking implementation
+ // create the cascade object using the default stack and tracking
+ // implementation
setup::Tracking tracking;
setup::Stack<EnvType> stack;
Cascade EAS(env, tracking, sequence, output, stack);
@@ -375,7 +390,8 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("Primary Energy: {}", E0);
CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates());
- CORSIKA_LOG_INFO("Shower Axis Length: {}", (showerCore - injectionPos).getNorm() * 1.2);
+ CORSIKA_LOG_INFO("Shower Axis Length: {}",
+ (showerCore - injectionPos).getNorm() * 1.2);
// trigger the output manager to open the library for writing
output.startOfLibrary();
@@ -387,8 +403,10 @@ int main(int argc, char** argv) {
// directory for outputs
string const outdir(app["--filename"]->as<std::string>());
- string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz";
- string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz";
+ string const labHist_file =
+ outdir + "/inthist_lab_" + to_string(i_shower) + ".npz";
+ string const cMSHist_file =
+ outdir + "/inthist_cms_" + to_string(i_shower) + ".npz";
// setup particle stack, and add primary particle
stack.clear();
@@ -410,14 +428,16 @@ int main(int argc, char** argv) {
HEPEnergyType const Efinal =
dEdX.getEnergyLost() + observationLevel.getEnergyGround();
- CORSIKA_LOG_INFO(
- "total energy budget (GeV): {} (dEdX={} ground={}), "
- "relative difference (%): {}",
- Efinal / 1_GeV, dEdX.getEnergyLost() / 1_GeV,
- observationLevel.getEnergyGround() / 1_GeV, (Efinal / E0 - 1) * 100);
+ CORSIKA_LOG_INFO("total energy budget (GeV): {} (dEdX={} ground={}), "
+ "relative difference (%): {}",
+ Efinal / 1_GeV, dEdX.getEnergyLost() / 1_GeV,
+ observationLevel.getEnergyGround() / 1_GeV,
+ (Efinal / E0 - 1) * 100);
- // auto const hists = heModelCounted.getHistogram() + urqmdCounted.getHistogram();
- auto const hists = sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
+ // auto const hists = heModelCounted.getHistogram() +
+ // urqmdCounted.getHistogram();
+ auto const hists =
+ sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
save_hist(hists.labHist(), labHist_file, true);
save_hist(hists.CMSHist(), cMSHist_file, true);
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 0c8e4c448..5fde7d39e 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -6,41 +6,41 @@
* the license.
*/
-#include <corsika/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/process/SwitchProcessSequence.hpp>
-#include <corsika/framework/process/InteractionCounter.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/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/geometry/PhysicalGeometry.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/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/SubWriter.hpp>
+#include <corsika/output/OutputManager.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
#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/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
-#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
-#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/Sibyll.hpp>
-#include <corsika/modules/PROPOSAL.hpp>
+#include <corsika/modules/Sophia.hpp>
+#include <corsika/modules/TrackWriter.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
@@ -68,6 +68,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("sibyll");
+ RNGManager<>::getInstance().registerRandomStream("sophia");
if (seed == 0) {
std::random_device rd;
seed = rd();
@@ -79,7 +80,7 @@ void registerRandomStreams(int seed) {
template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
-int main(int argc, char** argv) {
+int main(int argc, char **argv) {
logging::set_level(logging::level::warn);
@@ -91,15 +92,18 @@ int main(int argc, char** argv) {
feenableexcept(FE_INVALID);
int seed = 0;
- if (argc >= 3) { seed = std::stoi(std::string(argv[2])); }
+ if (argc >= 3) {
+ seed = std::stoi(std::string(argv[2]));
+ }
// initialize random number sequence(s)
registerRandomStreams(seed);
// setup environment, geometry
- using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+ using EnvironmentInterface =
+ IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env;
- CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+ CoordinateSystemPtr const &rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
// build a Linsley US Standard atmosphere into `env`
@@ -144,9 +148,11 @@ int main(int argc, char** argv) {
static_pow<2>(injectionHeight));
Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
Point const injectionPos =
- showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+ showerCore +
+ DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
- std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
+ std::cout << "point of injection: " << injectionPos.getCoordinates()
+ << std::endl;
stack.addParticle(std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
@@ -155,8 +161,8 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("shower axis length: {} ",
(showerCore - injectionPos).getNorm() * 1.02);
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
- false, 1000};
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02,
+ env, false, 1000};
OutputManager output("em_shower_outputs");
@@ -166,12 +172,15 @@ int main(int argc, char** argv) {
// setup processes, decays and interactions
- ParticleCut<SubWriter<decltype(dEdX)>> cut(2_MeV, 2_MeV, 100_GeV, 100_GeV, true, dEdX);
+ 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;
- corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
- heThresholdNN);
- corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
+ corsika::proposal::Interaction emCascade(
+ env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
+ corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(
+ env, dEdX);
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
// NOT possible right now, due to interface differenc in PROPOSAL
@@ -190,7 +199,8 @@ int main(int argc, char** argv) {
obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
output.add("particles", observationLevel);
- auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, observationLevel);
+ auto sequence =
+ make_sequence(emCascade, emContinuous, longprof, cut, observationLevel);
// define air shower object, run simulation
setup::Tracking tracking;
@@ -202,12 +212,12 @@ int main(int argc, char** argv) {
EAS.run();
- HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
+ HEPEnergyType const Efinal =
+ dEdX.getEnergyLost() + observationLevel.getEnergyGround();
- CORSIKA_LOG_INFO(
- "total energy budget (GeV): {}, "
- "relative difference (%): {}",
- Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
+ CORSIKA_LOG_INFO("total energy budget (GeV): {}, "
+ "relative difference (%): {}",
+ Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
output.endOfLibrary();
}
--
GitLab