adapted all examples

aed0552f · ralfulrich · ba0f72ac · aed0552f · aed0552f · aed0552f
Commit aed0552f authored 3 years ago by ralfulrich
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -24,8 +24,8 @@ namespace corsika {

  template <typename TOutput>
  inline ProcessReturn ObservationPlane<TOutput>::doContinuous(
-      corsika::setup::Stack::particle_type& particle,
-      corsika::setup::Trajectory&, bool const stepLimit) {
+      corsika::setup::Stack::particle_type& particle, corsika::setup::Trajectory&,
+      bool const stepLimit) {
    /*
       The current step did not yet reach the ObservationPlane, do nothing now and wait:
     */

--- a/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl
+++ b/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl
@@ -34,17 +34,23 @@ namespace corsika {

    // and build the streamer
    output_.buildStreamer();
+
+    setInit(true);
  }

  void ObservationPlaneWriterParquet::endOfShower() { ++shower_; }

  void ObservationPlaneWriterParquet::endOfLibrary() { output_.closeStreamer(); }

-  void ObservationPlaneWriterParquet::write(Code const& pid,
-                                            units::si::HEPEnergyType const& energy,
-                                            units::si::LengthType const& x,
-                                            units::si::LengthType const& y) {
-    using namespace units::si;
+  void ObservationPlaneWriterParquet::write(Code const& pid, HEPEnergyType const& energy,
+                                            LengthType const& x, LengthType const& y) {
+    if (!isInit()) {
+      std::runtime_error(
+          "ObservationPlaneWriterParquet not initialized. Either 1) add the "
+          "corresponding module to "
+          "the OutputManager, or 2) declare the module to write no output using "
+          "NoOutput.");
+    }

    // write the next row - we must write `shower_` first.
    *(output_.getWriter()) << shower_ << static_cast<int>(get_PDG(pid))

--- a/corsika/detail/modules/writers/TrackWriterParquet.inl
+++ b/corsika/detail/modules/writers/TrackWriterParquet.inl
@@ -38,16 +38,24 @@ namespace corsika {

    // and build the streamer
    output_.buildStreamer();
+
+    setInit(true);
  }

  void TrackWriterParquet::endOfShower() { ++shower_; }

  void TrackWriterParquet::endOfLibrary() { output_.closeStreamer(); }

-  void TrackWriterParquet::write(Code const& pid, units::si::HEPEnergyType const& energy,
+  void TrackWriterParquet::write(Code const& pid, HEPEnergyType const& energy,
                                 QuantityVector<length_d> const& start,
                                 QuantityVector<length_d> const& end) {
-    using namespace units::si;
+
+    if (!isInit()) {
+      std::runtime_error(
+          "TrackWriterParquet not initialized. Either 1) add the corresponding module to "
+          "the OutputManager, or 2) declare the module to write no output using "
+          "NoOutput.");
+    }

    // write the next row - we must write `shower_` first.
    // clang-format off

--- a/corsika/output/BaseOutput.hpp
+++ b/corsika/output/BaseOutput.hpp
@@ -20,8 +20,6 @@ namespace corsika {
  class BaseOutput {

  protected:
-    int shower_{0}; ///< The current event number.
-
    BaseOutput();

  public:
@@ -54,6 +52,21 @@ namespace corsika {
     * Get any summary information for the entire library.
     */
    virtual YAML::Node getSummary() { return YAML::Node(); };
+
+    /**
+     * Flag to indicate readiness.
+     */
+    bool isInit() const { return is_init_; }
+
+    /**
+     * Set init flag.
+     */
+    void setInit(bool const v) { is_init_ = v; }
+    
+  protected:
+    int shower_{0}; ///< The current event number.
+  private:
+    bool is_init_{false}; ///< flag to indicate readiness
  };

 } // namespace corsika

--- a/examples/boundary_example.cpp
+++ b/examples/boundary_example.cpp
@@ -13,7 +13,8 @@
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
 #include <corsika/framework/core/Logging.hpp>
-#include <corsika/output/DummyOutputManager.hpp>
+
+#include <corsika/output/OutputManager.hpp>

 #include <corsika/setup/SetupEnvironment.hpp>
 #include <corsika/setup/SetupStack.hpp>
@@ -114,6 +115,8 @@ int main() {
  world->addChild(std::move(target));
  universe.addChild(std::move(world));

+  OutputManager output("boundary_outputs");
+
  // setup processes, decays and interactions
  setup::Tracking tracking;

@@ -124,6 +127,8 @@ int main() {
  ParticleCut cut(50_GeV, true, true);

  TrackWriter trackWriter;
+  output.add("tracks", trackWriter); // register TrackWriter
+
  MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");

  // assemble all processes into an ordered process list
@@ -166,7 +171,6 @@ int main() {
  }

  // define air shower object, run simulation
-  DummyOutputManager output;
  Cascade EAS(env, tracking, sequence, output, stack);

  EAS.run();
@@ -177,4 +181,6 @@ int main() {
      (cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy());
  CORSIKA_LOG_INFO("Total energy (GeV): {} relative difference (%): {}", Efinal / 1_GeV,
                   (Efinal / E0 - 1.) * 100);
+
+  output.endOfLibrary();
 }
--- a/examples/cascade_example.cpp
+++ b/examples/cascade_example.cpp
@@ -11,10 +11,10 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
-
 #include <corsika/framework/utility/CorsikaFenv.hpp>
 #include <corsika/framework/core/Logging.hpp>
-#include <corsika/output/DummyOutputManager.hpp>
+
+#include <corsika/output/OutputManager.hpp>

 #include <corsika/media/Environment.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
@@ -55,8 +55,7 @@ using namespace std;
 //
 int main() {

-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-  logging::set_level(logging::level::trace);
+  logging::set_level(logging::level::info);

  std::cout << "cascade_example" << std::endl;

@@ -108,6 +107,8 @@ int main() {
      rootCS, 0_m, 0_m,
      height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe

+  OutputManager output("cascade_outputs");
+
  ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};

  {
@@ -142,6 +143,8 @@ int main() {
  ParticleCut cut(80_GeV, true, true);

  TrackWriter trackWriter;
+  output.add("tracks", trackWriter); // register TrackWriter
+  
  BetheBlochPDG eLoss{showerAxis};

  // assemble all processes into an ordered process list
@@ -149,7 +152,6 @@ int main() {
      make_sequence(stackInspect, sibyll, sibyllNuc, decay, eLoss, cut, trackWriter);

  // define air shower object, run simulation
-  DummyOutputManager output;
  Cascade EAS(env, tracking, sequence, output, stack);

  EAS.run();
@@ -164,4 +166,6 @@ int main() {
  cout << "total dEdX energy (GeV): " << eLoss.getTotal() / 1_GeV << endl
       << "relative difference (%): " << eLoss.getTotal() / E0 * 100 << endl;
  cut.reset();
+
+  output.endOfLibrary();
 }
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -11,10 +11,10 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
 #include <corsika/framework/geometry/Sphere.hpp>
-
 #include <corsika/framework/utility/CorsikaFenv.hpp>
 #include <corsika/framework/core/Logging.hpp>
-#include <corsika/output/DummyOutputManager.hpp>
+
+#include <corsika/output/OutputManager.hpp>

 #include <corsika/media/Environment.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
@@ -67,6 +67,8 @@ int main() {
  // initialize random number sequence(s)
  RNGManager::getInstance().registerRandomStream("cascade");

+  OutputManager output("cascade_proton_outputs");
+
  // setup environment, geometry
  using EnvType = setup::Environment;
  EnvType env;
@@ -116,7 +118,7 @@ int main() {

  // setup processes, decays and interactions
  setup::Tracking tracking;
-  StackInspector<setup::Stack> stackInspect(1, true, E0);
+  StackInspector<setup::Stack> stackInspect(1000, true, E0);

  RNGManager::getInstance().registerRandomStream("sibyll");
  RNGManager::getInstance().registerRandomStream("pythia");
@@ -132,6 +134,8 @@ int main() {
  // hadronicElastic(env);

  TrackWriter trackWriter;
+  output.add("tracks", trackWriter); // register TrackWriter
+
  ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
  BetheBlochPDG eLoss{showerAxis};

@@ -142,7 +146,6 @@ int main() {
  auto sequence = make_sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect);

  // define air shower object, run simulation
-  DummyOutputManager output;
  Cascade EAS(env, tracking, sequence, output, stack);
  EAS.run();

@@ -152,4 +155,6 @@ int main() {
      cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy();
  cout << "total energy (GeV): " << Efinal / 1_GeV << endl
       << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+
+  output.endOfLibrary();
 }
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -17,7 +17,8 @@
 #include <corsika/framework/utility/CorsikaFenv.hpp>
 #include <corsika/framework/process/InteractionCounter.hpp>
 #include <corsika/framework/core/Logging.hpp>
-#include <corsika/output/DummyOutputManager.hpp>
+
+#include <corsika/output/OutputManager.hpp>

 #include <corsika/media/Environment.hpp>
 #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
@@ -67,7 +68,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
 int main(int argc, char** argv) {

  logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");

  if (argc != 2) {
    std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
@@ -138,6 +138,7 @@ int main(int argc, char** argv) {
  std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
            << std::endl;

+  OutputManager output("em_shower_outputs");
  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};

  // setup processes, decays and interactions
@@ -148,6 +149,7 @@ int main(int argc, char** argv) {
  InteractionCounter emCascadeCounted(emCascade);

  TrackWriter trackWriter;
+  output.add("tracks", trackWriter); // register TrackWriter

  // long. profile; columns for gamma, e+, e- still need to be added
  LongitudinalProfile longprof{showerAxis};
@@ -155,12 +157,12 @@ int main(int argc, char** argv) {
  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
  ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
                                    "particles.dat");
+  output.add("obsplane", observationLevel);

  auto sequence = make_sequence(emCascadeCounted, emContinuous, longprof, cut,
                                observationLevel, trackWriter);
  // define air shower object, run simulation
  setup::Tracking tracking;
-  DummyOutputManager output;
  Cascade EAS(env, tracking, sequence, output, stack);

  // to fix the point of first interaction, uncomment the following two lines:
@@ -185,4 +187,6 @@ int main(int argc, char** argv) {
  save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
  save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
  longprof.save("longprof_emShower.txt");
+
+  output.endOfLibrary();
 }
--- a/examples/geometry_example.cpp
+++ b/examples/geometry_example.cpp
@@ -21,7 +21,6 @@ using namespace corsika;
 int main() {

  logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");

  CORSIKA_LOG_INFO("geometry_example");


--- a/examples/helix_example.cpp
+++ b/examples/helix_example.cpp
@@ -21,7 +21,6 @@ using namespace corsika;
 int main() {

  logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");

  CORSIKA_LOG_INFO("helix_example");


--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -23,9 +23,10 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/utility/CorsikaFenv.hpp>
 #include <corsika/framework/core/Cascade.hpp>
-#include <corsika/output/DummyOutputManager.hpp>
 #include <corsika/framework/geometry/PhysicalGeometry.hpp>

+#include <corsika/output/OutputManager.hpp>
+
 #include <corsika/media/Environment.hpp>
 #include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
@@ -87,7 +88,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
 int main(int argc, char** argv) {

  logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");

  CORSIKA_LOG_INFO("hybrid_MC");

@@ -172,6 +172,7 @@ int main(int argc, char** argv) {
  std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
            << std::endl;

+  OutputManager output("hybrid_MC_outputs");
  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};

  // setup processes, decays and interactions
@@ -220,7 +221,8 @@ int main(int argc, char** argv) {
  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
  ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
                                    "particles.dat");
-
+  output.add("obsplane", observationLevel);
+  
  corsika::urqmd::UrQMD urqmd_model;
  InteractionCounter urqmdCounted{urqmd_model};

@@ -241,7 +243,6 @@ int main(int argc, char** argv) {

  // define air shower object, run simulation
  setup::Tracking tracking;
-  DummyOutputManager output;
  Cascade EAS(env, tracking, sequence, output, stack);

  // to fix the point of first interaction, uncomment the following two lines:
@@ -274,4 +275,6 @@ int main(int argc, char** argv) {

  std::ofstream finish("finished");
  finish << "run completed without error" << std::endl;
+
+  output.endOfLibrary();
 }
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -99,7 +99,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;

 int main(int argc, char** argv) {

-  corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
  logging::set_level(logging::level::info);

  CORSIKA_LOG_INFO("vertical_EAS");