diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc
index c6647209bb13d4843020c2bce4d0354486f97753..93531a2a09f994c08ee5f926f15552f2fc1393f2 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/Documentation/Examples/hybrid_MC.cc
@@ -22,18 +22,19 @@
#include <corsika/geometry/Sphere.h>
#include <corsika/logging/Logging.h>
#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/SwitchProcessSequence.h>
#include <corsika/process/StackProcess.h>
-#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <corsika/process/on_shell_check/OnShellCheck.h>
+#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
+#include <corsika/process/switch_process/SwitchProcess.h>
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
@@ -73,16 +74,17 @@ void registerRandomStreams(const int seed) {
}
template <typename T>
-using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
+using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
+
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
- C8LOG_INFO("hybrid_MC");
+ C8LOG_INFO("vertical_EAS");
if (argc < 4) {
- std::cerr << "usage: hybrid_MC <A> <Z> <energy/GeV> [seed]" << std::endl;
+ std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed]" << std::endl;
std::cerr << " if no seed is given, a random seed is chosen" << std::endl;
return 1;
}
@@ -98,20 +100,26 @@ int main(int argc, char** argv) {
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
- auto builder = environment::make_layered_spherical_atmosphere_builder<
- setup::EnvironmentInterface,
- MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
- environment::Medium::AirDry1Atm,
- geometry::Vector{rootCS, 0_T, 0_T, 1_T});
+ environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{center};
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
- builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
- builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
- builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
- builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
- builder.addLinearLayer(1e9_cm, 112.8_km);
+ builder.addExponentialLayer<MEnv>(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km,
+ environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<MEnv>(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
+ environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<MEnv>(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
+ environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addExponentialLayer<MEnv>(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
+ environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+ builder.addLinearLayer<MEnv>(1e9_cm, 112.8_km, environment::EMediumType::eAir,
+ geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+
builder.assemble(env);
// setup particle stack, and add primary particle
@@ -223,23 +231,14 @@ int main(int argc, char** argv) {
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
// assemble all processes into an ordered process list
- struct EnergySwitch {
- HEPEnergyType cutE_;
- EnergySwitch(HEPEnergyType cutE)
- : cutE_(cutE) {}
- process::SwitchResult operator()(const setup::Stack::ParticleType& p) {
- if (p.GetEnergy() < cutE_)
- return process::SwitchResult::First;
- else
- return process::SwitchResult::Second;
- }
- };
- auto hadronSequence =
- process::select(urqmdCounted, process::sequence(sibyllNucCounted, sibyllCounted),
- EnergySwitch(55_GeV));
- auto decaySequence = process::sequence(decayPythia, decaySibyll);
- auto sequence = process::sequence(hadronSequence, reset_particle_mass, decaySequence,
- eLoss, cut, conex, longprof, observationLevel);
+
+ auto sibyllSequence = sibyllNucCounted << sibyllCounted;
+ process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
+ 55_GeV);
+ auto decaySequence = decayPythia << decaySibyll;
+
+ auto sequence = switchProcess << reset_particle_mass << decaySequence
+ << eLoss << cut << conex << longprof << observationLevel;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
@@ -252,21 +251,22 @@ int main(int argc, char** argv) {
EAS.Run();
cut.ShowResults();
- eLoss.showResults();
+ //eLoss.ShowResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
- cut.GetEmEnergy() + eLoss.energyLost() +
+ cut.GetEmEnergy() + //eLoss.GetEnergyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
- eLoss.reset();
+ //eLoss.Reset();
+ //conex.addParticle(0, E0, 0_eV, emPosition, momentum.normalized(), 0_s);
conex.SolveCE();
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
- urqmdCounted.GetHistogram();
+ urqmdCounted.GetHistogram();
hists.saveLab("inthist_lab.txt");
hists.saveCMS("inthist_cms.txt");
diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc
index 1078f447734d2e925b4b700029907bd6f3ffb324..84dad94cfea6e61825e98cedcf4bc9081cc62f38 100644
--- a/Processes/QGSJetII/testQGSJetII.cc
+++ b/Processes/QGSJetII/testQGSJetII.cc
@@ -137,7 +137,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
Interaction model;
- [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(*secViewPtr);
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
CHECK(length / (1_g / square(1_cm)) == Approx(93.04).margin(0.1));
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 68e9740987fac39c9ef80df4245e5d8371944a12..007be947e5d4c8103da17cabee6b39700120463a 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -115,8 +115,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
Interaction model;
- [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(*view);
- auto const pSum = sumMomentum(*view, cs);
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
+ auto const pSum = sumMomentum(view, cs);
/*
Interactions between hadrons (h) and nuclei (A) in Sibyll are treated in the
@@ -195,7 +195,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
Interaction hmodel;
NuclearInteraction model(hmodel, *env);
- [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(*view);
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
CHECK(length / 1_g * 1_cm * 1_cm == Approx(44.2).margin(.1));
CHECK(view.getSize() == 11);