thining: first try

corsika/detail/modules/thinning/EMThinning.inl

@@ -13,17 +13,62 @@
 
 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) {
+  template <typename TStackView>
+  void EMThinning::doSecondaries(TStackView& view) {
     if (view.size() != 2) return;
+
+    auto projectile = view.getProjecile();
+    if (!is_em(projectile.getCode())
+        return;
+        
     auto particle1 = view.begin();
     auto particle2 = std::next(particle1);
 
-    if (.getCode())
-	corsika::HEPEnergyType const E0 = 
+    if (!is_em(particle1.getCode()) || !is_em(particle2.getCode()))
+        return;
+
+	corsika::HEPEnergyType const Esum = particle1.getEnergy() + particle2.getEnergy();
+    
+    double const p1 = particle1.getEnergy() / Esum;
+    double const p2 = particle2.getEnergy() / Esum;
+    
+    double const parentWeight = projectile.getWeight();
+    
+    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_(d_) * w1prime <= parentWeight) { // keep 1st
+        particle1.setWeight(w1prime);
+      } else {
+        particle1.setWeight(0);
+      }
+      if (uniform_(d_) * w2prime <= parentWeight) { // keep 2nd
+        particle2.setWeight(w2prime);
+      } else {
+        particle2.setWeight(0);
+      }
+    }
     
-}
+    // TODO: skip this for multithinning
+    for (auto p: view) {
+      if (p.getWeight() == 0) p.erase();
+    }
+  }
 
-}
+} // namespace corsika

corsika/modules/thinning/EMThinning.hpp

@@ -9,35 +9,38 @@
 #pragma once
 
 #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 uniform_{};
+    HEPEnergyType const threshold_;
+    double const maxWeight_;
+  };
+} // namespace corsika
 
 #include <corsika/detail/process/thinning/EMThinning.inl>