From 21a72942500e33121cfa55b31793af48282e9531 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Mon, 5 Dec 2022 18:45:10 +0100
Subject: [PATCH] implemented EMThinning
---
.../detail/modules/thinning/EMThinning.inl | 72 ++++++++++++++++---
corsika/modules/thinning/EMThinning.hpp | 41 ++++++-----
2 files changed, 85 insertions(+), 28 deletions(-)
diff --git a/corsika/detail/modules/thinning/EMThinning.inl b/corsika/detail/modules/thinning/EMThinning.inl
index 9eea31d14..c0198da2c 100644
--- a/corsika/detail/modules/thinning/EMThinning.inl
+++ b/corsika/detail/modules/thinning/EMThinning.inl
@@ -13,17 +13,69 @@
namespace corsika {
-EMThinning::EMThinning(HEPEnergyType threshold, double maxWeight) : threshold_{threshold}, maxWeight_{maxWeight} {}
+ EMThinning::EMThinning(HEPEnergyType threshold, double maxWeight)
+ : threshold_{threshold}
+ , maxWeight_{maxWeight} {}
+
+ 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; }
-template <typename TStackView>
- void EMThinning::doSecondaries(TStackView& view) {
- if (view.size() != 2) return;
auto particle1 = view.begin();
- auto particle2 = std::next(particle1);
+ 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;
+
+ double const w1 = parentWeight / p1;
+ double const w2 = parentWeight / p2;
+
+ if (w1 <= maxWeight_ && w2 <= maxWeight_) { // apply Hillas thinning
+ if (uniform_(rng_) <= p1) { // keep 1st with probability p1
+ particle2.setWeight(0);
+ particle1.setWeight(w1);
+ } else { // keep 2nd
+ particle1.setWeight(0);
+ particle2.setWeight(w2);
+ }
+ } else { // weight limitation kicks in, do statistical thinning
+ double const w1prime = std::min(w1, maxWeight_);
+ double const w2prime = std::min(w2, maxWeight_);
+
+ if (uniform_(rng_) * w1prime <= parentWeight) { // keep 1st
+ particle1.setWeight(w1prime);
+ } else {
+ particle1.setWeight(0);
+ }
+ if (uniform_(rng_) * w2prime <= parentWeight) { // keep 2nd
+ particle2.setWeight(w2prime);
+ } else {
+ particle2.setWeight(0);
+ }
+ }
- if (.getCode())
- corsika::HEPEnergyType const E0 =
-
-}
+ // TODO: skip this for multithinning
+ for (auto& p : view) {
+ if (auto const w = p.getWeight(); w == 0) { p.erase(); }
+ }
+ }
-}
+} // namespace corsika
diff --git a/corsika/modules/thinning/EMThinning.hpp b/corsika/modules/thinning/EMThinning.hpp
index 017d0b770..b414be492 100644
--- a/corsika/modules/thinning/EMThinning.hpp
+++ b/corsika/modules/thinning/EMThinning.hpp
@@ -8,36 +8,41 @@
#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);
+ //! 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);
/**
- * Apply thinning to secondaries. Only EM primaries with two EM secondaries are considered.
+ * Apply thinning to secondaries. Only EM primaries with two EM secondaries are
+ * considered.
*
* @tparam TStackView
*/
template <typename TStackView>
void doSecondaries(TStackView&);
-private:
- HEPEnergyType const threshold_;
- double const maxWeight_;
-};
-}
+ private:
+ default_prng_type& rng_ = RNGManager<>::getInstance().getRandomStream("thinning");
+ std::uniform_real_distribution<double> uniform_{};
+ HEPEnergyType const threshold_;
+ double const maxWeight_;
+ };
+} // namespace corsika
-#include <corsika/detail/process/thinning/EMThinning.inl>
+#include <corsika/detail/modules/thinning/EMThinning.inl>
--
GitLab