diff --git a/conanfile.txt b/conanfile.txt
index b4471ff3dbf87295603111cdfc2f80534d46102a..78cda273ada5bf3c944cb45f6f22df53d0567a8c 100644
--- a/conanfile.txt
+++ b/conanfile.txt
@@ -1,9 +1,28 @@
[requires]
+readline/8.0 # this is only needed to fix a missing dependency in "bison" which is pulled-in by arrow
spdlog/1.8.5
catch2/2.13.3
eigen/3.3.8
-boost/1.74.0
+boost/1.76.0
zlib/1.2.11
+proposal/7.0.2
+yaml-cpp/0.6.3
+arrow/2.0.0
[generators]
cmake
+
+[options]
+readline:shared=True
+arrow:shared=False
+arrow:parquet=True
+arrow:fPIC=False
+arrow:with_re2=False
+arrow:with_protobuf=False
+arrow:with_openssl=False
+arrow:with_gflags=False
+arrow:with_glog=False
+arrow:with_grpc=False
+arrow:with_utf8proc=False
+arrow:with_zstd=False
+arrow:with_bz2=False
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 7b9ec0dfb060da6284f976e61fd6fcaaf50c27f0..5eb3d20ae785b7645101b362cb4697282a29ecdd 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -350,8 +350,7 @@ namespace corsika {
A_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
// if A_ did succeed, stop routine. Not checking other static branch B_.
if (ret != ProcessReturn::Ok) { return ret; }
- } else if constexpr (is_process_v<process1_type> &&
- is_interaction_process_v<process1_type>) {
+ } else if constexpr (is_interaction_process_v<process1_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
@@ -375,8 +374,7 @@ namespace corsika {
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
- } else if constexpr (is_process_v<process2_type> &&
- is_interaction_process_v<process2_type>) {
+ } else if constexpr (is_interaction_process_v<process2_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
// soon as SecondaryView::parent() is migrated!
diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
index 0b73a606f95da6324a5f76276718583e1af4deee..3a68a3d6ebcf39290807c3402c09fcfccfe8ab6e 100644
--- a/corsika/detail/framework/utility/CubicSolver.inl
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -8,6 +8,7 @@
#pragma once
+#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/CubicSolver.hpp>
#include <cmath>
@@ -16,13 +17,10 @@ namespace corsika {
namespace andre {
//---------------------------------------------------------------------------
- // solve cubic equation x^3 + a*x^2 + b*x + c = 0
- // x - array of size 3
- // In case 3 real roots: => x[0], x[1], x[2], return 3
- // 2 real roots: x[0], x[1], return 2
- // 1 real root : x[0], x[1] ± i*x[2], return 1
- inline std::vector<double> solveP3(long double A, long double B, long double C,
- long double D, double const epsilon) {
+ // solve cubic equation A x^3 + B*x^2 + C*x + D = 0
+ inline std::vector<double> solve_cubic_real_analytic(long double A, long double B,
+ long double C, long double D,
+ double const epsilon) {
if (std::abs(A) < epsilon) { return solve_quadratic_real(B, C, epsilon); }
@@ -66,17 +64,17 @@ namespace corsika {
double(q * std::cos((t + 2 * M_PI) / 3) - a),
double(q * std::cos((t - 2 * M_PI) / 3) - a)};
} else {
- long double A = -std::pow(std::fabs(r) + std::sqrt(-disc), 1. / 3);
- if (r < 0) A = -A;
- long double B = (0 == A ? 0 : q / A);
+ long double term1 = -cbrt(std::fabs(r) + std::sqrt(-disc));
+ if (r < 0) term1 = -term1;
+ long double term2 = (0 == term1 ? 0 : q / term1);
a /= 3;
- long double test = 0.5 * std::sqrt(3.) * (A - B);
+ long double test = 0.5 * std::sqrt(3.) * (term1 - term2);
if (std::fabs(test) < epsilon) {
- return {double((A + B) - 1), double(-0.5 * (A + B) - a)};
+ return {double((term1 + term2) - 1), double(-0.5 * (term1 + term2) - a)};
}
- return {double((A + B) - a)};
+ return {double((term1 + term2) - a)};
}
}
} // namespace andre
@@ -169,7 +167,7 @@ namespace corsika {
}
/**
- Analytical approach. Not very stable in all conditions.
+ * Analytical approach. Not very stable in all conditions.
*/
inline std::vector<double> solve_cubic_real_analytic(long double a, long double b,
long double c, long double d,
@@ -206,28 +204,31 @@ namespace corsika {
double const b0 = b1 * z + c / a;
std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
- std::pow(z, 3) * a + std::pow(z, 2) * b + std::pow(z, 1) * c + d);
+ static_pow<3>(z) * a + static_pow<2>(z) * b + z * c + d);
}
CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(zs, ", "));
return zs;
}
- inline long double cubic_function(long double x, long double a, long double b,
- long double c, long double d) {
- return std::pow(x, 3) * a + std::pow(x, 2) * b + x * c + d;
+ template <typename T> // T must be floating point type
+ inline T cubic_function(T x, T a, T b, T c, T d) {
+ T x2 = x * x;
+ return x2 * x * a + x2 * b + x * c + d;
}
- inline long double cubic_function_dfdx(long double x, long double a, long double b,
- long double c) {
- return std::pow(x, 2) * a * 3 + x * b * 2 + c;
+ template <typename T> // T must be floating point type
+ inline T cubic_function_dfdx(T x, T a, T b, T c) {
+ T x2 = x * x;
+ return x2 * a * 3 + x * b * 2 + c;
}
- inline long double cubic_function_d2fd2x(long double x, long double a, long double b) {
+ template <typename T> // T must be floating point type
+ inline T cubic_function_d2fd2x(T x, T a, T b) {
return x * a * 6 + b * 2;
}
/**
- Iterative approach.
+ * Iterative approach.
*/
inline std::vector<double> solve_cubic_real(long double a, long double b, long double c,
@@ -270,7 +271,8 @@ namespace corsika {
if (std::abs(f_prime_x1) < epsilon) {
x1 -= std::cbrt(f_x1);
} else {
- x1 -= f_x1 * f_prime_x1 / (std::pow(f_prime_x1, 2) - 0.5 * f_x1 * f_prime2_x1);
+ x1 -=
+ f_x1 * f_prime_x1 / (static_pow<2>(f_prime_x1) - 0.5 * f_x1 * f_prime2_x1);
}
f_x1 = cubic_function(x1, a, b, c, d);
CORSIKA_LOG_TRACE("niter={} x1={} f_x1={} f_prime={} f_prime2={} eps={}", niter,
diff --git a/corsika/detail/framework/utility/LinearSolver.inl b/corsika/detail/framework/utility/LinearSolver.inl
index b2fc96ee7ae9a2a6947f1858b7c097e001b920bc..6eccee91ce4614d5cf95cd0948c6497824741213 100644
--- a/corsika/detail/framework/utility/LinearSolver.inl
+++ b/corsika/detail/framework/utility/LinearSolver.inl
@@ -29,7 +29,7 @@ namespace corsika {
return {}; // no (b!=0), or infinite number (b==0) of solutions....
}
- return {{-b / a, 0}};
+ return {std::complex<double>(-b / a, 0)};
}
} // namespace corsika
diff --git a/corsika/detail/framework/utility/QuadraticSolver.inl b/corsika/detail/framework/utility/QuadraticSolver.inl
index 9e67ea031c11d6160602a520746c4e3aead33a48..9c9759e975ae6d322dee7a478f935c095bdf0ff9 100644
--- a/corsika/detail/framework/utility/QuadraticSolver.inl
+++ b/corsika/detail/framework/utility/QuadraticSolver.inl
@@ -6,6 +6,8 @@
* the license.
*/
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
namespace corsika {
inline std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
@@ -19,7 +21,7 @@ namespace corsika {
return lin_result;
}
- long double const radicant = std::pow(b, 2) - a * c * 4;
+ long double const radicant = static_pow<2>(b) - a * c * 4;
if (radicant < -epsilon) { // two complex solutions
double const rpart = -b / 2 * a;
diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
index fc613fc86f49e96695f244bf241c2eec426486b4..988bbe80b315d59344c8c71f8960ac5833e270aa 100644
--- a/corsika/detail/framework/utility/QuarticSolver.inl
+++ b/corsika/detail/framework/utility/QuarticSolver.inl
@@ -8,6 +8,7 @@
#pragma once
+#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/CubicSolver.hpp>
#include <cmath>
@@ -71,30 +72,6 @@ namespace corsika {
// solving quadratic eqs. x^2 + p1*x + q1 = 0
// x^2 + p2*x + q2 = 0
-#ifdef DEBUG
- {
- [[maybe_unused]] long double Det1 = p1 * p1 - 4 * q1;
- [[maybe_unused]] long double Det2 = p2 * p2 - 4 * q2;
- CORSIKA_LOG_TRACE("Det1={} Det2={}", Det1, Det2);
- if (Det1 >= 0 && Det2 >= 0) {
- [[maybe_unused]] long double sqDet1 = sqrt(Det1);
- [[maybe_unused]] long double sqDet2 = sqrt(Det2);
- CORSIKA_LOG_TRACE("check1 {} {} {} {}", double(-p1 + sqDet1) * 0.5,
- double(-p1 - sqDet1) * 0.5, double(-p2 + sqDet2) * 0.5,
- double(-p2 - sqDet2) * 0.5);
- }
- if (Det1 >= 0) {
- [[maybe_unused]] long double sqDet1 = sqrt(Det1);
- CORSIKA_LOG_TRACE("check2 {} {} ", double(-p1 + sqDet1) * 0.5,
- double(-p1 - sqDet1) * 0.5);
- }
- if (Det2 >= 0) {
- [[maybe_unused]] long double sqDet2 = sqrt(Det2);
- CORSIKA_LOG_TRACE("check3 {} {}", double(-p1 + sqDet2) * 0.5,
- double(-p1 - sqDet2) * 0.5);
- }
- }
-#endif
std::vector<double> quad1 = solve_quadratic_real(1, p1, q1);
std::vector<double> quad2 = solve_quadratic_real(1, p2, q2);
@@ -112,8 +89,8 @@ namespace corsika {
CORSIKA_LOG_TRACE("quartic-depressed: p={:f}, q={:f}, r={:f}, epsilon={}", p, q, r,
epsilon);
- long double const p2 = std::pow(p, 2);
- long double const q2 = std::pow(q, 2);
+ long double const p2 = static_pow<2>(p);
+ long double const q2 = static_pow<2>(q);
std::vector<double> const resolve_cubic =
solve_cubic_real(1, p, p2 / 4 - r, -q2 / 8, epsilon);
@@ -125,12 +102,12 @@ namespace corsika {
long double m = 0;
for (auto const& v : resolve_cubic) {
- CORSIKA_LOG_TRACE("check pol3(v)={}", (std::pow(v, 3) + std::pow(v, 2) * p +
+ CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
v * (p2 / 4 - r) - q2 / 8));
if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
}
CORSIKA_LOG_TRACE("check m={}", m);
- if (m == 0) { return {{0}}; }
+ if (m == 0) { return {0}; }
CORSIKA_LOG_TRACE("check m={}", m);
@@ -162,12 +139,12 @@ namespace corsika {
return solve_quartic_depressed_real(c, d, e, epsilon);
}
- long double const b2 = std::pow(b, 2);
- long double const b3 = std::pow(b, 3);
- long double const b4 = std::pow(b, 4);
- long double const a2 = std::pow(a, 2);
- long double const a3 = std::pow(a, 3);
- long double const a4 = std::pow(a, 4);
+ long double const b2 = static_pow<2>(b);
+ long double const b3 = static_pow<3>(b);
+ long double const b4 = static_pow<4>(b);
+ long double const a2 = static_pow<2>(a);
+ long double const a3 = static_pow<3>(a);
+ long double const a4 = static_pow<4>(a);
long double const p = (c * a * 8 - b2 * 3) / (a4 * 8);
long double const q = (b3 - b * c * a * 4 + d * a2 * 8) / (a4 * 8);
diff --git a/corsika/detail/modules/LongitudinalProfile.inl b/corsika/detail/modules/LongitudinalProfile.inl
index ce1c1352ae34b1d25c04da6eb2bb1d8346da7021..5b1f8438bc0438c45940e34c9502cd31ec637196 100644
--- a/corsika/detail/modules/LongitudinalProfile.inl
+++ b/corsika/detail/modules/LongitudinalProfile.inl
@@ -26,8 +26,6 @@ namespace corsika {
, shower_axis_{shower_axis}
, profiles_{static_cast<unsigned int>(shower_axis.getMaximumX() / dX_) + 1} {}
- inline LongitudinalProfile::~LongitudinalProfile() {}
-
template <typename TParticle, typename TTrack>
inline ProcessReturn LongitudinalProfile::doContinuous(TParticle const& vP,
TTrack const& vTrack,
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 96dd34ff1c7980bf32849e4448dd0c4f68995468..08d5a4c9629d5901f0964c2e52218d186e31ae49 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -184,8 +184,8 @@ namespace corsika {
energy * 0.99, // either 10% relative loss max., or
get_energy_threshold(vParticle.getPID()) // energy thresholds globally defined for
// individual particles
- * 0.99999 // need to go 1% below global e-cut to assure removal in
- // ParticleCut. The 1% does not matter since at cut-time the entire
+ * 0.99999 // need to go slightly below global e-cut to assure removal in
+ // ParticleCut. This does not matter since at cut-time the entire
// energy is removed.
);
auto const maxGrammage = (vParticle.getEnergy() - energy_lim) / dEdX;
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 35ebad050be46604a48a364ecefe8b9f7f35e449..a6553a71dd3646534645bd8a1ee90c9fa29b08ca 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -59,7 +59,7 @@ namespace corsika {
inline auto Tracking::getTrack(TParticle const& particle) {
VelocityVector const initialVelocity = particle.getVelocity();
- auto const position = particle.getPosition();
+ auto const& position = particle.getPosition();
CORSIKA_LOG_DEBUG(
"TrackingLeapfrog_Curved pid: {}"
" , E = {} GeV \n"
@@ -254,7 +254,7 @@ namespace corsika {
VelocityVector const velocity = particle.getVelocity();
auto const absVelocity = velocity.getNorm();
DirectionVector const direction = velocity.normalized();
- Point const position = particle.getPosition();
+ Point const& position = particle.getPosition();
auto const magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 2799cde00172eb5dacad4ba15b7ed8745e248597..cad8ac39d99ebf67af35c19eea7974897806b3f9 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -29,7 +29,7 @@ namespace corsika::tracking_line {
VelocityVector const initialVelocity =
particle.getMomentum() / particle.getEnergy() * constants::c;
- auto const initialPosition = particle.getPosition();
+ auto const& initialPosition = particle.getPosition();
CORSIKA_LOG_DEBUG(
"TrackingStraight pid: {}"
" , E = {} GeV \n"
@@ -51,7 +51,7 @@ namespace corsika::tracking_line {
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Sphere const& sphere) {
- auto const position = particle.getPosition();
+ auto const& position = particle.getPosition();
auto const delta = position - sphere.getCenter();
auto const velocity = particle.getMomentum() / particle.getEnergy() * constants::c;
auto const vSqNorm = velocity.getSquaredNorm();
diff --git a/corsika/framework/process/ContinuousProcess.hpp b/corsika/framework/process/ContinuousProcess.hpp
index c549bc714da8e4c6b655afbb8690568a2f4c4e24..51a7ea7109439401c23e39b722fdd731ff8b1ddb 100644
--- a/corsika/framework/process/ContinuousProcess.hpp
+++ b/corsika/framework/process/ContinuousProcess.hpp
@@ -18,47 +18,46 @@
namespace corsika {
/**
- @ingroup Processes
- @{
-
- Processes with continuous effects along a particle Trajectory
-
- Create a new ContinuousProcess, e.g. for XYModel, via
- @code{.cpp}
- class XYModel : public ContinuousProcess<XYModel> {};
- @endcode
-
- and provide two necessary interface methods:
- @code{.cpp}
- template <typename TParticle, typename TTrack>
- LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const;
- @endcode
-
- which allows any ContinuousProcess to tell to CORSIKA a maximum
- allowed step length. Such step-length limitation, if it turns out
- to be smaller/sooner than any other limit (decay length,
- interaction length, other continuous processes, geometry, etc.)
- will lead to a limited step length.
-
- @code{.cpp}
- template <typename TParticle, typename TTrack>
- ProcessReturn doContinuous(TParticle& p, TTrack const& t, bool const stepLimit)
- const;
- @endcode
-
- which applied any continuous effects on Particle p along
- Trajectory t. The particle in all typical scenarios will be
- altered by a doContinuous. The flag stepLimit will be true if the
- preious evaluation of getMaxStepLength resulted in this
- particular ContinuousProcess to be responsible for the step
- length limit on the current track t. This information can be
- expoited and avoid e.g. any uncessary calculations.
-
- Particle and Track are the valid classes to
- access particles and track (Trajectory) data on the Stack. Those two methods
- do not need to be templated, they could use the types
- e.g. corsika::setup::Stack::particle_type -- but by the cost of
- loosing all flexibility otherwise provided.
+ * @ingroup Processes
+ * @{
+ * Processes with continuous effects along a particle Trajectory.
+ *
+ * Create a new ContinuousProcess, e.g. for XYModel, via:
+ * @code{.cpp}
+ * class XYModel : public ContinuousProcess<XYModel> {};
+ * @endcode
+ *
+ * and provide two necessary interface methods:
+ * @code{.cpp}
+ * template <typename TParticle, typename TTrack>
+ * LengthType getMaxStepLength(TParticle const& p, TTrack const& track) const;
+ * @endcode
+
+ * which allows any ContinuousProcess to tell to CORSIKA a maximum
+ * allowed step length. Such step-length limitation, if it turns out
+ * to be smaller/sooner than any other limit (decay length,
+ * interaction length, other continuous processes, geometry, etc.)
+ * will lead to a limited step length.
+
+ * @code{.cpp}
+ * template <typename TParticle, typename TTrack>
+ * ProcessReturn doContinuous(TParticle& p, TTrack const& t, bool const stepLimit)
+ * const;
+ * @endcode
+
+ * which applied any continuous effects on Particle p along
+ * Trajectory t. The particle in all typical scenarios will be
+ * altered by a doContinuous. The flag stepLimit will be true if the
+ * preious evaluation of getMaxStepLength resulted in this
+ * particular ContinuousProcess to be responsible for the step
+ * length limit on the current track t. This information can be
+ * expoited and avoid e.g. any uncessary calculations.
+
+ * Particle and Track are the valid classes to
+ * access particles and track (Trajectory) data on the Stack. Those two methods
+ * do not need to be templated, they could use the types
+ * e.g. corsika::setup::Stack::particle_type -- but by the cost of
+ * loosing all flexibility otherwise provided.
*/
@@ -68,8 +67,8 @@ namespace corsika {
};
/**
- * ProcessTraits specialization to flag ContinuousProcess objects
- **/
+ * ProcessTraits specialization to flag ContinuousProcess objects.
+ */
template <typename TProcess>
struct is_continuous_process<
TProcess, std::enable_if_t<
diff --git a/corsika/framework/process/ContinuousProcessIndex.hpp b/corsika/framework/process/ContinuousProcessIndex.hpp
index e36bf2812a33c6f03e4058033585f846a871b897..7192b86bf968ffdc74b5420940ea549befa12c73 100644
--- a/corsika/framework/process/ContinuousProcessIndex.hpp
+++ b/corsika/framework/process/ContinuousProcessIndex.hpp
@@ -11,12 +11,10 @@
namespace corsika {
/**
- @ingroup Processes
-
- To index individual processes (continuous processes) inside a
- ProcessSequence.
-
- **/
+ * @ingroup Processes
+ * To index individual processes (continuous processes) inside a
+ * ProcessSequence.
+ */
class ContinuousProcessIndex {
public:
diff --git a/corsika/framework/process/ContinuousProcessStepLength.hpp b/corsika/framework/process/ContinuousProcessStepLength.hpp
index 575be9063d49a32b5853320dc9d25bac40a7b9d8..55c2f417a470b8e4a475de8269c5e19fe8247ff5 100644
--- a/corsika/framework/process/ContinuousProcessStepLength.hpp
+++ b/corsika/framework/process/ContinuousProcessStepLength.hpp
@@ -14,12 +14,10 @@
namespace corsika {
/**
- @ingroup Processes
-
- To store step length in LengthType and unique index in ProcessSequence of shortest
- step ContinuousProcess.
-
- **/
+ * @ingroup Processes
+ * To store step length in LengthType and unique index in ProcessSequence of shortest
+ * step ContinuousProcess.
+ */
class ContinuousProcessStepLength {
public:
diff --git a/corsika/framework/utility/CubicSolver.hpp b/corsika/framework/utility/CubicSolver.hpp
index 694b6e3bcfe849e33270532ea53446cce6d95b12..0959a4fa84c3bee6aad856ccb57c887d21e30c9c 100644
--- a/corsika/framework/utility/CubicSolver.hpp
+++ b/corsika/framework/utility/CubicSolver.hpp
@@ -20,14 +20,14 @@
#include <boost/multiprecision/cpp_bin_float.hpp>
/**
- @file CubicSolver.hpp
+ * @file CubicSolver.hpp
*/
namespace corsika {
/**
- @ingroup Utility
- @{
+ * @ingroup Utility
+ * @{
*/
namespace andre {
@@ -57,21 +57,26 @@ namespace corsika {
double const epsilon = 1e-12);
/**
- @return a x^3 + b x^2 + c x + d
- */
- long double cubic_function(long double x, long double a, long double b, long double c,
- long double d);
+ * Cubic function.
+ *
+ * T must be a floating point type.
+ *
+ * @return a x^3 + b x^2 + c x + d
+ */
+ template <typename T>
+ T cubic_function(T x, T a, T b, T c, T d);
/**
@return 3 a x^2 + 2 b x + c
*/
- long double cubic_function_dfdx(long double x, long double a, long double b,
- long double c);
+ template <typename T>
+ T cubic_function_dfdx(T x, T a, T b, T c);
/**
@return 6 a x + 2 b
*/
- long double cubic_function_d2fd2x(long double x, long double a, long double b);
+ template <typename T>
+ T cubic_function_d2fd2x(T x, T a, T b);
/**
Iterative approach to solve: a x^3 + b x^2 + c x + d = 0
diff --git a/corsika/modules/LongitudinalProfile.hpp b/corsika/modules/LongitudinalProfile.hpp
index 770d26ed88dbb4c70890a6a1f261129ab8343244..b06e2ef852541342f474881a0209a154cf85602e 100644
--- a/corsika/modules/LongitudinalProfile.hpp
+++ b/corsika/modules/LongitudinalProfile.hpp
@@ -40,8 +40,6 @@ namespace corsika {
LongitudinalProfile(ShowerAxis const&,
GrammageType dX = 10_g / square(1_cm)); // profile binning);
- ~LongitudinalProfile();
-
template <typename TParticle, typename TTrack>
ProcessReturn doContinuous(
TParticle const&, TTrack const&,
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 1459b709a8d7fb8f1882eb0c2144225bb6a179e4..e5f6455d3f6782cef10d69c128e38534c50a2d8d 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -134,7 +134,6 @@ int main(int argc, char** argv) {
{{Code::Nitrogen, 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, 2_km);
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);
@@ -146,7 +145,7 @@ int main(int argc, char** argv) {
"environment setup: universe={}, layer1={}, layer2={}, layer3={}, layer4={}, "
"layer5={}",
fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 130000_km, 0_m, 0_m}))),
+ Point(rootCS, {constants::EarthRadius::Mean + 130_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 110_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
@@ -290,46 +289,6 @@ int main(int argc, char** argv) {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
- unsigned short const A = std::stoi(std::string(argv[1]));
- Code beamCode;
- HEPEnergyType mass;
- unsigned short Z = 0;
- if (A > 0) {
- beamCode = Code::Nucleus;
- Z = std::stoi(std::string(argv[2]));
- mass = get_nucleus_mass(A, Z);
- } else {
- int pdg = std::stoi(std::string(argv[2]));
- beamCode = convert_from_PDG(PDGCode(pdg));
- mass = get_mass(beamCode);
- }
- HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
- double theta = 0.;
- auto const thetaRad = theta / 180. * M_PI;
-
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt((Elab - m) * (Elab + m));
- };
- HEPMomentumType P0 = elab2plab(E0, mass);
- auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
- return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
- };
-
- auto const [px, py, pz] = momentumComponents(thetaRad, P0);
- auto plab = MomentumVector(rootCS, {px, py, pz});
- cout << "input particle: " << beamCode << endl;
- cout << "input angles: theta=" << theta << endl;
- cout << "input momentum: " << plab.getComponents() / 1_GeV
- << ", norm = " << plab.getNorm() << endl;
-
- auto const observationHeight = 0_km + builder.getEarthRadius();
- auto const injectionHeight = 111.75_km + builder.getEarthRadius();
- auto const t = (injectionHeight - observationHeight) / cos(thetaRad);
- Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
- Point const injectionPos =
- showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
-
- std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
@@ -349,60 +308,7 @@ int main(int argc, char** argv) {
}
}
- // we make the axis much longer than the inj-core distance since the
- // profile will go beyond the core, depending on zenith angle
- std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
- << std::endl;
-
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env,
- false};
-
- // setup processes, decays and interactions
-
- // corsika::qgsjetII::Interaction qgsjet;
- corsika::sibyll::Interaction sibyll;
- InteractionCounter sibyllCounted(sibyll);
-
- corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
- InteractionCounter sibyllNucCounted(sibyllNuc);
-
- corsika::pythia8::Decay decayPythia;
-
- // use sibyll decay routine for decays of particles unknown to pythia
- corsika::sibyll::Decay decaySibyll{{
- Code::N1440Plus,
- Code::N1440MinusBar,
- Code::N1440_0,
- Code::N1440_0Bar,
- Code::N1710Plus,
- Code::N1710MinusBar,
- Code::N1710_0,
- Code::N1710_0Bar,
-
- Code::Pi1300Plus,
- Code::Pi1300Minus,
- Code::Pi1300_0,
-
- Code::KStar0_1430_0,
- Code::KStar0_1430_0Bar,
- Code::KStar0_1430_Plus,
- Code::KStar0_1430_MinusBar,
- }};
-
- decaySibyll.printDecayConfig();
-
- ParticleCut cut{60_GeV, true, true};
- // corsika::proposal::Interaction emCascade(env);
- // corsika::proposal::ContinuousProcess emContinuous(env);
- // InteractionCounter emCascadeCounted(emCascade);
- BetheBlochPDG emContinuous(showerAxis);
-
- OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
TrackWriter trackWriter(tracks_file);
-
- LongitudinalProfile longprof{showerAxis};
-
- Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
particles_file);
@@ -420,16 +326,16 @@ int main(int argc, char** argv) {
EAS.run();
cut.showResults();
- // emContinuous.showResults();
+ emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
- cut.getEmEnergy() + // emContinuous.getEnergyLost() +
+ cut.getEmEnergy() + emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
- // emContinuous.reset();
+ emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
diff --git a/tests/framework/testLogging.cpp b/tests/framework/testLogging.cpp
index 92ad036e5b03eb940e740d060e91b561a2d84e9c..bf92cf8df00e8f12d64c0d5779f59b7976c4e09d 100644
--- a/tests/framework/testLogging.cpp
+++ b/tests/framework/testLogging.cpp
@@ -14,8 +14,6 @@ using namespace corsika;
TEST_CASE("Logging", "[Logging]") {
- // corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
SECTION("top level functions using default corsika logger") {
logging::info("(1) This is an info message!");
logging::warn("(1) This is a warning message!");
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index e687f76e1596d7497e4e440b56af52dd0d62297e..3d68eb0cd4d14cd07c09bfd9bbc8afc0a2116350 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -216,6 +216,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Electron, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
GrammageType const length = model.getInteractionLength(particle);
@@ -225,7 +226,8 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Pi0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
- setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
CHECK(view.getSize() == Approx(18).margin(2)); // this is not physics validation
@@ -234,7 +236,8 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Rho0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
- setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
CHECK(view.getSize() == Approx(12).margin(8)); // this is not physics validation
@@ -243,7 +246,8 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Lambda0, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
- setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
CHECK(view.getSize() == Approx(15).margin(10)); // this is not physics validation
@@ -252,7 +256,8 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Lambda0Bar, 0, 0, 1000_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
- setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] setup::StackView& view = *(secViewPtr.get());
+ [[maybe_unused]] auto particle = stackPtr->first();
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
CHECK(view.getSize() == Approx(40).margin(20)); // this is not physics validation
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index f09297e7699483f6fb25b9a9d6514eefcb76d6d9..3a74b386165f6016948d9104a9a671c6070256c7 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -255,7 +255,7 @@ TEST_CASE("TrackingLeapFrogCurved") {
worldPtr->addChild(std::move(target));
auto [stack, viewPtr] = setup::testing::setup_stack(PID, 0, 0, P0, targetPtr, cs);
- { [[maybe_unused]] auto& viewPtr_dum = viewPtr; }
+ { [[maybe_unused]] auto& viewPtr_dum = viewPtr; } // prevent warning
auto particle = stack->first();
// Note: momentum in X-direction
// magnetic field in X-direction
@@ -264,8 +264,9 @@ TEST_CASE("TrackingLeapFrogCurved") {
particle.setPosition(Point(cs, 0_m, 0_m, 0_m));
auto [traj, nextVol] = tracking.getTrack(particle);
- // this must be a "linear trajectory" with no curvature
+ { [[maybe_unused]] auto const& dummy = nextVol; } // prevent warning
+ // this must be a "linear trajectory" with no curvature
CHECK(traj.getDirection(0).getComponents() == traj.getDirection(1).getComponents());
}
}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index fd3ba64fe76da93403e3653635e7faca1108ffa3..f55aa032948c22026b170776d8946dc074c6e688 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -98,6 +98,7 @@ TEST_CASE("UrQMD") {
{ [[maybe_unused]] auto const& env_dummy = env; }
auto [stack, view] =
setup::testing::setup_stack(Code::Proton, 0, 0, 100_GeV, nodePtr, cs);
+ [[maybe_unused]] setup::StackView& viewRef = *(view.get());
CHECK(urqmd.getInteractionLength(stack->getNextParticle()) / 1_g * square(1_cm) ==
Approx(105).margin(5));
}
@@ -108,6 +109,7 @@ TEST_CASE("UrQMD") {
{ [[maybe_unused]] auto const& env_dummy = env; }
auto [stack, view] =
setup::testing::setup_stack(Code::Neutron, 0, 0, 100_GeV, nodePtr, cs);
+ [[maybe_unused]] setup::StackView& viewRef = *(view.get());
CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle()));
}
@@ -120,6 +122,7 @@ TEST_CASE("UrQMD") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::Nucleus, A, Z, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
+ [[maybe_unused]] setup::StackView& viewRef = *(secViewPtr.get());
CHECK(stackPtr->getEntries() == 1);
CHECK(secViewPtr->getEntries() == 0);
@@ -171,6 +174,7 @@ TEST_CASE("UrQMD") {
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
*csPtr);
+ [[maybe_unused]] auto particle = stackPtr->first();
CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); // Code::Proton not a valid target
}