From b7deecfabe0b063b6e511f9d7bf72fd7f1fb360b Mon Sep 17 00:00:00 2001 From: ralfulrich Date: Tue, 15 Jun 2021 18:14:35 +0200 Subject: [PATCH 1/8] some density modelling refactoring --- corsika/detail/media/BaseExponential.inl | 47 ++++++++----- corsika/detail/media/FlatExponential.inl | 14 ++-- corsika/detail/media/HomogeneousMedium.inl | 6 +- corsika/detail/media/InhomogeneousMedium.inl | 4 +- .../LayeredSphericalAtmosphereBuilder.inl | 35 ++++++++++ .../media/LinearApproximationIntegrator.inl | 6 +- .../detail/media/SlidingPlanarExponential.inl | 14 ++-- .../modules/energy_loss/BetheBlochPDG.inl | 3 +- .../modules/proposal/ContinuousProcess.inl | 3 +- corsika/media/BaseExponential.hpp | 21 +++--- corsika/media/FlatExponential.hpp | 3 +- corsika/media/HomogeneousMedium.hpp | 3 +- corsika/media/IMediumModel.hpp | 24 +++++-- corsika/media/InhomogeneousMedium.hpp | 3 +- .../LayeredSphericalAtmosphereBuilder.hpp | 4 ++ .../media/LinearApproximationIntegrator.hpp | 2 +- corsika/media/SlidingPlanarExponential.hpp | 3 +- tests/media/testEnvironment.cpp | 69 +++++++++++-------- tests/media/testMagneticField.cpp | 5 -- tests/media/testMedium.cpp | 9 ++- tests/media/testRefractiveIndex.cpp | 29 ++++---- 21 files changed, 193 insertions(+), 114 deletions(-) diff --git a/corsika/detail/media/BaseExponential.inl b/corsika/detail/media/BaseExponential.inl index 3d7790c29..dbd418041 100644 --- a/corsika/detail/media/BaseExponential.inl +++ b/corsika/detail/media/BaseExponential.inl @@ -14,23 +14,43 @@ namespace corsika { + template + inline BaseExponential::BaseExponential(Point const& point, + LengthType const referenceHeight, + MassDensityType rho0, + LengthType lambda) + : rho0_(rho0) + , lambda_(lambda) + , invLambda_(1 / lambda) + , point_(point) + , referenceHeight_(referenceHeight) {} + template inline auto const& BaseExponential::getImplementation() const { return *static_cast(this); } + template + inline MassDensityType BaseExponential::getMassDensity( + LengthType const height) const { + return rho0_ * exp(invLambda_ * (height - referenceHeight_)); + } + template inline GrammageType BaseExponential::getIntegratedGrammage( - BaseTrajectory const& traj, LengthType vL, DirectionVector const& axis) const { - if (vL == LengthType::zero()) { return GrammageType::zero(); } + BaseTrajectory const& traj, DirectionVector const& axis) const { + LengthType const length = traj.getLength(); + if (length == LengthType::zero()) { return GrammageType::zero(); } - auto const uDotA = traj.getDirection(0).dot(axis).magnitude(); - auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0)); + // this corresponds to height: + double const uDotA = traj.getDirection(0).dot(axis).magnitude(); + MassDensityType const rhoStart = + getImplementation().getMassDensity(traj.getPosition(0)); if (uDotA == 0) { - return vL * rhoStart; + return length * rhoStart; } else { - return rhoStart * (lambda_ / uDotA) * (exp(uDotA * vL * invLambda_) - 1); + return rhoStart * (lambda_ / uDotA) * (exp(uDotA * length * invLambda_) - 1); } } @@ -38,8 +58,10 @@ namespace corsika { inline LengthType BaseExponential::getArclengthFromGrammage( BaseTrajectory const& traj, GrammageType grammage, DirectionVector const& axis) const { - auto const uDotA = traj.getDirection(0).dot(axis).magnitude(); - auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0)); + // this corresponds to height: + double const uDotA = traj.getDirection(0).dot(axis).magnitude(); + MassDensityType const rhoStart = + getImplementation().getMassDensity(traj.getPosition(0)); if (uDotA == 0) { return grammage / rhoStart; @@ -54,13 +76,4 @@ namespace corsika { } } - template - inline BaseExponential::BaseExponential(Point const& point, - MassDensityType rho0, - LengthType lambda) - : rho0_(rho0) - , lambda_(lambda) - , invLambda_(1 / lambda) - , point_(point) {} - } // namespace corsika diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl index b64b2ded2..f665dda99 100644 --- a/corsika/detail/media/FlatExponential.inl +++ b/corsika/detail/media/FlatExponential.inl @@ -18,19 +18,17 @@ namespace corsika { template inline FlatExponential::FlatExponential(Point const& point, - Vector const& axis, + DirectionVector const& axis, MassDensityType rho, LengthType lambda, NuclearComposition const& nuclComp) - : BaseExponential>(point, rho, lambda) + : BaseExponential>(point, 0_m, rho, lambda) , axis_(axis) , nuclComp_(nuclComp) {} template inline MassDensityType FlatExponential::getMassDensity(Point const& point) const { - return BaseExponential>::getRho0() * - exp(BaseExponential>::getInvLambda() * - (point - BaseExponential>::getAnchorPoint()) - .dot(axis_)); + return BaseExponential>::getMassDensity( + (point - BaseExponential>::getAnchorPoint()).getNorm()); } template @@ -40,8 +38,8 @@ namespace corsika { template inline GrammageType FlatExponential::getIntegratedGrammage( - BaseTrajectory const& line, LengthType to) const { - return BaseExponential>::getIntegratedGrammage(line, to, axis_); + BaseTrajectory const& line) const { + return BaseExponential>::getIntegratedGrammage(line, axis_); } template diff --git a/corsika/detail/media/HomogeneousMedium.inl b/corsika/detail/media/HomogeneousMedium.inl index 368073c9f..032209bd3 100644 --- a/corsika/detail/media/HomogeneousMedium.inl +++ b/corsika/detail/media/HomogeneousMedium.inl @@ -32,9 +32,9 @@ namespace corsika { } template - inline GrammageType HomogeneousMedium::getIntegratedGrammage(BaseTrajectory const&, - LengthType to) const { - return to * density_; + inline GrammageType HomogeneousMedium::getIntegratedGrammage( + BaseTrajectory const& track) const { + return track.getLength() * density_; } template diff --git a/corsika/detail/media/InhomogeneousMedium.inl b/corsika/detail/media/InhomogeneousMedium.inl index 8706765ad..f53f2c3be 100644 --- a/corsika/detail/media/InhomogeneousMedium.inl +++ b/corsika/detail/media/InhomogeneousMedium.inl @@ -36,8 +36,8 @@ namespace corsika { template inline GrammageType InhomogeneousMedium::getIntegratedGrammage( - BaseTrajectory const& line, LengthType to) const { - return densityFunction_.getIntegrateGrammage(line, to); + BaseTrajectory const& line) const { + return densityFunction_.getIntegrateGrammage(line); } template diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl index 74f00f101..791d57972 100644 --- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl +++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl @@ -13,6 +13,7 @@ #include #include #include +//#include namespace corsika { @@ -103,6 +104,40 @@ namespace corsika { layers_.push(std::move(node)); } + template typename TMediumModelExtra, + typename... TModelArgs> + inline void + LayeredSphericalAtmosphereBuilder:: + addTabularLayer(GrammageType b, + std::function const& funcRho, + LengthType upperBoundary) { + + auto const radius = earthRadius_ + upperBoundary; + checkRadius(radius); + previousRadius_ = radius; + + auto node = std::make_unique>( + std::make_unique(center_, radius)); + /* + if constexpr (detail::has_extra_models::value) { + // helper lambda in which the last 5 arguments to make_shared<...> are bound + auto lastBound = [&](auto... argPack) { + return std::make_shared< + TMediumModelExtra>>( + argPack..., center_, rho0, -c, *composition_, earthRadius_); + }; + + // now unpack the additional arguments + auto model = std::apply(lastBound, additionalModelArgs_); + node->setModelProperties(std::move(model)); + } else { + node->template setModelProperties>( + center_, rho0, -c, *composition_, earthRadius_); + } + */ + layers_.push(std::move(node)); + } + template typename TMediumModelExtra, typename... TModelArgs> inline Environment LayeredSphericalAtmosphereBuilder< diff --git a/corsika/detail/media/LinearApproximationIntegrator.inl b/corsika/detail/media/LinearApproximationIntegrator.inl index ced2dee38..8ee60d984 100644 --- a/corsika/detail/media/LinearApproximationIntegrator.inl +++ b/corsika/detail/media/LinearApproximationIntegrator.inl @@ -19,10 +19,14 @@ namespace corsika { template inline auto LinearApproximationIntegrator::getIntegrateGrammage( - BaseTrajectory const& line, LengthType length) const { + BaseTrajectory const& line) const { + LengthType const length = line.getLength(); auto const c0 = getImplementation().evaluateAt(line.getPosition(0)); auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0), line.getDirection(0)); + CORSIKA_LOG_INFO("length={} c0={} c1={} pos={} dir={} return={}", length, c0, c1, + line.getPosition(0), line.getDirection(0), + (c0 + 0.5 * c1 * length) * length); return (c0 + 0.5 * c1 * length) * length; } diff --git a/corsika/detail/media/SlidingPlanarExponential.inl b/corsika/detail/media/SlidingPlanarExponential.inl index 61911772a..80fc9bf6b 100644 --- a/corsika/detail/media/SlidingPlanarExponential.inl +++ b/corsika/detail/media/SlidingPlanarExponential.inl @@ -16,19 +16,17 @@ namespace corsika { inline SlidingPlanarExponential::SlidingPlanarExponential( Point const& p0, MassDensityType rho0, LengthType lambda, NuclearComposition const& nuclComp, LengthType referenceHeight) - : BaseExponential>(p0, rho0, lambda) + : BaseExponential>(p0, referenceHeight, rho0, lambda) , nuclComp_(nuclComp) , referenceHeight_(referenceHeight) {} template inline MassDensityType SlidingPlanarExponential::getMassDensity( Point const& point) const { - auto const height = + auto const heightFull = (point - BaseExponential>::getAnchorPoint()) - .getNorm() - - referenceHeight_; - return BaseExponential>::getRho0() * - exp(BaseExponential>::getInvLambda() * height); + .getNorm(); + return BaseExponential>::getMassDensity(heightFull); } template @@ -39,11 +37,11 @@ namespace corsika { template inline GrammageType SlidingPlanarExponential::getIntegratedGrammage( - BaseTrajectory const& traj, LengthType l) const { + BaseTrajectory const& traj) const { auto const axis = (traj.getPosition(0) - BaseExponential>::getAnchorPoint()) .normalized(); - return BaseExponential>::getIntegratedGrammage(traj, l, + return BaseExponential>::getIntegratedGrammage(traj, axis); } diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl index 7772f2660..537a3d85f 100644 --- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl +++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl @@ -154,8 +154,7 @@ namespace corsika { if (p.getChargeNumber() == 0) return ProcessReturn::Ok; - GrammageType const dX = - p.getNode()->getModelProperties().getIntegratedGrammage(t, t.getLength()); + GrammageType const dX = p.getNode()->getModelProperties().getIntegratedGrammage(t); CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(), p.getChargeNumber(), dX / 1_g * square(1_cm)); HEPEnergyType dE = getTotalEnergyLoss(p, dX); diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl index aa50f5866..30f0b10a1 100644 --- a/corsika/detail/modules/proposal/ContinuousProcess.inl +++ b/corsika/detail/modules/proposal/ContinuousProcess.inl @@ -97,8 +97,7 @@ namespace corsika::proposal { if (vT.getLength() == 0_m) return ProcessReturn::Ok; // calculate passed grammage - auto dX = - vP.getNode()->getModelProperties().getIntegratedGrammage(vT, vT.getLength()); + auto dX = vP.getNode()->getModelProperties().getIntegratedGrammage(vT); // get or build corresponding track integral calculator and solve the // integral diff --git a/corsika/media/BaseExponential.hpp b/corsika/media/BaseExponential.hpp index 66f907943..27a267df6 100644 --- a/corsika/media/BaseExponential.hpp +++ b/corsika/media/BaseExponential.hpp @@ -23,9 +23,20 @@ namespace corsika { */ template class BaseExponential { + + public: + BaseExponential(Point const& point, LengthType const referenceHeight, + MassDensityType rho0, LengthType lambda); + + Point const& getAnchorPoint() const { return point_; } + MassDensityType getRho0() const { return rho0_; } + InverseLengthType getInvLambda() const { return invLambda_; } + protected: auto const& getImplementation() const; + MassDensityType getMassDensity(LengthType const height) const; + // clang-format off /** * For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized) @@ -40,7 +51,7 @@ namespace corsika { * \f] */ // clang-format on - GrammageType getIntegratedGrammage(BaseTrajectory const& line, LengthType vL, + GrammageType getIntegratedGrammage(BaseTrajectory const& line, DirectionVector const& axis) const; // clang-format off @@ -64,18 +75,12 @@ namespace corsika { LengthType getArclengthFromGrammage(BaseTrajectory const& line, GrammageType grammage, DirectionVector const& axis) const; - public: - BaseExponential(Point const& point, MassDensityType rho0, LengthType lambda); - - Point const& getAnchorPoint() const { return point_; } - MassDensityType getRho0() const { return rho0_; } - InverseLengthType getInvLambda() const { return invLambda_; } - private: MassDensityType const rho0_; LengthType const lambda_; InverseLengthType const invLambda_; Point const point_; + LengthType const referenceHeight_; }; // class BaseExponential diff --git a/corsika/media/FlatExponential.hpp b/corsika/media/FlatExponential.hpp index 71767d705..6f64295e8 100644 --- a/corsika/media/FlatExponential.hpp +++ b/corsika/media/FlatExponential.hpp @@ -41,8 +41,7 @@ namespace corsika { NuclearComposition const& getNuclearComposition() const override; - GrammageType getIntegratedGrammage(BaseTrajectory const& line, - LengthType to) const override; + GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override; LengthType getArclengthFromGrammage(BaseTrajectory const& line, GrammageType grammage) const override; diff --git a/corsika/media/HomogeneousMedium.hpp b/corsika/media/HomogeneousMedium.hpp index 95a9fb0ae..16331957d 100644 --- a/corsika/media/HomogeneousMedium.hpp +++ b/corsika/media/HomogeneousMedium.hpp @@ -30,8 +30,7 @@ namespace corsika { NuclearComposition const& getNuclearComposition() const override; - GrammageType getIntegratedGrammage(BaseTrajectory const&, - LengthType to) const override; + GrammageType getIntegratedGrammage(BaseTrajectory const&) const override; LengthType getArclengthFromGrammage(BaseTrajectory const&, GrammageType grammage) const override; diff --git a/corsika/media/IMediumModel.hpp b/corsika/media/IMediumModel.hpp index 87e702c5d..e3e030392 100644 --- a/corsika/media/IMediumModel.hpp +++ b/corsika/media/IMediumModel.hpp @@ -21,11 +21,25 @@ namespace corsika { virtual MassDensityType getMassDensity(Point const&) const = 0; - // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory - // approach; for now, only lines are supported - virtual GrammageType getIntegratedGrammage(BaseTrajectory const&, - LengthType) const = 0; - + /** + * Integrate the matter density along trajectory. + * + * @return GrammageType as integrated matter density along the BaseTrajectory + * + * @todo think about the mixin inheritance of the trajectory vs the BaseTrajectory + * approach; for now, only lines are supported (?). + */ + virtual GrammageType getIntegratedGrammage(BaseTrajectory const&) const = 0; + + /** + * Calculates the length along the trajectory. + * + * The length along the trajectory is determined at which the integrated matter + * density is reached. If the specified matter density cannot be reached (is too + * large) the result becomes meaningless and could be "infinity" (discuss this). + * + * @return LengthType the length corresponding to grammage. + */ virtual LengthType getArclengthFromGrammage(BaseTrajectory const&, GrammageType) const = 0; diff --git a/corsika/media/InhomogeneousMedium.hpp b/corsika/media/InhomogeneousMedium.hpp index a257bc508..66cf00f73 100644 --- a/corsika/media/InhomogeneousMedium.hpp +++ b/corsika/media/InhomogeneousMedium.hpp @@ -32,8 +32,7 @@ namespace corsika { NuclearComposition const& getNuclearComposition() const override; - GrammageType getIntegratedGrammage(BaseTrajectory const& line, - LengthType to) const override; + GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override; LengthType getArclengthFromGrammage(BaseTrajectory const& pLine, GrammageType grammage) const override; diff --git a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp index 5de160d13..41ad15801 100644 --- a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp +++ b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp @@ -78,6 +78,10 @@ namespace corsika { void addExponentialLayer(GrammageType b, LengthType c, LengthType upperBoundary); void addLinearLayer(LengthType c, LengthType upperBoundary); + void addTabularLayer(GrammageType b, + std::function const& funcRho, + LengthType upperBoundary); + int getSize() const { return layers_.size(); } void assemble(Environment& env); diff --git a/corsika/media/LinearApproximationIntegrator.hpp b/corsika/media/LinearApproximationIntegrator.hpp index 0e2eb0a29..ef442ebdb 100644 --- a/corsika/media/LinearApproximationIntegrator.hpp +++ b/corsika/media/LinearApproximationIntegrator.hpp @@ -20,7 +20,7 @@ namespace corsika { auto const& getImplementation() const; public: - auto getIntegrateGrammage(BaseTrajectory const& line, LengthType length) const; + auto getIntegrateGrammage(BaseTrajectory const& line) const; auto getArclengthFromGrammage(BaseTrajectory const& line, GrammageType grammage) const; diff --git a/corsika/media/SlidingPlanarExponential.hpp b/corsika/media/SlidingPlanarExponential.hpp index a06c293a4..63f8238ba 100644 --- a/corsika/media/SlidingPlanarExponential.hpp +++ b/corsika/media/SlidingPlanarExponential.hpp @@ -46,8 +46,7 @@ namespace corsika { NuclearComposition const& getNuclearComposition() const override; - GrammageType getIntegratedGrammage(BaseTrajectory const& line, - LengthType l) const override; + GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override; LengthType getArclengthFromGrammage(BaseTrajectory const& line, GrammageType grammage) const override; diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp index b39fbb357..88e143b7b 100644 --- a/tests/media/testEnvironment.cpp +++ b/tests/media/testEnvironment.cpp @@ -97,33 +97,35 @@ TEST_CASE("FlatExponential") { auto const rho0 = 1_g / static_pow<3>(1_cm); FlatExponential const medium(gOrigin, axis, rho0, lambda, protonComposition); - auto const tEnd = 5_s; + SpeedType const speed = 20_m / second; + LengthType const length = 2_m; + TimeType const tEnd = length / speed; SECTION("horizontal") { Line const line(gOrigin, Vector( - gCS, {20_cm / second, 0_m / second, 0_m / second})); + gCS, {speed, 0_m / second, 0_m / second})); setup::Trajectory const trajectory = setup::testing::make_track(line, tEnd); - CHECK((medium.getIntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1)); - CHECK((medium.getArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1)); + CHECK((medium.getIntegratedGrammage(trajectory) / (rho0 * length)) == Approx(1)); + CHECK((medium.getArclengthFromGrammage(trajectory, rho0 * length) / length) == + Approx(1)); } SECTION("vertical") { Line const line(gOrigin, Vector( - gCS, {0_m / second, 0_m / second, 5_m / second})); + gCS, {0_m / second, 0_m / second, speed})); setup::Trajectory const trajectory = setup::testing::make_track(line, tEnd); - LengthType const length = 2 * lambda; GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1); - CHECK((medium.getIntegratedGrammage(trajectory, length) / exact) == Approx(1)); + CHECK((medium.getIntegratedGrammage(trajectory) / exact) == Approx(1)); CHECK((medium.getArclengthFromGrammage(trajectory, exact) / length) == Approx(1)); } SECTION("escape grammage") { Line const line(gOrigin, Vector( - gCS, {0_m / second, 0_m / second, -5_m / second})); + gCS, {SpeedType::zero(), SpeedType::zero(), -speed})); setup::Trajectory const trajectory = setup::testing::make_track(line, tEnd); GrammageType const escapeGrammage = rho0 * lambda; @@ -134,15 +136,15 @@ TEST_CASE("FlatExponential") { } SECTION("inclined") { - Line const line(gOrigin, Vector( - gCS, {0_m / second, 5_m / second, 5_m / second})); + Line const line(gOrigin, + Vector( + gCS, {0_m / second, speed / sqrt(2.), speed / sqrt(2.)})); setup::Trajectory const trajectory = setup::testing::make_track(line, tEnd); double const cosTheta = M_SQRT1_2; - LengthType const length = 2 * lambda; GrammageType const exact = rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta; - CHECK((medium.getIntegratedGrammage(trajectory, length) / exact) == Approx(1)); + CHECK((medium.getIntegratedGrammage(trajectory) / exact) == Approx(1)); CHECK((medium.getArclengthFromGrammage(trajectory, exact) / length) == Approx(1)); } } @@ -179,16 +181,16 @@ TEST_CASE("SlidingPlanarExponential") { CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() == flat.getMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude()); - CHECK(medium.getIntegratedGrammage(trajectory, 2_m).magnitude() == - flat.getIntegratedGrammage(trajectory, 2_m).magnitude()); + CHECK(medium.getIntegratedGrammage(trajectory).magnitude() == + flat.getIntegratedGrammage(trajectory).magnitude()); CHECK(medium.getArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude() == flat.getArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude()); } } -auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m; +MassDensityType constexpr rho0 = 1_kg / 1_m / 1_m / 1_m; -struct Exponential { +struct ExponentialTest { auto operator()(Point const& p) const { return exp(p.getCoordinates()[0] / 1_m) * rho0; } @@ -213,41 +215,50 @@ TEST_CASE("InhomogeneousMedium") { Vector direction(gCS, QuantityVector(1, 0, 0)); + SpeedType const speed = 20_m / second; Line line(gOrigin, Vector( - gCS, {20_m / second, 0_m / second, 0_m / second})); + gCS, {speed, SpeedType::zero(), SpeedType::zero()})); - auto const tEnd = 5_s; + // the tested LinearApproximationIntegrator really does a single step only. It is very + // poor for exponentials with a bit larger step-width. + TimeType const tEnd = 0.001_s; setup::Trajectory const trajectory = setup::testing::make_track(line, tEnd); - Exponential const e; - DensityFunction const rho(e); + ExponentialTest const expTest; + DensityFunction const rho(expTest); SECTION("DensityFunction") { - CHECK(e.getDerivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) == + CHECK(expTest.getDerivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) == Approx(1)); - CHECK(rho.evaluateAt(gOrigin) == e(gOrigin)); + CHECK(rho.evaluateAt(gOrigin) == expTest(gOrigin)); } auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); }; auto const exactLength = [](auto X) { return 1_m * log(1 + X / (rho0 * 1_m)); }; - auto constexpr l = 15_cm; + LengthType constexpr length = tEnd * speed; NuclearComposition const composition{{Code::Proton}, {1.f}}; InhomogeneousMedium const inhMedium(composition, rho); + CORSIKA_LOG_INFO("test={} l={} {} {}", rho.getIntegrateGrammage(trajectory), length, + exactGrammage(length), 1_m * rho0 * (exp(length / 1_m) - 1)); + SECTION("Integration") { - CHECK(rho.getIntegrateGrammage(trajectory, l) / exactGrammage(l) == + CORSIKA_LOG_INFO("test={} {} {}", rho.getIntegrateGrammage(trajectory), + exactGrammage(length), + rho.getIntegrateGrammage(trajectory) / exactGrammage(length)); + CHECK(rho.getIntegrateGrammage(trajectory) / exactGrammage(length) == Approx(1).epsilon(1e-2)); - CHECK(rho.getArclengthFromGrammage(trajectory, exactGrammage(l)) / - exactLength(exactGrammage(l)) == + CHECK(rho.getArclengthFromGrammage(trajectory, exactGrammage(length)) / + exactLength(exactGrammage(length)) == Approx(1).epsilon(1e-2)); CHECK(rho.getMaximumLength(trajectory, 1e-2) > - l); // todo: write reasonable test when implementation is working + length); // todo: write reasonable test when implementation is working - CHECK(rho.getIntegrateGrammage(trajectory, l) == - inhMedium.getIntegratedGrammage(trajectory, l)); + CHECK(rho.getIntegrateGrammage(trajectory) == + inhMedium.getIntegratedGrammage(trajectory)); CHECK(rho.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) == inhMedium.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm))); CHECK(inhMedium.getNuclearComposition() == composition); diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp index b3ccf4d75..6e9f347b8 100644 --- a/tests/media/testMagneticField.cpp +++ b/tests/media/testMagneticField.cpp @@ -25,7 +25,6 @@ using namespace corsika; TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { logging::set_level(logging::level::info); - corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v"); CoordinateSystemPtr const& gCS = get_root_CoordinateSystem(); Point const gOrigin(gCS, {0_m, 0_m, 0_m}); @@ -85,8 +84,4 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { // and the associated trajectory setup::Trajectory const track = setup::testing::make_track(line, tEnd); - - // and check the integrated grammage - CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); - CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); } diff --git a/tests/media/testMedium.cpp b/tests/media/testMedium.cpp index a10348717..bb5a9f0e8 100644 --- a/tests/media/testMedium.cpp +++ b/tests/media/testMedium.cpp @@ -92,18 +92,21 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") { CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m))); medium.getNuclearComposition(); + SpeedType const speed = 1_m / second; + // create a line of length 1 m - Line const line(gOrigin, - VelocityVector(gCS, {1_m / second, 0_m / second, 0_m / second})); + Line const line(gOrigin, VelocityVector(gCS, {speed, 0_m / second, 0_m / second})); // the end time of our line auto const tEnd = 1_s; + LengthType const length = tEnd * speed; + // and the associated trajectory setup::Trajectory const trajectory = corsika::setup::testing::make_track(line, tEnd); // and check the integrated grammage - CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1)); + CHECK((medium.getIntegratedGrammage(trajectory) / (density * length)) == Approx(1)); CHECK((medium.getArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1)); } diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp index acef32c4f..62422f7a8 100644 --- a/tests/media/testRefractiveIndex.cpp +++ b/tests/media/testRefractiveIndex.cpp @@ -76,19 +76,22 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") { CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m))); medium.getNuclearComposition(); + SpeedType const speed = 1_m / second; + // create a line of length 1 m - Line const line(gOrigin, - VelocityVector(gCS, {1_m / second, 0_m / second, 0_m / second})); + Line const line(gOrigin, VelocityVector(gCS, {speed, 0_m / second, 0_m / second})); // the end time of our line auto const tEnd = 1_s; + LengthType const length = tEnd * speed; + // and the associated trajectory setup::Trajectory const track = setup::testing::make_track(line, tEnd); // and check the integrated grammage - CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); + CHECK((medium.getIntegratedGrammage(track) / (density * length)) == Approx(1)); CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); } @@ -148,24 +151,26 @@ TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") { REQUIRE(density == medium__.getMassDensity(Point(gCS, 0_m, 0_m, 0_m))); medium__.getNuclearComposition(); - // create a line of length 1 m - Line const line(gOrigin, Vector( - gCS, {1_m / second, 0_m / second, 0_m / second})); + SpeedType const velocity = 1_m / second; // the end time of our line - auto const tEnd = 1_s; + TimeType const tEnd = 1_s; + + LengthType const length = tEnd * velocity; + + // create a line of length 1 m + Line const line(gOrigin, Vector( + gCS, {velocity, 0_m / second, 0_m / second})); // and the associated trajectory setup::Trajectory const track = setup::testing::make_track(line, tEnd); - // // and the associated trajectory - // Trajectory const trajectory(line, tEnd); // and check the integrated grammage - REQUIRE((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); + REQUIRE((medium.getIntegratedGrammage(track) / (density * length)) == Approx(1)); REQUIRE((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); - REQUIRE((medium_.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); + REQUIRE((medium_.getIntegratedGrammage(track) / (density * length)) == Approx(1)); REQUIRE((medium_.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); - REQUIRE((medium__.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1)); + REQUIRE((medium__.getIntegratedGrammage(track) / (density * length)) == Approx(1)); REQUIRE((medium__.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1)); } -- GitLab From 8dbeb7b80df07b8d455c32483b67f59e3550e026 Mon Sep 17 00:00:00 2001 From: ralfulrich Date: Fri, 18 Jun 2021 13:40:33 +0200 Subject: [PATCH 2/8] updates and fixes --- .../LayeredSphericalAtmosphereBuilder.inl | 59 ++-- .../detail/media/SlidingPlanarExponential.inl | 55 ++-- corsika/media/BaseExponential.hpp | 2 - .../LayeredSphericalAtmosphereBuilder.hpp | 17 +- corsika/media/SlidingPlanarExponential.hpp | 2 - examples/CMakeLists.txt | 2 + examples/corsika.cpp | 47 +++- examples/em_shower.cpp | 4 +- examples/hybrid_MC.cpp | 4 +- examples/vertical_EAS.cpp | 4 +- tests/media/testEnvironment.cpp | 258 +++++++++++++++++- 11 files changed, 368 insertions(+), 86 deletions(-) diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl index 791d57972..17371ff3b 100644 --- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl +++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl @@ -13,7 +13,7 @@ #include #include #include -//#include +#include namespace corsika { @@ -42,7 +42,7 @@ namespace corsika { TModelArgs...>::addExponentialLayer(GrammageType b, LengthType c, LengthType upperBoundary) { - auto const radius = earthRadius_ + upperBoundary; + auto const radius = planetRadius_ + upperBoundary; checkRadius(radius); previousRadius_ = radius; @@ -56,7 +56,7 @@ namespace corsika { auto lastBound = [&](auto... argPack) { return std::make_shared< TMediumModelExtra>>( - argPack..., center_, rho0, -c, *composition_, earthRadius_); + argPack..., center_, rho0, -c, *composition_, planetRadius_); }; // now unpack the additional arguments @@ -64,7 +64,7 @@ namespace corsika { node->setModelProperties(std::move(model)); } else { node->template setModelProperties>( - center_, rho0, -c, *composition_, earthRadius_); + center_, rho0, -c, *composition_, planetRadius_); } layers_.push(std::move(node)); @@ -75,14 +75,14 @@ namespace corsika { inline void LayeredSphericalAtmosphereBuilder< TMediumInterface, TMediumModelExtra, TModelArgs...>::addLinearLayer(LengthType c, LengthType upperBoundary) { - auto const radius = earthRadius_ + upperBoundary; + auto const radius = planetRadius_ + upperBoundary; checkRadius(radius); previousRadius_ = radius; auto node = std::make_unique>( std::make_unique(center_, radius)); - units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm); + units::si::GrammageType constexpr b = 1_g / (1_cm * 1_cm); auto const rho0 = b / c; if constexpr (detail::has_extra_models::value) { @@ -94,7 +94,6 @@ namespace corsika { // now unpack the additional arguments auto model = std::apply(lastBound, additionalModelArgs_); - node->setModelProperties(std::move(model)); } else { node->template setModelProperties>( @@ -108,33 +107,33 @@ namespace corsika { typename... TModelArgs> inline void LayeredSphericalAtmosphereBuilder:: - addTabularLayer(GrammageType b, - std::function const& funcRho, - LengthType upperBoundary) { + addTabularLayer(std::function const& funcRho, + unsigned int const nBins, LengthType const deltaHeight, + LengthType const upperBoundary) { - auto const radius = earthRadius_ + upperBoundary; + auto const radius = planetRadius_ + upperBoundary; checkRadius(radius); previousRadius_ = radius; auto node = std::make_unique>( std::make_unique(center_, radius)); - /* - if constexpr (detail::has_extra_models::value) { - // helper lambda in which the last 5 arguments to make_shared<...> are bound - auto lastBound = [&](auto... argPack) { - return std::make_shared< - TMediumModelExtra>>( - argPack..., center_, rho0, -c, *composition_, earthRadius_); - }; - - // now unpack the additional arguments - auto model = std::apply(lastBound, additionalModelArgs_); - node->setModelProperties(std::move(model)); - } else { - node->template setModelProperties>( - center_, rho0, -c, *composition_, earthRadius_); - } - */ + + if constexpr (detail::has_extra_models::value) { + // helper lambda in which the last 5 arguments to make_shared<...> are bound + auto lastBound = [&](auto... argPack) { + return std::make_shared< + TMediumModelExtra>>( + argPack..., center_, funcRho, nBins, deltaHeight, *composition_, + planetRadius_); + }; + + // now unpack the additional arguments + auto model = std::apply(lastBound, additionalModelArgs_); + node->setModelProperties(std::move(model)); + } else { + node->template setModelProperties>( + center_, funcRho, nBins, deltaHeight, *composition_, planetRadius_); + } layers_.push(std::move(node)); } @@ -167,10 +166,10 @@ namespace corsika { template typename MExtraEnvirnoment> struct make_layered_spherical_atmosphere_builder { template - static auto create(Point const& center, LengthType earthRadius, TArgs... args) { + static auto create(Point const& center, LengthType planetRadius, TArgs... args) { return LayeredSphericalAtmosphereBuilder{std::forward(args)..., - center, earthRadius}; + center, planetRadius}; } }; diff --git a/corsika/detail/media/SlidingPlanarExponential.inl b/corsika/detail/media/SlidingPlanarExponential.inl index 80fc9bf6b..9b9218728 100644 --- a/corsika/detail/media/SlidingPlanarExponential.inl +++ b/corsika/detail/media/SlidingPlanarExponential.inl @@ -8,50 +8,51 @@ #pragma once -#include - namespace corsika { - template - inline SlidingPlanarExponential::SlidingPlanarExponential( + template + inline SlidingPlanarExponential::SlidingPlanarExponential( Point const& p0, MassDensityType rho0, LengthType lambda, NuclearComposition const& nuclComp, LengthType referenceHeight) - : BaseExponential>(p0, referenceHeight, rho0, lambda) - , nuclComp_(nuclComp) - , referenceHeight_(referenceHeight) {} + : BaseExponential>(p0, referenceHeight, rho0, + lambda) + , nuclComp_(nuclComp) {} - template - inline MassDensityType SlidingPlanarExponential::getMassDensity( + template + inline MassDensityType SlidingPlanarExponential::getMassDensity( Point const& point) const { auto const heightFull = - (point - BaseExponential>::getAnchorPoint()) + (point - BaseExponential>::getAnchorPoint()) .getNorm(); - return BaseExponential>::getMassDensity(heightFull); + return BaseExponential>::getMassDensity( + heightFull); } - template - inline NuclearComposition const& SlidingPlanarExponential::getNuclearComposition() - const { + template + inline NuclearComposition const& + SlidingPlanarExponential::getNuclearComposition() const { return nuclComp_; } - template - inline GrammageType SlidingPlanarExponential::getIntegratedGrammage( + template + inline GrammageType SlidingPlanarExponential::getIntegratedGrammage( BaseTrajectory const& traj) const { - auto const axis = (traj.getPosition(0) - - BaseExponential>::getAnchorPoint()) - .normalized(); - return BaseExponential>::getIntegratedGrammage(traj, - axis); + auto const axis = + (traj.getPosition(0) - + BaseExponential>::getAnchorPoint()) + .normalized(); + return BaseExponential>::getIntegratedGrammage( + traj, axis); } - template - inline LengthType SlidingPlanarExponential::getArclengthFromGrammage( + template + inline LengthType SlidingPlanarExponential::getArclengthFromGrammage( BaseTrajectory const& traj, GrammageType const grammage) const { - auto const axis = (traj.getPosition(0) - - BaseExponential>::getAnchorPoint()) - .normalized(); - return BaseExponential>::getArclengthFromGrammage( + auto const axis = + (traj.getPosition(0) - + BaseExponential>::getAnchorPoint()) + .normalized(); + return BaseExponential>::getArclengthFromGrammage( traj, grammage, axis); } diff --git a/corsika/media/BaseExponential.hpp b/corsika/media/BaseExponential.hpp index 27a267df6..507ed640f 100644 --- a/corsika/media/BaseExponential.hpp +++ b/corsika/media/BaseExponential.hpp @@ -29,8 +29,6 @@ namespace corsika { MassDensityType rho0, LengthType lambda); Point const& getAnchorPoint() const { return point_; } - MassDensityType getRho0() const { return rho0_; } - InverseLengthType getInvLambda() const { return invLambda_; } protected: auto const& getImplementation() const; diff --git a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp index 41ad15801..6c178f458 100644 --- a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp +++ b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp @@ -22,6 +22,7 @@ #include #include #include +#include namespace corsika { @@ -68,9 +69,9 @@ namespace corsika { protected: LayeredSphericalAtmosphereBuilder(TModelArgs... args, Point const& center, - LengthType earthRadius) + LengthType planetRadius) : center_(center) - , earthRadius_(earthRadius) + , planetRadius_(planetRadius) , additionalModelArgs_{args...} {} public: @@ -78,9 +79,9 @@ namespace corsika { void addExponentialLayer(GrammageType b, LengthType c, LengthType upperBoundary); void addLinearLayer(LengthType c, LengthType upperBoundary); - void addTabularLayer(GrammageType b, - std::function const& funcRho, - LengthType upperBoundary); + void addTabularLayer(std::function const& funcRho, + unsigned int const nBins, LengthType const deltaHeight, + LengthType const upperBoundary); int getSize() const { return layers_.size(); } @@ -88,9 +89,9 @@ namespace corsika { Environment assemble(); /** - * Get the current Earth radius. + * Get the current planet radius. */ - LengthType getEarthRadius() const { return earthRadius_; } + LengthType getPlanetRadius() const { return planetRadius_; } private: void checkRadius(LengthType r) const; @@ -98,7 +99,7 @@ namespace corsika { std::unique_ptr composition_; Point center_; LengthType previousRadius_{LengthType::zero()}; - LengthType earthRadius_; + LengthType planetRadius_; std::tuple const additionalModelArgs_; std::stack::VTNUPtr> diff --git a/corsika/media/SlidingPlanarExponential.hpp b/corsika/media/SlidingPlanarExponential.hpp index 63f8238ba..4ef91bb4e 100644 --- a/corsika/media/SlidingPlanarExponential.hpp +++ b/corsika/media/SlidingPlanarExponential.hpp @@ -13,7 +13,6 @@ #include #include #include -#include #include #include @@ -53,7 +52,6 @@ namespace corsika { private: NuclearComposition const nuclComp_; - LengthType const referenceHeight_; }; } // namespace corsika diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index e3d056e22..70816c0d3 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -60,4 +60,6 @@ CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.) add_executable (corsika corsika.cpp) target_link_libraries (corsika CORSIKA8::CORSIKA8) +add_executable (mars mars.cpp) +target_link_libraries (mars CORSIKA8::CORSIKA8) diff --git a/examples/corsika.cpp b/examples/corsika.cpp index fd048470d..ee424c1de 100644 --- a/examples/corsika.cpp +++ b/examples/corsika.cpp @@ -103,8 +103,6 @@ using MyExtraEnv = MediumPropertyModel>; int main(int argc, char** argv) { - logging::set_level(logging::level::info); - // the main command line description CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."}; @@ -154,10 +152,31 @@ int main(int argc, char** argv) { ->group("Misc."); app.add_flag("--force-interaction", "Force the location of the first interaction.") ->group("Misc."); + app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.") + ->default_val("info") + ->check(CLI::IsMember({"warn", "info", "debug", "trace"})) + ->group("Misc."); // parse the command line options into the variables CLI11_PARSE(app, argc, argv); + if (app.count("--verbosity")) { + string const loglevel = app["verbosity"]->as(); + if (loglevel == "warn") { + logging::set_level(logging::level::warn); + } else if (loglevel == "info") { + logging::set_level(logging::level::info); + } else if (loglevel == "debug") { + logging::set_level(logging::level::debug); + } else if (loglevel == "trace") { +#ifndef DEBUG + CORSIKA_LOG_ERROR("trace log level requires a Debug build."); + return 1; +#endif + logging::set_level(logging::level::trace); + } + } + // check that we got either PDG or A/Z // this can be done with option_groups but the ordering // gets all messed up @@ -185,17 +204,29 @@ int main(int argc, char** argv) { 50_uT, 0_T}); builder.setNuclearComposition( {{Code::Nitrogen, Code::Oxygen}, - {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now + {0.7847, 1. - 0.7847}}); // 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); builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km); - builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean); + builder.addLinearLayer(1e9_cm, 112.8_km); builder.assemble(env); /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */ + ofstream atmout("earth.dat"); + for (LengthType h = 0_m; h < 110_km; h += 100_m) { + Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h}; + auto rho = + env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity( + ptest); + atmout << h / 1_m << " " << rho / 1_kg * cube(1_m) << "\n"; + } + atmout.close(); + + /* === START: CONSTRUCT PRIMARY PARTICLE === */ + /* === START: CONSTRUCT PRIMARY PARTICLE === */ // parse the primary ID as a PDG or A/Z code @@ -231,8 +262,8 @@ int main(int argc, char** argv) { /* === END: CONSTRUCT PRIMARY PARTICLE === */ /* === START: CONSTRUCT GEOMETRY === */ - auto const observationHeight = 0_km + builder.getEarthRadius(); - auto const injectionHeight = 111.75_km + builder.getEarthRadius(); + auto const observationHeight = 0_km + builder.getPlanetRadius(); + auto const injectionHeight = 111.75_km + builder.getPlanetRadius(); auto const t = -observationHeight * cos(thetaRad) + sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) + static_pow<2>(injectionHeight)); @@ -281,9 +312,9 @@ int main(int argc, char** argv) { Code::KStar0_1430_MinusBar, }}; - decaySibyll.printDecayConfig(); + // decaySibyll.printDecayConfig(); - ParticleCut cut{50_GeV, 50_GeV, 50_GeV, 50_GeV, false}; + ParticleCut cut{1_GeV, 1_GeV, 1_GeV, 1_GeV, false}; corsika::proposal::Interaction emCascade(env); corsika::proposal::ContinuousProcess emContinuous(env); InteractionCounter emCascadeCounted(emCascade); diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp index 0ff4d28b7..d8b8385f3 100644 --- a/examples/em_shower.cpp +++ b/examples/em_shower.cpp @@ -128,8 +128,8 @@ int main(int argc, char** argv) { cout << "input momentum: " << plab.getComponents() / 1_GeV << ", norm = " << plab.getNorm() << endl; - auto const observationHeight = 1.4_km + builder.getEarthRadius(); - auto const injectionHeight = 112.75_km + builder.getEarthRadius(); + auto const observationHeight = 1.4_km + builder.getPlanetRadius(); + auto const injectionHeight = 112.75_km + builder.getPlanetRadius(); auto const t = -observationHeight * cos(thetaRad) + sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) + static_pow<2>(injectionHeight)); diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp index e065c4ee8..4674512cc 100644 --- a/examples/hybrid_MC.cpp +++ b/examples/hybrid_MC.cpp @@ -178,8 +178,8 @@ int main(int argc, char** argv) { cout << "input momentum: " << plab.getComponents() / 1_GeV << ", norm = " << plab.getNorm() << endl; - auto const observationHeight = 0_km + builder.getEarthRadius(); - auto const injectionHeight = 112.75_km + builder.getEarthRadius(); + auto const observationHeight = 0_km + builder.getPlanetRadius(); + auto const injectionHeight = 112.75_km + builder.getPlanetRadius(); auto const t = -observationHeight * cos(thetaRad) + sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) + static_pow<2>(injectionHeight)); diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp index 6ae766bb1..6800a84be 100644 --- a/examples/vertical_EAS.cpp +++ b/examples/vertical_EAS.cpp @@ -177,8 +177,8 @@ int main(int argc, char** argv) { 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 observationHeight = 0_km + builder.getPlanetRadius(); + auto const injectionHeight = 111.75_km + builder.getPlanetRadius(); auto const t = -observationHeight * cos(thetaRad) + sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) + static_pow<2>(injectionHeight)); diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp index 88e143b7b..fffdff7fe 100644 --- a/tests/media/testEnvironment.cpp +++ b/tests/media/testEnvironment.cpp @@ -11,6 +11,7 @@ #include #include #include + #include #include #include @@ -26,6 +27,7 @@ #include #include #include +#include #include #include @@ -94,7 +96,7 @@ TEST_CASE("FlatExponential") { Vector const axis(gCS, QuantityVector(0, 0, 1)); LengthType const lambda = 3_m; - auto const rho0 = 1_g / static_pow<3>(1_cm); + auto const rho0 = 1_g / cube(1_cm); FlatExponential const medium(gOrigin, axis, rho0, lambda, protonComposition); SpeedType const speed = 20_m / second; @@ -188,6 +190,256 @@ TEST_CASE("SlidingPlanarExponential") { } } +struct RhoFuncConst { + MassDensityType operator()(LengthType) const { return 1_g / cube(1_cm); } + static GrammageType integrate(LengthType dL) { return dL * 1_g / cube(1_cm); } +}; + +struct RhoFuncExp { + MassDensityType operator()(LengthType height) const { + return 1_g / cube(1_cm) * exp(-height / 1000_m); + } + static GrammageType integrate(BaseTrajectory const& traj, Point const& origin, + LengthType const& refH) { + LengthType height1 = (traj.getPosition(0) - origin).getNorm() - refH; + LengthType height2 = (traj.getPosition(1) - origin).getNorm() - refH; + if (height1 > height2) { std::swap(height1, height2); } + + DirectionVector const axis( + (traj.getPosition(0) - origin).normalized()); // to gravity center + double const cosTheta = axis.dot(traj.getDirection(0)); + + CORSIKA_LOG_INFO("h1={} h2={} cT={} rho1={}, rho2={}", height1, height2, cosTheta, + 1_g / cube(1_cm) * exp(-height1 / 1000_m), + 1_g / cube(1_cm) * exp(-height2 / 1000_m)); + return (1_km * 1_g / cube(1_cm) * exp(-height1 / 1000_m) - + 1_km * 1_g / cube(1_cm) * exp(-height2 / 1000_m)) / + cosTheta; + } + static GrammageType integrate(BaseTrajectory const& traj, LengthType const& length, + Point const& origin, LengthType const& refH) { + LengthType height1 = (traj.getPosition(0) - origin).getNorm() - refH; + LengthType height2 = + (traj.getPosition(0) + traj.getDirection(0) * length - origin).getNorm() - refH; + if (height1 > height2) { std::swap(height1, height2); } + + DirectionVector const axis( + (traj.getPosition(0) - origin).normalized()); // to gravity center + double const cosTheta = axis.dot(traj.getDirection(0)); + + CORSIKA_LOG_INFO("h1={} h2={} cT={}", height1, height2, cosTheta); + return (1_km * 1_g / cube(1_cm) * exp(-height1 / 1000_m) - + 1_km * 1_g / cube(1_cm) * exp(-height2 / 1000_m)) / + cosTheta; + } +}; + +TEST_CASE("SlidingPlanarTabular") { + + logging::set_level(logging::level::info); + + NuclearComposition const protonComposition(std::vector{Code::Proton}, + std::vector{1.f}); + + RhoFuncConst rhoFunc; + SlidingPlanarTabular const medium(gOrigin, rhoFunc, 1000, 10_m, + protonComposition); + + SECTION("density") { + CHECK(medium.getMassDensity({gCS, {0_m, 0_m, 3_m}}) / + medium.getMassDensity({gCS, {0_m, 3_m, 0_m}}) == + Approx(1)); + CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 3_m}}) == 1_g / cube(1_cm)); + CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 300_m}}) == 1_g / cube(1_cm)); + } + + SECTION("vertical") { + SpeedType const speed = 5_m / second; + TimeType const tEnd1 = 1_s; + LengthType const length1 = speed * tEnd1; + TimeType const tEnd2 = 300_s; + LengthType const length2 = speed * tEnd2; + Line const line( + {gCS, {0_m, 0_m, 1_m}}, + Vector(gCS, {0_m / second, 0_m / second, speed})); + setup::Trajectory const trajectory1 = + setup::testing::make_track(line, tEnd1); + Line const line1Reverse( + trajectory1.getPosition(1), + Vector(gCS, {0_m / second, 0_m / second, -speed})); + setup::Trajectory const trajectory1Reverse = + setup::testing::make_track(line1Reverse, tEnd1); + + setup::Trajectory const trajectory2 = + setup::testing::make_track(line, tEnd2); + Line const line2Reverse( + trajectory2.getPosition(0), + Vector(gCS, {0_m / second, 0_m / second, -speed})); + setup::Trajectory const trajectory2Reverse = + setup::testing::make_track(line2Reverse, tEnd2); + + // failures + CHECK_THROWS(medium.getArclengthFromGrammage(trajectory1, -1_kg / square(1_cm))); + + MassDensityType const rho0 = 1_g / cube(1_cm); + + // short track + CHECK(medium.getIntegratedGrammage(trajectory1) == length1 * rho0); + LengthType const testD1 = length1 / 200; // within bin + CHECK(medium.getArclengthFromGrammage(trajectory1, rho0 * testD1) / testD1 == + Approx(1)); + // short track, reverse + CHECK(medium.getIntegratedGrammage(trajectory1Reverse) == length1 * rho0); + CHECK(medium.getArclengthFromGrammage(trajectory1Reverse, rho0 * testD1) / testD1 == + Approx(1)); + + // long track + CHECK(medium.getIntegratedGrammage(trajectory2) == length2 * 1_g / cube(1_cm)); + LengthType const testD2 = length2 / 25; // multi bin + CHECK(medium.getArclengthFromGrammage(trajectory2, rho0 * testD2) == testD2); + } + + SECTION("inclined") { + SpeedType const speed = 5_m / second; + TimeType const tEnd1 = 1_s; + LengthType const length1 = speed * tEnd1; + TimeType const tEnd2 = 300_s; + LengthType const length2 = speed * tEnd2; + Line const line({gCS, {0_m, 0_m, 1_m}}, + Vector( + gCS, {speed / sqrt(2.), 0_m / second, speed / sqrt(2.)})); + setup::Trajectory const trajectory1 = + setup::testing::make_track(line, tEnd1); + Line const line1Reverse( + trajectory1.getPosition(1), + Vector( + gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)})); + setup::Trajectory const trajectory1Reverse = + setup::testing::make_track(line1Reverse, tEnd1); + + setup::Trajectory const trajectory2 = + setup::testing::make_track(line, tEnd2); + Line const line2Reverse( + trajectory2.getPosition(1), + Vector( + gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)})); + setup::Trajectory const trajectory2Reverse = + setup::testing::make_track(line2Reverse, tEnd2); + + MassDensityType const rho0 = 1_g / cube(1_cm); + + // short track + CHECK(medium.getIntegratedGrammage(trajectory1) / (length1 * rho0) == Approx(1)); + LengthType const testD1 = length1 / 200; // within bin + CHECK(medium.getArclengthFromGrammage(trajectory1, RhoFuncConst::integrate(testD1)) / + testD1 == + Approx(1)); + // short track, reverse + CHECK(medium.getIntegratedGrammage(trajectory1Reverse) / (length1 * rho0) == + Approx(1)); + CHECK(medium.getArclengthFromGrammage(trajectory1Reverse, rho0 * testD1) / testD1 == + Approx(1)); + + // long track + CHECK(medium.getIntegratedGrammage(trajectory2) / (length2 * rho0) == + Approx(1).epsilon(0.01)); + LengthType const testD2 = length2 / 25; // multi bin + CHECK(medium.getArclengthFromGrammage(trajectory2, rho0 * testD2) / testD2 == + Approx(1).epsilon(0.01)); + // long track reverse + CORSIKA_LOG_INFO("length2={}", length2); + CHECK(medium.getIntegratedGrammage(trajectory2Reverse) / (length2 * rho0) == + Approx(1).epsilon(0.01)); + CHECK(medium.getArclengthFromGrammage(trajectory2Reverse, rho0 * testD2) / testD2 == + Approx(1).epsilon(0.01)); + } + + /*The exponential test is taken over phase-space where the exponential is not so steep + * and is samples in sufficient substeps. An reference-height offset of 1000_km is used. + * Thus, density is given from 1000 to 1010 km. And curvature effects are small. + */ + + RhoFuncExp rhoFuncExp; + SlidingPlanarTabular const mediumExp(gOrigin, rhoFuncExp, 1000, 10_m, + protonComposition, 1000_km); + + SECTION("exponential") { + + SpeedType const speed = 5_m / second; + TimeType const tEnd1 = 1_s; + LengthType const length1 = speed * tEnd1; + TimeType const tEnd2 = 300_s; + LengthType const length2 = speed * tEnd2; + Line const line({gCS, {0_m, 0_m, 1000.005_km}}, + Vector( + gCS, {speed / sqrt(2.), 0_m / second, speed / sqrt(2.)})); + setup::Trajectory const trajectory1 = + setup::testing::make_track(line, tEnd1); + Line const line1Reverse( + trajectory1.getPosition(1), + Vector( + gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)})); + setup::Trajectory const trajectory1Reverse = + setup::testing::make_track(line1Reverse, tEnd1); + + setup::Trajectory const trajectory2 = + setup::testing::make_track(line, tEnd2); + + CORSIKA_LOG_INFO("{} {}", RhoFuncExp::integrate(trajectory1, gOrigin, 1000_km), + length1); + + // short track + GrammageType const testShortX = RhoFuncExp::integrate(trajectory1, gOrigin, 1000_km); + CHECK(mediumExp.getIntegratedGrammage(trajectory1) / testShortX == + Approx(1).epsilon(0.01)); + LengthType const testD1 = length1 / 200; // within bin + GrammageType const testD1X = + RhoFuncExp::integrate(trajectory1, testD1, gOrigin, 1000_km); + CHECK(mediumExp.getArclengthFromGrammage(trajectory1, testD1X) / testD1 == + Approx(1).epsilon(0.01)); + // short track, reverse + CHECK(mediumExp.getIntegratedGrammage(trajectory1Reverse) / testShortX == + Approx(1).epsilon(0.01)); + CHECK(mediumExp.getArclengthFromGrammage(trajectory1Reverse, testD1X) / testD1 == + Approx(1).epsilon(0.01)); + + // long track + GrammageType const testLongX = RhoFuncExp::integrate(trajectory2, gOrigin, 1000_km); + CORSIKA_LOG_INFO("testLongX={}", testLongX); + CHECK(mediumExp.getIntegratedGrammage(trajectory2) / testLongX == + Approx(1).epsilon(0.01)); + LengthType const testD2 = length2 / 25; // multi bin + GrammageType const testD2X = + RhoFuncExp::integrate(trajectory2, testD2, gOrigin, 1000_km); + CHECK(mediumExp.getArclengthFromGrammage(trajectory2, testD2X) / testD2 == + Approx(1).epsilon(0.01)); + // long track, reverse + + // first full trajectory2 reverse + Line line2Reverse(trajectory2.getPosition(1), + Vector( + gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)})); + setup::Trajectory trajectory2Reverse = + setup::testing::make_track(line2Reverse, tEnd2); + + CHECK(mediumExp.getIntegratedGrammage(trajectory2Reverse) / testLongX == + Approx(1).epsilon(0.01)); + + // but now shorter trajectory2 reversed to correspond 100% to testD2 + + line2Reverse = Line(trajectory2.getPosition(0) + trajectory2.getDirection(0) * testD2, + Vector( + gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)})); + auto const trajectory2ReverseShort = + setup::testing::make_track(line2Reverse, testD2 / speed); + + CORSIKA_LOG_INFO("here {} {} {}", trajectory2ReverseShort.getLength(), testD2, + testD2X / 1_g * square(1_cm)); + CHECK(mediumExp.getArclengthFromGrammage(trajectory2ReverseShort, testD2X) / testD2 == + Approx(1).epsilon(0.01)); + } +} + MassDensityType constexpr rho0 = 1_kg / 1_m / 1_m / 1_m; struct ExponentialTest { @@ -289,7 +541,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") { CHECK(builder.getSize() == 0); - auto const R = builder.getEarthRadius(); + auto const R = builder.getPlanetRadius(); CHECK(univ->getChildNodes().size() == 1); @@ -334,7 +586,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") { CHECK(builder.getSize() == 0); CHECK(univ->getChildNodes().size() == 1); - auto const R = builder.getEarthRadius(); + auto const R = builder.getPlanetRadius(); // check magnetic field at several locations const Point pTest(gCS, -10_m, 4_m, R + 35_m); -- GitLab From 920af65e81b3232dd28d93b75df55c0d6e30ffac Mon Sep 17 00:00:00 2001 From: ralfulrich Date: Fri, 18 Jun 2021 13:41:15 +0200 Subject: [PATCH 3/8] added Mars --- examples/mars.cpp | 461 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 461 insertions(+) create mode 100644 examples/mars.cpp diff --git a/examples/mars.cpp b/examples/mars.cpp new file mode 100644 index 000000000..ac25be75c --- /dev/null +++ b/examples/mars.cpp @@ -0,0 +1,461 @@ +/* + * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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. + */ + +/* clang-format off */ +// InteractionCounter used boost/histogram, which +// fails if boost/type_traits have been included before. Thus, we have +// to include it first... +#include +/* clang-format on */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +/* + NOTE, WARNING, ATTENTION + + The .../Random.hpp implement the hooks of external modules to the C8 random + number generator. It has to occur excatly ONCE per linked + executable. If you include the header below multiple times and + link this togehter, it will fail. + */ +#include + +using namespace corsika; +using namespace std; + +using Particle = setup::Stack::particle_type; + +typedef decltype(1 * pascal) PressureType; +typedef decltype(1 * degree_celsius) TemperatureType; + +class MarsAtmModel { +public: + MarsAtmModel() = delete; + MarsAtmModel(PressureType a, InverseLengthType b, TemperatureType c, + decltype(1 * degree_celsius / 1_m) d) + : a_(a) + , b_(b) + , c_(c) + , d_(d) {} + + MassDensityType operator()(LengthType height) const { + PressureType const pressure = a_ * exp(-b_ * height); + TemperatureType const temperature = -c_ - d_ * height + 273.1_K; // in K + constexpr decltype(square(1_m) / (square(1_s) * 1_K)) constant = + 1000 * 0.1921 * square(1_m) / (square(1_s) * 1_K); + return pressure / (constant * temperature); + } + +private: + PressureType a_; + InverseLengthType b_; + TemperatureType c_; + decltype(1_K / 1_m) d_; +}; + +void registerRandomStreams(int seed) { + RNGManager<>::getInstance().registerRandomStream("cascade"); + RNGManager<>::getInstance().registerRandomStream("qgsjet"); + RNGManager<>::getInstance().registerRandomStream("sibyll"); + RNGManager<>::getInstance().registerRandomStream("pythia"); + RNGManager<>::getInstance().registerRandomStream("urqmd"); + RNGManager<>::getInstance().registerRandomStream("proposal"); + if (seed == 0) { + std::random_device rd; + seed = rd(); + cout << "new random seed (auto) " << seed << endl; + } + RNGManager<>::getInstance().setSeed(seed); +} + +template +using MyExtraEnv = MediumPropertyModel>; + +// argv : 1.number of nucleons, 2.number of protons, +// 3.total energy in GeV, 4.number of showers, +// 5.seed (0 by default to generate random values for all) + +int main(int argc, char** argv) { + + // the main command line description + CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."}; + + // some options that we want to fill in + int A, Z, nevent = 0; + + // the following section adds the options to the parser + + // we start by definining a sub-group for the primary ID + auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary") + ->check(CLI::Range(0, 26)) + ->group("Primary"); + auto opt_A = app.add_option("-A", A, "Atomic mass number for primary") + ->needs(opt_Z) + ->check(CLI::Range(1, 58)) + ->group("Primary"); + app.add_option("-p,--pdg", "PDG code for primary.") + ->excludes(opt_A) + ->excludes(opt_Z) + ->group("Primary"); + // the remainding options + app.add_option("-E,--energy", "Primary energy in GeV") + ->required() + ->check(CLI::PositiveNumber) + ->group("Primary"); + app.add_option("-z,--zenith", "Primary zenith angle (deg)") + ->required() + ->default_val(0.) + ->check(CLI::Range(0, 90)) + ->group("Primary"); + app.add_option("-a,--azimuth", "Primary azimuth angle (deg)") + ->default_val(0.) + ->check(CLI::Range(0, 360)) + ->group("Primary"); + app.add_option("-N,--nevent", nevent, "The number of events/showers to run.") + ->required() + ->check(CLI::PositiveNumber) + ->group("Library/Output"); + app.add_option("-f,--filename", "Filename for output library.") + ->required() + ->default_val("corsika_library") + ->check(CLI::NonexistentPath) + ->group("Library/Output"); + app.add_option("-s,--seed", "The random number seed.") + ->default_val(12351739) + ->check(CLI::NonNegativeNumber) + ->group("Misc."); + app.add_flag("--force-interaction", "Force the location of the first interaction.") + ->group("Misc."); + app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.") + ->default_val("info") + ->check(CLI::IsMember({"warn", "info", "debug", "trace"})) + ->group("Misc."); + + // parse the command line options into the variables + CLI11_PARSE(app, argc, argv); + + if (app.count("--verbosity")) { + string const loglevel = app["verbosity"]->as(); + if (loglevel == "warn") { + logging::set_level(logging::level::warn); + } else if (loglevel == "info") { + logging::set_level(logging::level::info); + } else if (loglevel == "debug") { + logging::set_level(logging::level::debug); + } else if (loglevel == "trace") { +#ifndef DEBUG + CORSIKA_LOG_ERROR("trace log level requires a Debug build."); + return 1; +#endif + logging::set_level(logging::level::trace); + } + } + + // check that we got either PDG or A/Z + // this can be done with option_groups but the ordering + // gets all messed up + if (app.count("--pdg") == 0) { + if ((app.count("-A") == 0) || (app.count("-Z") == 0)) { + std::cerr << "If --pdg is not provided, then both -A and -Z are required." + << std::endl; + return 1; + } + } + + // initialize random number sequence(s) + registerRandomStreams(app["--seed"]->as()); + + /* === START: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */ + using EnvType = setup::Environment; + EnvType env; + CoordinateSystemPtr const& rootCS = env.getCoordinateSystem(); + Point const center{rootCS, 0_m, 0_m, 0_m}; + LengthType const radiusMars = 3389.5_km; + auto builder = + make_layered_spherical_atmosphere_builder:: + create(center, + radiusMars, // Mars + Medium::AirDry1Atm, // Mars, close enough + MagneticFieldVector{rootCS, 0_T, 0_uT, 0_T}); // Mars + + builder.setNuclearComposition( // Mars + {{Code::Nitrogen, Code::Oxygen}, {1. / 3., 2. / 3.}}); // simplified + //{{Code::Carbon, Code::Oxygen, // 95.97 CO2 + // Code::Nitrogen}, // 1.89 N2 + 1.93 Argon + 0.146 O2 + // {0.9597 / 3, 0.9597 * 2 / 3, + // 1 - 0.9597}}); // values taken from AIRES manual, Ar removed for now + + MarsAtmModel layer1(0.699e3 * pascal, 0.00009 / 1_m, 31.0 * degree_celsius, + 0.000998 * 1 * degree_celsius / 1_m); + MarsAtmModel layer2(0.699e3 * pascal, 0.00009 / 1_m, 23.4 * degree_celsius, + 0.00222 * 1 * degree_celsius / 1_m); + + builder.addTabularLayer(layer1, 100, 100_m, 7_km); + builder.addTabularLayer(layer2, 300, 500_m, 100_km); + builder.addLinearLayer(1e9_cm, 112.8_km); + builder.assemble(env); + /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */ + + ofstream atmout("mars.dat"); + for (LengthType h = 0_m; h < 110_km; h += 100_m) { + Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h}; + auto rho = + env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity( + ptest); + atmout << h / 1_m << " " << rho / 1_kg * cube(1_m) << "\n"; + } + atmout.close(); + + /* === START: CONSTRUCT PRIMARY PARTICLE === */ + + // parse the primary ID as a PDG or A/Z code + Code beamCode; + HEPEnergyType mass; + + // check if we want to use a PDG code instead + if (app.count("--pdg") > 0) { + beamCode = convert_from_PDG(PDGCode(app["--pdg"]->as())); + mass = get_mass(beamCode); + } else { + // check manually for proton and neutrons + if ((A == 0) && (Z == 1)) beamCode = Code::Proton; + if ((A == 1) && (Z == 1)) beamCode = Code::Neutron; + mass = get_nucleus_mass(A, Z); + } + + // particle energy + HEPEnergyType const E0 = 1_GeV * app["--energy"]->as(); + + // direction of the shower in (theta, phi) space + auto const thetaRad = app["--zenith"]->as() / 180. * M_PI; + auto const phiRad = app["--azimuth"]->as() / 180. * M_PI; + + // convert Elab to Plab + HEPMomentumType P0 = sqrt((E0 - mass) * (E0 + mass)); + + // convert the momentum to the zenith and azimuth angle of the primary + auto const [px, py, pz] = + std::make_tuple(P0 * sin(thetaRad) * cos(phiRad), P0 * sin(thetaRad) * sin(phiRad), + -P0 * cos(thetaRad)); + auto plab = MomentumVector(rootCS, {px, py, pz}); + /* === END: CONSTRUCT PRIMARY PARTICLE === */ + + /* === START: CONSTRUCT GEOMETRY === */ + auto const observationHeight = 0_km + builder.getPlanetRadius(); + auto const injectionHeight = 111.75_km + builder.getPlanetRadius(); + auto const t = -observationHeight * cos(thetaRad) + + sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) + + static_pow<2>(injectionHeight)); + Point const showerCore{rootCS, 0_m, 0_m, observationHeight}; + Point const injectionPos = + showerCore + DirectionVector{rootCS, + {-sin(thetaRad) * cos(phiRad), + -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} * + t; + + // we make the axis much longer than the inj-core distance since the + // profile will go beyond the core, depending on zenith angle + ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env}; + /* === END: CONSTRUCT GEOMETRY === */ + + // create the output manager that we then register outputs with + OutputManager output(app["--filename"]->as()); + + /* === START: SETUP PROCESS LIST === */ + 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{1_GeV, 1_GeV, 1_GeV, 1_GeV, false}; + corsika::proposal::Interaction emCascade(env); + corsika::proposal::ContinuousProcess emContinuous(env); + InteractionCounter emCascadeCounted(emCascade); + + LongitudinalProfile longprof{showerAxis}; + + corsika::urqmd::UrQMD urqmd; + InteractionCounter urqmdCounted{urqmd}; + StackInspector stackInspect(5000, false, E0); + + // assemble all processes into an ordered process list + struct EnergySwitch { + HEPEnergyType cutE_; + EnergySwitch(HEPEnergyType cutE) + : cutE_(cutE) {} + bool operator()(const Particle& p) { return (p.getKineticEnergy() < cutE_); } + }; + auto hadronSequence = make_select(EnergySwitch(80_GeV), urqmdCounted, + make_sequence(sibyllNucCounted, sibyllCounted)); + auto decaySequence = make_sequence(decayPythia, decaySibyll); + + // track writer + TrackWriter trackWriter; + output.add("tracks", trackWriter); // register TrackWriter + + // observation plane + Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); + ObservationPlane observationLevel( + obsPlane, DirectionVector(rootCS, {1., 0., 0.})); + // register the observation plane with the output + output.add("particles", observationLevel); + + // assemble the final process sequence + auto sequence = + make_sequence(stackInspect, hadronSequence, decaySequence, emCascadeCounted, + emContinuous, cut, trackWriter, observationLevel, longprof); + /* === END: SETUP PROCESS LIST === */ + + // create the cascade object using the default stack and tracking implementation + setup::Tracking tracking; + setup::Stack stack; + Cascade EAS(env, tracking, sequence, output, stack); + + // print our primary parameters all in one place + if (app["--pdg"]->count() > 0) { + CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as()); + } else { + CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A); + } + CORSIKA_LOG_INFO("Primary Energy: {}", E0); + CORSIKA_LOG_INFO("Primary Momentum: {}", P0); + CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates()); + CORSIKA_LOG_INFO("Shower Axis Length: {}", (showerCore - injectionPos).getNorm() * 1.2); + + // trigger the output manager to open the library for writing + output.startOfLibrary(); + + // loop over each shower + for (int i_shower = 1; i_shower < nevent + 1; i_shower++) { + CORSIKA_LOG_INFO("Shower {} / {} ", i_shower, nevent); + + // trigger the start of the outputs for this shower + output.startOfShower(); + + // directory for outputs + string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz"; + string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz"; + string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt"; + + // setup particle stack, and add primary particle + stack.clear(); + + // add the desired particle to the stack + if (A > 1) { + stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z)); + } else { + stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns)); + } + + // if we want to fix the first location of the shower + if (app["--force-interaction"]) EAS.forceInteraction(); + + // run the shower + EAS.run(); + + cut.showResults(); + // emContinuous.showResults(); + observationLevel.showResults(); + const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() + + 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(); + + auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() + + urqmdCounted.getHistogram(); + + save_hist(hists.labHist(), labHist_file, true); + save_hist(hists.CMSHist(), cMSHist_file, true); + longprof.save(longprof_file); + + // trigger the output manager to save this shower to disk + output.endOfShower(); + } + + // and finalize the output on disk + output.endOfLibrary(); +} -- GitLab From b6f964afa2368247be80fcc8e41bf9832588a6ff Mon Sep 17 00:00:00 2001 From: ralfulrich Date: Fri, 18 Jun 2021 13:58:22 +0200 Subject: [PATCH 4/8] missing files --- corsika/detail/media/BaseTabular.inl | 201 ++++++++++++++++++ corsika/detail/media/SlidingPlanarTabular.inl | 47 ++++ corsika/media/BaseTabular.hpp | 89 ++++++++ corsika/media/SlidingPlanarTabular.hpp | 62 ++++++ 4 files changed, 399 insertions(+) create mode 100644 corsika/detail/media/BaseTabular.inl create mode 100644 corsika/detail/media/SlidingPlanarTabular.inl create mode 100644 corsika/media/BaseTabular.hpp create mode 100644 corsika/media/SlidingPlanarTabular.hpp diff --git a/corsika/detail/media/BaseTabular.inl b/corsika/detail/media/BaseTabular.inl new file mode 100644 index 000000000..916fb0320 --- /dev/null +++ b/corsika/detail/media/BaseTabular.inl @@ -0,0 +1,201 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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 +#include +#include + +#include + +namespace corsika { + + template + inline BaseTabular::BaseTabular( + Point const& point, LengthType const referenceHeight, + std::function const& rho, unsigned int const nBins, + LengthType const deltaHeight) + : nBins_(nBins) + , deltaHeight_(deltaHeight) + , point_(point) + , referenceHeight_(referenceHeight) { + density_.resize(nBins_); + for (unsigned int bin = 0; bin < nBins; ++bin) { + density_[bin] = rho(deltaHeight_ * bin); + CORSIKA_LOG_DEBUG("new tabulated atm bin={} h={} rho={}", bin, deltaHeight_ * bin, + density_[bin]); + } + } + + template + inline auto const& BaseTabular::getImplementation() const { + return *static_cast(this); + } + + template + inline MassDensityType BaseTabular::getMassDensity( + LengthType const height) const { + double const fbin = (height - referenceHeight_) / deltaHeight_; + int const bin = int(fbin); + if (bin < 0) return MassDensityType::zero(); + if (bin >= int(nBins_ - 1)) { + CORSIKA_LOG_ERROR( + "invalid height {} (corrected {}) in BaseTabular atmosphere. Min 0, max {}. If " + "max is too low: increase!", + height, height - referenceHeight_, nBins_ * deltaHeight_); + throw std::runtime_error("invalid height"); + } + return density_[bin] + (fbin - bin) * (density_[bin + 1] - density_[bin]); + } + + template + inline GrammageType BaseTabular::getIntegratedGrammage( + BaseTrajectory const& traj) const { + + Point pCurr = traj.getPosition(0); + DirectionVector dCurr = traj.getDirection(0); + LengthType height1 = (traj.getPosition(0) - point_).getNorm() - referenceHeight_; + LengthType height2 = (traj.getPosition(1) - point_).getNorm() - referenceHeight_; + + LengthType const fullLength = traj.getLength(1); + + int sign = +1; // normal direction + if (height1 > height2) { + std::swap(height1, height2); + pCurr = traj.getPosition(1); + dCurr = traj.getDirection(1); + sign = -1; // inverted direction + } + + double const fbin1 = height1 / deltaHeight_; + unsigned int const bin1 = int(fbin1); + + double const fbin2 = height2 / deltaHeight_; + unsigned int const bin2 = int(fbin2); + + if (fbin1 == fbin2) { return GrammageType::zero(); } + + if (bin1 >= nBins_ - 1 || bin2 >= nBins_ - 1) { + CORSIKA_LOG_ERROR("invalid height {} {} in BaseTabular atmosphere integration", + height1, height2); + throw std::runtime_error("invalid height"); + } + + // interpolated start/end densities + MassDensityType const rho1 = getMassDensity(height1 + referenceHeight_); + MassDensityType const rho2 = getMassDensity(height2 + referenceHeight_); + + // within first bin + if (bin1 == bin2) { return fullLength * (rho2 + rho1) / 2; } + + // inclination of trajectory (local) + DirectionVector axis((pCurr - point_).normalized()); // to gravity center + double cosTheta = axis.dot(dCurr); + + // distance to next height bin + unsigned int bin = bin1; + LengthType dD = (deltaHeight_ * (bin + 1) - height1) / cosTheta * sign; + LengthType distance = dD; + + GrammageType X = dD * (rho1 + density_[bin + 1]) / 2; + double frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength); + pCurr = traj.getPosition(frac); + dCurr = traj.getDirection(frac); + + for (++bin; bin < bin2; ++bin) { + // inclination of trajectory + axis = (pCurr - point_).normalized(); + cosTheta = axis.dot(dCurr); + // distance to next height bin + dD = deltaHeight_ / cosTheta * sign; + distance += dD; + GrammageType const dX = dD * (density_[bin] + density_[bin + 1]) / 2; + X += dX; + frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength); + pCurr = traj.getPosition(frac); + dCurr = traj.getDirection(frac); + } + + // inclination of trajectory + axis = ((pCurr - point_).normalized()); + cosTheta = axis.dot(dCurr); + // distance to next height bin + dD = (height2 - deltaHeight_ * bin2) / cosTheta * sign; + X += dD * (rho2 + density_[bin2]) / 2; + distance += dD; + return X; + } + + template + inline LengthType BaseTabular::getArclengthFromGrammage( + BaseTrajectory const& traj, GrammageType const grammage) const { + + if (grammage < GrammageType::zero()) { + CORSIKA_LOG_ERROR("cannot integrate negative grammage"); + throw std::runtime_error("negative grammage error"); + } + + LengthType const height = (traj.getPosition(0) - point_).getNorm() - referenceHeight_; + + double const fbin = height / deltaHeight_; + int bin = int(fbin); + + if (bin >= int(nBins_ - 1)) { + CORSIKA_LOG_ERROR("invalid height {} in BaseTabular atmosphere integration", + height); + throw std::runtime_error("invalid height"); + } + + // interpolated start/end densities + MassDensityType const rho = getMassDensity(height + referenceHeight_); + + // inclination of trajectory + Point pCurr = traj.getPosition(0); + DirectionVector dCurr = traj.getDirection(0); + DirectionVector axis((pCurr - point_).normalized()); + double cosTheta = axis.dot(dCurr); + int sign = +1; // height increasing along traj + if (cosTheta < 0) { + cosTheta = -cosTheta; // absolute value only + sign = -1; // height decreasing along traj + } + + // height -> distance + LengthType const deltaDistance = deltaHeight_ / cosTheta; + + // start with 0 g/cm2 + GrammageType X(GrammageType::zero()); + LengthType distance(LengthType::zero()); + + // within first bin + distance = + (sign > 0 ? deltaDistance * (bin + 1 - fbin) : deltaDistance * (fbin - bin)); + GrammageType binGrammage = (sign > 0 ? distance * (rho + density_[bin + 1]) / 2 + : distance * (rho + density_[bin]) / 2); + if (X + binGrammage > grammage) { + double const binFraction = (grammage - X) / binGrammage; + return distance * binFraction; + } + X += binGrammage; + + // the following bins (along trajectory) + for (bin += sign; bin < int(nBins_) && bin >= 0; bin += sign) { + + binGrammage = deltaDistance * (density_[bin] + density_[bin + 1]) / 2; + if (X + binGrammage > grammage) { + double const binFraction = (grammage - X) / binGrammage; + return distance + deltaDistance * binFraction; + } + X += binGrammage; + distance += deltaDistance; + } + return std::numeric_limits::infinity() * meter; + } + +} // namespace corsika diff --git a/corsika/detail/media/SlidingPlanarTabular.inl b/corsika/detail/media/SlidingPlanarTabular.inl new file mode 100644 index 000000000..c8e15792c --- /dev/null +++ b/corsika/detail/media/SlidingPlanarTabular.inl @@ -0,0 +1,47 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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 + +namespace corsika { + + template + inline SlidingPlanarTabular::SlidingPlanarTabular( + Point const& p0, std::function const& rho, + unsigned int const nBins, LengthType const deltaHeight, + NuclearComposition const& nuclComp, LengthType referenceHeight) + : BaseTabular>(p0, referenceHeight, rho, nBins, deltaHeight) + , nuclComp_(nuclComp) {} + + template + inline MassDensityType SlidingPlanarTabular::getMassDensity( + Point const& point) const { + auto const heightFull = + (point - BaseTabular>::getAnchorPoint()).getNorm(); + return BaseTabular>::getMassDensity(heightFull); + } + + template + inline NuclearComposition const& SlidingPlanarTabular::getNuclearComposition() + const { + return nuclComp_; + } + + template + inline GrammageType SlidingPlanarTabular::getIntegratedGrammage( + BaseTrajectory const& traj) const { + return BaseTabular>::getIntegratedGrammage(traj); + } + + template + inline LengthType SlidingPlanarTabular::getArclengthFromGrammage( + BaseTrajectory const& traj, GrammageType const grammage) const { + return BaseTabular>::getArclengthFromGrammage(traj, grammage); + } + +} // namespace corsika diff --git a/corsika/media/BaseTabular.hpp b/corsika/media/BaseTabular.hpp new file mode 100644 index 000000000..a552b6648 --- /dev/null +++ b/corsika/media/BaseTabular.hpp @@ -0,0 +1,89 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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 +#include +#include +#include +#include + +#include +#include + +namespace corsika { + + /** + * This class provides the grammage/length conversion functionality for + * (locally) flat exponential atmospheres. + */ + template + class BaseTabular { + + public: + BaseTabular(Point const& point, LengthType const referenceHeight, + std::function const& rho, + unsigned int const nBins, LengthType const deltaHeight); + + Point const& getAnchorPoint() const { return point_; } + + protected: + auto const& getImplementation() const; + + MassDensityType getMassDensity(LengthType const height) const; + + // clang-format off + /** + * For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized) + * direction \f$ \vec{u} \f$ is given by + * \f[ + * X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right) + * \f], where \f$ \varrho_0 \f$ is the density at the starting point. + * + * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density: + * \f[ + * X = \varrho_0 l; + * \f] + */ + // clang-format on + GrammageType getIntegratedGrammage(BaseTrajectory const& line) const; + + // clang-format off + /** + * For a (normalized) axis \f$ \vec{a} \f$, the length of a non-orthogonal line with (normalized) + * direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by + * \f[ + * l = \begin{cases} + * \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} Y := 0 > 1 + + * \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda} + * \infty & \text{else,} + * \end{cases} + * \f] where \f$ \varrho_0 \f$ is the density at the starting point. + * + * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density: + * \f[ + * l = \frac{X}{\varrho_0} + * \f] + */ + // clang-format on + LengthType getArclengthFromGrammage(BaseTrajectory const& line, + GrammageType const grammage) const; + + private: + unsigned int const nBins_; + LengthType deltaHeight_; + std::vector density_; + Point const point_; + LengthType const referenceHeight_; + + }; // class BaseTabular + +} // namespace corsika + +#include diff --git a/corsika/media/SlidingPlanarTabular.hpp b/corsika/media/SlidingPlanarTabular.hpp new file mode 100644 index 000000000..4550792aa --- /dev/null +++ b/corsika/media/SlidingPlanarTabular.hpp @@ -0,0 +1,62 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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 +#include +#include +#include +#include +#include +#include +#include + +namespace corsika { + + // clang-format off + /** + * The SlidingPlanarTabular models mass density as + * \f[ + * \varrho(r) = \varrho_0 \exp\left( \frac{|p_0 - r|}{\lambda} \right). + * \f] + * For grammage/length conversion, the density distribution is approximated as + * locally flat at the starting point \f$ r_0 \f$ of the trajectory with the axis pointing + * from \f$ p_0 \f$ to \f$ r_0 \f$. + */ + // clang-format on + + template + class SlidingPlanarTabular : public BaseTabular>, + public TDerived { + + using Base = BaseTabular>; + + public: + SlidingPlanarTabular(Point const& p0, + std::function const& rho, + unsigned int const nBins, LengthType const deltaHeight, + NuclearComposition const& nuclComp, + LengthType referenceHeight = LengthType::zero()); + + MassDensityType getMassDensity(Point const& point) const override; + + NuclearComposition const& getNuclearComposition() const override; + + GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override; + + LengthType getArclengthFromGrammage(BaseTrajectory const& line, + GrammageType grammage) const override; + + private: + NuclearComposition const nuclComp_; + }; + +} // namespace corsika + +#include -- GitLab From 5db5e0692a59f1c23c8d3af937f0ff53b5f17622 Mon Sep 17 00:00:00 2001 From: ralfulrich Date: Fri, 18 Jun 2021 15:06:37 +0200 Subject: [PATCH 5/8] fix CI --- examples/corsika.cpp | 39 ++++++++++++++------------------- examples/mars.cpp | 39 ++++++++++++++------------------- tests/media/testEnvironment.cpp | 2 +- 3 files changed, 35 insertions(+), 45 deletions(-) diff --git a/examples/corsika.cpp b/examples/corsika.cpp index ee424c1de..078b9d2be 100644 --- a/examples/corsika.cpp +++ b/examples/corsika.cpp @@ -97,10 +97,6 @@ void registerRandomStreams(int seed) { template using MyExtraEnv = MediumPropertyModel>; -// argv : 1.number of nucleons, 2.number of protons, -// 3.total energy in GeV, 4.number of showers, -// 5.seed (0 by default to generate random values for all) - int main(int argc, char** argv) { // the main command line description @@ -160,21 +156,19 @@ int main(int argc, char** argv) { // parse the command line options into the variables CLI11_PARSE(app, argc, argv); - if (app.count("--verbosity")) { - string const loglevel = app["verbosity"]->as(); - if (loglevel == "warn") { - logging::set_level(logging::level::warn); - } else if (loglevel == "info") { - logging::set_level(logging::level::info); - } else if (loglevel == "debug") { - logging::set_level(logging::level::debug); - } else if (loglevel == "trace") { + string const loglevel = app["verbosity"]->as(); + if (loglevel == "warn") { + logging::set_level(logging::level::warn); + } else if (loglevel == "info") { + logging::set_level(logging::level::info); + } else if (loglevel == "debug") { + logging::set_level(logging::level::debug); + } else if (loglevel == "trace") { #ifndef DEBUG - CORSIKA_LOG_ERROR("trace log level requires a Debug build."); - return 1; + CORSIKA_LOG_ERROR("trace log level requires a Debug build."); + return 1; #endif - logging::set_level(logging::level::trace); - } + logging::set_level(logging::level::trace); } // check that we got either PDG or A/Z @@ -216,7 +210,7 @@ int main(int argc, char** argv) { /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */ ofstream atmout("earth.dat"); - for (LengthType h = 0_m; h < 110_km; h += 100_m) { + for (LengthType h = 0_m; h < 110_km; h += 10_m) { Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h}; auto rho = env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity( @@ -342,7 +336,7 @@ int main(int argc, char** argv) { // observation plane Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); - ObservationPlane observationLevel( + ObservationPlane observationLevel( obsPlane, DirectionVector(rootCS, {1., 0., 0.})); // register the observation plane with the output output.add("particles", observationLevel); @@ -381,9 +375,10 @@ int main(int argc, char** argv) { output.startOfShower(); // directory for outputs - string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz"; - string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz"; - string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt"; + string const outdir(app["--filename"]->as()); + string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz"; + string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz"; + string const longprof_file = outdir + "/longprof_" + to_string(i_shower) + ".txt"; // setup particle stack, and add primary particle stack.clear(); diff --git a/examples/mars.cpp b/examples/mars.cpp index ac25be75c..0a9d7610b 100644 --- a/examples/mars.cpp +++ b/examples/mars.cpp @@ -125,10 +125,6 @@ void registerRandomStreams(int seed) { template using MyExtraEnv = MediumPropertyModel>; -// argv : 1.number of nucleons, 2.number of protons, -// 3.total energy in GeV, 4.number of showers, -// 5.seed (0 by default to generate random values for all) - int main(int argc, char** argv) { // the main command line description @@ -188,21 +184,19 @@ int main(int argc, char** argv) { // parse the command line options into the variables CLI11_PARSE(app, argc, argv); - if (app.count("--verbosity")) { - string const loglevel = app["verbosity"]->as(); - if (loglevel == "warn") { - logging::set_level(logging::level::warn); - } else if (loglevel == "info") { - logging::set_level(logging::level::info); - } else if (loglevel == "debug") { - logging::set_level(logging::level::debug); - } else if (loglevel == "trace") { + string const loglevel = app["verbosity"]->as(); + if (loglevel == "warn") { + logging::set_level(logging::level::warn); + } else if (loglevel == "info") { + logging::set_level(logging::level::info); + } else if (loglevel == "debug") { + logging::set_level(logging::level::debug); + } else if (loglevel == "trace") { #ifndef DEBUG - CORSIKA_LOG_ERROR("trace log level requires a Debug build."); - return 1; + CORSIKA_LOG_ERROR("trace log level requires a Debug build."); + return 1; #endif - logging::set_level(logging::level::trace); - } + logging::set_level(logging::level::trace); } // check that we got either PDG or A/Z @@ -251,7 +245,7 @@ int main(int argc, char** argv) { /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */ ofstream atmout("mars.dat"); - for (LengthType h = 0_m; h < 110_km; h += 100_m) { + for (LengthType h = 0_m; h < 110_km; h += 10_m) { Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h}; auto rho = env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity( @@ -375,7 +369,7 @@ int main(int argc, char** argv) { // observation plane Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); - ObservationPlane observationLevel( + ObservationPlane observationLevel( obsPlane, DirectionVector(rootCS, {1., 0., 0.})); // register the observation plane with the output output.add("particles", observationLevel); @@ -413,9 +407,10 @@ int main(int argc, char** argv) { output.startOfShower(); // directory for outputs - string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz"; - string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz"; - string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt"; + string const outdir(app["--filename"]->as()); + string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz"; + string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz"; + string const longprof_file = outdir + "/longprof" + to_string(i_shower) + ".txt"; // setup particle stack, and add primary particle stack.clear(); diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp index fffdff7fe..224bed806 100644 --- a/tests/media/testEnvironment.cpp +++ b/tests/media/testEnvironment.cpp @@ -489,7 +489,7 @@ TEST_CASE("InhomogeneousMedium") { auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); }; auto const exactLength = [](auto X) { return 1_m * log(1 + X / (rho0 * 1_m)); }; - LengthType constexpr length = tEnd * speed; + LengthType const length = tEnd * speed; NuclearComposition const composition{{Code::Proton}, {1.f}}; InhomogeneousMedium const inhMedium(composition, rho); -- GitLab From a995ef0d443f18462adaebeb8d3eaca1c3ccbae1 Mon Sep 17 00:00:00 2001 From: ralfulrich Date: Fri, 18 Jun 2021 15:23:02 +0200 Subject: [PATCH 6/8] cli default value handling --- examples/corsika.cpp | 3 ++- examples/mars.cpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/corsika.cpp b/examples/corsika.cpp index 078b9d2be..09e024287 100644 --- a/examples/corsika.cpp +++ b/examples/corsika.cpp @@ -156,7 +156,8 @@ int main(int argc, char** argv) { // parse the command line options into the variables CLI11_PARSE(app, argc, argv); - string const loglevel = app["verbosity"]->as(); + string const loglevel = + (app.count("--verbosity") ? app["verbosity"]->as() : "info"); if (loglevel == "warn") { logging::set_level(logging::level::warn); } else if (loglevel == "info") { diff --git a/examples/mars.cpp b/examples/mars.cpp index 0a9d7610b..73ca6fbf9 100644 --- a/examples/mars.cpp +++ b/examples/mars.cpp @@ -184,7 +184,8 @@ int main(int argc, char** argv) { // parse the command line options into the variables CLI11_PARSE(app, argc, argv); - string const loglevel = app["verbosity"]->as(); + string const loglevel = + (app.count("--verbosity") ? app["verbosity"]->as() : "info"); if (loglevel == "warn") { logging::set_level(logging::level::warn); } else if (loglevel == "info") { -- GitLab From d7917b08e59797e639ebad1c2c59b0a8bb746052 Mon Sep 17 00:00:00 2001 From: ralfulrich Date: Sat, 19 Jun 2021 10:56:16 +0200 Subject: [PATCH 7/8] fix verbosity flag --- examples/corsika.cpp | 6 +++--- examples/mars.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/corsika.cpp b/examples/corsika.cpp index 09e024287..275ca88e2 100644 --- a/examples/corsika.cpp +++ b/examples/corsika.cpp @@ -148,8 +148,8 @@ int main(int argc, char** argv) { ->group("Misc."); app.add_flag("--force-interaction", "Force the location of the first interaction.") ->group("Misc."); - app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.") - ->default_val("info") + app.add_option("-v,--verbosity", "Verbosity level") + ->default_str("info") ->check(CLI::IsMember({"warn", "info", "debug", "trace"})) ->group("Misc."); @@ -157,7 +157,7 @@ int main(int argc, char** argv) { CLI11_PARSE(app, argc, argv); string const loglevel = - (app.count("--verbosity") ? app["verbosity"]->as() : "info"); + (app.count("--verbosity") ? app["--verbosity"]->as() : "info"); if (loglevel == "warn") { logging::set_level(logging::level::warn); } else if (loglevel == "info") { diff --git a/examples/mars.cpp b/examples/mars.cpp index 73ca6fbf9..92ed62ad9 100644 --- a/examples/mars.cpp +++ b/examples/mars.cpp @@ -347,7 +347,7 @@ int main(int argc, char** argv) { corsika::proposal::ContinuousProcess emContinuous(env); InteractionCounter emCascadeCounted(emCascade); - LongitudinalProfile longprof{showerAxis}; + LongitudinalProfile longprof{showerAxis, 1_g / square(1_cm)}; corsika::urqmd::UrQMD urqmd; InteractionCounter urqmdCounted{urqmd}; -- GitLab From 9cf2ab6b4e4432df7b60da9a2be7699d9cfeee69 Mon Sep 17 00:00:00 2001 From: ralfulrich Date: Mon, 21 Jun 2021 23:26:47 +0200 Subject: [PATCH 8/8] coverage --- corsika/detail/media/BaseTabular.inl | 1 - corsika/detail/modules/urqmd/UrQMD.inl | 9 ++--- corsika/media/BaseTabular.hpp | 34 +------------------ corsika/media/SlidingPlanarTabular.hpp | 7 ++-- corsika/stack/DummyStack.hpp | 2 ++ tests/media/testEnvironment.cpp | 27 +++++++++++++++ tests/modules/testEpos.cpp | 1 + .../stack/testGeometryNodeStackExtension.cpp | 26 ++++++++++++++ 8 files changed, 65 insertions(+), 42 deletions(-) diff --git a/corsika/detail/media/BaseTabular.inl b/corsika/detail/media/BaseTabular.inl index 916fb0320..3a02f63ad 100644 --- a/corsika/detail/media/BaseTabular.inl +++ b/corsika/detail/media/BaseTabular.inl @@ -140,7 +140,6 @@ namespace corsika { CORSIKA_LOG_ERROR("cannot integrate negative grammage"); throw std::runtime_error("negative grammage error"); } - LengthType const height = (traj.getPosition(0) - point_).getNorm() - referenceHeight_; double const fbin = height / deltaHeight_; diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl index 4b40f19f4..9c9bcfe1a 100644 --- a/corsika/detail/modules/urqmd/UrQMD.inl +++ b/corsika/detail/modules/urqmd/UrQMD.inl @@ -88,7 +88,8 @@ namespace corsika::urqmd { break; default: { // LCOV_EXCL_START since this can never happen due to canInteract CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode); - return CrossSectionType::zero(); // LCOV_EXCL_STOP + return CrossSectionType::zero(); + // LCOV_EXCL_STOP } } @@ -395,9 +396,9 @@ namespace corsika::urqmd { boost::filesystem::ifstream file(filename, std::ios::in); if (!file.is_open()) { - throw std::runtime_error( - filename.native() + - " could not be opened."); // LCOV_EXCL_LINE since this is pointless to test + // LCOV_EXCL_START since this is pointless to test + throw std::runtime_error(filename.native() + " could not be opened."); + // LCOV_EXCL_STOP } std::string line; diff --git a/corsika/media/BaseTabular.hpp b/corsika/media/BaseTabular.hpp index a552b6648..cecf7cc8a 100644 --- a/corsika/media/BaseTabular.hpp +++ b/corsika/media/BaseTabular.hpp @@ -21,7 +21,7 @@ namespace corsika { /** * This class provides the grammage/length conversion functionality for - * (locally) flat exponential atmospheres. + * (locally) flat tabulated atmospheres. */ template class BaseTabular { @@ -38,40 +38,8 @@ namespace corsika { MassDensityType getMassDensity(LengthType const height) const; - // clang-format off - /** - * For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized) - * direction \f$ \vec{u} \f$ is given by - * \f[ - * X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right) - * \f], where \f$ \varrho_0 \f$ is the density at the starting point. - * - * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density: - * \f[ - * X = \varrho_0 l; - * \f] - */ - // clang-format on GrammageType getIntegratedGrammage(BaseTrajectory const& line) const; - // clang-format off - /** - * For a (normalized) axis \f$ \vec{a} \f$, the length of a non-orthogonal line with (normalized) - * direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by - * \f[ - * l = \begin{cases} - * \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} Y := 0 > 1 + - * \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda} - * \infty & \text{else,} - * \end{cases} - * \f] where \f$ \varrho_0 \f$ is the density at the starting point. - * - * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density: - * \f[ - * l = \frac{X}{\varrho_0} - * \f] - */ - // clang-format on LengthType getArclengthFromGrammage(BaseTrajectory const& line, GrammageType const grammage) const; diff --git a/corsika/media/SlidingPlanarTabular.hpp b/corsika/media/SlidingPlanarTabular.hpp index 4550792aa..afd0d88e2 100644 --- a/corsika/media/SlidingPlanarTabular.hpp +++ b/corsika/media/SlidingPlanarTabular.hpp @@ -12,7 +12,6 @@ #include #include #include -#include #include #include #include @@ -23,11 +22,11 @@ namespace corsika { /** * The SlidingPlanarTabular models mass density as * \f[ - * \varrho(r) = \varrho_0 \exp\left( \frac{|p_0 - r|}{\lambda} \right). + * \varrho(r) = \varrho_0 \rho\left( |p_0 - r| \right). * \f] * For grammage/length conversion, the density distribution is approximated as - * locally flat at the starting point \f$ r_0 \f$ of the trajectory with the axis pointing - * from \f$ p_0 \f$ to \f$ r_0 \f$. + * locally flat at the starting point \f$ r_0 \f$ of the trajectory with the + * axis pointing rom \f$ p_0 \f$ to \f$ r_0 \f$ defining the local height. */ // clang-format on diff --git a/corsika/stack/DummyStack.hpp b/corsika/stack/DummyStack.hpp index 9bbd51403..bcba30f1d 100644 --- a/corsika/stack/DummyStack.hpp +++ b/corsika/stack/DummyStack.hpp @@ -71,6 +71,8 @@ namespace corsika::dummy_stack { */ void copy(const int /*i1*/, const int /*i2*/) {} + void swap(const int, const int) {} + void incrementSize() { entries_++; } void decrementSize() { entries_--; } diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp index 224bed806..372cfdb9a 100644 --- a/tests/media/testEnvironment.cpp +++ b/tests/media/testEnvironment.cpp @@ -83,6 +83,9 @@ TEST_CASE("HomogeneousMedium") { std::vector{1.f}); HomogeneousMedium const medium(19.2_g / cube(1_cm), protonComposition); + CHECK(protonComposition.getFractions() == std::vector{1.}); + CHECK(protonComposition.getComponents() == std::vector{Code::Proton}); + CHECK_THROWS(NuclearComposition({Code::Proton}, {1.1})); CHECK_THROWS(NuclearComposition({Code::Proton}, {0.99})); } @@ -103,6 +106,10 @@ TEST_CASE("FlatExponential") { LengthType const length = 2_m; TimeType const tEnd = length / speed; + CHECK(medium.getNuclearComposition().getFractions() == std::vector{1.}); + CHECK(medium.getNuclearComposition().getComponents() == + std::vector{Code::Proton}); + SECTION("horizontal") { Line const line(gOrigin, Vector( gCS, {speed, 0_m / second, 0_m / second})); @@ -245,6 +252,26 @@ TEST_CASE("SlidingPlanarTabular") { SlidingPlanarTabular const medium(gOrigin, rhoFunc, 1000, 10_m, protonComposition); + SECTION("not possible") { + CHECK_THROWS(medium.getMassDensity({gCS, {0_m, 1e10_m, 0_m}})); + + SpeedType const speed = 5_m / second; + TimeType const tEnd = 1e10_s; + Line const line( + {gCS, {0_m, 0_m, 1_m}}, + Vector(gCS, {0_m / second, 0_m / second, speed})); + setup::Trajectory const trajectory = + setup::testing::make_track(line, tEnd); + CHECK_THROWS(medium.getIntegratedGrammage(trajectory)); + + Line const line2( + {gCS, {0_m, 0_m, 1e9_m}}, + Vector(gCS, {0_m / second, 0_m / second, speed})); + setup::Trajectory const trajectory2 = + setup::testing::make_track(line2, tEnd); + CHECK_THROWS(medium.getArclengthFromGrammage(trajectory2, 1e3_g / square(1_cm))); + } + SECTION("density") { CHECK(medium.getMassDensity({gCS, {0_m, 0_m, 3_m}}) / medium.getMassDensity({gCS, {0_m, 3_m, 0_m}}) == diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp index a4c4ef593..aa2a006f2 100644 --- a/tests/modules/testEpos.cpp +++ b/tests/modules/testEpos.cpp @@ -38,6 +38,7 @@ TEST_CASE("epos", "modules") { corsika::epos::convertFromEpos(corsika::epos::EposCode::Electron)); CHECK(Code::Proton == corsika::epos::convertFromEpos(corsika::epos::EposCode::Proton)); + CHECK_THROWS(corsika::epos::convertFromEpos(corsika::epos::EposCode::Unknown)); } SECTION("corsika -> epos") { diff --git a/tests/stack/testGeometryNodeStackExtension.cpp b/tests/stack/testGeometryNodeStackExtension.cpp index 9cc29edcc..2bc2e1160 100644 --- a/tests/stack/testGeometryNodeStackExtension.cpp +++ b/tests/stack/testGeometryNodeStackExtension.cpp @@ -86,4 +86,30 @@ TEST_CASE("GeometryNodeStackExtension", "stack") { CHECK(v == 99 * data); CHECK(s.getEntries() == 0); } + + SECTION("copy and swap") { + + const int data1 = 16; + const int data2 = 17; + + TestStack s; + // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2! + for (int i = 0; i < 4; ++i) { + auto p = s.addParticle(std::tuple{noData}); + p.setNode(i % 2 == 0 ? &data1 : &data2); + } + + CHECK(*((s.first() + 0)).getNode() == 16); + CHECK(*((s.first() + 1)).getNode() == 17); + CHECK(*((s.first() + 2)).getNode() == 16); + CHECK(*((s.first() + 3)).getNode() == 17); + + s.copy(s.first() + 0, s.first() + 1); + CHECK(*((s.first() + 0)).getNode() == 16); + CHECK(*((s.first() + 1)).getNode() == 16); + + s.swap(s.first() + 3, s.first() + 1); + CHECK(*((s.first() + 1)).getNode() == 17); + CHECK(*((s.first() + 3)).getNode() == 16); + } } -- GitLab