diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 91308390ba4f9b5b8e87bb0d94e79646e9a6f101..165f0dd60973bfcf03c989594ee322cd75b645d3 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -1,3 +1,4 @@
+
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -132,9 +133,9 @@ namespace corsika::proposal {
std::max(energy * 0.9, // either 10% relative loss max., or
get_kinetic_energy_threshold(
code) // energy thresholds globally defined for individual particles
- * 0.9999 // need to go slightly below global e-cut to assure removal
- // in ParticleCut. This does not matter since at cut-time
- // the entire energy is removed.
+ * 0.9999 // need to go slightly below global e-cut to assure removal in
+ // ParticleCut. This does not matter since at cut-time the
+ // entire energy is removed.
);
// solving the track integral for giving energy lim
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index 019d31270abeb1da6fe2a99121081363762dcac6..44e7bc5a158dfaca830abef68540473302b5acf7 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -98,7 +98,8 @@ namespace corsika::proposal {
auto vec = QuantityVector(vecProposal.GetX() * E, vecProposal.GetY() * E,
vecProposal.GetZ() * E);
auto p = MomentumVector(labCS, vec);
- auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
+ auto sec_code =
+ convert_from_PDG(static_cast<PDGCode>(s.type));
view.addSecondary(
std::make_tuple(sec_code, p, projectile.getPosition(), projectile.getTime()));
}
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 3240a632cb15c390454e99e0515c52f7945336e7..2d20738b31793611db793fa90bf1cc1bca1b1dd5 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -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. 1)
+CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
add_executable (stopping_power stopping_power.cpp)
target_link_libraries (stopping_power CORSIKA8::CORSIKA8)
@@ -56,3 +56,7 @@ CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS "100.")
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 (corsika corsika.cpp)
+target_link_libraries (corsika CORSIKA8::CORSIKA8)
+CORSIKA_REGISTER_EXAMPLE (corsika RUN_OPTIONS 4 2 10000. 1)
diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp
index 3ce8df868981a77af9d4df509f414264a2c3138d..d9522ce15535b4edf56c3f99c4f1857862d1ceff 100644
--- a/examples/boundary_example.cpp
+++ b/examples/boundary_example.cpp
@@ -26,7 +26,6 @@
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
-#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp>
@@ -38,9 +37,6 @@
executable. If you include the header below multiple times and
link this togehter, it will fail.
*/
-#include <corsika/modules/sibyll/Random.hpp>
-#include <corsika/modules/urqmd/Random.hpp>
-
#include <iostream>
#include <limits>
#include <typeinfo>
@@ -83,7 +79,6 @@ private:
int main() {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("boundary_example");
@@ -120,10 +115,6 @@ int main() {
// setup processes, decays and interactions
setup::Tracking tracking;
- RNGManager::getInstance().registerRandomStream("sibyll");
- corsika::sibyll::Interaction sibyll;
- corsika::sibyll::Decay decay;
-
ParticleCut cut(50_GeV, true, true);
TrackWriter trackWriter;
@@ -132,7 +123,7 @@ int main() {
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
- auto sequence = make_sequence(sibyll, decay, cut, boundaryCrossing, trackWriter);
+ auto sequence = make_sequence(cut, boundaryCrossing, trackWriter);
// setup particle stack, and add primary particles
setup::Stack stack;
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index 71be6f4f9edc102cbc5d5d190955abe757d07191..e6ab7d76fff01d3b0081e7e5a1bb784a01fb59a3 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -59,7 +59,6 @@ using namespace std;
int main() {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::cout << "cascade_proton_example" << std::endl;
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 12400723bf9c29550f8c424c9b55d117c8b7a82c..2d67cf52d2141ab7be08fbddfae179c04fb95536 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -139,7 +139,8 @@ int main(int argc, char** argv) {
<< std::endl;
OutputManager output("em_shower_outputs");
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
+ false, 1000.};
// setup processes, decays and interactions
@@ -189,4 +190,4 @@ int main(int argc, char** argv) {
longprof.save("longprof_emShower.txt");
output.endOfLibrary();
-}
+}
\ No newline at end of file
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index 25603ad8af4c3f797d5cde7ae47d0190cf960a29..3a4f6757708a4a2ff081d8240b2bdf537a90a6f6 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -173,7 +173,8 @@ int main(int argc, char** argv) {
<< std::endl;
OutputManager output("hybrid_MC_outputs");
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env, true,
+ 1000};
// setup processes, decays and interactions
diff --git a/examples/stack_example.cpp b/examples/stack_example.cpp
index 41beda4f9c5c6f220411f675fe8e3d8c76e7eda3..da97201afad3001f6a4b9e5cec6f2701df227091 100644
--- a/examples/stack_example.cpp
+++ b/examples/stack_example.cpp
@@ -44,11 +44,11 @@ void read(VectorStack& s) {
int main() {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
- std::cout << "stack_example" << std::endl;
+ CORSIKA_LOG_INFO("stack_example");
VectorStack s;
fill(s);
read(s);
+ CORSIKA_LOG_INFO("done");
return 0;
}
diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp
index 4c9e43efc1883e933433aa83240896442491bb9b..e971a279d27edd940098b228d335efdca9d27594 100644
--- a/examples/stopping_power.cpp
+++ b/examples/stopping_power.cpp
@@ -30,7 +30,6 @@ using namespace std;
int main() {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("stopping_power");
@@ -48,15 +47,16 @@ int main() {
rootCS, 0_m, 0_m,
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
- ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env};
+ ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env, false,
+ 100};
BetheBlochPDG eLoss{showerAxis};
setup::Stack stack;
std::ofstream file("dEdX.dat");
- file << "# beta*gamma, dE/dX / eV/(g/cm²)" << std::endl;
+ file << "# beta*gamma, dE/dX / MeV/(g/cm²)" << std::endl;
- for (HEPEnergyType E0 = 300_MeV; E0 < 1_PeV; E0 *= 1.05) {
+ for (HEPEnergyType E0 = 200_MeV; E0 < 1_PeV; E0 *= 1.05) {
stack.clear();
const Code beamCode = Code::MuPlus;
const HEPMassType mass = get_mass(beamCode);
@@ -72,16 +72,14 @@ int main() {
-ptot * cos(theta));
};
auto const [px, py, pz] =
- momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
+ momentumComponents(theta / 180. * constants::pi, phi / 180. * constants::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, plab, injectionPos, 0_ns));
auto const p = stack.getNextParticle();
HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm));
- file << P0 / mass << "\t" << -dE / 1_eV << std::endl;
+ file << P0 / mass << "\t" << -dE / 1_MeV << std::endl;
}
+ CORSIKA_LOG_INFO("finished writing dEdX.dat");
}
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 21354adcfc31d8fec2ec28381a7aeb4ceb39e39e..2fc5ad8f0d855c5ed8e31c9fb38006d801c7bb21 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -43,15 +43,12 @@
#include <corsika/modules/BetheBlochPDG.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/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.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>
@@ -71,7 +68,6 @@
*/
#include <corsika/modules/sibyll/Random.hpp>
#include <corsika/modules/urqmd/Random.hpp>
-#include <corsika/modules/qgsjetII/Random.hpp>
using namespace corsika;
using namespace std;
@@ -80,7 +76,6 @@ using Particle = setup::Stack::particle_type;
void registerRandomStreams(const int seed) {
RNGManager::getInstance().registerRandomStream("cascade");
- RNGManager::getInstance().registerRandomStream("qgsjet");
RNGManager::getInstance().registerRandomStream("sibyll");
RNGManager::getInstance().registerRandomStream("pythia");
RNGManager::getInstance().registerRandomStream("urqmd");
@@ -105,8 +100,8 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("vertical_EAS");
- if (argc < 5) {
- std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> <Nevt> [seed] \n"
+ 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;
@@ -115,9 +110,8 @@ int main(int argc, char** argv) {
feenableexcept(FE_INVALID);
int seed = 0;
- int number_showers = std::stoi(std::string(argv[4]));
- if (argc > 5) { seed = std::stoi(std::string(argv[5])); }
+ if (argc > 4) { seed = std::stoi(std::string(argv[4])); }
// initialize random number sequence(s)
registerRandomStreams(seed);
@@ -162,8 +156,8 @@ int main(int argc, char** argv) {
HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
double phi = 180.;
- auto const thetaRad = theta / 180. * M_PI;
- auto const phiRad = phi / 180. * M_PI;
+ auto const thetaRad = theta / 180. * constants::pi;
+ auto const phiRad = phi / 180. * constants::pi;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
@@ -200,14 +194,14 @@ int main(int argc, char** argv) {
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
<< std::endl;
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
+ 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");
// setup processes, decays and interactions
- // corsika::qgsjetII::Interaction qgsjet;
corsika::sibyll::Interaction sibyll;
InteractionCounter sibyllCounted(sibyll);
@@ -240,15 +234,6 @@ int main(int argc, char** argv) {
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);
-
- LongitudinalProfile longprof{showerAxis};
-
- Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
@@ -265,170 +250,71 @@ int main(int argc, char** argv) {
make_sequence(sibyllNucCounted, sibyllCounted));
auto decaySequence = make_sequence(decayPythia, decaySibyll);
- for (int i_shower = 1; i_shower < number_showers + 1; i_shower++) {
-
- // directory for outputs
- string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
- string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
- string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
-
- 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 {
- 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 = (injectionHeight - observationHeight) / cos(thetaRad);
- 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, plab, injectionPos, 0_ns, A, Z));
+ // directory for outputs
+ string const labHist_file = "inthist_lab_verticalEAS.npz";
+ string const cMSHist_file = "inthist_cms_verticalEAS.npz";
+ string const longprof_file = "longprof_verticalEAS.txt";
- } else {
- if (A == 1) {
- if (Z == 1) {
- stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
- } else if (Z == 0) {
- stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
- } else {
- std::cerr << "illegal parameters" << std::endl;
- return EXIT_FAILURE;
- }
+ // setup particle stack, and add primary particle
+ setup::Stack stack;
+ stack.clear();
+
+ if (A > 1) {
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
+
+ } else {
+ if (A == 1) {
+ if (Z == 1) {
+ stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
+ } else if (Z == 0) {
+ stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
} else {
- stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+ std::cerr << "illegal parameters" << std::endl;
+ return EXIT_FAILURE;
}
+ } else {
+ stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
}
+ }
- // 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;
+ BetheBlochPDG emContinuous(showerAxis);
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env,
- false};
-
- // 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, true, true};
- // corsika::proposal::Interaction emCascade(env);
- // corsika::proposal::ContinuousProcess emContinuous(env);
- // InteractionCounter emCascadeCounted(emCascade);
- BetheBlochPDG emContinuous(showerAxis);
-
- OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
- TrackWriter trackWriter;
- output.add("tracks", trackWriter); // register TrackWriter
-
- LongitudinalProfile longprof{showerAxis};
-
- Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
- // register the observation plane with the output
- output.add("particles", observationLevel);
-
- auto sequence = make_sequence( // emCascadeCounted,
- stackInspect, hadronSequence, reset_particle_mass, decaySequence,
- // emContinuous,
- BetheBlochPDG(showerAxis), cut, trackWriter, observationLevel, longprof);
-
- // define air shower object, run simulation
- setup::Tracking tracking;
- Cascade EAS(env, tracking, sequence, output, 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_file, true);
- save_hist(hists.CMSHist(), cMSHist_file, true);
- longprof.save(longprof_file);
-
- output.endOfLibrary();
- }
+ TrackWriter trackWriter;
+ output.add("tracks", trackWriter); // register TrackWriter
+
+ LongitudinalProfile longprof{showerAxis};
+
+ Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+ ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
+ // register the observation plane with the output
+ output.add("particles", observationLevel);
+
+ auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, emContinuous,
+ cut, trackWriter, observationLevel, longprof);
+
+ // define air shower object, run simulation
+ setup::Tracking tracking;
+ Cascade EAS(env, tracking, sequence, output, stack);
+ 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_file, true);
+ save_hist(hists.CMSHist(), cMSHist_file, true);
+ longprof.save(longprof_file);
+
+ output.endOfLibrary();
}