From 00d47c87a3b5c4fe537e31f34703f51edcf7a9c1 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Mon, 15 Jun 2020 11:59:54 +0200
Subject: [PATCH] style, explanations
---
Environment/ShowerAxis.cc | 17 +++++++++++------
Environment/ShowerAxis.h | 25 +++++++++++++------------
2 files changed, 24 insertions(+), 18 deletions(-)
diff --git a/Environment/ShowerAxis.cc b/Environment/ShowerAxis.cc
index e7dd237d5..3ae28638a 100644
--- a/Environment/ShowerAxis.cc
+++ b/Environment/ShowerAxis.cc
@@ -1,5 +1,5 @@
/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
@@ -9,6 +9,7 @@
*/
#include <corsika/environment/ShowerAxis.h>
+#include <sstream>
using namespace corsika::environment;
using namespace corsika::units::si;
@@ -19,21 +20,25 @@ GrammageType ShowerAxis::X(LengthType l) const {
auto const lambda = fractionalBin - lower;
int const upper = lower + 1;
- if (upper >= steps) {
- throw std::runtime_error("shower axis too short, cannot extrapolate");
+ if (upper >= X_.size()) {
+ std::stringstream errormsg;
+ errormsg << "shower axis too short, cannot extrapolate (l / max_length_ = "
+ << l / max_length_ << ")";
+ throw std::runtime_error(errormsg.str().c_str());
} else if (lower < 0) {
throw std::runtime_error("cannot extrapolate to points behind point of injection");
}
assert(0 <= lambda && lambda <= 1.);
- return (*X_)[lower] * lambda + (*X_)[upper] * (1 - lambda);
+ // linear interpolation between X[lower] and X[upper]
+ return X_[lower] * lambda + X_[upper] * (1 - lambda);
}
LengthType ShowerAxis::steplength() const { return steplength_; }
-GrammageType ShowerAxis::maximumX() const { return *(X_->rbegin()); }
+GrammageType ShowerAxis::maximumX() const { return *X_.rbegin(); }
-GrammageType ShowerAxis::minimumX() const { return *(X_->cbegin()); }
+GrammageType ShowerAxis::minimumX() const { return *X_.cbegin(); }
GrammageType ShowerAxis::projectedX(geometry::Point const& p) const {
auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
diff --git a/Environment/ShowerAxis.h b/Environment/ShowerAxis.h
index 888f53398..3d31a0c7c 100644
--- a/Environment/ShowerAxis.h
+++ b/Environment/ShowerAxis.h
@@ -31,12 +31,13 @@ namespace corsika::environment {
template <typename TEnvModel>
ShowerAxis(geometry::Point const& pStart,
geometry::Vector<units::si::length_d> length,
- environment::Environment<TEnvModel> const& env)
+ environment::Environment<TEnvModel> const& env, int steps = 10'000)
: pointStart_(pStart)
, length_(length)
, max_length_(length_.norm())
, steplength_(max_length_ / steps)
- , axis_normalized_(length / max_length_) {
+ , axis_normalized_(length / max_length_)
+ , X_(steps) {
auto const* const universe = env.GetUniverse().get();
auto rho = [pStart, length, universe](double x) {
@@ -56,15 +57,15 @@ namespace corsika::environment {
auto const result =
units::si::MassDensityType(phys::units::detail::magnitude_tag, r) *
max_length_;
- auto const resultTotal = result + (*X_)[i - 1];
- (*X_)[i] = resultTotal;
+ auto const resultTotal = result + X_[i - 1];
+ X_[i] = resultTotal;
for (; resultTotal > k * X_binning_; ++k) {
d_.emplace_back(d_prev + k * X_binning_ * steplength_ / result);
}
}
- assert(std::is_sorted(X_->cbegin(), X_->cend()));
+ assert(std::is_sorted(X_.cbegin(), X_.cend()));
assert(std::is_sorted(d_.cbegin(), d_.cend()));
}
@@ -79,17 +80,17 @@ namespace corsika::environment {
units::si::GrammageType X(units::si::LengthType) const;
private:
- static constexpr int steps = 10'000;
- units::si::GrammageType const X_binning_ = std::invoke([]() {
- using namespace units::si;
- return 1_g / 1_cm / 1_cm;
- });
- std::unique_ptr<std::array<units::si::GrammageType, steps>> const X_{
- new std::array<units::si::GrammageType, steps>};
geometry::Point const pointStart_;
geometry::Vector<units::si::length_d> const length_;
units::si::LengthType const max_length_, steplength_;
geometry::Vector<units::si::dimensionless_d> const axis_normalized_;
+ std::vector<units::si::GrammageType> X_;
+
+ // for storing the lengths corresponding to equidistant X values
+ units::si::GrammageType const X_binning_ = std::invoke([]() {
+ using namespace units::si;
+ return 1_g / 1_cm / 1_cm;
+ });
std::vector<units::si::LengthType> d_;
};
} // namespace corsika::environment
--
GitLab