diff --git a/corsika/detail/modules/thinning/EMThinning.inl b/corsika/detail/modules/thinning/EMThinning.inl
new file mode 100644
index 0000000000000000000000000000000000000000..29ed5ad46d2268700b505fde6f386782cf06e8bc
--- /dev/null
+++ b/corsika/detail/modules/thinning/EMThinning.inl
@@ -0,0 +1,90 @@
+/*
+ * (c) Copyright 2022 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/PhysicalUnits.hpp>
+#include <corsika/framework/process/SecondariesProcess.hpp>
+
+namespace corsika {
+
+ EMThinning::EMThinning(HEPEnergyType threshold, double maxWeight, bool const eraseParticles)
+ : threshold_{threshold}
+ , maxWeight_{maxWeight}
+ , eraseParticles_ {eraseParticles} {}
+
+ template <typename TStackView>
+ void EMThinning::doSecondaries(TStackView& view) {
+ if (view.getSize() != 2) return;
+
+ auto projectile = view.getProjectile();
+ if (!is_em(projectile.getPID())) return;
+
+ double const parentWeight = projectile.getWeight();
+ if (parentWeight >= maxWeight_) { return; }
+
+ auto particle1 = view.begin();
+ auto particle2 = ++(view.begin());
+
+ // skip non-EM interactions
+ if (!is_em(particle1.getPID()) || !is_em(particle2.getPID())) { return; }
+
+ /* skip high-energy interactions
+ * We consider only the primary energy. A more sophisticated algorithm might
+ * sort out above-threshold secondaries and apply thinning only on the subset of
+ * those below the threshold.
+ */
+ if (projectile.getEnergy() > threshold_) { return; }
+
+ corsika::HEPEnergyType const Esum = particle1.getEnergy() + particle2.getEnergy();
+
+ double const p1 = particle1.getEnergy() / Esum;
+ double const p2 = particle2.getEnergy() / Esum;
+
+ // weight factors
+ double const w1 = parentWeight / p1;
+ double const w2 = parentWeight / p2;
+
+ // max. allowed weight factor
+ double const maxWeightFactor = maxWeight_ / parentWeight;
+
+ if (w1 <= maxWeightFactor && w2 <= maxWeightFactor) { // apply Hillas thinning
+ if (uniform_(rng_) <= p1) { // keep 1st with probability p1
+ particle2.setWeight(0);
+ particle1.setWeight(w1 * parentWeight);
+ } else { // keep 2nd
+ particle1.setWeight(0);
+ particle2.setWeight(w2 * parentWeight);
+ }
+ } else { // weight limitation kicks in, do statistical thinning
+ double const w1prime = std::min(w1, maxWeightFactor);
+ double const w2prime = std::min(w2, maxWeightFactor);
+
+ if (uniform_(rng_) * w1prime <= parentWeight) { // keep 1st
+ particle1.setWeight(w1prime * parentWeight);
+ } else {
+ particle1.setWeight(0);
+ }
+ if (uniform_(rng_) * w2prime <= parentWeight) { // keep 2nd
+ particle2.setWeight(w2prime * parentWeight);
+ } else {
+ particle2.setWeight(0);
+ }
+ }
+
+ // erase discared particles in case of multithinning
+ if (eraseParticles_) {
+ for (auto& p : view) {
+ if (auto const w = p.getWeight(); w == 0) { p.erase(); }
+ }
+ } else {
+ return;
+ }
+ }
+
+} // namespace corsika
diff --git a/corsika/modules/thinning/EMThinning.hpp b/corsika/modules/thinning/EMThinning.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a32c05a6345d533217c0b1653b8ecfb6886d5308
--- /dev/null
+++ b/corsika/modules/thinning/EMThinning.hpp
@@ -0,0 +1,57 @@
+/*
+ * (c) Copyright 2022 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 <random>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/process/SecondariesProcess.hpp>
+
+namespace corsika {
+
+ /**
+ * This process implements thinning for EM splitting processes (1 -> 2).
+ */
+
+ class EMThinning : public SecondariesProcess<EMThinning> {
+ public:
+ /**
+ * Construct a new EMThinning process.
+ *
+ * @param threshold: thinning applied below this energy
+ * @param maxWeight: maximum allowed weight
+ */
+ EMThinning(HEPEnergyType threshold, double maxWeight, bool const eraseParticles=true);
+
+ /**
+ * Apply thinning to secondaries. Only EM primaries with two EM secondaries are
+ * considered.
+ *
+ * If the maximum weight is still out of reach, Hillas thinning is applied (i.e.
+ * one of the two secondaries is kept, the other one discarded). If the acceptance
+ * probabilities would lead to a weight factor exceeding the maximum weight, we resort
+ * to statistical thinning (i.e. the secondaries are kept/discared randomly each on is
+ * own). In that case, acceptance probabilities can be assigned without constraints
+ * (sum does not need to be 1) and we increase the acceptance probability such that
+ * the maximum weight is not exceeded.
+ */
+ template <typename TStackView>
+ void doSecondaries(TStackView&);
+
+ private:
+ default_prng_type& rng_ = RNGManager<>::getInstance().getRandomStream("thinning");
+ std::uniform_real_distribution<double> uniform_{};
+ HEPEnergyType const threshold_;
+ double const maxWeight_;
+ bool const eraseParticles_;
+ };
+} // namespace corsika
+
+#include <corsika/detail/modules/thinning/EMThinning.inl>
diff --git a/tests/common/SetupStack.hpp b/tests/common/SetupStack.hpp
index e2bd85f2aeaa903f65854357d5a51834f43f0c57..4a7ceff3ff4590f1bd20a8736aec0ddc4581d704 100644
--- a/tests/common/SetupStack.hpp
+++ b/tests/common/SetupStack.hpp
@@ -31,11 +31,11 @@ namespace corsika::test {
/*
* the version without history
*/
- using Stack = detail::StackWithGeometry;
+ using Stack = detail::StackWithWeight;
#endif
// the correct secondary stack view
using StackView = typename Stack::stack_view_type;
-} // namespace corsika::test
\ No newline at end of file
+} // namespace corsika::test
diff --git a/tests/common/TestStack.hpp b/tests/common/TestStack.hpp
index 2488fc1f059cfa4b4df6ed27a931feff60a62606..934575cb1deef376109b8dcd7c55d1c06e4b1533 100644
--- a/tests/common/TestStack.hpp
+++ b/tests/common/TestStack.hpp
@@ -48,20 +48,33 @@ namespace corsika {
node::GeometryData<DummyEnvironment>, StackWithGeometryInterface,
DefaultSecondaryProducer>;
+ template <typename TStackIter>
+ using SetupWeightDataInterface =
+ typename weights::MakeWeightDataInterface<TStackIter>::type;
+
+ template <typename TStackIter>
+ using StackWithWeightInterface =
+ CombinedParticleInterface<StackWithGeometry::pi_type, SetupWeightDataInterface,
+ TStackIter>;
+
+ using StackWithWeight =
+ CombinedStack<typename StackWithGeometry::stack_data_type, weights::WeightData,
+ StackWithWeightInterface, DefaultSecondaryProducer>;
+
// ------------------------------------------
// Add [optional] history data to stack, too:
// combine dummy stack with geometry information for tracking
template <typename TStackIter>
using StackWithHistoryInterface =
- CombinedParticleInterface<StackWithGeometry::pi_type,
+ CombinedParticleInterface<StackWithWeight::pi_type,
history::HistoryEventDataInterface, TStackIter>;
using StackWithHistory =
- CombinedStack<typename StackWithGeometry::stack_data_type,
+ CombinedStack<typename StackWithWeight::stack_data_type,
history::HistoryEventData, StackWithHistoryInterface,
history::HistorySecondaryProducer>;
} // namespace test::detail
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 3bb6cf8db10c42d58b4862d01dde1e297b07ea2b..015093588f7ec42d0496ba5bb0d323badc425f4d 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -13,6 +13,7 @@ set (test_modules_sources
testSibyll.cpp
testEpos.cpp
testRadio.cpp
+ testEMThinning.cpp
testSophia.cpp
)
diff --git a/tests/modules/testEMThinning.cpp b/tests/modules/testEMThinning.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b23d0c129d94fcce3c63873b7af8254ffe6eae5
--- /dev/null
+++ b/tests/modules/testEMThinning.cpp
@@ -0,0 +1,223 @@
+/*
+ * (c) Copyright 2022 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.
+ */
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/media/Environment.hpp>
+#include <corsika/modules/thinning/EMThinning.hpp>
+
+#include <SetupStack.hpp>
+#include <SetupTestTrajectory.hpp>
+#include <SetupTestEnvironment.hpp>
+
+#include <catch2/catch.hpp>
+
+#include <boost/accumulators/accumulators.hpp>
+#include <boost/accumulators/statistics/stats.hpp>
+#include <boost/accumulators/statistics/mean.hpp>
+
+using namespace corsika;
+
+using DummyEnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
+using DummyEnvironment = Environment<DummyEnvironmentInterface>;
+
+TEST_CASE("EMThinning", "process,secondary") {
+ RNGManager<>::getInstance().registerRandomStream("thinning");
+ logging::set_level(logging::level::info);
+
+ feenableexcept(FE_INVALID);
+ using EnvType = DummyEnvironment;
+
+ EnvType env;
+ CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+
+ // setup empty particle stack
+ test::Stack stack;
+ stack.clear();
+
+ SECTION("non-EM primary") {
+ auto prim = stack.addParticle(std::make_tuple(Code::Proton, 10_GeV,
+ DirectionVector(rootCS, {1, 0, 0}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+ test::StackView view{prim};
+ auto projectile = view.getProjectile();
+
+ auto sec1 = projectile.addSecondary(
+ std::make_tuple(Code::PiPlus, 2_GeV, DirectionVector(rootCS, {1, 0, 0})));
+ auto sec2 = projectile.addSecondary(
+ std::make_tuple(Code::KMinus, 8_GeV, DirectionVector(rootCS, {1, 0, 0})));
+
+ EMThinning thinning{100_GeV, 1e20};
+ thinning.doSecondaries(view);
+ REQUIRE(view.getEntries() == 2);
+ REQUIRE(sec1.getWeight() == 1.);
+ REQUIRE(sec2.getWeight() == 1.);
+ }
+
+ SECTION("non-EM interaction") {
+ auto prim = stack.addParticle(std::make_tuple(Code::Photon, 10_GeV,
+ DirectionVector(rootCS, {1, 0, 0}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+ test::StackView view{prim};
+ auto projectile = view.getProjectile();
+
+ // mimic photohadronic interaction, more than 2 secondaries
+ auto sec1 = projectile.addSecondary(
+ std::make_tuple(Code::PiPlus, 2_GeV, DirectionVector(rootCS, {1, 0, 0})));
+ auto sec2 = projectile.addSecondary(
+ std::make_tuple(Code::KMinus, 3_GeV, DirectionVector(rootCS, {1, 0, 0})));
+ auto sec3 = projectile.addSecondary(
+ std::make_tuple(Code::Proton, 5_GeV, DirectionVector(rootCS, {1, 0, 0})));
+
+ EMThinning thinning{100_GeV, 1e20};
+ thinning.doSecondaries(view);
+ REQUIRE(view.getEntries() == 3);
+ REQUIRE(sec1.getWeight() == 1.);
+ REQUIRE(sec2.getWeight() == 1.);
+ REQUIRE(sec3.getWeight() == 1.);
+ }
+
+ SECTION("above threshold") {
+ auto prim = stack.addParticle(std::make_tuple(Code::Photon, 10_GeV,
+ DirectionVector(rootCS, {1, 0, 0}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+
+ test::StackView view{prim};
+ auto projectile = view.getProjectile();
+
+ auto sec1 = projectile.addSecondary(
+ std::make_tuple(Code::Electron, 2_GeV, DirectionVector(rootCS, {1, 0, 0})));
+ auto sec2 = projectile.addSecondary(
+ std::make_tuple(Code::Positron, 8_GeV, DirectionVector(rootCS, {1, 0, 0})));
+
+ EMThinning thinning{1_MeV, 1e20};
+ thinning.doSecondaries(view);
+ REQUIRE(view.getEntries() == 2);
+ REQUIRE(sec1.getWeight() == 1.);
+ REQUIRE(sec2.getWeight() == 1.);
+ }
+
+ SECTION("below threshold") {
+ EMThinning thinning{100_GeV, 1e20};
+
+ boost::accumulators::accumulator_set<
+ double, boost::accumulators::stats<boost::accumulators::tag::mean>>
+ acc;
+ boost::accumulators::accumulator_set<
+ HEPEnergyType, boost::accumulators::stats<boost::accumulators::tag::mean>>
+ accE;
+
+ // statistical test, 10000 iterations
+ for (int i = 0; i < 10'000; ++i) {
+ stack.clear();
+ auto prim = stack.addParticle(std::make_tuple(Code::Photon, 10_GeV,
+ DirectionVector(rootCS, {1, 0, 0}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+
+ test::StackView view{prim};
+ auto projectile = view.getProjectile();
+
+ projectile.addSecondary(
+ std::make_tuple(Code::Electron, 2_GeV, DirectionVector(rootCS, {1, 0, 0})));
+ projectile.addSecondary(
+ std::make_tuple(Code::Positron, 8_GeV, DirectionVector(rootCS, {1, 0, 0})));
+
+ thinning.doSecondaries(view);
+ REQUIRE(view.getEntries() == 1);
+
+ double sumWeight{};
+ HEPEnergyType sumE{};
+
+ for (auto& p : view) {
+ sumWeight += p.getWeight();
+ sumE += p.getWeight() * p.getEnergy();
+ }
+
+ acc(sumWeight);
+ accE(sumE);
+ }
+
+ REQUIRE(boost::accumulators::mean(acc) == Approx(2).epsilon(2e-2));
+ REQUIRE(boost::accumulators::mean(accE) / 1_GeV == Approx(10).epsilon(2e-2));
+ }
+
+ SECTION("weight limitation") {
+ double constexpr maxWeight = 3; // min acceptance prob. >= 0.3333
+ EMThinning thinning{100_GeV, maxWeight};
+
+ boost::accumulators::accumulator_set<
+ double, boost::accumulators::stats<boost::accumulators::tag::mean>>
+ acc;
+ boost::accumulators::accumulator_set<
+ HEPEnergyType, boost::accumulators::stats<boost::accumulators::tag::mean>>
+ accE;
+
+ // statistical test, 10000 iterations
+ for (int i = 0; i < 10'000; ++i) {
+ stack.clear();
+ auto prim = stack.addParticle(std::make_tuple(Code::Photon, 10_GeV,
+ DirectionVector(rootCS, {1, 0, 0}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+
+ test::StackView view{prim};
+ auto projectile = view.getProjectile();
+
+ projectile.addSecondary(
+ std::make_tuple(Code::Electron, 2_GeV, DirectionVector(rootCS, {1, 0, 0})));
+ projectile.addSecondary(
+ std::make_tuple(Code::Positron, 8_GeV, DirectionVector(rootCS, {1, 0, 0})));
+
+ thinning.doSecondaries(view);
+
+ double sumWeight{};
+ HEPEnergyType sumE{};
+
+ for (auto& p : view) {
+ auto const w = p.getWeight();
+ REQUIRE(w <= maxWeight);
+ sumWeight += w;
+ sumE += p.getEnergy() * w;
+ }
+
+ acc(sumWeight);
+ accE(sumE);
+ }
+
+ REQUIRE(boost::accumulators::mean(acc) == Approx(2).epsilon(2e-2));
+ REQUIRE(boost::accumulators::mean(accE) / 1_GeV == Approx(10).epsilon(2e-2));
+ }
+
+ SECTION("max. weight reached") {
+ double constexpr maxWeight = 3; // min acceptance prob. >= 0.3333
+ EMThinning thinning{100_GeV, maxWeight};
+
+ for (int i = 0; i < 10'000; ++i) {
+ stack.clear();
+ auto prim = stack.addParticle(std::make_tuple(Code::Photon, 10_GeV,
+ DirectionVector(rootCS, {1, 0, 0}),
+ Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+ prim.setWeight(3.5);
+
+ test::StackView view{prim};
+ auto projectile = view.getProjectile();
+
+ auto sec1 = projectile.addSecondary(
+ std::make_tuple(Code::Electron, 2_GeV, DirectionVector(rootCS, {1, 0, 0})));
+ auto sec2 = projectile.addSecondary(
+ std::make_tuple(Code::Positron, 8_GeV, DirectionVector(rootCS, {1, 0, 0})));
+
+ // cannot perform thinning anymore, always keep all secondaries, weight unchanged
+ thinning.doSecondaries(view);
+ REQUIRE(view.getEntries() == 2);
+ REQUIRE(sec1.getWeight() == 3.5);
+ REQUIRE(sec2.getWeight() == 3.5);
+ }
+ }
+}