diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc
index 515d9ac95237473cd7a5010ed9beed9561af5f3f..33dec55b01d5a4d9d122abda6705f1236b6d878b 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/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};
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 53089979ab1ec77efa55d02093a7f2e7d2cc4e92..5cf491d29abc901763189aeab63c1adfe9875734 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/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;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 2fd8de8da5a8e97460ede156aa24d5eeb1e263a7..4596b238a53f5bea8c62616937f11186eec111a9 100644
--- a/Framework/Cascade/Cascade.h
+++ b/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(
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 7cf7b9224327fe3422a86b32e3515f1a44f4af48..f0514598e8d40195a9d356599bf62e29056a4a25 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/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>
diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h
index 55740c5601694b2c982c54607bcd4182f519171f..8c78f740c1eb0338ee05eb36766c736bf1b5cf7b 100644
--- a/Framework/Cascade/testCascade.h
+++ b/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 =
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index c77ce9b6c0528abfe6008545611ba68a78fe0963..27bed61acc8582b54f411051137ad748706128c4 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/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;
diff --git a/Processes/ObservationPlane/testObservationPlane.cc b/Processes/ObservationPlane/testObservationPlane.cc
index 7f3f183a6def750818aa1eb05be4bec81937d587..42837e8b2504d27ceead6eb28ce4ca70b7a61bdc 100644
--- a/Processes/ObservationPlane/testObservationPlane.cc
+++ b/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);
diff --git a/Processes/TrackingLine/dump_bh.hpp b/Processes/TrackingLine/dump_bh.hpp
index 01e0dac20726b2d2b4787fad3bb3a9f0fff27b09..97f8460d7fa38f99d8fece4c3a227cb43fb38d22 100644
--- a/Processes/TrackingLine/dump_bh.hpp
+++ b/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();
}
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index a72e870b762fc10ed61cf58604c5eb535490b2f3..7634758e90b959c46d16700f3ef95d6c6c9182ea 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/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));
diff --git a/Processes/TrackingLine/testTrackingLineStack.h b/Processes/TrackingLine/testTrackingLineStack.h
index 5c714048a5eb1e3978f39f6815c80eaaebcdc8a0..16b07e0cba8d944472e46660319dc4aed847c62b 100644
--- a/Processes/TrackingLine/testTrackingLineStack.h
+++ b/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 =
diff --git a/do-copyright.py b/do-copyright.py
index 4c2b29cc7a2c79418a5a8517018806e6c0aa1a9e..716b2c4bc811fb431aa441c40368abbe85ef9b1f 100755
--- a/do-copyright.py
+++ b/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"]