IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 8883a66e authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

Merge branch '82-handling-of-boundary-crossings-in-geometry-tree' of... 

Merge branch '82-handling-of-boundary-crossings-in-geometry-tree' of gitlab.ikp.kit.edu:AirShowerPhysics/corsika into 82-handling-of-boundary-crossings-in-geometry-tree
parents eb49cb5b 9bcb3367
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/StackInterface/ParticleBase.h
+ 2
− 1
View file @ 8883a66e
...@@ -117,13 +117,14 @@ namespace corsika::stack { ...@@ -117,13 +117,14 @@ namespace corsika::stack {
using T::GetStackData; using T::GetStackData;
public: public:
/*
template <typename... Args1, typename... Args2> template <typename... Args1, typename... Args2>
void SetParticleData(Args1... args1, Args2... args2) { void SetParticleData(Args1... args1, Args2... args2) {
T::SetParticleData(args1...); T::SetParticleData(args1...);
} }
template <typename... Args1, typename... Args2> template <typename... Args1, typename... Args2>
void SetParticleData(T& p, Args1... args1, Args2... args2) {} void SetParticleData(T& p, Args1... args1, Args2... args2) {}
*/
}; };
} // namespace corsika::stack } // namespace corsika::stack
......
Framework/StackInterface/Stack.h
+ 13
− 7
View file @ 8883a66e
...@@ -67,11 +67,17 @@ namespace corsika::stack { ...@@ -67,11 +67,17 @@ namespace corsika::stack {
delete; ///< since Stack can be very big, we don't want to copy it delete; ///< since Stack can be very big, we don't want to copy it
public: public:
// Stack() { Init(); }
/** /**
* if StackDataType is a reference member we *have* to initialize * if StackDataType is a reference member we *HAVE* to initialize
* it in the constructor, this is typically needed for SecondaryView * it in the constructor, this is typically needed for SecondaryView
*/ */
template <typename = std::enable_if_t<std::is_reference<StackDataType>{}>> template <
typename _StackDataType = StackDataType,
typename = std::enable_if<std::is_same<StackDataType, _StackDataType>::value &&
std::is_reference<_StackDataType>::value,
void>>
Stack(StackDataType vD) Stack(StackDataType vD)
: fData(vD) {} : fData(vD) {}
...@@ -80,8 +86,11 @@ namespace corsika::stack { ...@@ -80,8 +86,11 @@ namespace corsika::stack {
* StackDataType user class. If the user did not provide a suited * StackDataType user class. If the user did not provide a suited
* constructor this will fail with an error message. * constructor this will fail with an error message.
*/ */
template <typename... Args, template <
typename = std::enable_if_t<!std::is_reference<StackDataType>{}>> typename... Args, typename _StackDataType = StackDataType,
typename = std::enable_if<std::is_same<StackDataType, _StackDataType>::value &&
!std::is_reference<_StackDataType>::value,
void>>
Stack(Args... args) Stack(Args... args)
: fData(args...) {} : fData(args...) {}
...@@ -223,9 +232,6 @@ namespace corsika::stack { ...@@ -223,9 +232,6 @@ namespace corsika::stack {
StackIterator GetNextParticle() { return last(); } StackIterator GetNextParticle() { return last(); }
protected: protected:
// typename std::enable_if<HasGetIndexFromIterator<T>::value, unsigned int>::type
// typename std::enable_if<std::is_base_of<decltype(*this)>,
// SecondaryView<StackDataType, ParticleInterface>>::value, unsigned int>::type
/** /**
* Function to perform eventual transformation from * Function to perform eventual transformation from
* StackIterator::GetIndex() to index in data stored in * StackIterator::GetIndex() to index in data stored in
......
Framework/StackInterface/testStackInterface.cc
+ 2
− 3
View file @ 8883a66e
...@@ -33,7 +33,6 @@ using namespace corsika::stack; ...@@ -33,7 +33,6 @@ using namespace corsika::stack;
using namespace std; using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest; typedef Stack<TestStackData, TestParticleInterface> StackTest;
typedef StackTest::ParticleInterfaceType Particle;
TEST_CASE("Stack", "[Stack]") { TEST_CASE("Stack", "[Stack]") {
...@@ -48,7 +47,6 @@ TEST_CASE("Stack", "[Stack]") { ...@@ -48,7 +47,6 @@ TEST_CASE("Stack", "[Stack]") {
// construct a valid Stack object // construct a valid Stack object
StackTest s; StackTest s;
s.Init();
s.Clear(); s.Clear();
s.AddParticle(std::tuple{0.}); s.AddParticle(std::tuple{0.});
s.Copy(s.cbegin(), s.begin()); s.Copy(s.cbegin(), s.begin());
...@@ -98,7 +96,8 @@ TEST_CASE("Stack", "[Stack]") { ...@@ -98,7 +96,8 @@ TEST_CASE("Stack", "[Stack]") {
StackTest s; StackTest s;
REQUIRE(s.GetSize() == 0); REQUIRE(s.GetSize() == 0);
auto iter = s.AddParticle(std::tuple{9.9}); auto iter = s.AddParticle(std::tuple{9.9});
Particle& p = *iter; // also this is valid to access particle data StackTest::ParticleInterfaceType& p =
*iter; // also this is valid to access particle data
REQUIRE(s.GetSize() == 1); REQUIRE(s.GetSize() == 1);
p.AddSecondary(std::tuple{4.4}); p.AddSecondary(std::tuple{4.4});
REQUIRE(s.GetSize() == 2); REQUIRE(s.GetSize() == 2);
......
Framework/Utilities/testCOMBoost.cc
+ 87
− 88
View file @ 8883a66e
...@@ -32,99 +32,98 @@ using corsika::units::constants::cSquared; ...@@ -32,99 +32,98 @@ using corsika::units::constants::cSquared;
double constexpr absMargin = 1e-6; double constexpr absMargin = 1e-6;
CoordinateSystem const& rootCS = CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// helper function for energy-momentum
// relativistic energy
auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
return sqrt(m * m + p.squaredNorm());
};
// helper function for mandelstam-s // helper function for energy-momentum
auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) { // relativistic energy
return E * E - p.squaredNorm(); auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
}; return sqrt(m * m + p.squaredNorm());
};
// helper function for mandelstam-s
auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
return E * E - p.squaredNorm();
};
TEST_CASE("rotation") { TEST_CASE("rotation") {
// define projectile kinematics in lab frame // define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV; HEPMassType const projectileMass = 1_GeV;
HEPMassType const targetMass = 1.0e300_eV; HEPMassType const targetMass = 1.0e300_eV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, 1_GeV}}; Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, 1_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab); HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab); const FourVector PprojLab(eProjectileLab, pProjectileLab);
Eigen::Vector3d e1, e2, e3; Eigen::Vector3d e1, e2, e3;
e1 << 1, 0, 0; e1 << 1, 0, 0;
e2 << 0, 1, 0; e2 << 0, 1, 0;
e3 << 0, 0, 1; e3 << 0, 0, 1;
// define boost to com frame
SECTION("pos. z-axis") {
COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 0_GeV, 1_GeV}}}, targetMass);
auto const& rot = boost.GetRotationMatrix();
CHECK((rot * e3 - e3).norm() == Approx(0).margin(absMargin));
CHECK((rot * e1).norm() == Approx(1));
CHECK((rot * e2).norm() == Approx(1));
CHECK((rot * e3).norm() == Approx(1));
CHECK(rot.determinant() == Approx(1));
}
// define boost to com frame SECTION("y-axis in upper half") {
SECTION("pos. z-axis") { COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 1_GeV, 1_meV}}}, targetMass);
COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 0_GeV, 1_GeV}}}, targetMass); auto const& rot = boost.GetRotationMatrix();
auto const& rot = boost.GetRotationMatrix();
CHECK((rot * e2 - e3).norm() == Approx(0).margin(absMargin));
CHECK((rot * e3 - e3).norm() == Approx(0).margin(absMargin)); CHECK((rot * e1).norm() == Approx(1));
CHECK((rot * e1).norm() == Approx(1)); CHECK((rot * e2).norm() == Approx(1));
CHECK((rot * e2).norm() == Approx(1)); CHECK((rot * e3).norm() == Approx(1));
CHECK((rot * e3).norm() == Approx(1)); CHECK(rot.determinant() == Approx(1));
CHECK(rot.determinant() == Approx(1)); }
}
SECTION("x-axis in upper half") {
SECTION("y-axis in upper half") { COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, 1_meV}}}, targetMass);
COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 1_GeV, 1_meV}}}, targetMass); auto const& rot = boost.GetRotationMatrix();
auto const& rot = boost.GetRotationMatrix();
CHECK((rot * e1 - e3).norm() == Approx(0).margin(absMargin));
CHECK((rot * e2 - e3).norm() == Approx(0).margin(absMargin)); CHECK((rot * e1).norm() == Approx(1));
CHECK((rot * e1).norm() == Approx(1)); CHECK((rot * e2).norm() == Approx(1));
CHECK((rot * e2).norm() == Approx(1)); CHECK((rot * e3).norm() == Approx(1));
CHECK((rot * e3).norm() == Approx(1)); CHECK(rot.determinant() == Approx(1));
CHECK(rot.determinant() == Approx(1)); }
}
SECTION("neg. z-axis") {
SECTION("x-axis in upper half") { COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 0_GeV, -1_GeV}}}, targetMass);
COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, 1_meV}}}, targetMass); auto const& rot = boost.GetRotationMatrix();
auto const& rot = boost.GetRotationMatrix();
CHECK((rot * (-e3) - e3).norm() == Approx(0).margin(absMargin));
CHECK((rot * e1 - e3).norm() == Approx(0).margin(absMargin)); CHECK((rot * e1).norm() == Approx(1));
CHECK((rot * e1).norm() == Approx(1)); CHECK((rot * e2).norm() == Approx(1));
CHECK((rot * e2).norm() == Approx(1)); CHECK((rot * e3).norm() == Approx(1));
CHECK((rot * e3).norm() == Approx(1)); CHECK(rot.determinant() == Approx(1));
CHECK(rot.determinant() == Approx(1)); }
}
SECTION("x-axis lower half") {
SECTION("neg. z-axis") { COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, -1_meV}}}, targetMass);
COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 0_GeV, -1_GeV}}}, targetMass); auto const& rot = boost.GetRotationMatrix();
auto const& rot = boost.GetRotationMatrix();
CHECK((rot * e1 - e3).norm() == Approx(0).margin(absMargin));
CHECK((rot * (-e3) - e3).norm() == Approx(0).margin(absMargin)); CHECK((rot * e1).norm() == Approx(1));
CHECK((rot * e1).norm() == Approx(1)); CHECK((rot * e2).norm() == Approx(1));
CHECK((rot * e2).norm() == Approx(1)); CHECK((rot * e3).norm() == Approx(1));
CHECK((rot * e3).norm() == Approx(1)); CHECK(rot.determinant() == Approx(1));
CHECK(rot.determinant() == Approx(1)); }
}
SECTION("y-axis lower half") {
SECTION("x-axis lower half") { COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 1_GeV, -1_meV}}}, targetMass);
COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, -1_meV}}}, targetMass); auto const& rot = boost.GetRotationMatrix();
auto const& rot = boost.GetRotationMatrix();
CHECK((rot * e2 - e3).norm() == Approx(0).margin(absMargin));
CHECK((rot * e1 - e3).norm() == Approx(0).margin(absMargin)); CHECK((rot * e1).norm() == Approx(1));
CHECK((rot * e1).norm() == Approx(1)); CHECK((rot * e2).norm() == Approx(1));
CHECK((rot * e2).norm() == Approx(1)); CHECK((rot * e3).norm() == Approx(1));
CHECK((rot * e3).norm() == Approx(1)); CHECK(rot.determinant() == Approx(1));
CHECK(rot.determinant() == Approx(1)); }
}
SECTION("y-axis lower half") {
COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 1_GeV, -1_meV}}}, targetMass);
auto const& rot = boost.GetRotationMatrix();
CHECK((rot * e2 - e3).norm() == Approx(0).margin(absMargin));
CHECK((rot * e1).norm() == Approx(1));
CHECK((rot * e2).norm() == Approx(1));
CHECK((rot * e3).norm() == Approx(1));
CHECK(rot.determinant() == Approx(1));
}
} }
TEST_CASE("boosts") { TEST_CASE("boosts") {
......
Processes/Sibyll/Decay.cc
+ 2
− 2
View file @ 8883a66e
...@@ -157,8 +157,6 @@ namespace corsika::process::sibyll { ...@@ -157,8 +157,6 @@ namespace corsika::process::sibyll {
// remember position // remember position
Point const decayPoint = p.GetPosition(); Point const decayPoint = p.GetPosition();
TimeType const t0 = p.GetTime(); TimeType const t0 = p.GetTime();
// remove original particle from corsika stack
p.Delete();
// set all particles/hadrons unstable // set all particles/hadrons unstable
// setHadronsUnstable(); // setHadronsUnstable();
setUnstable(pCode); setUnstable(pCode);
...@@ -184,6 +182,8 @@ namespace corsika::process::sibyll { ...@@ -184,6 +182,8 @@ namespace corsika::process::sibyll {
} }
// empty sibyll stack // empty sibyll stack
ss.Clear(); ss.Clear();
// remove original particle from corsika stack
p.Delete();
} }
} // namespace corsika::process::sibyll } // namespace corsika::process::sibyll
Processes/Sibyll/Interaction.cc
+ 2
− 3
View file @ 8883a66e
...@@ -320,9 +320,6 @@ namespace corsika::process::sibyll { ...@@ -320,9 +320,6 @@ namespace corsika::process::sibyll {
sib_list_(print_unit); sib_list_(print_unit);
fNucCount += get_nwounded() - 1; fNucCount += get_nwounded() - 1;
// delete current particle
p.Delete();
// add particles from sibyll to stack // add particles from sibyll to stack
// link to sibyll stack // link to sibyll stack
SibStack ss; SibStack ss;
...@@ -356,6 +353,8 @@ namespace corsika::process::sibyll { ...@@ -356,6 +353,8 @@ namespace corsika::process::sibyll {
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl; << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
} }
} }
// delete current particle
p.Delete();
return process::EProcessReturn::eOk; return process::EProcessReturn::eOk;
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum