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") {