Merge remote-tracking branch 'origin/64-random-number-access-in-different-modules' into sibyll

CMakeLists.txt

@@ -15,7 +15,7 @@ set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/CMakeModules)
 include (CorsikaUtilities) # a few cmake function
 
 # enable warnings and disallow non-standard language
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -pedantic -Wextra")
+add_definitions(-Wall -pedantic -Wextra)
 
 set (CMAKE_CXX_STANDARD 17)
 enable_testing ()
@@ -33,16 +33,16 @@ ENABLE_LANGUAGE (Fortran)
 #  #-j ${PROCESSOR_COUNT}
 #  #     DEPENDENCIES corsika
 #  )
-   
+
 #add_custom_target (corsika_pre_build)
 #add_custom_command (TARGET corsika_pre_build PRE_BUILD COMMAND "${PROJECT_SOURCE_DIR}/pre_compile.py")
 
 # dependencies
-#find_package (Boost 1.40 COMPONENTS program_options REQUIRED)
+find_package (Boost 1.60 REQUIRED)
 find_package (Eigen3 REQUIRED)
 #find_package (HDF5) # not yet needed
 
-# order of subdirectories 
+# order of subdirectories
 add_subdirectory (ThirdParty)
 #add_subdirectory (Utilities)
 add_subdirectory (Framework)

Documentation/Doxygen/CMakeLists.txt

 
-find_package (Doxygen REQUIRED dot OPTIONAL_COMPONENTS mscgen dia)
+find_package (Doxygen OPTIONAL_COMPONENTS dot mscgen dia)
 
 if (DOXYGEN_FOUND)
