From 2ffa9f2d0f368adff38862003b83d49e164cabee Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Thu, 9 Dec 2021 18:45:17 +0100
Subject: [PATCH] fixed test and various things around it
---
corsika/detail/framework/core/Cascade.inl | 2 +-
.../tracking/TrackingLeapFrogStraight.inl | 20 +-
.../detail/stack/NuclearStackExtension.inl | 384 ------------------
corsika/stack/DummyStack.hpp | 15 +-
corsika/stack/NuclearStackExtension.hpp | 307 --------------
corsika/stack/VectorStack.hpp | 22 +-
.../stack/history/HistoryStackExtension.hpp | 6 +-
tests/modules/testObservationPlane.cpp | 2 +-
tests/modules/testTracking.cpp | 55 ++-
tests/stack/CMakeLists.txt | 1 +
tests/stack/testHistoryView.cpp | 79 ++--
11 files changed, 80 insertions(+), 813 deletions(-)
delete mode 100644 corsika/detail/stack/NuclearStackExtension.inl
delete mode 100644 corsika/stack/NuclearStackExtension.hpp
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index f11427148..82b6887a3 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -211,7 +211,7 @@ namespace corsika {
}
particle.setTime(particle.getTime() + step.getDuration());
particle.setPosition(step.getPosition(1));
- particle.setMomentum(step.getDirection(1) * particle.getMomentum().getNorm());
+ particle.setDirection(step.getDirection(1));
if (isContinuous) {
return; // there is nothing further, step is finished
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
index 675636a17..a39e8f59c 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
@@ -100,9 +100,7 @@ namespace corsika {
CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
// calculate first halve step for "steplimit"
- auto const initialMomentum = particle.getMomentum();
- auto const absMomentum = initialMomentum.getNorm();
- DirectionVector const direction = initialVelocity.normalized();
+ DirectionVector const initialDirection = particle.getDirection();
// avoid any intersections within first halve steplength
// aim for intersection in second halve step
@@ -112,26 +110,28 @@ namespace corsika {
CORSIKA_LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
firstHalveSteplength, steplimit, initialTrackLength);
// perform the first halve-step
- Point const position_mid = initialPosition + direction * firstHalveSteplength;
+ Point const position_mid =
+ initialPosition + initialDirection * firstHalveSteplength;
auto const k = charge / (constants::c * convert_HEP_to_SI<MassType::dimension_type>(
particle.getMomentum().getNorm()));
DirectionVector const new_direction =
- direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
+ initialDirection +
+ initialDirection.cross(magneticfield) * firstHalveSteplength * 2 * k;
auto const new_direction_norm = new_direction.getNorm(); // by design this is >1
CORSIKA_LOG_DEBUG(
"position_mid={}, new_direction={}, (new_direction_norm)={}, deflection={}",
position_mid.getCoordinates(), new_direction.getComponents(),
new_direction_norm,
- acos(std::min(1.0, direction.dot(new_direction) / new_direction_norm)) * 180 /
- M_PI);
+ acos(std::min(1.0, initialDirection.dot(new_direction) / new_direction_norm)) *
+ 180 / M_PI);
// check, where the second halve-step direction has geometric intersections
particle.setPosition(position_mid);
- particle.setMomentum(new_direction * absMomentum);
+ particle.setDirection(new_direction);
auto const [finalTrack, finalTrackNextVolume] =
tracking_line::Tracking::getTrack(particle);
- particle.setPosition(initialPosition); // this is not nice...
- particle.setMomentum(initialMomentum); // this is not nice...
+ particle.setPosition(initialPosition); // this is not nice...
+ particle.setDirection(initialDirection); // this is not nice...
LengthType const finalTrackLength = finalTrack.getLength(1);
LengthType const secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
diff --git a/corsika/detail/stack/NuclearStackExtension.inl b/corsika/detail/stack/NuclearStackExtension.inl
deleted file mode 100644
index cdd377d12..000000000
--- a/corsika/detail/stack/NuclearStackExtension.inl
+++ /dev/null
@@ -1,384 +0,0 @@
-/*
- * (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.
- */
-
-#pragma once
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/stack/Stack.hpp>
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/stack/VectorStack.hpp>
-
-#include <algorithm>
-#include <tuple>
-#include <vector>
-
-namespace corsika::nuclear_stack {
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(typename super_type::particle_data_type const& v) {
-
- if (std::get<0>(v) == Code::Nucleus) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- super_type::setParticleData(v);
- setNucleusRef(-1); // this is not a nucleus
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(nuclear_particle_data_type const& v) {
- unsigned short const A = std::get<5>(v);
- unsigned short const Z = std::get<6>(v);
- Code const PID = std::get<0>(v);
- if (PID != Code::Nucleus || A == 0 || Z == 0) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
- setNucleusRef(
- super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
- setNuclearA(A);
- setNuclearZ(Z);
- super_type::setParticleData(typename super_type::particle_data_type{
- PID, std::get<1>(v), std::get<2>(v), std::get<3>(v), std::get<4>(v)});
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(super_type& p, typename super_type::particle_data_type const& v) {
- if (std::get<0>(v) == Code::Nucleus) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- super_type::setParticleData(p, v);
- setNucleusRef(-1); // this is not a nucleus
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(super_type& p, nuclear_particle_data_type const& v) {
-
- unsigned short const A = std::get<5>(v);
- unsigned short const Z = std::get<6>(v);
-
- if (std::get<0>(v) != Code::Nucleus || A == 0 || Z == 0) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- setNucleusRef(
- super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
- setNuclearA(A);
- setNuclearZ(Z);
- super_type::setParticleData(p, typename super_type::particle_data_type{
- std::get<0>(v), std::get<1>(v), std::get<2>(v),
- std::get<3>(v), std::get<4>(v)});
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(typename super_type::particle_data_momentum_type const& v) {
-
- if (std::get<0>(v) == Code::Nucleus) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- super_type::setParticleData(v);
- setNucleusRef(-1); // this is not a nucleus
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(nuclear_particle_data_momentum_type const& v) {
- unsigned short const A = std::get<4>(v);
- unsigned short const Z = std::get<5>(v);
- Code const PID = std::get<0>(v);
- if (PID != Code::Nucleus || A == 0 || Z == 0) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
- setNucleusRef(
- super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
- setNuclearA(A);
- setNuclearZ(Z);
- HEPMassType m = 0_GeV;
- if (PID == Code::Nucleus) {
- m = get_nucleus_mass(get_nucleus_code(A, Z));
- } else {
- m = get_mass(PID);
- }
- MomentumVector const& p = std::get<1>(v);
- auto const P2 = p.getSquaredNorm();
- HEPEnergyType Ekin = sqrt(P2 + square(m)) - m;
- super_type::setParticleData(typename super_type::particle_data_type{
- PID, Ekin, p / sqrt(P2), std::get<2>(v), std::get<3>(v)});
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(super_type& p,
- typename super_type::particle_data_momentum_type const& v) {
- if (std::get<0>(v) == Code::Nucleus) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- super_type::setParticleData(p, v);
- setNucleusRef(-1); // this is not a nucleus
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::
- setParticleData(super_type& parent, nuclear_particle_data_momentum_type const& v) {
-
- unsigned short const A = std::get<4>(v);
- unsigned short const Z = std::get<5>(v);
- Code const PID = std::get<0>(v);
- if (PID != Code::Nucleus || A == 0 || Z == 0) {
- std::ostringstream err;
- err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
- throw std::runtime_error(err.str());
- }
-
- setNucleusRef(
- super_type::getStackData().getNucleusNextRef()); // store this nucleus data ref
- setNuclearA(A);
- setNuclearZ(Z);
- HEPMassType m = 0_GeV;
- if (PID == Code::Nucleus) {
- m = get_nucleus_mass(get_nucleus_code(A, Z));
- } else {
- m = get_mass(PID);
- }
- MomentumVector const& p = std::get<1>(v);
- auto const P2 = p.getSquaredNorm();
- HEPEnergyType Ekin = sqrt(P2 + square(m)) - m;
- super_type::setParticleData(
- parent, typename super_type::particle_data_type{PID, Ekin, p / sqrt(P2),
- std::get<2>(v), std::get<3>(v)});
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline std::string NuclearParticleInterface<InnerParticleInterface,
- StackIteratorInterface>::asString() const {
- return fmt::format(
- "{}, nuc({})", super_type::asString(),
- (isNucleus() ? fmt::format("A={}, Z={}", getNuclearA(), getNuclearZ()) : "n/a"));
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline PDGCode NuclearParticleInterface<InnerParticleInterface,
- StackIteratorInterface>::getPDG() const {
- return (isNucleus() ? PDGCode(1000000000 + getNuclearZ() * 10000 + getNuclearA() * 10)
- : super_type::getPDG());
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void
- NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::setMomentum(
- MomentumVector const& v) {
- HEPMomentumType const P = v.getNorm();
- if (P == 0_eV) {
- super_type::getStackData().setKineticEnergy(super_type::getIndex(), 0_eV);
- super_type::getStackData().setDirection(
- super_type::getIndex(), DirectionVector(v.getCoordinateSystem(), {0, 0, 0}));
- } else {
- super_type::getStackData().setKineticEnergy(
- super_type::getIndex(),
- sqrt(square(this->getMass()) + square(P)) - this->getMass());
- super_type::getStackData().setDirection(super_type::getIndex(), v / P);
- }
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline MomentumVector NuclearParticleInterface<
- InnerParticleInterface, StackIteratorInterface>::getMomentum() const {
- auto const P = sqrt(square(this->getKineticEnergy() + this->getMass()) -
- square(this->getMass()));
- return super_type::getStackData().getDirection(super_type::getIndex()) * P;
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline VelocityVector NuclearParticleInterface<
- InnerParticleInterface, StackIteratorInterface>::getVelocity() const {
- return this->getMomentum() / this->getEnergy() * constants::c;
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline HEPMassType NuclearParticleInterface<InnerParticleInterface,
- StackIteratorInterface>::getMass() const {
- if (super_type::getPID() == Code::Nucleus)
- return get_nucleus_mass(get_nucleus_code(getNuclearA(), getNuclearZ()));
- return super_type::getMass();
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline ElectricChargeType NuclearParticleInterface<
- InnerParticleInterface, StackIteratorInterface>::getCharge() const {
- if (super_type::getPID() == Code::Nucleus) return getNuclearZ() * constants::e;
- return super_type::getCharge();
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline HEPEnergyType NuclearParticleInterface<
- InnerParticleInterface, StackIteratorInterface>::getEnergy() const {
- return this->getKineticEnergy() + this->getMass();
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline void
- NuclearParticleInterface<InnerParticleInterface, StackIteratorInterface>::setEnergy(
- HEPEnergyType const& e) {
- super_type::getStackData().setKineticEnergy(super_type::getIndex(),
- e - this->getMass());
- }
-
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- inline int16_t NuclearParticleInterface<
- InnerParticleInterface, StackIteratorInterface>::getChargeNumber() const {
- if (super_type::getPID() == Code::Nucleus) return getNuclearZ();
- return super_type::getChargeNumber();
- }
-
- template <typename InnerStackImpl>
- inline int NuclearStackExtensionImpl<InnerStackImpl>::getNucleusNextRef() {
- nuclearA_.push_back(0);
- nuclearZ_.push_back(0);
- return nuclearA_.size() - 1;
- }
-
- template <typename InnerStackImpl>
- inline int NuclearStackExtensionImpl<InnerStackImpl>::getNucleusRef(
- const unsigned int i) const {
- if (nucleusRef_[i] >= 0) return nucleusRef_[i];
- std::ostringstream err;
- err << "NuclearStackExtension: no nucleus at ref=" << i;
- throw std::runtime_error(err.str());
- }
-
- template <typename InnerStackImpl>
- inline void NuclearStackExtensionImpl<InnerStackImpl>::copy(const unsigned int i1,
- const unsigned int i2) {
- // index range check
- if (i1 >= getSize() || i2 >= getSize()) {
- std::ostringstream err;
- err << "NuclearStackExtension: trying to access data beyond size of stack!";
- throw std::runtime_error(err.str());
- }
- // copy internal particle data p[i2] = p[i1]
- super_type::copy(i1, i2);
- // check if any of p[i1] or p[i2] was a Code::Nucleus
- const int ref1 = nucleusRef_[i1];
- const int ref2 = nucleusRef_[i2];
- if (ref2 < 0) {
- if (ref1 >= 0) {
- // i1 is nucleus, i2 is not
- nucleusRef_[i2] = getNucleusNextRef();
- nuclearA_[nucleusRef_[i2]] = nuclearA_[ref1];
- nuclearZ_[nucleusRef_[i2]] = nuclearZ_[ref1];
- } else {
- // neither i1 nor i2 are nuclei
- }
- } else {
- if (ref1 >= 0) {
- // both are nuclei, i2 is overwritten with nucleus i1
- // fNucleusRef stays the same, but A and Z data is overwritten
- nuclearA_[ref2] = nuclearA_[ref1];
- nuclearZ_[ref2] = nuclearZ_[ref1];
- } else {
- // i2 is overwritten with non-nucleus i1
- nucleusRef_[i2] = -1; // flag as non-nucleus
- nuclearA_.erase(nuclearA_.cbegin() + ref2); // remove data for i2
- nuclearZ_.erase(nuclearZ_.cbegin() + ref2); // remove data for i2
- const int n = nucleusRef_.size(); // update fNucleusRef: indices above ref2
- // must be decremented by 1
- for (int i = 0; i < n; ++i) {
- if (nucleusRef_[i] > ref2) { nucleusRef_[i] -= 1; }
- }
- }
- }
- }
-
- template <typename InnerStackImpl>
- inline void NuclearStackExtensionImpl<InnerStackImpl>::clear() {
- super_type::clear();
- nucleusRef_.clear();
- nuclearA_.clear();
- nuclearZ_.clear();
- }
-
- template <typename InnerStackImpl>
- inline void NuclearStackExtensionImpl<InnerStackImpl>::swap(const unsigned int i1,
- const unsigned int i2) {
- // index range check
- if (i1 >= getSize() || i2 >= getSize()) {
- std::ostringstream err;
- err << "NuclearStackExtension: trying to access data beyond size of stack!";
- throw std::runtime_error(err.str());
- }
- // swap original particle data
- super_type::swap(i1, i2);
- // swap corresponding nuclear reference data
- std::swap(nucleusRef_[i2], nucleusRef_[i1]);
- }
-
- template <typename InnerStackImpl>
- inline void NuclearStackExtensionImpl<InnerStackImpl>::incrementSize() {
- super_type::incrementSize();
- nucleusRef_.push_back(-1);
- }
-
- template <typename InnerStackImpl>
- inline void NuclearStackExtensionImpl<InnerStackImpl>::decrementSize() {
- super_type::decrementSize();
- if (nucleusRef_.size() > 0) {
- const int ref = nucleusRef_.back();
- nucleusRef_.pop_back();
- if (ref >= 0) {
- nuclearA_.erase(nuclearA_.begin() + ref);
- nuclearZ_.erase(nuclearZ_.begin() + ref);
- const int n = nucleusRef_.size();
- for (int i = 0; i < n; ++i) {
- if (nucleusRef_[i] >= ref) { nucleusRef_[i] -= 1; }
- }
- }
- }
- }
-
-} // namespace corsika::nuclear_stack
diff --git a/corsika/stack/DummyStack.hpp b/corsika/stack/DummyStack.hpp
index bcba30f1d..a37fe53d7 100644
--- a/corsika/stack/DummyStack.hpp
+++ b/corsika/stack/DummyStack.hpp
@@ -10,6 +10,7 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/stack/Stack.hpp>
#include <string>
@@ -22,8 +23,8 @@ namespace corsika::dummy_stack {
*/
/**
- however, conceptually we need to provide fake data. A stack without data does not
- work...
+ * However, conceptually we need to provide fake data. A stack without data does not
+ * work...
*/
struct NoData { /* nothing */
@@ -40,6 +41,14 @@ namespace corsika::dummy_stack {
void setParticleData(super_type& /*parent*/, const std::tuple<NoData>& /*v*/) {}
std::string asString() const { return "dummy-data"; }
+
+ // unfortunately we need those dummy getter
+ // for some more complex tests with "history"
+ HEPEnergyType getEnergy() const { return 0_GeV; }
+ MomentumVector getMomentum() const {
+ return MomentumVector(get_root_CoordinateSystem(), {0_GeV, 0_GeV, 0_GeV});
+ }
+ Code getPID() const { return Code::Unknown; }
};
/**
@@ -67,7 +76,7 @@ namespace corsika::dummy_stack {
int getCapacity() const { return entries_; }
/**
- * Function to copy particle at location i2 in stack to i1
+ * Function to copy particle at location i2 in stack to i1.
*/
void copy(const int /*i1*/, const int /*i2*/) {}
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
deleted file mode 100644
index 03870ccf5..000000000
--- a/corsika/stack/NuclearStackExtension.hpp
+++ /dev/null
@@ -1,307 +0,0 @@
-/*
- * (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.
- */
-
-#pragma once
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/stack/Stack.hpp>
-#include <corsika/framework/geometry/Point.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
-#include <corsika/stack/VectorStack.hpp>
-
-#include <algorithm>
-#include <tuple>
-#include <vector>
-
-namespace corsika::nuclear_stack {
-
- /**
- *
- * Define ParticleInterface for NuclearStackExtension Stack derived from
- * ParticleInterface of Inner stack class
- *
- * Add A and Z data to existing stack (currently VectorStack) of particle
- * properties. This is done via inheritance, not via CombinedStack since the nuclear
- * data is stored ONLY when needed (for nuclei) and not for all particles. Thus, this is
- * a new, derived Stack object.
- *
- * Only for Code::Nucleus particles A and Z are stored, not for all
- * normal elementary particles.
- *
- * Thus in your code, make sure to always check <code>
- * particle.getPID()==Code::Nucleus </code> before attempting to
- * read any nuclear information.
- */
- template <template <typename> class InnerParticleInterface,
- typename StackIteratorInterface>
- struct NuclearParticleInterface
- : public InnerParticleInterface<StackIteratorInterface> {
-
- typedef InnerParticleInterface<StackIteratorInterface> super_type;
-
- public:
- typedef std::tuple<Code, HEPEnergyType, DirectionVector, Point, TimeType,
- unsigned short, unsigned short>
- nuclear_particle_data_type;
-
- typedef std::tuple<Code, MomentumVector, Point, TimeType, unsigned short,
- unsigned short>
- nuclear_particle_data_momentum_type;
-
- /**
- *
- * @param v which is a tuple containing: PID, kinetic Energy, DirectionVector,
- * Position, Time
- */
- void setParticleData(typename super_type::particle_data_type const& v);
-
- /**
- *
- * @param v which is a tuple containing: PID, kinetic Energy, DirectionVector,
- * Position, Time, A, Z
- */
- void setParticleData(nuclear_particle_data_type const& v);
-
- /**
- *
- * @param p the parent particle
- * @param v which is a tuple containing: PID, Momentum Vector, Position,
- * Time
- */
- void setParticleData(super_type& p, typename super_type::particle_data_type const& v);
-
- /**
- *
- * @param p the parent particle
- * @param v which is a tuple containing: PID, Momentum Vector, Position,
- * Time, A, Z
- */
- void setParticleData(super_type& p, nuclear_particle_data_type const& v);
-
- /**
- *
- * @param v which is a tuple containing: PID, Total Energy, MomentumVector, Position,
- * Time
- */
- void setParticleData(typename super_type::particle_data_momentum_type const& v);
-
- /**
- *
- * @param v which is a tuple containing: PID, Total Energy, MomentumVector, Position,
- * Time, A, Z
- */
- void setParticleData(nuclear_particle_data_momentum_type const& v);
-
- /**
- *
- * @param p parent particle
- * @param v which is a tuple containing: PID, Total Energy, MomentumVector, Position,
- * Time
- */
- void setParticleData(super_type& p,
- typename super_type::particle_data_momentum_type const& v);
-
- /**
- *
- * @param p parent particle
- * @param v which is a tuple containing: PID, Total Energy, MomentumVector, Position,
- * Time, A, Z
- */
- void setParticleData(super_type& p, nuclear_particle_data_momentum_type const& v);
-
- std::string asString() const;
-
- /**
- * @name individual setters
- * @{
- */
- void setNuclearA(const unsigned short vA) {
- super_type::getStackData().setNuclearA(super_type::getIndex(), vA);
- }
- void setNuclearZ(const unsigned short vZ) {
- super_type::getStackData().setNuclearZ(super_type::getIndex(), vZ);
- }
- /// @}
-
- /**
- * @name individual getters
- * @{
- */
- int getNuclearA() const {
- return super_type::getStackData().getNuclearA(super_type::getIndex());
- }
- int getNuclearZ() const {
- return super_type::getStackData().getNuclearZ(super_type::getIndex());
- }
- /// @}
-
- /**
- * Overwrite normal getPDG function with nuclear version
- */
- PDGCode getPDG() const;
-
- /**
- * Overwrite normal setMomentum function with nuclear version
- */
- void setMomentum(MomentumVector const& v);
-
- /**
- * Overwrite normal getMomentum function with nuclear version
- */
- MomentumVector getMomentum() const;
-
- /**
- * Overwrite normal getEnergy function with nuclear version
- */
- void setEnergy(HEPEnergyType const& e);
-
- /**
- * Overwrite normal getVelocity function with nuclear version
- */
- VelocityVector getVelocity() const;
-
- /**
- * Overwrite normal getMass function with nuclear version
- */
- HEPMassType getMass() const;
-
- /**
- * Overwrite normal getParticleCharge function with nuclear version
- */
- ElectricChargeType getCharge() const;
-
- /**
- * Overwrite normal getEnergy function with nuclear version
- */
- HEPEnergyType getEnergy() const;
-
- /**
- * Overwirte normal getChargeNumber function with nuclear version
- **/
- int16_t getChargeNumber() const;
-
- int getNucleusRef() const {
- return super_type::getStackData().getNucleusRef(super_type::getIndex());
- } // LCOV_EXCL_LINE
-
- protected:
- void setNucleusRef(const int vR) {
- super_type::getStackData().setNucleusRef(super_type::getIndex(), vR);
- }
-
- bool isNucleus() const {
- return super_type::getStackData().isNucleus(super_type::getIndex());
- }
- };
-
- /**
- * @class NuclearStackExtension
- *
- * Memory implementation of adding nuclear inforamtion to the
- * existing particle stack defined in class InnerStackImpl.
- *
- * Inside the NuclearStackExtension class there is a dedicated
- * fNucleusRef index, where fNucleusRef[i] is referring to the
- * correct A and Z for a specific particle index i. fNucleusRef[i]
- * == -1 means that this is not a nucleus, and a subsequent call to
- * getNucleusA would produce an exception.
- */
- template <typename InnerStackImpl>
- class NuclearStackExtensionImpl : public InnerStackImpl {
-
- typedef InnerStackImpl super_type;
-
- public:
- typedef std::vector<int> nucleus_ref_type;
- typedef std::vector<unsigned short> nuclear_a_type;
- typedef std::vector<unsigned short> nuclear_z_type;
-
- NuclearStackExtensionImpl() = default;
-
- NuclearStackExtensionImpl(NuclearStackExtensionImpl<InnerStackImpl> const&) = default;
-
- NuclearStackExtensionImpl(NuclearStackExtensionImpl<InnerStackImpl>&&) = default;
-
- NuclearStackExtensionImpl<InnerStackImpl>& operator=(
- NuclearStackExtensionImpl<InnerStackImpl> const&) = default;
-
- NuclearStackExtensionImpl<InnerStackImpl>& operator=(
- NuclearStackExtensionImpl<InnerStackImpl>&&) = default;
-
- void init() { super_type::init(); }
-
- void dump() { super_type::dump(); }
-
- void clear();
-
- unsigned int getSize() const { return nucleusRef_.size(); }
-
- unsigned int getCapacity() const { return nucleusRef_.capacity(); }
-
- void setNuclearA(const unsigned int i, const unsigned short vA) {
- nuclearA_[getNucleusRef(i)] = vA;
- }
-
- void setNuclearZ(const unsigned int i, const unsigned short vZ) {
- nuclearZ_[getNucleusRef(i)] = vZ;
- }
-
- void setNucleusRef(const unsigned int i, const int v) { nucleusRef_[i] = v; }
-
- int getNuclearA(const unsigned int i) const { return nuclearA_[getNucleusRef(i)]; }
-
- int getNuclearZ(const unsigned int i) const { return nuclearZ_[getNucleusRef(i)]; }
-
- // this function will create new storage for Nuclear Properties, and return the
- // reference to it
- int getNucleusNextRef();
-
- int getNucleusRef(const unsigned int i) const;
-
- bool isNucleus(const unsigned int i) const { return nucleusRef_[i] >= 0; }
-
- /**
- * Function to copy particle at location i1 in stack to i2
- */
- void copy(const unsigned int i1, const unsigned int i2);
- /**
- * Function to copy particle at location i2 in stack to i1
- */
- void swap(const unsigned int i1, const unsigned int i2);
-
- void incrementSize();
-
- void decrementSize();
-
- private:
- /// the actual memory to store particle data
-
- nucleus_ref_type nucleusRef_;
- nuclear_a_type nuclearA_;
- nuclear_z_type nuclearZ_;
-
- }; // end class NuclearStackExtensionImpl
-
- template <typename TInnerStack, template <typename> typename PI_>
- using NuclearStackExtension =
- Stack<NuclearStackExtensionImpl<typename TInnerStack::stack_data_type>, PI_>;
-
- //
- template <typename TStackIter>
- using ExtendedParticleInterfaceType =
- NuclearParticleInterface<VectorStack::pi_type, TStackIter>;
-
- // the particle data stack with extra nuclear information:
- using ParticleDataStack =
- NuclearStackExtension<VectorStack, ExtendedParticleInterfaceType>;
-
-} // namespace corsika::nuclear_stack
-
-#include <corsika/detail/stack/NuclearStackExtension.inl>
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index 32dc0af17..f58f6cd56 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -10,6 +10,7 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/stack/Stack.hpp>
#include <corsika/framework/geometry/Point.hpp>
@@ -63,8 +64,6 @@ namespace corsika {
*
* @param v tuple containing of type particle_data_type
*
- * MomentumVector is only used to determine the DirectionVector, the normalization
- * is lost.
*/
void setParticleData(particle_data_type const& v);
@@ -102,23 +101,6 @@ namespace corsika {
super_type::getStackData().setKineticEnergy(super_type::getIndex(), ekin);
}
- /**
- * The MomentumVector v is used to determine the DirectionVector, and to update the
- * particle energy.
- */
- void setMomentum(MomentumVector const& v) {
- HEPMomentumType const P = v.getNorm();
- if (P == 0_eV) {
- super_type::getStackData().setKineticEnergy(super_type::getIndex(), 0_eV);
- super_type::getStackData().setDirection(
- super_type::getIndex(), DirectionVector(v.getCoordinateSystem(), {0, 0, 0}));
- } else {
- super_type::getStackData().setKineticEnergy(
- super_type::getIndex(),
- sqrt(square(getMass()) + square(P)) - this->getMass());
- super_type::getStackData().setDirection(super_type::getIndex(), v / P);
- }
- }
//! Set direction
void setDirection(DirectionVector const& v) {
super_type::getStackData().setDirection(super_type::getIndex(), v);
@@ -165,7 +147,7 @@ namespace corsika {
}
//! Get momentum
MomentumVector getMomentum() const {
- auto const P = sqrt(square(getEnergy()) - square(this->getMass()));
+ auto const P = calculate_momentum(this->getEnergy(), this->getMass());
return super_type::getStackData().getDirection(super_type::getIndex()) * P;
}
//! Get mass of particle
diff --git a/corsika/stack/history/HistoryStackExtension.hpp b/corsika/stack/history/HistoryStackExtension.hpp
index 9d517b0ac..4dc5f1df9 100644
--- a/corsika/stack/history/HistoryStackExtension.hpp
+++ b/corsika/stack/history/HistoryStackExtension.hpp
@@ -67,10 +67,8 @@ namespace corsika::history {
/**
* @class HistoryDataInterface
*
- * corresponding defintion of a stack-readout object, the iteractor
- * dereference operator will deliver access to these function
- // defintion of a stack-readout object, the iteractor dereference
- // operator will deliver access to these function
+ * corresponding definition of a stack-readout object, the iteractor
+ * dereference operator will deliver access to these function.
*/
template <typename T, typename TEvent>
class HistoryDataInterface : public T {
diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp
index e350c2fb2..6411d9738 100644
--- a/tests/modules/testObservationPlane.cpp
+++ b/tests/modules/testObservationPlane.cpp
@@ -100,7 +100,7 @@ TEST_CASE("ObservationPlane", "interface") {
MomentumVector const pnew = MomentumVector(cs, {1_GeV, 0.5_GeV, -0.4_GeV});
HEPEnergyType const enew = sqrt(pnew.dot(pnew));
- particle.setMomentum(pnew);
+ particle.setDirection(pnew.normalized());
particle.setEnergy(enew);
Plane const obsPlane(Point(cs, {10_m, 5_m, 5_m}),
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index 1710844b1..4906bc58c 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -35,31 +35,30 @@ int sgn(T val) {
}
/**
- @file testTracking.cpp
-
- This is the unified and common unit test for all Tracking algorithms:
-
- - tracking_leapfrog_curved::Tracking
- - tracking_leapfrog_straight::Tracking
- - tracking_line::Tracking
-
-
- The main part of tests are to inject particles at 10GeV momentum at
- (-Rg,0,0) in +x direction into a sphere of radius Rg, where Rg is
- the gyroradius (or 10m for neutral particles). Then it is checked
- where the particles leaves the sphere for different charges
- (-1,0,+1) and field strength (-50uT, 0T, +50uT).
-
- Each test is perfromed once, with the particles starting logically
- outside of the Rg sphere (thus it first has to enter insides) and a
- second time with the particle already logically inside the sphere.
-
- There is a second smaller, internal sphere at +z displacement. Only
- neutral particles are allowed and expected to hit this.
-
- All those tests are parameterized, thus, they can be easily extended
- or applied to new algorithms.
-
+ * @file testTracking.cpp
+ *
+ * This is the unified and common unit test for all Tracking algorithms:
+ *
+ * - tracking_leapfrog_curved::Tracking
+ * - tracking_leapfrog_straight::Tracking
+ * - tracking_line::Tracking
+ *
+ *
+ * The main part of tests are to inject particles at 10GeV momentum at
+ * (-Rg,0,0) in +x direction into a sphere of radius Rg, where Rg is
+ * the gyroradius (or 10m for neutral particles). Then it is checked
+ * where the particles leaves the sphere for different charges
+ * (-1,0,+1) and field strength (-50uT, 0T, +50uT).
+ *
+ * Each test is perfromed once, with the particles starting logically
+ * outside of the Rg sphere (thus it first has to enter insides) and a
+ * second time with the particle already logically inside the sphere.
+ *
+ * There is a second smaller, internal sphere at +z displacement. Only
+ * neutral particles are allowed and expected to hit this.
+ *
+ * All those tests are parameterized, thus, they can be easily extended
+ * or applied to new algorithms.
*/
TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
@@ -185,7 +184,7 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
auto [traj, nextVol] = tracking.getTrack(particle);
particle.setNode(nextVol);
particle.setPosition(traj.getPosition(1));
- particle.setMomentum(traj.getDirection(1) * particle.getMomentum().getNorm());
+ particle.setDirection(traj.getDirection(1));
SpeedType const speed_0 = particle.getVelocity().getNorm();
if (outer) {
// now we know we are in target volume, depending on "outer"
@@ -206,7 +205,7 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
nextVol = nextVol2;
particle.setNode(nextVol);
particle.setPosition(traj2.getPosition(1));
- particle.setMomentum(traj2.getDirection(1) * particle.getMomentum().getNorm());
+ particle.setDirection(traj2.getDirection(1));
CORSIKA_LOG_TRACE("pos={}, p={}, |p|={} |v|={}, delta-l={}, delta-t={}",
particle.getPosition(), particle.getMomentum(),
particle.getMomentum().getNorm(),
@@ -234,7 +233,7 @@ TEMPLATE_TEST_CASE("Tracking", "tracking", tracking_leapfrog_curved::Tracking,
}
}
-/** specifc test for curved leap-frog algorithm **/
+/** specifc test for curved leap-frog algorithm. */
TEST_CASE("TrackingLeapFrogCurved") {
diff --git a/tests/stack/CMakeLists.txt b/tests/stack/CMakeLists.txt
index 1314fd435..149aa0533 100644
--- a/tests/stack/CMakeLists.txt
+++ b/tests/stack/CMakeLists.txt
@@ -1,6 +1,7 @@
set (test_stack_sources
TestMain.cpp
testHistoryStack.cpp
+ testHistoryView.cpp
testGeometryNodeStackExtension.cpp
testWeightStackExtension.cpp
testDummyStack.cpp
diff --git a/tests/stack/testHistoryView.cpp b/tests/stack/testHistoryView.cpp
index 1633b8212..198e98993 100644
--- a/tests/stack/testHistoryView.cpp
+++ b/tests/stack/testHistoryView.cpp
@@ -9,52 +9,34 @@
#include <corsika/stack/history/Event.hpp>
#include <corsika/stack/history/HistorySecondaryProducer.hpp>
#include <corsika/stack/history/HistoryStackExtension.hpp>
-
-#include <corsika/framework/stack/CombinedStack.hpp>
+#include <corsika/stack/history/HistorySecondaryProducer.hpp>
#include <corsika/stack/DummyStack.hpp>
-#include <corsika/stack/NuclearStackExtension.hpp>
-
+#include <corsika/framework/stack/CombinedStack.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <catch2/catch.hpp>
using namespace corsika;
-/**
- Need to replicate setup::SetupStack in a maximally simplified
- way, but with real particle data
- */
+// the GeometryNode stack needs to know the type of geometry-nodes from the DummyEnv:
+template <typename TStackIter>
+using DummyHistoryDataInterface =
+ typename history::MakeHistoryDataInterface<TStackIter, history::Event>::type;
// combine dummy stack with geometry information for tracking
template <typename TStackIter>
using StackWithHistoryInterface =
- CombinedParticleInterface<nuclear_stack::ParticleDataStack::pi_type,
- history::HistoryEventDataInterface, TStackIter>;
+ CombinedParticleInterface<dummy_stack::DummyStack::pi_type, DummyHistoryDataInterface,
+ TStackIter>;
using TestStack =
- CombinedStack<typename nuclear_stack::ParticleDataStack::stack_data_type,
- history::HistoryEventData, StackWithHistoryInterface>;
-
-/*
- See Issue 161
-
- unfortunately clang does not support this in the same way (yet) as
- gcc, so we have to distinguish here. If clang cataches up, we
- could remove the clang branch here and also in
- corsika::Cascade. The gcc code is much more generic and
- universal. If we could do the gcc version, we won't had to define
- StackView globally, we could do it with MakeView whereever it is
- actually needed. Keep an eye on this!
- */
-#if defined(__clang__)
-using TheTestStackView =
- SecondaryView<typename TestStack::stack_data_type, StackWithHistoryInterface,
+ CombinedStack<typename dummy_stack::DummyStack::stack_data_type,
+ history::HistoryData<history::Event>, StackWithHistoryInterface,
history::HistorySecondaryProducer>;
-#elif defined(__GNUC__) || defined(__GNUG__)
-using TheTestStackView = MakeView<TestStack, history::HistorySecondaryProducer>::type;
-#endif
-using TestStackView = TheTestStackView;
+// the correct secondary stack view
+using TestStackView = typename TestStack::stack_view_type;
+using EvtPtr = std::shared_ptr<history::Event>;
template <typename Event>
int count_generations(Event const* event) {
@@ -70,27 +52,24 @@ int count_generations(Event const* event) {
TEST_CASE("HistoryStackExtensionView", "[stack]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
- CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
// in this test we only use one singel stack !
+ const dummy_stack::NoData noData;
TestStack stack;
// add primary particle
- auto p0 = stack.addParticle(
- std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto p0 = stack.addParticle(std::make_tuple(noData));
CHECK(stack.getEntries() == 1);
- corsika::history::EventPtr evt = p0.getEvent();
+ EvtPtr evt = p0.getEvent();
CHECK(evt == nullptr);
CHECK(count_generations(evt.get()) == 0);
SECTION("interface test, view") {
// add secondaries, 1st generation
- TestStackView hview0(p0);
+ auto pnext = stack.getNextParticle();
+ TestStackView hview0(pnext);
auto const ev0 = p0.getEvent();
CHECK(ev0 == nullptr);
@@ -99,12 +78,10 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
// add 5 secondaries
for (int i = 0; i < 5; ++i) {
- auto sec = hview0.addSecondary(
- std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = hview0.addSecondary(std::make_tuple(noData));
CHECK(sec.getParentEventIndex() == i);
- CHECK(sec.getEvent() != nullptr);
+ CHECK(sec.getEvent().get() != nullptr);
CHECK(sec.getEvent()->parentEvent() == nullptr);
CHECK(count_generations(sec.getEvent().get()) == 1);
}
@@ -120,9 +97,7 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
// add second generation of secondaries
// add 10 secondaries
for (int i = 0; i < 10; ++i) {
- auto sec = hview1.addSecondary(
- std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = hview1.addSecondary(std::make_tuple(noData));
CHECK(sec.getParentEventIndex() == i);
CHECK(sec.getEvent()->parentEvent() == ev1);
@@ -146,9 +121,7 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
for (int i = 0; i < 15; ++i) {
CORSIKA_LOG_TRACE("loop, view: " + std::to_string(i));
- auto sec = hview2.addSecondary(
- std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = hview2.addSecondary(std::make_tuple(noData));
CORSIKA_LOG_TRACE("loop, ---- ");
CHECK(sec.getParentEventIndex() == i);
@@ -175,9 +148,7 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
// add 5 secondaries
for (int i = 0; i < 5; ++i) {
CORSIKA_LOG_TRACE("loop " + std::to_string(i));
- auto sec = proj0.addSecondary(
- std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = proj0.addSecondary(std::make_tuple(noData));
CHECK(sec.getParentEventIndex() == i);
CHECK(sec.getEvent() != nullptr);
@@ -195,9 +166,7 @@ TEST_CASE("HistoryStackExtensionView", "[stack]") {
// add second generation of secondaries
// add 10 secondaries
for (unsigned int i = 0; i < 10; ++i) {
- auto sec = proj1.addSecondary(
- std::make_tuple(Code::Electron, MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
- Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s));
+ auto sec = proj1.addSecondary(std::make_tuple(noData));
CHECK(sec.getParentEventIndex() == int(i));
CHECK(sec.getEvent()->parentEvent() == ev1);
--
GitLab