Merge branch '385-add-optional-weight-for-particles-on-stack-thinning' into 'master'

Resolve "Add optional "weight" for particles on stack (thinning)"

Closes #385

See merge request !319

corsika/detail/setup/SetupStack.hpp

@@ -11,6 +11,7 @@
 #include <corsika/framework/stack/CombinedStack.hpp>
 #include <corsika/stack/GeometryNodeStackExtension.hpp>
 #include <corsika/stack/NuclearStackExtension.hpp>
+#include <corsika/stack/WeightStackExtension.hpp>
 #include <corsika/stack/history/HistorySecondaryProducer.hpp>
 #include <corsika/stack/history/HistoryStackExtension.hpp>
 

corsika/detail/stack/WeightStackExtension.inl

+/*
+ * (c) Copyright 2018 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/Logging.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <tuple>
+#include <utility>
+#include <vector>
+
+namespace corsika::weights {
+
+  // default version for particle-creation from input data
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setParticleData(
+      std::tuple<double> const v) {
+    setWeight(std::get<0>(v));
+  }
+
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setParticleData(
+      WeightDataInterface const& parent, std::tuple<double> const) {
+    setWeight(parent.getWeight()); // copy Weight from parent particle!
+  }
+
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setParticleData() {
+    setWeight(1);
+  } // default weight
+
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setParticleData(
+      WeightDataInterface const& parent) {
+    setWeight(parent.getWeight()); // copy Weight from parent particle!
+  }
+
+  template <typename TParentStack>
+  inline std::string WeightDataInterface<TParentStack>::asString() const {
+    return fmt::format("weight={}", getWeight());
+  }
+
+  template <typename TParentStack>
+  inline void WeightDataInterface<TParentStack>::setWeight(double const v) {
+    super_type::getStackData().setWeight(super_type::getIndex(), v);
+  }
+
+  template <typename TParentStack>
+  inline double WeightDataInterface<TParentStack>::getWeight() const {
+    return super_type::getStackData().getWeight(super_type::getIndex());
+  }
+
+  // definition of stack-data object to store geometry information
+
+  // these functions are needed for the Stack interface
+  inline void WeightData::clear() { weight_vector_.clear(); }
+
+  inline unsigned int WeightData::getSize() const { return weight_vector_.size(); }
+
+  inline unsigned int WeightData::getCapacity() const { return weight_vector_.size(); }
+
+  inline void WeightData::copy(int const i1, int const i2) {
+    weight_vector_[i2] = weight_vector_[i1];
+  }
+
+  inline void WeightData::swap(int const i1, int const i2) {
+    std::swap(weight_vector_[i1], weight_vector_[i2]);
+  }
+
+  // custom data access function
+  inline void WeightData::setWeight(int const i, double const v) {
+    weight_vector_[i] = v;
+  }
+
+  inline double WeightData::getWeight(int const i) const { return weight_vector_[i]; }
+
+  // these functions are also needed by the Stack interface
+  inline void WeightData::incrementSize() {
+    weight_vector_.push_back(1);
+  } // default weight
+
+  inline void WeightData::decrementSize() {
+    if (weight_vector_.size() > 0) { weight_vector_.pop_back(); }
+  }
+
+} // namespace corsika::weights

corsika/framework/process/NullModel.hpp

@@ -18,7 +18,7 @@ namespace corsika {
      @{
 
      Process that does nothing. It is not even derived from
-     BaseProcess
+     BaseProcess. But it can be added to a ProcessSequence.
    */
 
   class NullModel {
@@ -30,7 +30,7 @@ namespace corsika {
     static bool const is_process_sequence = false;
     static bool const is_switch_process_sequence = false;
 
-    //! Default number of processes is just one, obviously
+    //! Default number of processes is zero, obviously
     static unsigned int constexpr getNumberOfProcesses() { return 0; }
   };
 

corsika/framework/utility/COMBoost.hpp

@@ -25,7 +25,7 @@ namespace corsika {
    **/
 
   /**
-     @class
+     @class COMBoost
      @ingroup Utilities
 
      This utility class handles Lorentz boost between different

corsika/framework/utility/SaveBoostHistogram.hpp

@@ -13,6 +13,7 @@
 namespace corsika {
 
   /**
+   * @file SaveBoostHistogram.hpp
    * @ingroup Utilities
    *
    * This functions saves a boost::histogram into a numpy file. Only rather basic axis

corsika/stack/WeightStackExtension.hpp

+/*
+ * (c) Copyright 2018 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/Logging.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <tuple>
+#include <utility>
+#include <vector>
+
+namespace corsika::weights {
+
+  /**
+   * Describe "particle weights" on a Stack.
+   *
+   * Corresponding defintion of a stack-readout object, the iteractor
+   * dereference operator will deliver access to these function
+   * defintion of a stack-readout object, the iteractor dereference
+   * operator will deliver access to these function
+   */
+
+  /**
+   * \tparam TParentStack The stack to be exteneded here with weight data
+   */
+  template <typename TParentStack>
+  struct WeightDataInterface : public TParentStack {
+
+    typedef TParentStack super_type;
+
+  public:
+    // default version for particle-creation from input data
+    void setParticleData(std::tuple<double> const v);
+    void setParticleData(WeightDataInterface const& parent, std::tuple<double> const);
+    void setParticleData();
+    void setParticleData(WeightDataInterface const& parent);
+
+    std::string asString() const;
+
+    void setWeight(double const v);
+
+    double getWeight() const;
+  };
+
+  // definition of stack-data object to store geometry information
+
+  /**
+   * @class WeightData
+   *
+   * definition of stack-data object to store geometry information
+   */
+  class WeightData {
+
+  public:
+    typedef std::vector<double> weight_vector_type;
+
+    WeightData() = default;
+    WeightData(WeightData const&) = default;
+    WeightData(WeightData&&) = default;
+    WeightData& operator=(WeightData const&) = default;
+    WeightData& operator=(WeightData&&) = default;
+
+    // these functions are needed for the Stack interface
+    void clear();
+
+    unsigned int getSize() const;
+
+    unsigned int getCapacity() const;
+
+    void copy(int const i1, int const i2);
+
+    void swap(int const i1, int const i2);
+
+    // custom data access function
+    void setWeight(int const i, double const v);
+
+    double getWeight(int const i) const;
+
+    // these functions are also needed by the Stack interface
+    void incrementSize();
+
+    void decrementSize();
+
+    // custom private data section
+  private:
+    weight_vector_type weight_vector_;
+  };
+
+  template <typename TParentStack>
+  struct MakeWeightDataInterface {
+    typedef WeightDataInterface<TParentStack> type;
+  };
+
+} // namespace corsika::weights
+
+#include <corsika/detail/stack/WeightStackExtension.inl>

tests/stack/CMakeLists.txt

@@ -3,6 +3,7 @@ set (test_stack_sources
   testHistoryStack.cpp
   testHistoryView.cpp
   testGeometryNodeStackExtension.cpp
+  testWeightStackExtension.cpp
   testDummyStack.cpp
   testVectorStack.cpp
   testNuclearStackExtension.cpp

tests/stack/testWeightStackExtension.cpp

+/*
+ * (c) Copyright 2018 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/stack/CombinedStack.hpp>
+#include <corsika/stack/DummyStack.hpp>
+#include <corsika/stack/WeightStackExtension.hpp>
+
+using namespace corsika;
+
+#include <catch2/catch.hpp>
+
+#include <iostream>
+using namespace std;
+
+// the Weight stack:
+template <typename TStackIter>
+using DummyWeightDataInterface =
+    typename weights::MakeWeightDataInterface<TStackIter>::type;
+
+// combine dummy stack with geometry information for tracking
+template <typename TStackIter>
+using StackWithGeometryInterface =
+    CombinedParticleInterface<dummy_stack::DummyStack::pi_type, DummyWeightDataInterface,
+                              TStackIter>;
+
+using TestStack =
+    CombinedStack<typename dummy_stack::DummyStack::stack_implementation_type,
+                  weights::WeightData, StackWithGeometryInterface>;
+
+TEST_CASE("WeightStackExtension", "[stack]") {
+
+  logging::set_level(logging::level::info);
+  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+
+  dummy_stack::NoData noData;
+
+  SECTION("write weights") {
+
+    double const weight = 5.1;
+
+    TestStack s;
+    s.addParticle(std::make_tuple(noData), std::tuple<double>{weight});
+
+    CHECK(s.getEntries() == 1);
+  }
+
+  SECTION("write/read weights") {
+    double const weight = 15;
+
+    TestStack s;
+    auto p = s.addParticle(std::make_tuple(noData));
+    p.setWeight(weight);
+    CHECK(s.getEntries() == 1);
+
+    const auto pout = s.getNextParticle();
+    CHECK(pout.getWeight() == weight);
+  }
+
+  SECTION("stack fill and cleanup") {
+
+    double const weight = 16;
+
+    TestStack s;
+    // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
+    for (int i = 0; i < 99; ++i) {
+      auto p = s.addParticle(std::tuple<dummy_stack::NoData>{noData});
+      p.setWeight(weight);
+    }
+
+    CHECK(s.getEntries() == 99);
+    double v = 0;
+    for (int i = 0; i < 99; ++i) {
+      auto p = s.getNextParticle();
+      v += p.getWeight();
+      p.erase();
+    }
+    CHECK(v == Approx(99 * weight));
+    CHECK(s.getEntries() == 0);
+  }
+}