fixed failing tests

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");

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));

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);