IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 234f6123 authored by Felix Riehn's avatar Felix Riehn Committed by Maximilian Reininghaus
Browse files

use LE photon model in examples

parent d05e531b
No related branches found
No related tags found
1 merge request!465Resolve "SOPHIA for low energy photo-hadronic interaction"
Hide whitespace changes
Inline Side-by-side
corsika/modules/Random.hpp
+ 1
− 0
View file @ 234f6123
......@@ -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>
......
examples/corsika.cpp
+ 80
− 60
View file @ 234f6123
......@@ -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);
......
examples/em_shower.cpp
+ 42
− 32
View file @ 234f6123
......@@ -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();
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum