diff --git a/.gitignore b/.gitignore
index 8fd6d88b5ed5509ee059aa5605929ee9e78450d3..6924790a75ed639a638e4e0fc65bb0ac53487f3b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@ build/
Framework/Particles/GeneratedParticleProperties.inc
flymd.html
flymd.md
+Testing
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index cb746cecd7265a9c2dba0d5141644cc42468a4d6..a16c2cb1df60382498c09ba989f1bc8bffdcd2a5 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -237,6 +237,7 @@ test-clang-8:
- ctest -j4 -VV | gzip -v -9 > test.log.gz
rules:
- if: '$CI_MERGE_REQUEST_ID && $CI_MERGE_REQUEST_TITLE =~ /^WIP:/'
+ allow_failure: false
- if: $CI_MERGE_REQUEST_ID
when: manual
allow_failure: true
@@ -342,7 +343,7 @@ example-clang-8:
- cd build
- cmake --build . -- -j4
- set -o pipefail
- - ctest -j4 -VV | gzip -v -9 > test.log.gz
+ - ctest -j4 -VV | gzip -v -9 > test.log.gz
- make -j4 run_examples | gzip -v -9 > examples.log.gz
rules:
- if: '$CI_MERGE_REQUEST_ID && $CI_MERGE_REQUEST_TITLE =~ /^WIP:/'
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 86acd1300bd6f0609e0670d59e3269a55c3dba75..9b2a52509c41ab385cb38dda0a7bffd834d5dbda 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -115,9 +115,12 @@ if (CMAKE_BUILD_TYPE STREQUAL Coverage)
add_custom_target (coverage DEPENDS coverage-report)
endif ()
-# add call to ./do-copyright.py to run as unit-test-case
-add_test (NAME copyright_notices COMMAND ./do-copyright.py WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
-
+# include this test only if NOT run on gitlab-ci, since there we have a dedicated job for it:
+if (NOT DEFINED ENV{CI})
+ # add call to ./do-copyright.py to run as unit-test-case
+ add_test (NAME copyright_notices COMMAND ./do-copyright.py WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+endif (NOT DEFINED ENV{CI})
+
if (DEFINED ENV{CORSIKA_DATA})
message ("Found corsika-data in $ENV{CORSIKA_DATA}")
set (CORSIKA_DATA $ENV{CORSIKA_DATA})
diff --git a/Data b/Data
index b65a9961760359531f4ed8668e61d4e1350c7437..221d945be7d2e08b9c1432b0bcc2707e8687eee9 160000
--- a/Data
+++ b/Data
@@ -1 +1 @@
-Subproject commit b65a9961760359531f4ed8668e61d4e1350c7437
+Subproject commit 221d945be7d2e08b9c1432b0bcc2707e8687eee9
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 0118d7c87713be4f5225ec1e8de13d4a0a32c1f7..0871e4add1393c013c84a38ae509d6f2077b6e97 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -93,6 +93,7 @@ if (Pythia8_FOUND)
ProcessUrQMD
ProcessSwitch
CORSIKAcascade
+ ProcessCONEXSourceCut
ProcessEnergyLoss
ProcessObservationPlane
ProcessInteractionCounter
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index c4fd63e69fa697c5f3602ad73339d851c1e2e600..f96362addfba2a3b5f3fd30e6070be1d33c4c4fb 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -18,6 +18,7 @@
#include <corsika/geometry/Sphere.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/StackProcess.h>
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/interaction_counter/InteractionCounter.h>
#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
@@ -118,7 +119,7 @@ int main(int argc, char** argv) {
cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
<< endl;
- auto const observationHeight = 1.4_km + builder.getEarthRadius();
+ auto const observationHeight = 0_km + builder.getEarthRadius();
auto const injectionHeight = 112.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
@@ -193,17 +194,24 @@ int main(int argc, char** argv) {
process::observation_plane::ObservationPlane observationLevel(obsPlane,
"particles.dat");
- // assemble all processes into an ordered process list
-
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
+ sibyllNuc.Init();
+ sibyll.Init();
+
+ process::conex_source_cut::CONEXSourceCut conexSource(
+ center, showerAxis, t, injectionHeight, E0,
+ particles::GetPDG(particles::Code::Proton));
+
+ // assemble all processes into an ordered process list
+
auto sibyllSequence = sibyllNucCounted << sibyllCounted;
process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
55_GeV);
auto decaySequence = decayPythia << decaySibyll;
- auto sequence = switchProcess << reset_particle_mass << decaySequence << longprof
+ auto sequence = switchProcess << reset_particle_mass << decaySequence << conexSource << longprof
<< eLoss << cut << observationLevel;
// define air shower object, run simulation
@@ -218,6 +226,7 @@ int main(int argc, char** argv) {
EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile
+ conexSource.SolveCE();
cut.ShowResults();
const HEPEnergyType Efinal =
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index e90160a67a2c74c55c4a40fbad479cfe86d05416..b062a66c5aa6dcb7e76377da8ac8b0365520f9cc 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -75,3 +75,11 @@ target_link_libraries (
CORSIKAenvironment
CORSIKAtesting
)
+
+CORSIKA_ADD_TEST(testShowerAxis)
+target_link_libraries (
+ testShowerAxis
+ CORSIKAsetup
+ CORSIKAenvironment
+ CORSIKAtesting
+ )
diff --git a/Environment/ShowerAxis.cc b/Environment/ShowerAxis.cc
index f67cfb42f8473729a9d679338f8a264150240d0e..c13fdd2972fecc9324d4aa07b165b5f77bba9916 100644
--- a/Environment/ShowerAxis.cc
+++ b/Environment/ShowerAxis.cc
@@ -13,6 +13,7 @@
using namespace corsika::environment;
using namespace corsika::units::si;
+using namespace corsika;
GrammageType ShowerAxis::X(LengthType l) const {
auto const fractionalBin = l / steplength_;
@@ -33,17 +34,25 @@ GrammageType ShowerAxis::X(LengthType l) const {
assert(0 <= lambda && lambda <= 1.);
+ std::cout << l << ": " << lower << " " << lambda << " " << upper << std::endl;
+
// linear interpolation between X[lower] and X[upper]
- return X_[lower] * lambda + X_[upper] * (1 - lambda);
+ return X_[upper] * lambda + X_[lower] * (1 - lambda);
}
LengthType ShowerAxis::steplength() const { return steplength_; }
GrammageType ShowerAxis::maximumX() const { return *X_.rbegin(); }
-GrammageType ShowerAxis::minimumX() const { return *X_.cbegin(); }
+GrammageType ShowerAxis::minimumX() const { return GrammageType::zero(); }
GrammageType ShowerAxis::projectedX(geometry::Point const& p) const {
auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
return X(projectedLength);
}
+
+geometry::Vector<units::si::dimensionless_d> const& ShowerAxis::GetDirection() const {
+ return axis_normalized_;
+}
+
+geometry::Point const& ShowerAxis::GetStart() const { return pointStart_; }
diff --git a/Environment/ShowerAxis.h b/Environment/ShowerAxis.h
index 3d31a0c7c9eb3a73b92c052c67ce6ce2d85febd8..be5c605b6e38ac8febc570345cecad23154ef224 100644
--- a/Environment/ShowerAxis.h
+++ b/Environment/ShowerAxis.h
@@ -9,6 +9,7 @@
*/
#pragma once
+
#include <corsika/environment/Environment.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
@@ -23,6 +24,8 @@
#include <stdexcept>
#include <vector>
+#include <iostream>
+
#include <boost/math/quadrature/gauss_kronrod.hpp>
namespace corsika::environment {
@@ -37,7 +40,7 @@ namespace corsika::environment {
, max_length_(length_.norm())
, steplength_(max_length_ / steps)
, axis_normalized_(length / max_length_)
- , X_(steps) {
+ , X_(steps + 1) {
auto const* const universe = env.GetUniverse().get();
auto rho = [pStart, length, universe](double x) {
@@ -48,19 +51,23 @@ namespace corsika::environment {
double error;
int k = 0;
- for (int i = 1; i < steps; ++i) {
- auto const x_prev = (i - 1) / (steps - 1.);
+ X_[0] = units::si::GrammageType::zero();
+ auto sum = units::si::GrammageType::zero();
+
+ for (int i = 1; i <= steps; ++i) {
+ auto const x_prev = (i - 1.) / steps;
auto const d_prev = max_length_ * x_prev;
- auto const x = i / (steps - 1.);
+ auto const x = double(i) / steps;
auto const r = boost::math::quadrature::gauss_kronrod<double, 15>::integrate(
rho, x_prev, x, 15, 1e-9, &error);
auto const result =
units::si::MassDensityType(phys::units::detail::magnitude_tag, r) *
max_length_;
- auto const resultTotal = result + X_[i - 1];
- X_[i] = resultTotal;
- for (; resultTotal > k * X_binning_; ++k) {
+ sum += result;
+ X_[i] = sum;
+
+ for (; sum > k * X_binning_; ++k) {
d_.emplace_back(d_prev + k * X_binning_ * steplength_ / result);
}
}
@@ -79,6 +86,10 @@ namespace corsika::environment {
units::si::GrammageType X(units::si::LengthType) const;
+ geometry::Vector<units::si::dimensionless_d> const& GetDirection() const;
+
+ geometry::Point const& GetStart() const;
+
private:
geometry::Point const pointStart_;
geometry::Vector<units::si::length_d> const length_;
diff --git a/Environment/testShowerAxis.cc b/Environment/testShowerAxis.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8c2441b98fa72800a8b6394b6b30abf0b3561e2e
--- /dev/null
+++ b/Environment/testShowerAxis.cc
@@ -0,0 +1,86 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/environment/DensityFunction.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/ShowerAxis.h>
+#include <corsika/environment/VolumeTreeNode.h>
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika::geometry;
+using namespace corsika::environment;
+using namespace corsika::particles;
+using namespace corsika::units;
+using namespace corsika::units::si;
+using namespace corsika;
+
+const auto density = 1_kg / (1_m * 1_m * 1_m);
+
+auto setupEnvironment(particles::Code vTargetCode) {
+ // setup environment, geometry
+ auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
+ auto& universe = *(env->GetUniverse());
+ const geometry::CoordinateSystem& cs = env->GetCoordinateSystem();
+
+ auto theMedium =
+ environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ geometry::Point{cs, 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ density, environment::NuclearComposition(std::vector<particles::Code>{vTargetCode},
+ std::vector<float>{1.}));
+
+ auto const* nodePtr = theMedium.get();
+ universe.AddChild(std::move(theMedium));
+
+ return std::make_tuple(std::move(env), &cs, nodePtr);
+}
+
+TEST_CASE("Homogeneous Density") {
+ auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Nitrogen);
+ auto const& cs = *csPtr;
+ [[maybe_unused]] auto const& env_dummy = env;
+ [[maybe_unused]] auto const& node_dummy = nodePtr;
+
+ auto const observationHeight = 0_km;
+ auto const injectionHeight = 10_km;
+ auto const t = -observationHeight + injectionHeight;
+ Point const showerCore{cs, 0_m, 0_m, observationHeight};
+ Point const injectionPos = showerCore + Vector<dimensionless_d>{cs, {0, 0, 1}} * t;
+
+ environment::ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos),
+ *env, 20};
+
+ CHECK(showerAxis.steplength() == 500_m);
+
+ CHECK(showerAxis.maximumX() / (10_km * density) == Approx(1).epsilon(1e-8));
+
+ CHECK(showerAxis.minimumX() == 0_g / square(1_cm));
+
+ const Point p{cs, 10_km, 20_km, 8.3_km};
+ CHECK(showerAxis.projectedX(p) / (1.7_km * density) == Approx(1).epsilon(1e-8));
+
+ const units::si::LengthType d = 6.789_km;
+ CHECK(showerAxis.X(d) / (d * density) == Approx(1).epsilon(1e-8));
+
+ const Vector<dimensionless_d> dir{cs, {0, 0, -1}};
+ CHECK(showerAxis.GetDirection().GetComponents(cs) == dir.GetComponents(cs));
+ CHECK(showerAxis.GetStart().GetCoordinates() == injectionPos.GetCoordinates());
+}
diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index b0257c39ed5ca880242387c9d843382e8cea11a6..0f4e007eb937a778e890921fb486ffd00b104a4a 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/Framework/Geometry/CoordinateSystem.h
@@ -31,10 +31,6 @@ namespace corsika::geometry {
CoordinateSystem const* reference = nullptr;
EigenTransform transf;
- CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
- : reference(&reference)
- , transf(transf) {}
-
CoordinateSystem()
: // for creating the root CS
transf(EigenTransform::Identity()) {}
@@ -48,6 +44,10 @@ namespace corsika::geometry {
static EigenTransform GetTransformation(CoordinateSystem const& c1,
CoordinateSystem const& c2);
+ CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
+ : reference(&reference)
+ , transf(transf) {}
+
auto& operator=(const CoordinateSystem& pCS) {
reference = pCS.reference;
transf = pCS.transf;
@@ -60,6 +60,9 @@ namespace corsika::geometry {
return CoordinateSystem(*this, translation);
}
+ /**
+ * creates a new CS in which vVec points in direction of the new z-axis
+ */
template <typename TDim>
auto RotateToZ(Vector<TDim> vVec) const {
auto const a = vVec.normalized().GetComponents(*this).eVector;
@@ -117,6 +120,12 @@ namespace corsika::geometry {
auto const* GetReference() const { return reference; }
auto const& GetTransform() const { return transf; }
+
+ bool operator==(CoordinateSystem const& cs) const {
+ return reference == cs.reference && transf.matrix() == cs.transf.matrix();
+ }
+
+ bool operator!=(CoordinateSystem const& cs) const { return !(cs == *this); }
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 6e054d336b1a2a66e6e615b506079eae4f4aebbf..002eb6d5e478d136f7395e5ede98396d67ae3a48 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -131,7 +131,10 @@ TEST_CASE("transformations between CoordinateSystems") {
Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
- CHECK(zPrime.dot(v).magnitude() > 0);
+ CHECK(xPrime.dot(v).magnitude() == Approx(0).margin(absMargin));
+ CHECK(yPrime.dot(v).magnitude() == Approx(0).margin(absMargin));
+ CHECK((zPrime.dot(v) / 1_m).magnitude() == Approx(5 * sqrt(2)));
+
CHECK(zPrime.GetComponents(rootCS)[1].magnitude() ==
Approx(zPrime.GetComponents(rootCS)[2].magnitude()));
CHECK(zPrime.GetComponents(rootCS)[0].magnitude() == Approx(0));
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index 572839296b07634b219f518eac5192e086cda5e4..2f5c7c62f25735756de4c103cd6ffc4895168b84 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -35,6 +35,8 @@ auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
return sqrt(m * m + p.squaredNorm());
};
+auto const momentum = [](HEPEnergyType E, HEPMassType m) { return sqrt(E * E - m * m); };
+
// helper function for mandelstam-s
auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
return E * E - p.squaredNorm();
@@ -54,7 +56,7 @@ TEST_CASE("boosts") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1._GeV;
- Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 1_PeV, 0_GeV}};
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 20_GeV, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
@@ -101,18 +103,26 @@ TEST_CASE("boosts") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
- Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -1_PeV}};
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_GeV, -20_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
+ auto const sqrt_s_lab =
+ sqrt(s(eProjectileLab + targetMass, pProjectileLab.GetComponents(rootCS)));
+
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
+ auto const a = PprojCoM.GetSpaceLikeComponents().GetComponents(boost.GetRotatedCS());
+ CHECK(a.GetX() / 1_GeV == Approx(0));
+ CHECK(a.GetY() / 1_GeV == Approx(0));
+ CHECK(a.GetZ() / (momentum(sqrt_s_lab / 2, projectileMass)) == Approx(1));
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
+ CHECK(PtargCoM.GetTimeLikeComponent() / sqrt_s_lab == Approx(.5));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
@@ -183,31 +193,30 @@ TEST_CASE("boosts") {
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
}
+}
- SECTION("rest frame") {
- HEPMassType const projectileMass = 1_GeV;
- HEPMomentumType const P0 = 1_TeV;
- Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
- HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
-
- COMBoost boostRest(pProjectileLab, projectileMass);
- auto const& csPrime = boostRest.GetRotatedCS();
- FourVector const rest4Mom = boostRest.toCoM(PprojLab);
-
- CHECK(rest4Mom.GetTimeLikeComponent() / 1_GeV == Approx(projectileMass / 1_GeV));
- CHECK(rest4Mom.GetSpaceLikeComponents().norm() / 1_GeV ==
- Approx(0).margin(absMargin));
-
- FourVector const a{0_eV, Vector{csPrime, 0_eV, 5_GeV, 0_eV}};
- FourVector const b{0_eV, Vector{rootCS, 3_GeV, 0_eV, 0_eV}};
- auto const aLab = boostRest.fromCoM(a);
- auto const bLab = boostRest.fromCoM(b);
-
- CHECK(aLab.GetNorm() / a.GetNorm() == Approx(1));
- CHECK(aLab.GetSpaceLikeComponents().GetComponents(csPrime)[1].magnitude() ==
- Approx((5_GeV).magnitude()));
- CHECK(bLab.GetSpaceLikeComponents().GetComponents(rootCS)[0].magnitude() ==
- Approx((3_GeV).magnitude()));
- }
+TEST_CASE("rest frame") {
+ HEPMassType const projectileMass = 1_GeV;
+ HEPMomentumType const P0 = 1_TeV;
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
+ HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ COMBoost boostRest(pProjectileLab, projectileMass);
+ auto const& csPrime = boostRest.GetRotatedCS();
+ FourVector const rest4Mom = boostRest.toCoM(PprojLab);
+
+ CHECK(rest4Mom.GetTimeLikeComponent() / 1_GeV == Approx(projectileMass / 1_GeV));
+ CHECK(rest4Mom.GetSpaceLikeComponents().norm() / 1_GeV == Approx(0).margin(absMargin));
+
+ FourVector const a{0_eV, Vector{csPrime, 0_eV, 5_GeV, 0_eV}};
+ FourVector const b{0_eV, Vector{rootCS, 3_GeV, 0_eV, 0_eV}};
+ auto const aLab = boostRest.fromCoM(a);
+ auto const bLab = boostRest.fromCoM(b);
+
+ CHECK(aLab.GetNorm() / a.GetNorm() == Approx(1));
+ CHECK(aLab.GetSpaceLikeComponents().GetComponents(csPrime)[1].magnitude() ==
+ Approx((5_GeV).magnitude()));
+ CHECK(bLab.GetSpaceLikeComponents().GetComponents(rootCS)[0].magnitude() ==
+ Approx((3_GeV).magnitude()));
}
diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt
index f27d027a10fafc1aabbc5bb2790fe4b8262ca62d..08c5dca484b96237be9efa5b00181cfe898fd069 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -9,7 +9,7 @@ if (Pythia8_FOUND)
add_subdirectory (Pythia)
endif (Pythia8_FOUND)
if (CONEX_FOUND)
- add_subdirectory (CONEX)
+ add_subdirectory (CONEXSourceCut)
endif (CONEX_FOUND)
add_subdirectory (HadronicElasticModel)
add_subdirectory (UrQMD)
diff --git a/Processes/CONEX/CONEX.cc b/Processes/CONEX/CONEX.cc
deleted file mode 100644
index 111bbd3674ebd38bb728d93d3b04d0bb217cf8d7..0000000000000000000000000000000000000000
--- a/Processes/CONEX/CONEX.cc
+++ /dev/null
@@ -1,15 +0,0 @@
-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/process/conex/CONEX.h>
-
-using namespace corsika::process::CONEX;
-
-conex::conex() {}
diff --git a/Processes/CONEX/CONEX.h b/Processes/CONEX/CONEX.h
deleted file mode 100644
index bbc914f6e09639f054d191dcd5c1ca5ac25c9172..0000000000000000000000000000000000000000
--- a/Processes/CONEX/CONEX.h
+++ /dev/null
@@ -1,22 +0,0 @@
-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <ConexDynamicInterface.h>
-
-namespace corsika::process::CONEX {
-
- class conex {
- public:
- conex();
- };
-
-} // namespace corsika::process::CONEX
diff --git a/Processes/CONEX/testConex.cc b/Processes/CONEX/testConex.cc
deleted file mode 100644
index 1f0bb20591062d0f25a6dd43c61971ee17c6022d..0000000000000000000000000000000000000000
--- a/Processes/CONEX/testConex.cc
+++ /dev/null
@@ -1,51 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/process/conex/CONEX.h>
-
-#include <corsika/random/RNGManager.h>
-
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/utl/CorsikaFenv.h>
-#include <catch2/catch.hpp>
-
-TEST_CASE("CONEX", "[processes]") {
-
- SECTION("linking conex") {
- using std::cout;
- using std::endl;
-
- std::string parameterPathName = "";
- auto cxModel = eSibyll23;
- ConexDynamicInterface cx(cxModel);
-
- int randomSeeds[3];
- randomSeeds[0] = 1234;
- randomSeeds[1] = 0;
- randomSeeds[2] = 0;
-
- int nShower = 1; // large to avoid final stats.
- int maxDetail = 0;
- int particleListMode = 0;
- cx.Init(nShower, randomSeeds, maxDetail, particleListMode, parameterPathName);
-
- double energyInGeV = 100.;
- double zenith = 60;
- double azimuth = 0;
- double impactParameter = 0;
- int particleType = 100;
-
- cx.RunConex(randomSeeds, energyInGeV, zenith, azimuth, impactParameter, particleType);
- }
-}
diff --git a/Processes/CONEX/CMakeLists.txt b/Processes/CONEXSourceCut/CMakeLists.txt
similarity index 53%
rename from Processes/CONEX/CMakeLists.txt
rename to Processes/CONEXSourceCut/CMakeLists.txt
index 3f2ef6fd818db02c65c9a776c8b22f2f07a7fe31..1ef0f7c3f388a9fdad658b2c61a8e90d129e2ae7 100644
--- a/Processes/CONEX/CMakeLists.txt
+++ b/Processes/CONEXSourceCut/CMakeLists.txt
@@ -1,24 +1,24 @@
set (
MODEL_SOURCES
- CONEX.cc
- )
+ CONEXSourceCut.cc
+)
set (
MODEL_HEADERS
- CONEX.h
+ CONEXSourceCut.h
+ CONEX_f.h
)
set (
MODEL_NAMESPACE
- corsika/process/conex
+ corsika/process/conex_source_cut
)
-add_library (ProcessConex STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessConex ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
+add_library (ProcessCONEXSourceCut STATIC ${MODEL_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessCONEXSourceCut ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
- ProcessConex
+ ProcessCONEXSourceCut
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
@@ -26,44 +26,44 @@ set_target_properties (
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
- ProcessConex
- CORSIKAprocesssequence
- CORSIKAparticles
- CORSIKAutilities
+ ProcessCONEXSourceCut
CORSIKAunits
- CORSIKAgeometry
- CORSIKAenvironment
+ CORSIKAparticles
+ CORSIKAprocesssequence
CORSIKAsetup
- CORSIKArandom
- CorsikaData
C8::ext::conex
+ ProcessSibyll
+ ProcessUrQMD
)
target_include_directories (
- ProcessConex
- INTERFACE
+ ProcessCONEXSourceCut
+ INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
- TARGETS ProcessConex
+ TARGETS ProcessCONEXSourceCut
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
)
-
# --------------------
# code unit testing
-
-CORSIKA_ADD_TEST(testConex
- SOURCES
- testConex.cc
+CORSIKA_ADD_TEST(testCONEXSourceCut SOURCES
+ testCONEXSourceCut.cc
${MODEL_HEADERS}
)
+#
target_link_libraries (
- testConex
- ProcessConex
+ testCONEXSourceCut
+ ProcessCONEXSourceCut
+ CORSIKAunits
+ CORSIKAstackinterface
+ CORSIKAprocesssequence
+ CORSIKAsetup
+ CORSIKAgeometry
+ CORSIKAenvironment
CORSIKAtesting
- )
-
+)
diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.cc b/Processes/CONEXSourceCut/CONEXSourceCut.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5d0f68f8cb07dfc226a28b125c5d600dfd6061a5
--- /dev/null
+++ b/Processes/CONEXSourceCut/CONEXSourceCut.cc
@@ -0,0 +1,276 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
+#include <corsika/process/conex_source_cut/CONEX_f.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalConstants.h>
+#include <algorithm>
+#include <fstream>
+#include <iomanip>
+#include <utility>
+
+using namespace corsika::process::conex_source_cut;
+using namespace corsika::units::si;
+using namespace corsika::particles;
+using namespace corsika::setup;
+
+corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
+ corsika::setup::StackView& vS) {
+ auto p = vS.begin();
+
+ while (p != vS.end()) {
+ Code const pid = p.GetPID();
+
+ auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
+ [=](auto const& p) { return pid == p.first; });
+ if (it == egs_em_codes_.cend()) {
+ ++p;
+ continue; // no EM particle
+ }
+
+ int egs_pid = it->second;
+
+ addParticle(egs_pid, p.GetEnergy(), p.GetPosition(), p.GetMomentum().normalized(),
+ p.GetTime());
+
+ p.Delete();
+ }
+
+ return corsika::process::EProcessReturn::eOk;
+}
+
+void CONEXSourceCut::addParticle(int egs_pid, HEPEnergyType energy,
+ geometry::Point const& position,
+ geometry::Vector<dimensionless_d> const& direction,
+ TimeType t) {
+ std::cout << "position conexObs: " << position.GetCoordinates(conexObservationCS_)
+ << std::endl;
+
+ auto coords = position.GetCoordinates(conexObservationCS_) / 1_m;
+ double x = coords[0].magnitude();
+ double y = coords[1].magnitude();
+
+ double altitude = ((position - center_).norm() - conex::earthRadius) / 1_m;
+ auto const d = position - showerCore_;
+
+ // distance from core to particle projected along shower axis
+ double slantDistance = -d.dot(showerAxis_.GetDirection()) / 1_m;
+
+ // lateral coordinates in CONEX shower frame
+ auto const dShowerPlane = d - d.parallelProjectionOnto(showerAxis_.GetDirection());
+ double lateralX = dShowerPlane.dot(x_sf_) / 1_m;
+ double lateralY = dShowerPlane.dot(y_sf_) / 1_m;
+
+ double slantX = showerAxis_.projectedX(position) * (1_cm * 1_cm / 1_g);
+
+ double time = (t * units::constants::c - groundDist_) / 1_m;
+
+ // fill u,v,w momentum direction in EGS frame
+ double u = direction.dot(y_sf_).magnitude();
+ double v = direction.dot(x_sf_).magnitude();
+ double w = direction.dot(showerAxis_.GetDirection()).magnitude();
+
+ double weight = 1;
+
+ double E = energy / 1_GeV;
+
+ std::cout << "CONEXSourceCut: removing " << egs_pid << " " << std::scientific << energy
+ << std::endl;
+
+ std::cout << "#### parameters to show_() ####" << std::endl;
+ std::cout << "egs_pid = " << egs_pid << std::endl;
+ std::cout << "E = " << E << std::endl;
+ std::cout << "x = " << x << std::endl;
+ std::cout << "y = " << y << std::endl;
+ std::cout << "altitude = " << altitude << std::endl;
+ std::cout << "slantDistance = " << slantDistance << std::endl;
+ std::cout << "lateralX = " << lateralX << std::endl;
+ std::cout << "lateralY = " << lateralY << std::endl;
+ std::cout << "slantX = " << slantX << std::endl;
+ std::cout << "time = " << time << std::endl;
+ std::cout << "u = " << u << std::endl;
+ std::cout << "v = " << v << std::endl;
+ std::cout << "w = " << w << std::endl;
+
+ conex::cxoptl_.dptl[10 - 1] = egs_pid;
+ conex::cxoptl_.dptl[4 - 1] = E;
+ conex::cxoptl_.dptl[6 - 1] = x;
+ conex::cxoptl_.dptl[7 - 1] = y;
+ conex::cxoptl_.dptl[8 - 1] = altitude;
+ conex::cxoptl_.dptl[9 - 1] = time;
+ conex::cxoptl_.dptl[11 - 1] = weight;
+ conex::cxoptl_.dptl[13 - 1] = slantX;
+ conex::cxoptl_.dptl[14 - 1] = lateralX;
+ conex::cxoptl_.dptl[15 - 1] = lateralY;
+ conex::cxoptl_.dptl[16 - 1] = slantDistance;
+ conex::cxoptl_.dptl[2 - 1] = u;
+ conex::cxoptl_.dptl[1 - 1] = v;
+ conex::cxoptl_.dptl[3 - 1] = w;
+
+ int n = 1, i = 1;
+ conex::cegs4_(n, i);
+}
+
+void CONEXSourceCut::Init() {}
+
+void CONEXSourceCut::SolveCE() {
+
+ conex::conexcascade_();
+
+ int nX = conex::get_number_of_depth_bins_(); // make sure this works!
+
+ int icut = 1;
+ int icutg = 2;
+ int icute = 3;
+ int icutm = 2;
+ int icuth = 3;
+ int iSec = 0;
+
+ const int maxX = nX;
+
+ auto X = std::make_unique<float[]>(maxX);
+ auto H = std::make_unique<float[]>(maxX);
+ auto D = std::make_unique<float[]>(maxX);
+ auto N = std::make_unique<float[]>(maxX);
+ auto dEdX = std::make_unique<float[]>(maxX);
+ auto Mu = std::make_unique<float[]>(maxX);
+ auto dMu = std::make_unique<float[]>(maxX);
+ auto Gamma = std::make_unique<float[]>(maxX);
+ auto Electrons = std::make_unique<float[]>(maxX);
+ auto Hadrons = std::make_unique<float[]>(maxX);
+
+ float EGround[3], fitpars[13];
+
+ conex::get_shower_data_(icut, iSec, nX, X[0], N[0], fitpars[0], H[0], D[0]);
+ conex::get_shower_edep_(icut, nX, dEdX[0], EGround[0]);
+ conex::get_shower_muon_(icutm, nX, Mu[0], dMu[0]);
+ conex::get_shower_gamma_(icutg, nX, Gamma[0]);
+ conex::get_shower_electron_(icute, nX, Electrons[0]);
+ conex::get_shower_hadron_(icuth, nX, Hadrons[0]);
+
+ std::ofstream file{"conex_output.txt"};
+ for (int i = 0; i < nX; ++i) {
+ file << X[i] << " " << N[i] << " " << dEdX[i] << " " << Mu[i] << " " << dMu[i] << " "
+ << Gamma[i] << " " << Electrons[i] << " " << Hadrons[i] << std::endl;
+ }
+}
+
+
+CONEXSourceCut::CONEXSourceCut(geometry::Point center,
+ environment::ShowerAxis const& showerAxis,
+ units::si::LengthType groundDist,
+ units::si::LengthType injectionHeight,
+ units::si::HEPEnergyType primaryEnergy,
+ particles::PDGCode primaryID)
+ : center_{center}
+ , showerAxis_{showerAxis}
+ , groundDist_{groundDist}
+ , showerCore_{showerAxis_.GetStart() + showerAxis_.GetDirection() * groundDist_}
+ , conexObservationCS_{std::invoke([&]() {
+ auto const& c8cs = center.GetCoordinateSystem();
+ auto const translation = showerCore_ - center;
+ auto const intermediateCS = c8cs.translate(translation.GetComponents(c8cs));
+ auto const intermediateCS2 = intermediateCS.RotateToZ(translation);
+
+ std::cout << "translation C8/CONEX obs: " << translation.GetComponents()
+ << std::endl;
+
+ auto const transform = geometry::CoordinateSystem::GetTransformation(
+ intermediateCS2, c8cs); // either this way or vice versa... TODO: test this!
+ std::cout << transform.matrix() << std::endl << std::endl;
+ std::cout
+ << geometry::CoordinateSystem::GetTransformation(intermediateCS, c8cs).matrix()
+ << std::endl
+ << std::endl;
+ std::cout << geometry::CoordinateSystem::GetTransformation(intermediateCS2,
+ intermediateCS)
+ .matrix()
+ << std::endl;
+
+ return geometry::CoordinateSystem(c8cs, transform);
+ })}
+ , x_sf_{std::invoke([&]() {
+ geometry::Vector<length_d> const a{conexObservationCS_, 0._m, 0._m, 1._m};
+ auto b = a.cross(showerAxis_.GetDirection());
+ auto const lengthB = b.norm();
+ if (lengthB < 1e-10_m) {
+ b = geometry::Vector<length_d>{conexObservationCS_, 1_m, 0_m, 0_m};
+ }
+
+ return b.normalized();
+ })}
+ , y_sf_{showerAxis_.GetDirection().cross(x_sf_)} {
+
+ std::cout << "x_sf (conexObservationCS): " << x_sf_.GetComponents(conexObservationCS_)
+ << std::endl;
+ std::cout << "x_sf (C8): " << x_sf_.GetComponents(center.GetCoordinateSystem())
+ << std::endl;
+
+ std::cout << "y_sf (conexObservationCS): " << y_sf_.GetComponents(conexObservationCS_)
+ << std::endl;
+ std::cout << "y_sf (C8): " << y_sf_.GetComponents(center.GetCoordinateSystem())
+ << std::endl;
+
+ std::cout << "showerAxisDirection (conexObservationCS): "
+ << showerAxis_.GetDirection().GetComponents(conexObservationCS_) << std::endl;
+ std::cout << "showerAxisDirection (C8): "
+ << showerAxis_.GetDirection().GetComponents(center.GetCoordinateSystem())
+ << std::endl;
+
+ std::cout << "showerCore (conexObservationCS): "
+ << showerCore_.GetCoordinates(conexObservationCS_) << std::endl;
+ std::cout << "showerCore (C8): "
+ << showerCore_.GetCoordinates(center.GetCoordinateSystem()) << std::endl;
+
+ int randomSeeds[3] = {1234, 0, 0}; // will be overwritten later??
+ int heModel = eSibyll23;
+
+ int nShower = 1; // large to avoid final stats.
+ int maxDetail = 0;
+#ifdef CONEX_EXTENSIONS
+ int particleListMode = 0;
+#endif
+
+ std::string configPath = CONEX_CONFIG_PATH;
+ conex::initconex_(nShower, randomSeeds, heModel, maxDetail,
+#ifdef CONEX_EXTENSIONS
+ particleListMode,
+#endif
+ configPath.c_str(), configPath.size());
+
+ double eprima = primaryEnergy / 1_GeV;
+
+ // set phi, theta
+ geometry::Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, -1_m}};
+ auto const c = showerAxis_.GetDirection().dot(ez) / 1_m;
+ double theta = std::acos(c) * 180 / M_PI;
+
+ auto const showerAxisConex =
+ showerAxis_.GetDirection().GetComponents(conexObservationCS_);
+ double phi = std::atan2(-showerAxisConex.GetY().magnitude(),
+ showerAxisConex.GetX().magnitude()) *
+ 180 / M_PI;
+
+ std::cout << "theta (deg) = " << theta << "; phi (deg) = " << phi << std::endl;
+
+ int ipart = static_cast<int>(primaryID);
+ auto rng = corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
+
+ double dimpact = 0.; // valid only if shower core is fixed on the observation plane; for
+ // skimming showers an offset is needed like in CONEX
+
+ std::array<int, 3> ioseed{static_cast<int>(rng()), static_cast<int>(rng()),
+ static_cast<int>(rng())};
+
+ double xminp = injectionHeight / 1_m;
+
+ conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data());
+}
diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.h b/Processes/CONEXSourceCut/CONEXSourceCut.h
new file mode 100644
index 0000000000000000000000000000000000000000..9ea6e3d68c0038c00fe2a5d879a2fc9e49b8015f
--- /dev/null
+++ b/Processes/CONEXSourceCut/CONEXSourceCut.h
@@ -0,0 +1,69 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/environment/ShowerAxis.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/SecondariesProcess.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/process/conex_source_cut/CONEX_f.h>
+
+namespace conex {
+ corsika::units::si::LengthType constexpr earthRadius{6371315 *
+ corsika::units::si::meter};
+} // namespace conex
+
+namespace corsika::process {
+ namespace conex_source_cut {
+ class CONEXSourceCut : public process::SecondariesProcess<CONEXSourceCut> {
+
+ public:
+ CONEXSourceCut(geometry::Point center, environment::ShowerAxis const& showerAxis,
+ units::si::LengthType groundDist,
+ units::si::LengthType injectionHeight,
+ units::si::HEPEnergyType primaryEnergy,
+ particles::PDGCode primaryID);
+ corsika::process::EProcessReturn DoSecondaries(corsika::setup::StackView&);
+
+ void Init();
+
+ void SolveCE();
+
+ void addParticle(int egs_pid, units::si::HEPEnergyType energy,
+ geometry::Point const& position,
+ geometry::Vector<units::si::dimensionless_d> const& direction,
+ units::si::TimeType t);
+
+ auto const& GetObserverCS() const { return conexObservationCS_; }
+
+ private:
+
+ //! CONEX e.m. particle codes
+ static std::array<std::pair<particles::Code, int>, 3> constexpr egs_em_codes_{
+ {{particles::Code::Gamma, 0},
+ {particles::Code::Electron, -1},
+ {particles::Code::Positron, -1}}};
+
+ geometry::Point const center_; //!< center of CONEX Earth
+ environment::ShowerAxis const& showerAxis_;
+ units::si::LengthType groundDist_; //!< length from injection point to shower core
+ geometry::Point const showerCore_; //!< shower core
+ geometry::CoordinateSystem const conexObservationCS_; //!< CONEX observation frame
+ geometry::Vector<units::si::dimensionless_d> const x_sf_,
+ y_sf_; //!< unit vectors of CONEX shower frame, z_sf is shower axis direction
+ };
+ } // namespace conex_source_cut
+} // namespace corsika::process
+
diff --git a/Processes/CONEXSourceCut/CONEX_f.h b/Processes/CONEXSourceCut/CONEX_f.h
new file mode 100644
index 0000000000000000000000000000000000000000..a3076f5b63882d0cc992ab025cd1332c2bde24b5
--- /dev/null
+++ b/Processes/CONEXSourceCut/CONEX_f.h
@@ -0,0 +1,53 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+//----------------------------------------------
+// C++ interface for the CONEX fortran code
+//----------------------------------------------
+// wrapper
+
+#include <conexConfig.h>
+#include <conexHEModels.h>
+#include <array>
+
+namespace conex {
+ extern "C" {
+ extern struct { std::array<double, 16> dptl; } cxoptl_;
+
+ void cegs4_(int&, int&);
+
+ void initconex_(int&, int*, int&, int&,
+#ifdef CONEX_EXTENSIONS
+ int&,
+#endif
+ const char*, int);
+ void conexrun_(int& ipart, double& energy, double& theta, double& phi,
+ double& injectionHeight, double& dimpact, int ioseed[3]);
+ void conexcascade_();
+ void hadroncascade_(int&, int&, int&, int&);
+ void solvemomentequations_(int&);
+ void show_(int& iqi, double& ei, double& xmi, double& ymi, double& zmi, double& dmi,
+ double& xi, double& yi, double& zi, double& tmi, double& ui, double& vi,
+ double& wi, int& iri, double& wti, int& latchi);
+
+ int get_number_of_depth_bins_();
+
+ void get_shower_data_(const int&, const int&, const int&, float&, float&, float&,
+ float&, float&);
+ void get_shower_edep_(const int&, const int&, float&, float&);
+ void get_shower_muon_(const int&, const int&, float&, float&);
+ void get_shower_gamma_(const int&, const int&, float&);
+ void get_shower_electron_(const int&, const int&, float&);
+ void get_shower_hadron_(const int&, const int&, float&);
+ }
+
+} // namespace conex
diff --git a/Processes/CONEXSourceCut/testCONEXSourceCut.cc b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
new file mode 100644
index 0000000000000000000000000000000000000000..058714e9556c2c02b66c7a576de7601194743f46
--- /dev/null
+++ b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
@@ -0,0 +1,97 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
+#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::environment;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+TEST_CASE("CONEXSourceCut") {
+ random::RNGManager::GetInstance().RegisterRandomStream("cascade");
+ random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
+
+ feenableexcept(FE_INVALID);
+
+ // setup environment, geometry
+ using EnvType = Environment<setup::IEnvironmentModel>;
+ EnvType env;
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+ Point const center{rootCS, 0_m, 0_m, 0_m};
+ environment::LayeredSphericalAtmosphereBuilder builder{center, conex::earthRadius};
+ 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.assemble(env);
+
+ const HEPEnergyType E0 = 1_PeV;
+ double thetaDeg = 60.;
+ auto const thetaRad = thetaDeg / 180. * M_PI;
+
+ auto const observationHeight = 1.4_km + conex::earthRadius;
+ auto const injectionHeight = 112.75_km + conex::earthRadius;
+ auto const t =
+ -observationHeight * cos(thetaRad) +
+ sqrt(-units::si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
+ units::si::detail::static_pow<2>(injectionHeight));
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore +
+ Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ environment::ShowerAxis const showerAxis{injectionPos,
+ (showerCore - injectionPos) * 1.02, env};
+
+ corsika::process::sibyll::Interaction sibyll;
+ process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+ sibyll.Init();
+ sibyllNuc.Init();
+
+ corsika::process::conex_source_cut::CONEXSourceCut conex(
+ center, showerAxis, t, injectionHeight, E0,
+ particles::GetPDG(particles::Code::Proton));
+
+ HEPEnergyType const Eem{1_PeV};
+ auto const momentum = showerAxis.GetDirection() * Eem;
+
+ auto const emPosition = showerCore + showerAxis.GetDirection() * (-20_km);
+
+ std::cout << "position injection: "
+ << injectionPos.GetCoordinates(conex.GetObserverCS()) << " "
+ << injectionPos.GetCoordinates(rootCS) << std::endl;
+ std::cout << "position core: " << showerCore.GetCoordinates(conex.GetObserverCS())
+ << " " << showerCore.GetCoordinates(rootCS) << std::endl;
+ std::cout << "position EM: " << emPosition.GetCoordinates(conex.GetObserverCS()) << " "
+ << emPosition.GetCoordinates(rootCS) << std::endl;
+
+ conex.addParticle(0, Eem, emPosition, momentum.normalized(), 0_s);
+
+ conex.SolveCE();
+}
diff --git a/Processes/InteractionCounter/testInteractionCounter.cc b/Processes/InteractionCounter/testInteractionCounter.cc
index 797d0f87d1205fb0db28ae66b1ca21c50e709d3b..24518416aece74d2547b46d31ba6bf45e54f6b4d 100644
--- a/Processes/InteractionCounter/testInteractionCounter.cc
+++ b/Processes/InteractionCounter/testInteractionCounter.cc
@@ -119,9 +119,13 @@ TEST_CASE("InteractionCounter") {
}
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ [[maybe_unused]] auto& env_dummy = env;
+
SECTION("DoInteraction nucleus") {
unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setupStack(A, Z, 105_TeV, nodePtr, *csPtr);
+ REQUIRE(stackPtr->GetSize() == 1);
+ REQUIRE(secViewPtr->GetSize() == 0);
auto projectile = secViewPtr->GetProjectile();
auto const ret = countedProcess.DoInteraction(projectile);
@@ -140,6 +144,8 @@ TEST_CASE("InteractionCounter") {
auto constexpr code = particles::Code::Lambda0;
auto constexpr codeInt = static_cast<particles::CodeIntType>(code);
auto [stackPtr, secViewPtr] = setupStack(code, 105_TeV, nodePtr, *csPtr);
+ REQUIRE(stackPtr->GetSize() == 1);
+ REQUIRE(secViewPtr->GetSize() == 0);
auto projectile = secViewPtr->GetProjectile();
auto const ret = countedProcess.DoInteraction(projectile);
diff --git a/Processes/QGSJetII/qgsjet-II-04.f b/Processes/QGSJetII/qgsjet-II-04.f
index 6a935f62cbed7d30d41e737039507db99c0bcd55..b66cdfac2b49ae78057db05f664a7284997220fe 100644
--- a/Processes/QGSJetII/qgsjet-II-04.f
+++ b/Processes/QGSJetII/qgsjet-II-04.f
@@ -415,7 +415,9 @@ c reading cross sections from the file
if(ifIIdat.ne.1)then
inquire(file=DATDIR(1:INDEX(DATDIR,' ')-1)//'qgsdat-II-04.bz2'
* ,exist=lcalc)
- if (lcalc.and.CorDataCanDeCompress().ne.0) then
+ lcanRead=0
+ call CorDataCanDeCompress(lcanRead)
+ if (lcalc.and.lcanRead.ne.0) then
luseCompress=1
else
inquire(file=DATDIR(1:INDEX(DATDIR,' ')-1)//'qgsdat-II-04'
@@ -425,11 +427,12 @@ c reading cross sections from the file
luseCompress = (index(fnIIdat(1:nfnIIdat), ".bz2").eq.nfnIIdat-3)
inquire(file=fnIIdat(1:nfnIIdat), exist=lcalc) !used to link with nexus
endif
+
if(lcalc)then
if(ifIIdat.ne.1)then
if (luseCompress.ne.0) then
- call CorDataOpenFile(DATDIR(1:INDEX(DATDIR,' ')-1)
- * //'qgsdat-II-04.bz2')
+ call CorDataOpenFile(
+ * DATDIR(1:INDEX(DATDIR,' ')-1)//'qgsdat-II-04.bz2')
else
open(1,file=DATDIR(1:INDEX(DATDIR,' ')-1)//'qgsdat-II-04'
* ,status='old')
@@ -1970,11 +1973,13 @@ c nuclear cross sections
if(ifIIncs.ne.2)then
inquire(file=DATDIR(1:INDEX(DATDIR,' ')-1)//'sectnu-II-04.bz2'
* ,exist=lcalc)
- if (lcalc.and.CorDataCanDeCompress().ne.0) then
+ lcanRead = 0
+ call CorDataCanDeCompress(lcanRead)
+ if (lcalc.and.lcanRead.ne.0) then
luseCompress=1
else
- inquire(file=DATDIR(1:INDEX(DATDIR,' ')-1)//'sectnu-II-04'
- * ,exist=lcalc)
+ inquire(file=DATDIR(1:INDEX(DATDIR,' ')-1)//'sectnu-II-04'
+ * ,exist=lcalc)
endif
else !ctp
luseCompress = (index(fnIIncs(1:nfnIIncs), ".bz2").eq.nfnIIncs-3)
@@ -1982,12 +1987,13 @@ c nuclear cross sections
endif
if(lcalc)then
- if(debug.ge.0)write (moniou,207)
if(ifIIncs.ne.2)then
if (luseCompress.ne.0) then
+ if(debug.ge.0)write (moniou,207) 'sectnu-II-04.bz2'
call CorDataOpenFile(
* DATDIR(1:INDEX(DATDIR,' ')-1)//'sectnu-II-04.bz2')
else
+ if(debug.ge.0)write (moniou,207) 'sectnu-II-04'
open(2,file=DATDIR(1:INDEX(DATDIR,' ')-1)//'sectnu-II-04'
* ,status='old')
endif
@@ -2052,8 +2058,8 @@ c nuclear cross sections
205 format(2x,'qgaini: pretabulation of the interaction eikonals')
206 format(2x,'qgaini: initial particle energy:',e10.3,2x
*,'its type:',a7,2x,'target mass number:',i2)
-207 format(2x,'qgaini: nuclear cross sections readout from the file'
- *,' sectnu-II-04')
+ 207 format(2x,'qgaini: nuclear cross sections readout from the file:'
+ * ,A,2x)
208 format(2x,'qgaini: initial nucleus energy:',e10.3,2x
*,'projectile mass:',i2,2x,'target mass:',i2)
209 format(2x,'gtot',d10.3,' gprod',d10.3,' gabs',d10.3
@@ -2062,7 +2068,6 @@ c nuclear cross sections
212 format(2x,'qgaini: integrated Pomeron leg eikonals')
213 format(2x,'qgaini: integrated fan contributions')
214 format(2x,'qgaini: cross sections readout from the file: ', A,2x)
-c *,' qgsdat-II-04')
215 format(2x,'qgaini: integrated cut fan contributions')
c216 format(2x,'qgaini: integrated cut Pomeron eikonals')
218 format(2x,'qgaini - end')
diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc
index 8258b3c5eae4c52ab94633adae74f48a673a3d4b..707d6aea6d15cef8f2956cc5ced18ca0614c45c0 100644
--- a/Processes/QGSJetII/testQGSJetII.cc
+++ b/Processes/QGSJetII/testQGSJetII.cc
@@ -22,41 +22,60 @@
using namespace corsika;
using namespace corsika::process::qgsjetII;
+using namespace corsika::units::si;
+
+template <typename TStackView>
+auto sumCharge(TStackView const& view) {
+ int totalCharge = 0;
+
+ for (auto const& p : view) { totalCharge += particles::GetChargeNumber(p.GetPID()); }
+
+ return totalCharge;
+}
+
+template <typename TStackView>
+auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
+ geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
+
+ for (auto const& p : view) { sum += p.GetMomentum(); }
+
+ return sum;
+}
TEST_CASE("QgsjetII", "[processes]") {
SECTION("QgsjetII -> Corsika") {
- REQUIRE(particles::PiPlus::GetCode() == process::qgsjetII::ConvertFromQgsjetII(
- process::qgsjetII::QgsjetIICode::PiPlus));
+ CHECK(particles::PiPlus::GetCode() == process::qgsjetII::ConvertFromQgsjetII(
+ process::qgsjetII::QgsjetIICode::PiPlus));
}
SECTION("Corsika -> QgsjetII") {
- REQUIRE(process::qgsjetII::ConvertToQgsjetII(particles::PiMinus::GetCode()) ==
- process::qgsjetII::QgsjetIICode::PiMinus);
- REQUIRE(process::qgsjetII::ConvertToQgsjetIIRaw(particles::Proton::GetCode()) == 2);
+ CHECK(process::qgsjetII::ConvertToQgsjetII(particles::PiMinus::GetCode()) ==
+ process::qgsjetII::QgsjetIICode::PiMinus);
+ CHECK(process::qgsjetII::ConvertToQgsjetIIRaw(particles::Proton::GetCode()) == 2);
}
SECTION("canInteractInQgsjetII") {
- REQUIRE(process::qgsjetII::CanInteract(particles::Proton::GetCode()));
- REQUIRE(process::qgsjetII::CanInteract(particles::Code::KPlus));
- REQUIRE(process::qgsjetII::CanInteract(particles::Nucleus::GetCode()));
- // REQUIRE(process::qgsjetII::CanInteract(particles::Helium::GetCode()));
+ CHECK(process::qgsjetII::CanInteract(particles::Proton::GetCode()));
+ CHECK(process::qgsjetII::CanInteract(particles::Code::KPlus));
+ CHECK(process::qgsjetII::CanInteract(particles::Nucleus::GetCode()));
+ // CHECK(process::qgsjetII::CanInteract(particles::Helium::GetCode()));
- REQUIRE_FALSE(process::qgsjetII::CanInteract(particles::EtaC::GetCode()));
- REQUIRE_FALSE(process::qgsjetII::CanInteract(particles::SigmaC0::GetCode()));
+ CHECK_FALSE(process::qgsjetII::CanInteract(particles::EtaC::GetCode()));
+ CHECK_FALSE(process::qgsjetII::CanInteract(particles::SigmaC0::GetCode()));
}
SECTION("cross-section type") {
- REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) ==
- process::qgsjetII::QgsjetIIXSClass::Baryons);
- REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) ==
- process::qgsjetII::QgsjetIIXSClass::Kaons);
- REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Proton) ==
- process::qgsjetII::QgsjetIIXSClass::Baryons);
- REQUIRE(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::PiMinus) ==
- process::qgsjetII::QgsjetIIXSClass::LightMesons);
+ CHECK(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) ==
+ process::qgsjetII::QgsjetIIXSClass::Baryons);
+ CHECK(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) ==
+ process::qgsjetII::QgsjetIIXSClass::Kaons);
+ CHECK(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Proton) ==
+ process::qgsjetII::QgsjetIIXSClass::Baryons);
+ CHECK(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::PiMinus) ==
+ process::qgsjetII::QgsjetIIXSClass::LightMesons);
}
}
@@ -101,6 +120,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
random::RNGManager::GetInstance().RegisterRandomStream("qgran");
+ random::RNGManager::GetInstance().SeedAll(111);
SECTION("InteractionInterface") {
@@ -110,21 +130,26 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
- auto particle =
- stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType, unsigned int, unsigned int>{
- particles::Code::Nucleus, E0, plab, pos, 0_ns, 16, 8});
- // corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- // particles::Code::PiPlus, E0, plab, pos, 0_ns});
+ auto particle = stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Proton, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
+ auto const projectileMomentum = projectile.GetMomentum();
Interaction model;
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+
+ CHECK(length / (1_g / square(1_cm)) == Approx(93.47).margin(0.1));
+ CHECK(view.GetSize() == 14);
+ CHECK(sumCharge(view) == 1);
+ auto const secMomSum = sumMomentum(view, projectileMomentum.GetCoordinateSystem());
+ CHECK((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
+ Approx(0).margin(1e-2));
}
}
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 1a6e655982d437e0d302698c4e0e5e72ba94fb01..7332a2b0da8839f0d9989b1f854d37f9c7532866 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -164,9 +164,8 @@ namespace corsika::process::sibyll {
template <>
process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
-
- using namespace units;
using namespace utl;
+ using namespace units;
using namespace units::si;
using namespace geometry;
@@ -181,24 +180,22 @@ namespace corsika::process::sibyll {
}
if (process::sibyll::CanInteract(corsikaBeamId)) {
- const CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
// position and time of interaction, not used in Sibyll
- Point pOrig = vP.GetPosition();
- TimeType tOrig = vP.GetTime();
+ Point const pOrig = vP.GetPosition();
+ TimeType const tOrig = vP.GetTime();
+
+ // define projectile
+ HEPEnergyType const eProjectileLab = vP.GetEnergy();
+ auto const pProjectileLab = vP.GetMomentum();
+ const CoordinateSystem& originalCS = pProjectileLab.GetCoordinateSystem();
// define target
// for Sibyll is always a single nucleon
// FOR NOW: target is always at rest
const auto eTargetLab = 0_GeV + constants::nucleonMass;
- const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
+ const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
- // define projectile
- HEPEnergyType const eProjectileLab = vP.GetEnergy();
- auto const pProjectileLab = vP.GetMomentum();
-
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
<< endl;
@@ -211,6 +208,7 @@ namespace corsika::process::sibyll {
// define boost to and from CoM frame
// CoM frame definition in Sibyll projectile: +z
COMBoost const boost(PprojLab, constants::nucleonMass);
+ auto const& csPrime = boost.GetRotatedCS();
// just for show:
// boost projecticle
@@ -222,11 +220,11 @@ namespace corsika::process::sibyll {
cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: pbeam CoM: "
- << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ << PprojCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV << endl;
cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: ptarget CoM: "
- << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ << PtargCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV << endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
@@ -247,7 +245,7 @@ namespace corsika::process::sibyll {
*/
//#warning reading interaction cross section again, should not be necessary
auto const& compVec = mediumComposition.GetComponents();
- std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
+ std::vector<CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
@@ -303,7 +301,7 @@ namespace corsika::process::sibyll {
// link to sibyll stack
SibStack ss;
- MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
@@ -311,8 +309,10 @@ namespace corsika::process::sibyll {
if (psib.HasDecayed())
throw std::runtime_error("found particle that decayed in SIBYLL!");
- // transform energy to lab. frame
- auto const pCoM = psib.GetMomentum();
+ // transform 4-momentum to lab. frame
+ // note that the momentum needs to be rotated back
+ auto const tmp = psib.GetMomentum().GetComponents();
+ auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp);
HEPEnergyType const eCoM = psib.GetEnergy();
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
@@ -322,8 +322,8 @@ namespace corsika::process::sibyll {
auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
- pid, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
- tOrig});
+ pid, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
+ tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 779529e765e42aa97e10d4ec64ea7c992979f38d..a34405348a56910e52dd5b7f3f099642c7b5889c 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -30,41 +30,41 @@ using namespace corsika::units::si;
TEST_CASE("Sibyll", "[processes]") {
SECTION("Sibyll -> Corsika") {
- REQUIRE(particles::Electron::GetCode() ==
- process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
+ CHECK(particles::Electron::GetCode() ==
+ process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
}
SECTION("Corsika -> Sibyll") {
- REQUIRE(process::sibyll::ConvertToSibyll(particles::Electron::GetCode()) ==
- process::sibyll::SibyllCode::Electron);
- REQUIRE(process::sibyll::ConvertToSibyllRaw(particles::Proton::GetCode()) == 13);
- REQUIRE(process::sibyll::ConvertToSibyll(particles::XiStarC0::GetCode()) ==
- process::sibyll::SibyllCode::XiStarC0);
+ CHECK(process::sibyll::ConvertToSibyll(particles::Electron::GetCode()) ==
+ process::sibyll::SibyllCode::Electron);
+ CHECK(process::sibyll::ConvertToSibyllRaw(particles::Proton::GetCode()) == 13);
+ CHECK(process::sibyll::ConvertToSibyll(particles::XiStarC0::GetCode()) ==
+ process::sibyll::SibyllCode::XiStarC0);
}
SECTION("canInteractInSibyll") {
- REQUIRE(process::sibyll::CanInteract(particles::Proton::GetCode()));
- REQUIRE(process::sibyll::CanInteract(particles::Code::XiCPlus));
+ CHECK(process::sibyll::CanInteract(particles::Proton::GetCode()));
+ CHECK(process::sibyll::CanInteract(particles::Code::XiCPlus));
- REQUIRE_FALSE(process::sibyll::CanInteract(particles::Electron::GetCode()));
- REQUIRE_FALSE(process::sibyll::CanInteract(particles::SigmaC0::GetCode()));
+ CHECK_FALSE(process::sibyll::CanInteract(particles::Electron::GetCode()));
+ CHECK_FALSE(process::sibyll::CanInteract(particles::SigmaC0::GetCode()));
- REQUIRE_FALSE(process::sibyll::CanInteract(particles::Nucleus::GetCode()));
- REQUIRE_FALSE(process::sibyll::CanInteract(particles::Helium::GetCode()));
+ CHECK_FALSE(process::sibyll::CanInteract(particles::Nucleus::GetCode()));
+ CHECK_FALSE(process::sibyll::CanInteract(particles::Helium::GetCode()));
}
SECTION("cross-section type") {
- REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::Electron) == 0);
- REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::K0Long) == 3);
- REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::SigmaPlus) == 1);
- REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::PiMinus) == 2);
+ CHECK(process::sibyll::GetSibyllXSCode(particles::Code::Electron) == 0);
+ CHECK(process::sibyll::GetSibyllXSCode(particles::Code::K0Long) == 3);
+ CHECK(process::sibyll::GetSibyllXSCode(particles::Code::SigmaPlus) == 1);
+ CHECK(process::sibyll::GetSibyllXSCode(particles::Code::PiMinus) == 2);
}
SECTION("sibyll mass") {
- REQUIRE_FALSE(process::sibyll::GetSibyllMass(particles::Code::Electron) == 0_GeV);
+ CHECK_FALSE(process::sibyll::GetSibyllMass(particles::Code::Electron) == 0_GeV);
}
}
@@ -86,6 +86,15 @@ TEST_CASE("Sibyll", "[processes]") {
using namespace corsika::units::si;
using namespace corsika::units;
+template <typename TStackView>
+auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
+ geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
+
+ for (auto const& p : view) { sum += p.GetMomentum(); }
+
+ return sum;
+}
+
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
@@ -113,10 +122,10 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("InteractionInterface") {
setup::Stack stack;
- const HEPEnergyType E0 = 100_GeV;
+ const HEPEnergyType E0 = 60_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
@@ -130,6 +139,8 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
+ [[maybe_unused]] auto const pSum = sumMomentum(view, cs);
+
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
@@ -186,15 +197,15 @@ TEST_CASE("SibyllInterface", "[processes]") {
// run checks
// lambda decays into proton and pi- or neutron and pi+
- REQUIRE(stack.GetSize() == 3);
+ CHECK(stack.GetSize() == 3);
}
SECTION("DecayConfiguration") {
Decay model({particles::Code::PiPlus, particles::Code::PiMinus});
- REQUIRE(model.IsDecayHandled(particles::Code::PiPlus));
- REQUIRE(model.IsDecayHandled(particles::Code::PiMinus));
- REQUIRE_FALSE(model.IsDecayHandled(particles::Code::KPlus));
+ CHECK(model.IsDecayHandled(particles::Code::PiPlus));
+ CHECK(model.IsDecayHandled(particles::Code::PiMinus));
+ CHECK_FALSE(model.IsDecayHandled(particles::Code::KPlus));
const std::vector<particles::Code> particleTestList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
@@ -202,15 +213,15 @@ TEST_CASE("SibyllInterface", "[processes]") {
// setup decays
model.SetHandleDecay(particleTestList);
- for (auto& pCode : particleTestList) REQUIRE(model.IsDecayHandled(pCode));
+ for (auto& pCode : particleTestList) CHECK(model.IsDecayHandled(pCode));
// individually
model.SetHandleDecay(particles::Code::KMinus);
// possible decays
- REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Proton));
- REQUIRE_FALSE(model.CanHandleDecay(particles::Code::Electron));
- REQUIRE(model.CanHandleDecay(particles::Code::PiPlus));
- REQUIRE(model.CanHandleDecay(particles::Code::MuPlus));
+ CHECK_FALSE(model.CanHandleDecay(particles::Code::Proton));
+ CHECK_FALSE(model.CanHandleDecay(particles::Code::Electron));
+ CHECK(model.CanHandleDecay(particles::Code::PiPlus));
+ CHECK(model.CanHandleDecay(particles::Code::MuPlus));
}
}
diff --git a/Processes/SwitchProcess/SwitchProcess.h b/Processes/SwitchProcess/SwitchProcess.h
index 786f189c1a9dce98f13891ad8e937e8baa585960..26f32dd2c755e4a2cd0d7fcfd43bb9715a8177c7 100644
--- a/Processes/SwitchProcess/SwitchProcess.h
+++ b/Processes/SwitchProcess/SwitchProcess.h
@@ -74,7 +74,8 @@ namespace corsika::process::switch_process {
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if (vP.GetEnergy() < fThresholdEnergy) {
- if constexpr (is_process_sequence_v<TLowEProcess>) {
+ if constexpr (is_process_sequence_v<TLowEProcess> ||
+ is_switch_process_v<TLowEProcess>) {
return fLowEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
} else {
lambda_inv_count += fLowEProcess.GetInverseInteractionLength(vP);
@@ -87,7 +88,8 @@ namespace corsika::process::switch_process {
}
}
} else {
- if constexpr (is_process_sequence_v<THighEProcess>) {
+ if constexpr (is_process_sequence_v<THighEProcess> ||
+ is_switch_process_v<THighEProcess>) {
return fHighEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
} else {
lambda_inv_count += fHighEProcess.GetInverseInteractionLength(vP);
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index eeacb926f8e951ab45a3cd048996dbaffef78f5a..ecc75e185ec36d08bc530d7bceb1481965514380 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -33,7 +33,7 @@ using SetupProjectile = corsika::setup::StackView::StackIterator;
UrQMD::UrQMD(std::string const& xs_file) {
readXSFile(xs_file);
- iniurqmd_();
+ iniurqmdc8_();
}
CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode,
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index c477df119b536ed1bd54885dbdf0771563c4d8b3..ac24beb6f0a6198da4664702905745ff14711165 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -78,7 +78,7 @@ namespace corsika::process::UrQMD {
extern "C" {
// FORTRAN functions defined in UrQMD
- void iniurqmd_();
+ void iniurqmdc8_();
double ranf_(int&);
void cascinit_(int const&, int const&, int const&);
double nucrad_(int const&);
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index e8584ed5edf595ef7d919fd4f546469060aaff21..9127ef38167d520602fe44693002da80632529bd 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -138,6 +138,8 @@ TEST_CASE("UrQMD") {
SECTION("interaction length") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Nitrogen);
auto const& cs = *csPtr;
+ [[maybe_unused]] auto const& env_dummy = env;
+ [[maybe_unused]] auto const& node_dummy = nodePtr;
particles::Code validProjectileCodes[] = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Proton,
@@ -147,6 +149,7 @@ TEST_CASE("UrQMD") {
for (auto code : validProjectileCodes) {
auto [stack, view] = setupStack(code, 100_GeV, nodePtr, cs);
REQUIRE(stack->GetSize() == 1);
+ REQUIRE(view->GetSize() == 0);
// simple check whether the cross-section is non-vanishing
// only nuclei with available tabluated data so far
@@ -156,8 +159,13 @@ TEST_CASE("UrQMD") {
SECTION("nucleus projectile") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setupStack(A, Z, 400_GeV, nodePtr, *csPtr);
+ REQUIRE(stackPtr->GetSize() == 1);
+ REQUIRE(secViewPtr->GetSize() == 0);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
@@ -175,8 +183,13 @@ TEST_CASE("UrQMD") {
SECTION("\"special\" projectile") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
auto [stackPtr, secViewPtr] =
setupStack(particles::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
+ REQUIRE(stackPtr->GetSize() == 1);
+ REQUIRE(secViewPtr->GetSize() == 0);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
@@ -196,8 +209,13 @@ TEST_CASE("UrQMD") {
SECTION("K0Long projectile") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ [[maybe_unused]] auto const& env_dummy = env; // against warnings
+ [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
auto [stackPtr, secViewPtr] =
setupStack(particles::Code::K0Long, 400_GeV, nodePtr, *csPtr);
+ REQUIRE(stackPtr->GetSize() == 1);
+ REQUIRE(secViewPtr->GetSize() == 0);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
diff --git a/Processes/UrQMD/urqmdInterface.F b/Processes/UrQMD/urqmdInterface.F
index 259e3fd4739344c2fe5f67a8fa14bf4cd7ac17dc..ec2bd092a8386393368a6789df73b3efa3253160 100644
--- a/Processes/UrQMD/urqmdInterface.F
+++ b/Processes/UrQMD/urqmdInterface.F
@@ -11,7 +11,7 @@ c adapted by M. Reininghaus for linking UrQMD to CORSIKA 8 (Apr 2019)
c-----------------------------------------------------------------------
- subroutine IniUrQMD
+ subroutine IniUrQMDC8
c-----------------------------------------------------------------------
c Primary initialization for UrQMD 1.31
c-----------------------------------------------------------------------
diff --git a/ThirdParty/CMakeLists.txt b/ThirdParty/CMakeLists.txt
index 1ce5c22bc89df8127f5d8ff60f47736d11c01cbe..f6817bf7f8cca23ad769751575987f92cfb705ae 100644
--- a/ThirdParty/CMakeLists.txt
+++ b/ThirdParty/CMakeLists.txt
@@ -227,21 +227,18 @@ endif (NOT (${USE_CONEX_C8} IN_LIST ThirdPartyChoiceValues))
message (STATUS "USE_CONEX_C8='${USE_CONEX_C8}'")
add_library (C8::ext::conex STATIC IMPORTED GLOBAL)
-get_directory_property (Boost_iostreams_FOUND DIRECTORY ${CORSIKA_DATA}/readLib DEFINITION Boost_iostreams_FOUND)
+# this is from the corsika-data/readLib submodule:
+get_directory_property (Boost_IOSTREAMS_FOUND DIRECTORY ${CORSIKA_DATA}/readLib DEFINITION Boost_IOSTREAMS_FOUND)
set (conex_boost "")
-if (NOT Boost_iostreams_FOUND)
+set (boost_lib "boost_iostreams")
+if (NOT Boost_IOSTREAMS_FOUND)
set (conex_boost "NO_BOOST=1")
-endif (NOT Boost_iostreams_FOUND)
+ set (boost_lib "")
+endif (NOT Boost_IOSTREAMS_FOUND)
if ("x_${USE_CONEX_C8}" STREQUAL "x_SYSTEM")
find_package (CONEX REQUIRED)
message (STATUS "Using system-level CONEX version at ${CONEX_INCLUDE_DIR}")
- set_target_properties (
- C8::ext::conex PROPERTIES
- IMPORTED_LOCATION ${CONEX_LIBRARY}
- IMPORTED_LINK_INTERFACE_LIBRARIES dl
- INTERFACE_INCLUDE_DIRECTORIES ${CONEX_INCLUDE_DIR}
- )
set (CONEX_FOUND 1 PARENT_SCOPE)
else ()
@@ -257,33 +254,46 @@ else ()
GIT_REPOSITORY https://gitlab.ikp.kit.edu/AirShowerPhysics/cxroot.git
GIT_SUBMODULES ""
GIT_TAG origin/master
- GIT_SHALLOW 1
+ GIT_SHALLOW 5
GIT_PROGRESS 1
SOURCE_DIR ${CMAKE_CURRENT_BINARY_DIR}/cxroot/source
INSTALL_DIR ${CMAKE_CURRENT_BINARY_DIR}/cxroot/source
CONFIGURE_COMMAND ""
- BUILD_COMMAND make CX_NO_ROOT=1 ${conex_boost} CORSIKA_DATA=${CORSIKA_DATA} all
+ BUILD_COMMAND make CX_NO_ROOT=1 CORSIKA_8=1 ${conex_boost} CORSIKA_DATA=${CORSIKA_DATA} all
INSTALL_COMMAND ""
BUILD_IN_SOURCE ON
EXCLUDE_FROM_ALL TRUE
)
set (CONEX_FOUND 1 PARENT_SCOPE)
ExternalProject_Get_Property (cxroot INSTALL_DIR)
- set (CONEX_PREFIX ${INSTALL_DIR})
+ get_filename_component(INSTALL_DIR_ABS ${INSTALL_DIR} ABSOLUTE)
+ set (CONEX_PREFIX ${INSTALL_DIR_ABS})
set (CONEX_INCLUDE_DIR ${CONEX_PREFIX}/src)
set (CONEX_INCLUDE_DIR ${CONEX_PREFIX}/src PARENT_SCOPE)
add_dependencies (C8::ext::conex cxroot)
# create include directory at config time
file (MAKE_DIRECTORY ${CONEX_INCLUDE_DIR})
+
+endif ()
+if (Boost_IOSTREAMS_FOUND)
set_target_properties (
C8::ext::conex PROPERTIES
- IMPORTED_LOCATION ${CONEX_PREFIX}/lib/${CMAKE_SYSTEM_NAME}/libCONEXdynamic.a
- IMPORTED_LINK_INTERFACE_LIBRARIES dl
+ IMPORTED_LOCATION ${CONEX_PREFIX}/lib/${CMAKE_SYSTEM_NAME}/libCONEXsibyll.a
+ IMPORTED_NO_SONAME 1
+ SKIP_BUILD_RPATH FALSE
+ IMPORTED_LINK_INTERFACE_LIBRARIES boost_iostreams
INTERFACE_INCLUDE_DIRECTORIES
- $<BUILD_INTERFACE:${CONEX_INCLUDE_DIR}>
+ $<BUILD_INTERFACE:${CONEX_INCLUDE_DIR}>
)
-
-endif ()
-
+else (Boost_IOSTREAMS_FOUND)
+ set_target_properties (
+ C8::ext::conex PROPERTIES
+ IMPORTED_LOCATION ${CONEX_PREFIX}/lib/${CMAKE_SYSTEM_NAME}/libCONEXsibyll.a
+ IMPORTED_NO_SONAME 1
+ SKIP_BUILD_RPATH FALSE
+ INTERFACE_INCLUDE_DIRECTORIES
+ $<BUILD_INTERFACE:${CONEX_INCLUDE_DIR}>
+ )
+endif (Boost_IOSTREAMS_FOUND)