+  if (NOT DOXYGEN_DOT_EXECUTABLE)
+    message (FATAL_ERROR "Found doxygen but not 'dot' command, please install graphviz or set DOXYGEN_DOT_EXECUTABLE")
+  endif()
+
   set (DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
   set (DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile)
   set (DOXYGEN_GENERATE_HTML YES)
@@ -28,4 +32,3 @@ else (DOXYGEN_FOUND)
   message ("Doxygen need to be installed to generate the doxygen documentation")
   
 endif (DOXYGEN_FOUND)
-

Framework/CMakeLists.txt

 
+add_subdirectory (Utilities)
 add_subdirectory (Units)
 add_subdirectory (Geometry)
 add_subdirectory (Particles)

Framework/ProcessSequence/ProcessSequence.h

@@ -189,7 +189,7 @@ namespace corsika::process {
   };
 
   /// the +operator assembles many BaseProcess, ContinuousProcess, and
-  /// DiscreteProcess objects into a ProcessSequence, all combinatoris
+  /// DiscreteProcess objects into a ProcessSequence, all combinatorics
   /// must be allowed, this is why we define a macro to define all
   /// combinations here:
 

Framework/Random/RNGManager.h

@@ -12,6 +12,8 @@
 #ifndef _include_RNGManager_h_
 #define _include_RNGManager_h_
 
+#include <corsika/utl/Singleton.h>
+
 #include <map>
 #include <random>
 #include <sstream>
@@ -25,10 +27,16 @@ namespace corsika::random {
 
   using RNG = std::mt19937; //!< the actual RNG type that will be used
 
-  class RNGManager {
+  class RNGManager : public corsika::utl::Singleton<RNGManager> {
+
+    friend class corsika::utl::Singleton<RNGManager>;
+    
     std::map<std::string, RNG> rngs;
-    std::map<std::string, std::seed_seq> seeds;
+    std::map<std::string, std::seed_seq> seeds;    
 
+  protected:
+    RNGManager() {}
+    
   public:
     /*!
      * This function is to be called by a module requiring a random-number

Framework/Random/testRandom.cc

@@ -20,7 +20,7 @@ using namespace corsika::random;
 
 SCENARIO("random-number streams can be registered and retrieved") {
   GIVEN("a RNGManager") {
-    RNGManager rngManager;
+    RNGManager& rngManager = RNGManager::GetInstance();
 
     WHEN("a sequence is registered by name") {
       rngManager.RegisterRandomStream("stream_A");

Framework/Utilities/Bit.h

+#ifndef _corsika_utl_Bit_h_
+#define _corsika_utl_Bit_h_
+
+/**
+  \author Hans Dembinski
+  \author Lukas Nellen
+  \author Darko Veberic
+  \date 27 Jan 2014
+
+  \version $Id: Bit.h 25126 2014-02-03 22:13:10Z darko $
+*/
+
+#include <exception>
+
+// #include <utl/AugerException.h>
+
+
+namespace corsika::utl {
+
+  namespace Bit {
+
+    template<typename T>
+    class Array {
+    public:
+      Array(T& target) : fTarget(target) { }
+
+      class Bit {
+      public:
+        Bit(T& target, T mask) : fTarget(target), fMask(mask) { }
+
+        operator bool() const { return fTarget & fMask; }
+
+        bool operator~() const { return !bool(*this); }
+
+        Bit&
+        operator=(const bool value)
+        {
+          if (value)
+            fTarget |= fMask;
+          else
+            fTarget &= ~fMask;
+          return *this;
+        }
+
+        Bit& Flip() { return *this = ~(*this); }
+
+      private:
+        T& fTarget;
+        T fMask;
+      };
+
+      Bit operator[](unsigned int position)
+      { return Bit(fTarget, T(1) << position); }
+
+      Bit
+      At(unsigned int position)
+      {
+        if (position >= 8*sizeof(T))
+          //throw std::exceptionOutOfBoundException("Running out of bits.");
+	  throw std::exception("Running out of bits.");
+        return (*this)[position];
+      }
+
+      template<typename M>
+      Array& Mask(const M mask, const bool value)
+      { Bit(fTarget, mask) = value; return *this; }
+
+      template<typename M>
+      T Get(const M mask) { return fTarget & T(mask); }
+
+    private:
+      T& fTarget;
+    };
+
+  }
+
+  // helper
+  template<typename T>
+  inline
+  Bit::Array<T>
+  AsBitArray(T& target)
+  {
+    return Bit::Array<T>(target);
+  }
+
+}
+
+
+#endif

Framework/Utilities/CMakeLists.txt

+
+set (
+  UTILITIES_SOURCES  
+  Dummy.cc
+  )
+
+set (
+  UTILITIES_HEADERS
+  Dummy.h
+  Bit.h
+  Test.h
+  Singleton.h
+  )
+
+set (
+  UTILITIES_NAMESPACE
+  corsika/utl
+  )
+
+add_library (CORSIKAutilities STATIC ${UTILITIES_SOURCES})
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAutilities ${UTILITIES_NAMESPACE} ${UTILITIES_HEADERS})
+
+set_target_properties (
+  CORSIKAutilities
+  PROPERTIES
+  VERSION ${PROJECT_VERSION}
+  SOVERSION 1
+  PUBLIC_HEADER "${UTILITIES_HEADERS}"
+  )
+
+# target dependencies on other libraries (also the header onlys)
+#target_link_libraries (
+#  )
+
+target_include_directories (
+  CORSIKAutilities
+  PUBLIC
+  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+  $<INSTALL_INTERFACE:include/include>
+  )
+
+install (
+  TARGETS CORSIKAutilities
+  LIBRARY DESTINATION lib
+  ARCHIVE DESTINATION lib
+  PUBLIC_HEADER DESTINATION include/${UTILITIES_NAMESPACE}
+  )
+
+
+# --------------------
+# code unit testing
+# add_executable (testBit testBit.cc)
+
+# target_link_libraries (
+#   testBit
+#   CORSIKAutilities
+#   CORSIKAthirdparty # for catch2
+#   )
+
+# add_test (NAME testBit COMMAND testBit)
+

Framework/Utilities/Dummy.cc

+#include <corsika/utl/Dummy.h>
+
+using namespace corsika::utl;
+
+// big void...

Framework/Utilities/Dummy.h

+#ifndef _include_corsika_utilties_dummy_h_
+#define _include_corsika_utilties_dummy_h_
+
+namespace corsika::utl {
+
+  // void....
+  
+}
+
+#endif

Framework/Utilities/Singleton.h

+#ifndef _corsika_utl_Singleton_h_
+#define _corsika_utl_Singleton_h_
+
+//#define OFFLINE_USE_GAMMA_SINGLETON
+
+namespace corsika::utl {
+
+#ifndef OFFLINE_USE_GAMMA_SINGLETON
+
+  /**
+   * \class Singleton Singleton.h utl/Singleton.h
+   *
+   * \brief Curiously Recurring Template Pattern (CRTP) for Meyers singleton
+   *
+   * The singleton class is implemented as follows
+   * \code
+   * #include <utl/Singleton.h>
+   *
+   * class SomeClass : public utl::Singleton<SomeClass> {
+   *   ...
+   * private:
+   *   // prevent creation, destruction
+   *   SomeClass() { }
+   *   ~SomeClass() { }
+   *
+   *   friend class utl::Singleton<SomeClass>;
+   * };
+   * \endcode
+   * Singleton automatically prevents copying of the derived class.
+   *
+   * \author Darko Veberic
+   * \date 9 Aug 2006
+   * \version $Id: Singleton.h 25091 2014-01-30 09:49:57Z darko $
+   * \ingroup stl
+   */
+
+  template<typename T>
+  class Singleton {
+  public:
+    static
+    T&
+    GetInstance()
+#  ifdef __MAKECINT__
+    ;
+#  else
+    {
+      static T instance;
+      return instance;
+    }
+#  endif
+
+  protected:
+    // derived class can call ctor and dtor
+    Singleton() { }
+    ~Singleton() { }
+
+  private:
+    // no one should do copies
+    Singleton(const Singleton&);
+    Singleton& operator=(const Singleton&);
+  };
+
+#else
+
+  /// classical Gamma singleton
+  template<typename T>
+  class Singleton {
+  public:
+    static
+    T&
+    GetInstance()
+    {
+      if (!fgInstance)
+        fgInstance = new T;
+      return *fgInstance;
+    }
+
+  protected:
+    // derived class can call ctor and dtor
+    Singleton() { }
+    ~Singleton() { }
+
+  private:
+    // no one should do copies
+    Singleton(const Singleton&);
+    Singleton& operator=(const Singleton&);
+
+    static T* fgInstance = 0;
+  };
+
+#endif
+
+
+  /**
+   * \class LeakingSingleton Singleton.h utl/Singleton.h
+   *
+   * \brief CRTP for leaking singleton
+   *
+   * This type of creation (Gamma singleton) leaks the object at
+   * the end of the run, i.e. class T destructor does not get called
+   * in at_exit().
+   *
+   * This singleton can be implemented as follows
+   * \code
+   * #include <utl/Singleton.h>
+   *
+   * class SomeClass : public utl::LeakingSingleton<SomeClass> {
+   *   ...
+   * private:
+   *   // prevent creation, destruction
+   *   SomeClass() { }
+   *   ~SomeClass() { }
+   *
+   *   friend class utl::LeakingSingleton<SomeClass>;
+   * };
+   * \endcode
+   * LeakingSingleton automatically prevents copying of the derived
+   * class.
+   *
+   * \author Darko Veberic
+   * \date 9 Aug 2006
+   * \version $Id: Singleton.h 25091 2014-01-30 09:49:57Z darko $
+   * \ingroup stl
+   */
+
+  template<class T>
+  class LeakingSingleton {
+  public:
+    static
+    T&
+    GetInstance()
+    {
+      static T* const instance = new T;
+      return *instance;
+    }
+
+  protected:
+    // derived class can call ctor and dtor
+    LeakingSingleton() { }
+    ~LeakingSingleton() { }
+
+  private:
+    // no one should do copies
+    LeakingSingleton(const LeakingSingleton&);
+    LeakingSingleton& operator=(const LeakingSingleton&);
+  };
+
+}
+
+
+#endif

Framework/Utilities/Test.h

+#ifndef _utl_Test_h_
+#define _utl_Test_h_
+
+/**
+   \file
+   Tools to do simple testing in a readable way
+
+   \author Lukas Nellen
+   \author Darko Veberic
+   \version $Id: Test.h 31925 2018-09-25 16:02:12Z darko $
+   \date 08 Feb 2004
+
+   \ingroup testing
+*/
+
+#include <cmath>
+#include <boost/format.hpp>
+#include <boost/tuple/tuple.hpp>
+#include <boost/tuple/tuple_comparison.hpp>
+#include <boost/tuple/tuple_io.hpp>
+#include <utl/Triple.h>
+
+
+namespace utl {
+
+  //! Predicate used in STL for searching for whitespace
+  struct IsSpace {
+    bool operator()(const char x) const
+    { return x == ' ' || x == '\r' || x == '\n' || x == '\t'; }
+  };
+
+
+  /// Predicate for equality
+  class Equal {
+  public:
+    template<typename T>
+    bool operator()(const T& lhs, const T& rhs) const
+    { return lhs == rhs; }
+
+    static const char* Name() { return "equal"; }
+  };
+
+
+  /// Predicate for less
+  class Less {
+  public:
+    template<typename T>
+    bool operator()(const T& lhs, const T& rhs) const
+    { return lhs < rhs; }
+
+    static const char* Name() { return "less"; }
+  };
+
+
+  /// Predicate for less or equal
+  class LessOrEqual {
+  public:
+    template<typename T>
+    bool operator()(const T& lhs, const T& rhs) const
+    { return lhs <= rhs; }
+
+    static const char* Name() { return "less or equal"; }
+  };
+
+
+  /// Predicate for greater
+  class Greater {
+  public:
+    template<typename T>
+    bool operator()(const T& lhs, const T& rhs) const
+    { return lhs > rhs; }
+
+    static const char* Name() { return "greater"; }
+  };
+
+
+  /// Predicate for greater or equal
+  class GreaterOrEqual {
+  public:
+    template<typename T>
+    bool operator()(const T& lhs, const T& rhs) const
+    { return lhs >= rhs; }
+
+    static const char* Name() { return "greater or equal"; }
+  };
+
+  /// Predicate for approximate equality (for floating point)
+  /** The default precision is 1e-6, but it can be changed at
+      construction time.
+   */
+  class CloseTo {
+  public:
+    CloseTo(const double eps = 1e-6) : fEpsilon(eps) { }
+
+    template<typename T>
+    bool operator()(const T& lhs, const T& rhs) const
+    { return IsCloseTo(lhs, rhs); }
+
+    boost::format Name() const
+    { return boost::format("close (@%g) to") % fEpsilon; }
+
+  protected:
+    template<typename T>
+    bool IsCloseAbs(const T& lhs, const T& rhs) const
+    { return std::abs(double(lhs) - double(rhs)) < fEpsilon; }
+
+    bool IsCloseAbs(const utl::Triple& lhs, const utl::Triple& rhs) const
+    { return std::sqrt(TupleDist2(lhs, rhs)) < fEpsilon; }
+
+    template<typename T>
+    bool IsCloseRel(const T& lhs, const T& rhs) const
+    { return 2*std::abs(double(lhs) - double(rhs)) / (std::abs(double(lhs)) + std::abs(double(rhs))) < fEpsilon; }
+
+    bool
+    IsCloseRel(const utl::Triple& lhs, const utl::Triple& rhs)
+      const
+    {
+      return (2*sqrt(TupleDist2(lhs, rhs)) /
+        (sqrt(TupleDist2(lhs, utl::Triple(0,0,0))) +
+         sqrt(TupleDist2(rhs, utl::Triple(0,0,0))))) < fEpsilon;
+    }
+
+    template<typename T>
+    bool
+    IsCloseTo(const T& lhs, const T& rhs)
+      const
+    {
+      if (IsCloseAbs(lhs, rhs))
+        return true;
+      else
+        return IsCloseRel(lhs, rhs);
+    }
+
+    // tuple distance
+    template<typename Head, typename Tail>
+    static
+    double
+    TupleDist2(const boost::tuples::cons<Head, Tail>& lhs,
+               const boost::tuples::cons<Head, Tail>& rhs)
+    {
+      const double t = lhs.get_head() - rhs.get_head();
+      return t*t + TupleDist2(lhs.get_tail(), rhs.get_tail());
+    }
+
+    static double TupleDist2(const boost::tuples::null_type& /*lhs*/,
+                             const boost::tuples::null_type& /*rhs*/)
+    { return 0; }
+
+    double fEpsilon;
+  };
+
+
+  class CloseAbs : public CloseTo {
+  public:
+    CloseAbs(const double eps = 1e-6) : CloseTo(eps) { }
+
+    template<typename T>
+    bool operator()(const T& lhs, const T& rhs) const
+    { return IsCloseAbs(lhs, rhs); }
+
+    boost::format Name() const
+    { return boost::format("absolutely close (@%g) to") % fEpsilon; }
+  };
+
+
+  class CloseRel : public CloseTo {
+  public:
+    CloseRel(const double eps = 1e-6) : CloseTo(eps) { }
+
+    template<typename T>
+    bool operator()(const T& lhs, const T& rhs) const
+    { return IsCloseRel(lhs, rhs); }
+
+    boost::format Name() const
+    { return boost::format("relatively close (@%g) to") % fEpsilon; }
+  };
+
+
+  template<typename Predicate>
+  class Not : public Predicate {
+  public:
+    Not() : Predicate() { }
+    
+    Not(const double eps) : Predicate(eps) { }
+    
+    template<typename T>
+    bool operator()(const T& x) const
+    { return !Predicate::operator()(x); }
+
+    template<typename T, typename U>
+    bool operator()(const T& x, const U& y) const
+    { return !Predicate::operator()(x, y); }
+
+    template<typename T, typename U, typename W>
+    bool operator()(const T& x, const U& y, const W& z) const
+    { return !Predicate::operator()(x, y, z); }
+
+    static boost::format Name()
+    { return boost::format("not-%s") % Predicate().Name(); }
+  };
+
+
+  inline
+  utl::Triple
+  Diff(const utl::Triple& lhs, const utl::Triple& rhs)
+  {
+    return utl::Triple(lhs.get<0>() - rhs.get<0>(),
+                       lhs.get<1>() - rhs.get<1>(),
+                       lhs.get<2>() - rhs.get<2>());
+  }
+
+
+  /// Test condition by evaluating a predicate
+  /** If the predicate evaluates to false, we print a failure message
+      with the values of the left- and right-hand side and the name of
+      the predicate. This information is normally not available when
+      using the CPPUNIT_ASSERT macro.
+   */
+  template<class Predicate, typename T>
+  inline bool Test(const Predicate& pred, const T& lhs, const T& rhs)
+  { return pred(lhs, rhs); }
+
+
+  /// Main test function
+  template<class Predicate, typename T>
+  inline bool Test(const T& lhs, const T& rhs)
+  { return Test(Predicate(), lhs, rhs); }
+
+
+  /// Test function for predicates that take an option
+  template<class Predicate, typename T, typename U>
+  inline bool Test(const T& lhs, const T& rhs, const U& eps)
+  { return Test(Predicate(eps), lhs, rhs); }
+
+}
+
+
+#endif

Framework/Utilities/testBit.cc

+/**
+  \file
+  Test Bit functions
+
+  \author Hans Dembinski
+  \version $Id: testBit.cc 25126 2014-02-03 22:13:10Z darko $
+  \date 27 Jan 2014
+
+  \ingroup testing
+*/
+
+#include <corsika/utl/Test.h>
+#include <tst/Verify.h>
+#include <utl/Bit.h>
+#include <cppunit/extensions/HelperMacros.h>
+#include <cstdio>
+#include <iostream>
+#include <bitset>
+
+using namespace tst;
+using namespace utl;
+using namespace std;
+
+
+/**
+  \ingroup testing
+*/
+class TestBit : public CppUnit::TestFixture {
+
+  CPPUNIT_TEST_SUITE(TestBit);
+  CPPUNIT_TEST(TestGet);
+  CPPUNIT_TEST(TestSet);
+  CPPUNIT_TEST(TestMask);
+  CPPUNIT_TEST_SUITE_END();
+
+public:
+  void setUp() { }
+
+  void tearDown() { }
+
+  void
+  TestGet()
+  {
+    const int size = sizeof(int)*8;
+    const int bc2 = 12345;
+    int b2 = bc2;
+    bitset<size> b1(bc2);
+
+    ostringstream out1;
+    ostringstream out2;
+    ostringstream out3;
+    for (int i = 0; i < size; ++i) {
+      out1 << (b1[i] ? '^' : '.');
+      out2 << (AsBitArray(bc2)[i] ? '^' : '.');
+      out3 << (AsBitArray(b2)[i] ? '^' : '.');
+    }
+
+    CPPUNIT_ASSERT(Verify<Equal>(out1.str(), out2.str()));
+    CPPUNIT_ASSERT(Verify<Equal>(out1.str(), out3.str()));
+  }
+
+  void
+  TestSet()
+  {
+    const int size = sizeof(int)*8;
+    const int number = 12345;
+    bitset<size> b1(number);
+    int b2 = 11111;
+
+    for (int i = 0; i < size; ++i)
+      AsBitArray(b2)[i] = b1[i];
+
+    CPPUNIT_ASSERT(Verify<Equal>(b2, number));
+  }
+
+  void
+  TestMask()
+  {
+    const int n = (1 << 18) | (1 << 5);
+    int m = 0;
+
+    AsBitArray(m)[18] = true;
+    AsBitArray(m)[5] = true;
+    CPPUNIT_ASSERT(Verify<Equal>(n, m));
+
+    for (unsigned int i = 0; i < 8*sizeof(int); ++i)
+      AsBitArray(m)[i] = 0;
+    CPPUNIT_ASSERT(Verify<Equal>(m, 0));
+
+    m = 1;
+    AsBitArray(m).Mask(n, true);
+    CPPUNIT_ASSERT(Verify<Equal>(m, n+1));
+
+    AsBitArray(m).Mask(n, false);
+    CPPUNIT_ASSERT(Verify<Equal>(m, 1));
+  }
+
+};
+
+
+CPPUNIT_TEST_SUITE_REGISTRATION(TestBit);

Processes/StackInspector/StackInspector.cc

@@ -12,6 +12,8 @@
 #include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/logging/Logger.h>
+
 #include <iostream>
 using namespace std;
 
@@ -31,11 +33,8 @@ process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle
                                                                         Trajectory&,
                                                                         Stack& s) const {
 
-  // using namespace corsika::particles::io;
-
   static int countStep = 0;
   if (!fReport) return EProcessReturn::eOk;
-  // std::cout << "generation  " << countStep << std::endl;
   [[maybe_unused]] int i = 0;
   EnergyType Etot = 0_GeV;
 
@@ -55,9 +54,6 @@ double StackInspector<Stack, Trajectory>::MinStepLength(Particle&) const {
   return 0;
 }
 
-template <typename Stack, typename Trajectory>
-void StackInspector<Stack, Trajectory>::DoDiscrete(Particle&, Stack&) const {}
-
 template <typename Stack, typename Trajectory>
 void StackInspector<Stack, Trajectory>::Init() {}
 

Processes/StackInspector/StackInspector.h

@@ -36,9 +36,6 @@ namespace corsika::process {
       // template <typename Particle>
       double MinStepLength(Particle&) const;
 
-      // template <typename Particle, typename Stack>
-      void DoDiscrete(Particle&, Stack&) const;
-
     private:
       bool fReport;
     };

README.md

@@ -22,9 +22,9 @@ We also want to point you to the [MCnet guidelines](https://gitlab.ikp.kit.edu/A
 
 ## Installation
 
-Prerequisites: eigen3, cmake, g++, git. On Ubuntu 18.04, just do:
+Prerequisites: eigen3, boost, cmake, g++, git. On Ubuntu 18.04, just do:
 ```
-sudo apt-get install libeigen3-dev cmake g++ git
+sudo apt-get install libeigen3-dev libboost-dev cmake g++ git
 ```
 
 Follow these steps to download and install CORSIKA8-milestone1
@@ -40,3 +40,14 @@ make install
 make test
 ```
 and if you want to see how the Heitler model works and is implemented, see `Framework/Cascade/testCascade.cc` for a starting point. 
+
+### Generating doxygen documentation
+
+To generate the documentation, you need doxygen and graphviz. On Ubuntu 18.04, do:
+```
+sudo apt-get install doxygen graphviz
+```
+Switch to the corsika build directory and do
+```
+make doxygen
+```