made setup of tests mor robust and flexible

Documentation/Examples/hybrid_MC.cc

@@ -229,8 +229,8 @@ int main(int argc, char** argv) {
   process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
 
   Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
-  process::observation_plane::ObservationPlane observationLevel(obsPlane,
-                                                                "particles.dat");
+  process::observation_plane::ObservationPlane observationLevel(
+      obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
 
   process::UrQMD::UrQMD urqmd;
   process::interaction_counter::InteractionCounter urqmdCounted{urqmd};

Documentation/Examples/vertical_EAS.cc

@@ -17,6 +17,7 @@
 #include <corsika/environment/FlatExponential.h>
 #include <corsika/environment/HomogeneousMedium.h>
 #include <corsika/environment/IMagneticFieldModel.h>
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
 #include <corsika/environment/NuclearComposition.h>
 #include <corsika/environment/ShowerAxis.h>
 #include <corsika/environment/SlidingPlanarExponential.h>
@@ -34,15 +35,12 @@
 #include <corsika/process/particle_cut/ParticleCut.h>
 #include <corsika/process/proposal/ContinuousProcess.h>
 #include <corsika/process/proposal/Interaction.h>
-#include <corsika/process/track_writer/TrackWriter.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/tracking_line/TrackingLine.h>
 #include <corsika/process/urqmd/UrQMD.h>
-#include <corsika/process/proposal/ContinuousProcess.h>
-#include <corsika/process/proposal/Interaction.h>
 #include <corsika/random/RNGManager.h>
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
@@ -191,12 +189,9 @@ int main(int argc, char** argv) {
 
   // setup processes, decays and interactions
 
-  PROPOSAL::InterpolationDef::path_to_tables = "~/.local/share/PROPOSAL/tables/";
-  PROPOSAL::InterpolationDef::path_to_tables_readonly = "~/.local/share/PROPOSAL/tables/";
-
-  process::particle_cut::ParticleCut cut{3_GeV, true, true};
-  process::proposal::Interaction proposal(env, cut);
-  process::proposal::ContinuousProcess em_continuous(env, cut);
+  process::particle_cut::ParticleCut cut{60_GeV, true, true};
+  process::proposal::Interaction proposal(env, cut.GetECut());
+  process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
   process::interaction_counter::InteractionCounter proposalCounted(proposal);
 
   process::sibyll::Interaction sibyll;

Framework/Cascade/Cascade.h

@@ -285,15 +285,6 @@ namespace corsika::cascade {
       assert(currentLogicalNode != &*environment_.GetUniverse() ||
              environment_.GetUniverse()->HasModelProperties());
 
-      // convert next_step from grammage to length
-      LengthType const distance_interact =
-          currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
-                                                                         next_interact);
-
-      // determine the maximum geometric step length from continuous processes
-      LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step);
-      C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
-
       // determine combined total inverse decay time
       InverseTimeType const total_inv_lifetime =
           process_sequence_.GetInverseLifetime(vParticle);
@@ -321,7 +312,7 @@ namespace corsika::cascade {
       
       // determine the maximum geometric step length
       LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, stepWithoutB);
-      std::cout << "distance_max=" << distance_max << std::endl;
+      C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
 
       // take minimum of geometry, interaction, decay for next step
       auto min_distance = std::min(

Framework/Cascade/testCascade.cc

@@ -35,17 +35,24 @@ using namespace corsika::geometry;
 #include <limits>
 using namespace std;
 
+/*
+  The dummy env must support GetMagneticField(), and a density model
+
+ */
 auto MakeDummyEnv() {
   TestEnvironmentType env; // dummy environment
   auto& universe = *(env.GetUniverse());
+  const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
 
   auto theMedium = TestEnvironmentType::CreateNode<Sphere>(
       Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
-      100_km * std::numeric_limits<double>::infinity());
+      1_m * std::numeric_limits<double>::infinity());
+
+  using MyHomogeneousModel = environment::UniformMagneticField<
+      environment::HomogeneousMedium<TestEnvironmentInterface>>;
 
-  using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
   theMedium->SetModelProperties<MyHomogeneousModel>(
-      1_g / (1_cm * 1_cm * 1_cm),
+      geometry::Vector(cs, 0_T, 0_T, 0_T), 1_g / (1_cm * 1_cm * 1_cm),
       environment::NuclearComposition(
           std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
 
@@ -73,16 +80,12 @@ public:
     fCalls++;
     auto const projectile = view.GetProjectile();
     const HEPEnergyType E = projectile.GetEnergy();
-    view.AddSecondary(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            projectile.GetPID(), E / 2, projectile.GetMomentum(),
-            projectile.GetPosition(), projectile.GetTime()});
-    view.AddSecondary(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            projectile.GetPID(), E / 2, projectile.GetMomentum(),
-            projectile.GetPosition(), projectile.GetTime()});
+    view.AddSecondary(std::make_tuple(projectile.GetPID(), E / 2,
+                                      projectile.GetMomentum(), projectile.GetPosition(),
+                                      projectile.GetTime()));
+    view.AddSecondary(std::make_tuple(projectile.GetPID(), E / 2,
+                                      projectile.GetMomentum(), projectile.GetPosition(),
+                                      projectile.GetTime()));
     return EProcessReturn::eInteracted;
   }
 
@@ -141,12 +144,13 @@ TEST_CASE("Cascade", "[Cascade]") {
   auto sequence = process::sequence(nullModel, stackInspect, split, cut);
   TestCascadeStack stack;
   stack.Clear();
-  stack.AddParticle(
-      std::tuple<particles::Code, units::si::HEPEnergyType,
-                 corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-          particles::Code::Electron, E0,
-          corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
-          Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
+  stack.AddParticle(std::make_tuple(
+      particles::Code::Electron, E0,
+      corsika::stack::MomentumVector(
+          rootCS, {0_GeV, 0_GeV,
+                   -sqrt(E0 * E0 - units::static_pow<2>(
+                                       particles::GetMass(particles::Code::Electron)))}),
+      Point(rootCS, {0_m, 0_m, 10_km}), 0_ns));
 
   cascade::Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack,
                    TestCascadeStackView>

Framework/Cascade/testCascade.h

@@ -9,14 +9,15 @@
 #pragma once
 
 #include <corsika/environment/Environment.h>
+#include <corsika/environment/IMagneticFieldModel.h>
 
 #include <corsika/stack/CombinedStack.h>
 #include <corsika/stack/SecondaryView.h>
 #include <corsika/stack/node/GeometryNodeStackExtension.h>
 #include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
 
-using TestEnvironmentType =
-    corsika::environment::Environment<corsika::environment::IMediumModel>;
+using TestEnvironmentInterface = corsika::environment::IMagneticFieldModel<corsika::environment::IMediumModel>;
+using TestEnvironmentType = corsika::environment::Environment<TestEnvironmentInterface>;
 
 template <typename T>
 using SetupGeometryDataInterface =

Processes/ObservationPlane/ObservationPlane.cc

@@ -22,7 +22,7 @@ ObservationPlane::ObservationPlane(
     , outputStream_(filename)
     , deleteOnHit_(deleteOnHit)
     , energy_ground_(0_GeV)
-    , count_ground_(0) {
+    , count_ground_(0)
     , xAxis_(x_axis.normalized())
     , yAxis_(obsPlane.GetNormal().cross(xAxis_)) {
   outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" << std::endl;

Processes/ObservationPlane/testObservationPlane.cc

@@ -19,6 +19,10 @@
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/units/PhysicalUnits.h>
 
+//#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
+//#include <corsika/setup/SetupTrajectory.h>
+
 using namespace corsika::units::si;
 using namespace corsika::process::observation_plane;
 using namespace corsika;
@@ -27,63 +31,51 @@ using namespace corsika::particles;
 
 TEST_CASE("ContinuousProcess interface", "[proccesses][observation_plane]") {
 
-  auto const& rootCS = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+  auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
+  auto const& cs = *csPtr;
+  [[maybe_unused]] auto const& env_dummy = env;
+  [[maybe_unused]] auto const& node_dummy = nodePtr;
 
   /*
-    Test with downward going 1_GeV neutrino, starting at 0,1_m,10m
+    Test with downward going 1_GeV neutrino, starting at 0, 1_m, 10m
 
     ObservationPlane has origin at 0,0,0
    */
 
-  Point const start(rootCS, {0_m, 1_m, 10_m});
-  Vector<units::si::SpeedType::dimension_type> vec(rootCS, 0_m / second, 0_m / second,
+  auto [stack, viewPtr] =
+        setup::testing::setupStack(particles::Code::NuE, 0, 0, 1_GeV, nodePtr, cs);
+  [[maybe_unused]] setup::StackView& view = *viewPtr;
+  auto particle  = stack->GetNextParticle();
+  
+  Point const start(cs, {0_m, 1_m, 10_m});
+  Vector<units::si::SpeedType::dimension_type> vec(cs, 0_m / second, 0_m / second,
                                                    -units::constants::c);
   Line line(start, vec);
   Trajectory<Line> track(line, 12_m / units::constants::c);
 
-  // setup particle stack, and add primary particle
-  setup::Stack stack;
-  stack.Clear();
-  {
-    auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
-      return sqrt((Elab - m) * (Elab + m));
-    };
-    stack.AddParticle(
-        std::tuple<Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, Point,
-                   units::si::TimeType>{
-            Code::NuMu, 1_GeV,
-            corsika::stack::MomentumVector(
-                rootCS, {0_GeV, 0_GeV, -elab2plab(1_GeV, NuMu::GetMass())}),
-            Point(rootCS, {1_m, 1_m, 10_m}), 0_ns});
-  }
-  auto particle = stack.GetNextParticle();
+  particle.SetPosition(Point(cs, {1_m, 1_m, 10_m})); // moving already along -z
 
   SECTION("horizontal plane") {
 
-    Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
-                         Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
-    ObservationPlane obs(obsPlane, "particles.dat", true);
+    Plane const obsPlane(Point(cs, {0_m, 0_m, 0_m}),
+                         Vector<dimensionless_d>(cs, {0., 0., 1.}));
+    ObservationPlane obs(obsPlane, Vector<dimensionless_d>(cs, {1., 0., 0.}),
+                         "particles.dat", true);
 
     const LengthType length = obs.MaxStepLength(particle, track);
     const process::EProcessReturn ret = obs.DoContinuous(particle, track);
 
     REQUIRE(length / 10_m == Approx(1).margin(1e-4));
     REQUIRE(ret == process::EProcessReturn::eParticleAbsorbed);
-
-    /*
-    SECTION("horizontal plane") {
-      REQUIRE(true); // todo: we have to check content of output file...
-
-    }
-    */
   }
 
   SECTION("inclined plane") {}
 
   SECTION("transparent plane") {
-    Plane const obsPlane(Point(rootCS, {0_m, 0_m, 0_m}),
-                         Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
-    ObservationPlane obs(obsPlane, "particles.dat", false);
+    Plane const obsPlane(Point(cs, {0_m, 0_m, 0_m}),
+                         Vector<dimensionless_d>(cs, {0., 0., 1.}));
+    ObservationPlane obs(obsPlane, Vector<dimensionless_d>(cs, {1., 0., 0.}),
+                         "particles.dat", false);
 
     const LengthType length = obs.MaxStepLength(particle, track);
     const process::EProcessReturn ret = obs.DoContinuous(particle, track);

Processes/TrackingLine/dump_bh.hpp

@@ -75,7 +75,7 @@ void dump_bh_2d(std::ostream& os, T& h)
   // bins-x
   os << "  \"xbins\" : [";
   double lastx = 0;
-  for (unsigned int i=0; i<h.axis(0).size(); ++i) {
+  for (int i=0; i<h.axis(0).size(); ++i) {
     os << h.axis(0).bin(i).lower() << ", ";
     lastx = h.axis(0).bin(i).upper();
   }
@@ -84,7 +84,7 @@ void dump_bh_2d(std::ostream& os, T& h)
   // bins-y
   os << "  \"ybins\" : [";
   double lasty = 0;
-  for (unsigned int i=0; i<h.axis(1).size(); ++i) {
+  for (int i=0; i<h.axis(1).size(); ++i) {
     os << h.axis(1).bin(i).lower() << ", ";
     lasty = h.axis(1).bin(i).upper();
   }

Processes/TrackingLine/testTrackingLine.cc

@@ -30,8 +30,9 @@ using namespace corsika::geometry;
 using namespace std;
 using namespace corsika::units::si;
 
+
 TEST_CASE("TrackingLine") {
-  environment::Environment<environment::Empty> env; // dummy environment
+  environment::Environment<TestMagneticField> env; // dummy environment
   auto const& cs = env.GetCoordinateSystem();
 
   tracking_line::TrackingLine tracking;
@@ -64,7 +65,7 @@ TEST_CASE("TrackingLine") {
 
     auto const radius = 20_m;
 
-    auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>(
+    auto theMedium = environment::Environment<TestMagneticField>::CreateNode<Sphere>(
         Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
     auto const* theMediumPtr = theMedium.get();
     universe.AddChild(std::move(theMedium));
@@ -86,11 +87,15 @@ TEST_CASE("TrackingLine") {
                                                             0_m / second, 1_m / second);
     Line line(origin, v);
 
-    auto const [traj, geomMaxLength, nextVol, magMaxLength, beforeDirection, afterDirection] = tracking.GetTrack(p);
-    [[maybe_unused]] auto dummy_geomMaxLength = geomMaxLength;
-    [[maybe_unused]] auto dummy_nextVol = nextVol;
+    const auto [stepWithoutB, stepWithB, geomMaxLength, magMaxLength, nextVol] = tracking.GetTrack(p);
+    //auto const [traj, geomMaxLength, nextVol, magMaxLength, beforeDirection,
+    //          afterDirection] = tracking.GetTrack(p);
+    [[maybe_unused]] auto& dummy_1 = stepWithB;
+    [[maybe_unused]] auto& dummy_2 = magMaxLength;
+    [[maybe_unused]] auto& dummy_geomMaxLength = geomMaxLength;
+    [[maybe_unused]] auto& dummy_nextVol = nextVol;
 
-    REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
+    REQUIRE((stepWithoutB.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
                 .GetComponents(cs)
                 .norm()
                 .magnitude() == Approx(0).margin(1e-4));

Processes/TrackingLine/testTrackingLineStack.h

@@ -21,8 +21,21 @@
 
 #include <corsika/units/PhysicalUnits.h>
 
+
+class TestMagneticField {
+    using MagneticFieldVector =
+        corsika::geometry::Vector<corsika::units::si::magnetic_flux_density_d>;
+
+public:
+  MagneticFieldVector GetMagneticField(corsika::geometry::Point const& p) const {
+    using namespace corsika::units::si;
+    return MagneticFieldVector(p.GetCoordinateSystem(), 0_T, 0_T, 1_T);
+  }
+};
+
+
 using TestEnvironmentType =
-    corsika::environment::Environment<corsika::environment::Empty>;
+    corsika::environment::Environment<TestMagneticField>;
 
 template <typename T>
 using SetupGeometryDataInterface =

do-copyright.py

@@ -22,7 +22,7 @@ Debug settings are 0: nothing, 1: checking, 2: filesystem
 Debug = 0 
 
 excludeDirs = ["ThirdParty", "git", "build", "install", "PROPOSAL"]
-excludeFiles = ['PhysicalConstants.h','CorsikaFenvOSX.cc', 'sgn.h']
+excludeFiles = ['PhysicalConstants.h','CorsikaFenvOSX.cc', 'sgn.h', 'quartic.h']
 
 extensions = [".cc", ".h", ".test"]