diff --git a/Framework/StackInterface/ParticleBase.h b/Framework/StackInterface/ParticleBase.h
index ef85de7af7f61c829e66f04058b89d48706ddc03..d145c1e88610579c8b8ef2ddc219a261b48b695d 100644
--- a/Framework/StackInterface/ParticleBase.h
+++ b/Framework/StackInterface/ParticleBase.h
@@ -117,13 +117,14 @@ namespace corsika::stack {
using T::GetStackData;
public:
+ /*
template <typename... Args1, typename... Args2>
void SetParticleData(Args1... args1, Args2... args2) {
T::SetParticleData(args1...);
}
-
template <typename... Args1, typename... Args2>
void SetParticleData(T& p, Args1... args1, Args2... args2) {}
+ */
};
} // namespace corsika::stack
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index e917f0a721913666da4eea3df37e5c3162c9ec7d..9d259e1c282b278ab9bcec4f37a337685e29d641 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -67,11 +67,17 @@ namespace corsika::stack {
delete; ///< since Stack can be very big, we don't want to copy it
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
*/
- 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)
: fData(vD) {}
@@ -80,8 +86,11 @@ namespace corsika::stack {
* StackDataType user class. If the user did not provide a suited
* constructor this will fail with an error message.
*/
- template <typename... Args,
- typename = std::enable_if_t<!std::is_reference<StackDataType>{}>>
+ template <
+ typename... Args, typename _StackDataType = StackDataType,
+ typename = std::enable_if<std::is_same<StackDataType, _StackDataType>::value &&
+ !std::is_reference<_StackDataType>::value,
+ void>>
Stack(Args... args)
: fData(args...) {}
@@ -223,9 +232,6 @@ namespace corsika::stack {
StackIterator GetNextParticle() { return last(); }
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
* StackIterator::GetIndex() to index in data stored in
diff --git a/Framework/StackInterface/testStackInterface.cc b/Framework/StackInterface/testStackInterface.cc
index 5d4de9dae7b1cb49e6d4915eb565e2a9bf1c8416..1fa07872f3b695dda729335f2b6d52bfc69acecb 100644
--- a/Framework/StackInterface/testStackInterface.cc
+++ b/Framework/StackInterface/testStackInterface.cc
@@ -33,7 +33,6 @@ using namespace corsika::stack;
using namespace std;
typedef Stack<TestStackData, TestParticleInterface> StackTest;
-typedef StackTest::ParticleInterfaceType Particle;
TEST_CASE("Stack", "[Stack]") {
@@ -48,7 +47,6 @@ TEST_CASE("Stack", "[Stack]") {
// construct a valid Stack object
StackTest s;
- s.Init();
s.Clear();
s.AddParticle(std::tuple{0.});
s.Copy(s.cbegin(), s.begin());
@@ -98,7 +96,8 @@ TEST_CASE("Stack", "[Stack]") {
StackTest s;
REQUIRE(s.GetSize() == 0);
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);
p.AddSecondary(std::tuple{4.4});
REQUIRE(s.GetSize() == 2);
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index 159f97a7e598fe269d566fa11f53e36bab3cb24a..da7cffefe07da7cccbc6b350f80cf62d173c555a 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -32,99 +32,98 @@ using corsika::units::constants::cSquared;
double constexpr absMargin = 1e-6;
CoordinateSystem const& rootCS =
- 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());
- };
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- // helper function for mandelstam-s
- auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
- return E * E - p.squaredNorm();
- };
+// 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
+auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
+ return E * E - p.squaredNorm();
+};
TEST_CASE("rotation") {
- // define projectile kinematics in lab frame
- HEPMassType const projectileMass = 1_GeV;
- HEPMassType const targetMass = 1.0e300_eV;
- Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, 1_GeV}};
- HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
-
- Eigen::Vector3d e1, e2, e3;
- e1 << 1, 0, 0;
- e2 << 0, 1, 0;
- e3 << 0, 0, 1;
+ // define projectile kinematics in lab frame
+ HEPMassType const projectileMass = 1_GeV;
+ HEPMassType const targetMass = 1.0e300_eV;
+ Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, 1_GeV}};
+ HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+ Eigen::Vector3d e1, e2, e3;
+ e1 << 1, 0, 0;
+ e2 << 0, 1, 0;
+ 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("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));
- }
-
- SECTION("y-axis in upper 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));
- }
-
- SECTION("x-axis in upper half") {
- COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, 1_meV}}}, targetMass);
- auto const& rot = boost.GetRotationMatrix();
-
- CHECK((rot * e1 - 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));
- }
-
- SECTION("neg. 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));
- }
-
- SECTION("x-axis lower half") {
- COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, -1_meV}}}, targetMass);
- auto const& rot = boost.GetRotationMatrix();
-
- CHECK((rot * e1 - 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));
- }
-
- 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));
- }
+ SECTION("y-axis in upper 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));
+ }
+
+ SECTION("x-axis in upper half") {
+ COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, 1_meV}}}, targetMass);
+ auto const& rot = boost.GetRotationMatrix();
+
+ CHECK((rot * e1 - 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));
+ }
+
+ SECTION("neg. 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));
+ }
+
+ SECTION("x-axis lower half") {
+ COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, -1_meV}}}, targetMass);
+ auto const& rot = boost.GetRotationMatrix();
+
+ CHECK((rot * e1 - 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));
+ }
+
+ 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") {
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 1e0d1ca2b64be8db9444c94894fc6983a01e62ba..eec1f27f2239f6ba01bc4b77588313a62c6d518a 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -157,8 +157,6 @@ namespace corsika::process::sibyll {
// remember position
Point const decayPoint = p.GetPosition();
TimeType const t0 = p.GetTime();
- // remove original particle from corsika stack
- p.Delete();
// set all particles/hadrons unstable
// setHadronsUnstable();
setUnstable(pCode);
@@ -184,6 +182,8 @@ namespace corsika::process::sibyll {
}
// empty sibyll stack
ss.Clear();
+ // remove original particle from corsika stack
+ p.Delete();
}
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 5ccc1b743ce35f1e72a1ebd1e5cccba95a8d25e9..4c5f6044828cbeb651734158b7a58c81c4aaf477 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -320,9 +320,6 @@ namespace corsika::process::sibyll {
sib_list_(print_unit);
fNucCount += get_nwounded() - 1;
- // delete current particle
- p.Delete();
-
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
@@ -356,6 +353,8 @@ namespace corsika::process::sibyll {
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
}
+ // delete current particle
+ p.Delete();
return process::EProcessReturn::eOk;
}