diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 4b9c6d55fc9e91202568b46ad79194fecda4c590..8cd42babe67261c0f17e52805e67bfceb8a09ab7 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -7,7 +7,7 @@ set (CMAKE_VERBOSE_MAKEFILE OFF) # this can be changed with `make VERBOSE=1`
find_package (corsika CONFIG REQUIRED)
# this is only for CORSIKA_REGISTER_EXAMPLE
-include ("${CMAKE_CURRENT_SOURCE_DIR}/CMakeHelper.cmake")
+include ("${CMAKE_CURRENT_SOURCE_DIR}/CMakeHelper.cmake")
add_executable (helix_example helix_example.cpp)
target_link_libraries (helix_example CORSIKA8::CORSIKA8)
@@ -35,7 +35,7 @@ CORSIKA_REGISTER_EXAMPLE (cascade_proton_example)
add_executable (vertical_EAS vertical_EAS.cpp)
target_link_libraries (vertical_EAS CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
+CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000. 1)
add_executable (stopping_power stopping_power.cpp)
target_link_libraries (stopping_power CORSIKA8::CORSIKA8)
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index bdad125b05e7557a4db27da54ee2eecae27d7592..9f816ce19975a7be024ec08bfd4a8f82308e894f 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -91,6 +91,11 @@ void registerRandomStreams(const int 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.output directory
+// 6.seed (0 by default to generate random values for all)
+
int main(int argc, char** argv) {
corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
@@ -98,7 +103,7 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("vertical_EAS");
- if (argc < 4) {
+ if (argc < 5) {
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"
@@ -107,8 +112,16 @@ int main(int argc, char** argv) {
}
feenableexcept(FE_INVALID);
+ string output_dir = "";
int seed = 0;
- if (argc > 4) seed = std::stoi(std::string(argv[4]));
+ int number_showers = std::stoi(std::string(argv[4]));
+ if (argc > 5) {
+ output_dir = argv[5];
+ if (output_dir.back() != '/' && output_dir != "") output_dir += '/';
+ }
+
+ if (argc > 6) { seed = std::stoi(std::string(argv[6])); }
+
// initialize random number sequence(s)
registerRandomStreams(seed);
@@ -150,173 +163,190 @@ int main(int argc, char** argv) {
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m}))));
- // setup particle stack, and add primary particle
- setup::Stack stack;
- stack.clear();
- unsigned short const A = std::stoi(std::string(argv[1]));
- Code beamCode;
- HEPEnergyType mass;
- unsigned short Z = 0;
- if (A > 0) {
- beamCode = Code::Nucleus;
- Z = std::stoi(std::string(argv[2]));
- mass = get_nucleus_mass(A, Z);
- } else {
- int pdg = std::stoi(std::string(argv[2]));
- beamCode = convert_from_PDG(PDGCode(pdg));
- mass = get_mass(beamCode);
- }
- HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
- double theta = 0.;
- auto const thetaRad = theta / 180. * M_PI;
-
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt((Elab - m) * (Elab + m));
- };
- HEPMomentumType P0 = elab2plab(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 observationHeight = 0_km + builder.getEarthRadius();
- auto const injectionHeight = 111.75_km + builder.getEarthRadius();
- 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), 0, cos(thetaRad)}} * t;
-
- std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
-
- if (A > 1) {
- stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
-
- } else {
- if (A == 1) {
- if (Z == 1) {
- stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
- } else if (Z == 0) {
- stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
- } else {
- std::cerr << "illegal parameters" << std::endl;
- return EXIT_FAILURE;
- }
+ for (int i_shower = 1; i_shower < number_showers + 1; i_shower++) {
+
+ // directory for outputs
+ string labHist_dir =
+ output_dir + "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
+ string cMSHist_dir =
+ output_dir + "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
+ string longprof_dir =
+ output_dir + "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
+ string tracks_dir = output_dir + "tracks_" + to_string(i_shower) + ".dat";
+ string particles_dir = output_dir + "particles_" + to_string(i_shower) + ".dat";
+
+ std::cout << std::endl;
+ std::cout << "Shower " << i_shower << "/" << number_showers << std::endl;
+
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.clear();
+ unsigned short const A = std::stoi(std::string(argv[1]));
+ Code beamCode;
+ HEPEnergyType mass;
+ unsigned short Z = 0;
+ if (A > 0) {
+ beamCode = Code::Nucleus;
+ Z = std::stoi(std::string(argv[2]));
+ mass = get_nucleus_mass(A, Z);
} else {
- stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
+ int pdg = std::stoi(std::string(argv[2]));
+ beamCode = convert_from_PDG(PDGCode(pdg));
+ mass = get_mass(beamCode);
}
- }
+ HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(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 observationHeight = 0_km + builder.getEarthRadius();
+ auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+ 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), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
+
+ if (A > 1) {
+ stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
- // 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};
-
- // setup processes, decays and interactions
-
- // corsika::qgsjetII::Interaction qgsjet;
- corsika::sibyll::Interaction sibyll;
- InteractionCounter sibyllCounted(sibyll);
-
- corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
- InteractionCounter sibyllNucCounted(sibyllNuc);
-
- corsika::pythia8::Decay decayPythia;
-
- // use sibyll decay routine for decays of particles unknown to pythia
- corsika::sibyll::Decay decaySibyll{{
- Code::N1440Plus,
- Code::N1440MinusBar,
- Code::N1440_0,
- Code::N1440_0Bar,
- Code::N1710Plus,
- Code::N1710MinusBar,
- Code::N1710_0,
- Code::N1710_0Bar,
-
- Code::Pi1300Plus,
- Code::Pi1300Minus,
- Code::Pi1300_0,
-
- Code::KStar0_1430_0,
- Code::KStar0_1430_0Bar,
- Code::KStar0_1430_Plus,
- Code::KStar0_1430_MinusBar,
- }};
-
- decaySibyll.printDecayConfig();
-
- ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
- corsika::proposal::Interaction emCascade(env);
- corsika::proposal::ContinuousProcess emContinuous(env);
- InteractionCounter emCascadeCounted(emCascade);
-
- OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
- TrackWriter trackWriter("tracks.dat");
-
- LongitudinalProfile longprof{showerAxis};
-
- Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
- "particles.dat");
-
- corsika::urqmd::UrQMD urqmd;
- InteractionCounter urqmdCounted{urqmd};
- StackInspector<setup::Stack> stackInspect(50000, false, E0);
-
- // assemble all processes into an ordered process list
- struct EnergySwitch {
- HEPEnergyType cutE_;
- EnergySwitch(HEPEnergyType cutE)
- : cutE_(cutE) {}
- SwitchResult operator()(const Particle& p) {
- if (p.getEnergy() < cutE_)
- return SwitchResult::First;
- else
- return SwitchResult::Second;
+ } else {
+ if (A == 1) {
+ if (Z == 1) {
+ stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
+ } else if (Z == 0) {
+ stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
+ } else {
+ std::cerr << "illegal parameters" << std::endl;
+ return EXIT_FAILURE;
+ }
+ } else {
+ stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
+ }
}
- };
- auto hadronSequence = make_select(
- urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted), EnergySwitch(55_GeV));
- auto decaySequence = make_sequence(decayPythia, decaySibyll);
- auto sequence =
- make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
- emContinuous, cut, trackWriter, observationLevel, longprof);
-
- // define air shower object, run simulation
- setup::Tracking tracking;
- Cascade EAS(env, tracking, sequence, stack);
-
- // to fix the point of first interaction, uncomment the following two lines:
- // EAS.forceInteraction();
-
- EAS.run();
-
- cut.showResults();
- emContinuous.showResults();
- observationLevel.showResults();
- const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
- cut.getEmEnergy() + emContinuous.getEnergyLost() +
- observationLevel.getEnergyGround();
- cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
- << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
- observationLevel.reset();
- cut.reset();
- emContinuous.reset();
-
- auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
- urqmdCounted.getHistogram();
-
- save_hist(hists.labHist(), "inthist_lab_verticalEAS.npz", true);
- save_hist(hists.CMSHist(), "inthist_cms_verticalEAS.npz", true);
- longprof.save("longprof_verticalEAS.txt");
+
+ // 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};
+
+ // setup processes, decays and interactions
+
+ // corsika::qgsjetII::Interaction qgsjet;
+ corsika::sibyll::Interaction sibyll;
+ InteractionCounter sibyllCounted(sibyll);
+
+ corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+ InteractionCounter sibyllNucCounted(sibyllNuc);
+
+ corsika::pythia8::Decay decayPythia;
+
+ // use sibyll decay routine for decays of particles unknown to pythia
+ corsika::sibyll::Decay decaySibyll{{
+ Code::N1440Plus,
+ Code::N1440MinusBar,
+ Code::N1440_0,
+ Code::N1440_0Bar,
+ Code::N1710Plus,
+ Code::N1710MinusBar,
+ Code::N1710_0,
+ Code::N1710_0Bar,
+
+ Code::Pi1300Plus,
+ Code::Pi1300Minus,
+ Code::Pi1300_0,
+
+ Code::KStar0_1430_0,
+ Code::KStar0_1430_0Bar,
+ Code::KStar0_1430_Plus,
+ Code::KStar0_1430_MinusBar,
+ }};
+
+ decaySibyll.printDecayConfig();
+
+ ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
+ corsika::proposal::Interaction emCascade(env);
+ corsika::proposal::ContinuousProcess emContinuous(env);
+ InteractionCounter emCascadeCounted(emCascade);
+
+ OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
+ TrackWriter trackWriter(tracks_dir);
+
+ LongitudinalProfile longprof{showerAxis};
+
+ Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+ ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
+ particles_dir);
+
+ corsika::urqmd::UrQMD urqmd;
+ InteractionCounter urqmdCounted{urqmd};
+ StackInspector<setup::Stack> stackInspect(50000, false, E0);
+
+ // assemble all processes into an ordered process list
+ struct EnergySwitch {
+ HEPEnergyType cutE_;
+ EnergySwitch(HEPEnergyType cutE)
+ : cutE_(cutE) {}
+ SwitchResult operator()(const Particle& p) {
+ if (p.getEnergy() < cutE_)
+ return SwitchResult::First;
+ else
+ return SwitchResult::Second;
+ }
+ };
+ auto hadronSequence =
+ make_select(urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted),
+ EnergySwitch(55_GeV));
+ auto decaySequence = make_sequence(decayPythia, decaySibyll);
+ auto sequence =
+ make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
+ emContinuous, cut, trackWriter, observationLevel, longprof);
+
+ // define air shower object, run simulation
+ setup::Tracking tracking;
+ Cascade EAS(env, tracking, sequence, stack);
+
+ // to fix the point of first interaction, uncomment the following two lines:
+ // EAS.forceInteraction();
+
+ EAS.run();
+
+ cut.showResults();
+ emContinuous.showResults();
+ observationLevel.showResults();
+ const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
+ cut.getEmEnergy() + emContinuous.getEnergyLost() +
+ observationLevel.getEnergyGround();
+ cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+ observationLevel.reset();
+ cut.reset();
+ emContinuous.reset();
+
+ auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
+ urqmdCounted.getHistogram();
+
+ save_hist(hists.labHist(), labHist_dir, true);
+ save_hist(hists.CMSHist(), cMSHist_dir, true);
+ longprof.save(longprof_dir);
+ }
}