also re-enable all tests and examples

corsika/detail/modules/TrackWriter.inl

@@ -15,13 +15,14 @@
 
 namespace corsika {
 
-  template <typename TOutput>
-  inline TrackWriter<TOutput>::TrackWriter() {}
+  template <typename TOutputWriter>
+  inline TrackWriter<TOutputWriter>::TrackWriter() {}
 
-  template <typename TOutput>
+  template <typename TOutputWriter>
   template <typename TParticle, typename TTrack>
-  inline ProcessReturn TrackWriter<TOutput>::doContinuous(TParticle const& vP,
-                                                          TTrack const& vT, bool const) {
+  inline ProcessReturn TrackWriter<TOutputWriter>::doContinuous(TParticle const& vP,
+                                                                TTrack const& vT,
+                                                                bool const) {
 
     auto const start = vT.getPosition(0).getCoordinates();
     auto const end = vT.getPosition(1).getCoordinates();
@@ -32,15 +33,15 @@ namespace corsika {
     return ProcessReturn::Ok;
   }
 
-  template <typename TOutput>
+  template <typename TOutputWriter>
   template <typename TParticle, typename TTrack>
-  inline LengthType TrackWriter<TOutput>::getMaxStepLength(TParticle const&,
-                                                           TTrack const&) {
+  inline LengthType TrackWriter<TOutputWriter>::getMaxStepLength(TParticle const&,
+                                                                 TTrack const&) {
     return meter * std::numeric_limits<double>::infinity();
   }
 
-  template <typename TOutput>
-  YAML::Node TrackWriter<TOutput>::getConfig() const {
+  template <typename TOutputWriter>
+  YAML::Node TrackWriter<TOutputWriter>::getConfig() const {
     using namespace units::si;
 
     YAML::Node node;

corsika/detail/output/ParquetStreamer.inl

@@ -30,9 +30,9 @@ namespace corsika {
     fields_.push_back(parquet::schema::PrimitiveNode::Make(args...));
   }
 
-  void ParquetStreamer::enableCompression(int const level) {
-    builder_.compression(parquet::Compression::ZSTD);
-    builder_.compression_level(level);
+  void ParquetStreamer::enableCompression(int const /*level*/) {
+    // builder_.compression(parquet::Compression::ZSTD);
+    // builder_.compression_level(level);
   }
 
   void ParquetStreamer::buildStreamer() {

corsika/output/NoOutput.hpp

+/*
+ * (c) Copyright 2021 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
+
+namespace corsika {
+
+  /**
+   * This is the base class for all outputs so that they
+   * can be stored in homogeneous containers.
+   */
+  class NoOutput : public BaseOutput {
+
+  protected:
+    NoOutput() {}
+
+  public:
+    /**
+     * Called at the start of each run.
+     */
+    void startOfLibrary(boost::filesystem::path const&) final override {}
+
+    /**
+     * Called at the start of each event/shower.
+     */
+    void startOfShower() final override {}
+
+    /**
+     * Called at the end of each event/shower.
+     */
+    void endOfShower() final override {}
+
+    /**
+     * Called at the end of each run.
+     */
+    void endOfLibrary() final override {}
+
+    /**
+     * Get the configuration of this output.
+     */
+    YAML::Node getConfig() const { return YAML::Node(); };
+
+    /**
+     * Get any summary information for the entire library.
+     */
+    YAML::Node getSummary() final override { return YAML::Node(); };
+
+  protected:
+    void write(Code const&, units::si::HEPEnergyType const&, units::si::LengthType const&,
+               units::si::LengthType const&) {}
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/output/BaseOutput.inl>

examples/vertical_EAS.cpp

@@ -136,6 +136,7 @@ int main(int argc, char** argv) {
       {{Code::Nitrogen, Code::Oxygen},
        {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
 
+  builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
   builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
   builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
   builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
@@ -143,22 +144,6 @@ int main(int argc, char** argv) {
   builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
   builder.assemble(env);
 
-  CORSIKA_LOG_DEBUG(
-      "environment setup: universe={}, layer1={}, layer2={}, layer3={}, layer4={}, "
-      "layer5={}",
-      fmt::ptr(env.getUniverse()->getContainingNode(
-          Point(rootCS, {constants::EarthRadius::Mean + 130_km, 0_m, 0_m}))),
-      fmt::ptr(env.getUniverse()->getContainingNode(
-          Point(rootCS, {constants::EarthRadius::Mean + 110_km, 0_m, 0_m}))),
-      fmt::ptr(env.getUniverse()->getContainingNode(
-          Point(rootCS, {constants::EarthRadius::Mean + 50_km, 0_m, 0_m}))),
-      fmt::ptr(env.getUniverse()->getContainingNode(
-          Point(rootCS, {constants::EarthRadius::Mean + 20_km, 0_m, 0_m}))),
-      fmt::ptr(env.getUniverse()->getContainingNode(
-          Point(rootCS, {constants::EarthRadius::Mean + 5_km, 0_m, 0_m}))),
-      fmt::ptr(env.getUniverse()->getContainingNode(
-          Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m}))));
-
   // pre-setup particle stack
   unsigned short const A = std::stoi(std::string(argv[1]));
   Code beamCode;
@@ -216,6 +201,9 @@ int main(int argc, char** argv) {
 
   ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
 
+  // create the output manager that we then register outputs with
+  OutputManager output("vertical_EAS_outputs");
+
   // setup processes, decays and interactions
 
   // corsika::qgsjetII::Interaction qgsjet;
@@ -282,8 +270,6 @@ int main(int argc, char** argv) {
     string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
     string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
     string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
-    string const tracks_file = "tracks_" + to_string(i_shower) + ".dat";
-    string const particles_file = "particles_" + to_string(i_shower) + ".dat";
 
     std::cout << std::endl;
     std::cout << "Shower " << i_shower << "/" << number_showers << std::endl;
@@ -291,9 +277,50 @@ int main(int argc, char** argv) {
     // setup particle stack, and add primary particle
     setup::Stack stack;
     stack.clear();
+    unsigned short const A = std::stoi(std::string(argv[1]));
+    Code beamCode;
+    HEPEnergyType mass;
+    unsigned short Z = 0;
+    if (A > 0) {
+      beamCode = Code::Nucleus;
+      Z = std::stoi(std::string(argv[2]));
+      mass = get_nucleus_mass(A, Z);
+    } else {
+      int pdg = std::stoi(std::string(argv[2]));
+      beamCode = convert_from_PDG(PDGCode(pdg));
+      mass = get_mass(beamCode);
+    }
+    HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
+    double theta = 0.;
+    auto const thetaRad = theta / 180. * M_PI;
+
+    auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+      return sqrt((Elab - m) * (Elab + m));
+    };
+    HEPMomentumType P0 = elab2plab(E0, mass);
+    auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+      return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+    };
+
+    auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+    auto plab = MomentumVector(rootCS, {px, py, pz});
+    cout << "input particle: " << beamCode << endl;
+    cout << "input angles: theta=" << theta << endl;
+    cout << "input momentum: " << plab.getComponents() / 1_GeV
+         << ", norm = " << plab.getNorm() << endl;
+
+    auto const observationHeight = 0_km + builder.getEarthRadius();
+    auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+    auto const t = (injectionHeight - observationHeight) / cos(thetaRad);
+    Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+    Point const injectionPos =
+        showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+    std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
 
     if (A > 1) {
       stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
+
     } else {
       if (A == 1) {
         if (Z == 1) {
@@ -314,10 +341,8 @@ int main(int argc, char** argv) {
     std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
               << std::endl;
 
-    ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
-
-    // create the output manager that we then register outputs with
-    OutputManager output("vertical_EAS_outputs");
+    ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env,
+                                false};
 
     // setup processes, decays and interactions
 
@@ -353,10 +378,11 @@ int main(int argc, char** argv) {
 
     decaySibyll.printDecayConfig();
 
-    ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
-    corsika::proposal::Interaction emCascade(env);
-    corsika::proposal::ContinuousProcess emContinuous(env);
-    InteractionCounter emCascadeCounted(emCascade);
+    ParticleCut cut{60_GeV, true, true};
+    // corsika::proposal::Interaction emCascade(env);
+    // corsika::proposal::ContinuousProcess emContinuous(env);
+    // InteractionCounter emCascadeCounted(emCascade);
+    BetheBlochPDG emContinuous(showerAxis);
 
     OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
     TrackWriter trackWriter;
@@ -369,26 +395,6 @@ int main(int argc, char** argv) {
     // register the observation plane with the output
     output.add("obsplane", observationLevel);
 
-    corsika::urqmd::UrQMD urqmd;
-    InteractionCounter urqmdCounted{urqmd};
-    StackInspector<setup::Stack> stackInspect(50000, false, E0);
-
-    // assemble all processes into an ordered process list
-    struct EnergySwitch {
-      HEPEnergyType cutE_;
-      EnergySwitch(HEPEnergyType cutE)
-          : cutE_(cutE) {}
-      SwitchResult operator()(const Particle& p) {
-        if (p.getEnergy() < cutE_)
-          return SwitchResult::First;
-        else
-          return SwitchResult::Second;
-      }
-    };
-    auto hadronSequence =
-        make_select(urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted),
-                    EnergySwitch(55_GeV));
-    auto decaySequence = make_sequence(decayPythia, decaySibyll);
     auto sequence =
         make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
                       emContinuous, cut, trackWriter, observationLevel, longprof);
@@ -403,16 +409,16 @@ int main(int argc, char** argv) {
     EAS.run();
 
     cut.showResults();
-    emContinuous.showResults();
+    // emContinuous.showResults();
     observationLevel.showResults();
     const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
-                                 cut.getEmEnergy() + emContinuous.getEnergyLost() +
+                                 cut.getEmEnergy() + // emContinuous.getEnergyLost() +
                                  observationLevel.getEnergyGround();
     cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
          << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
     observationLevel.reset();
     cut.reset();
-    emContinuous.reset();
+    // emContinuous.reset();
 
     auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
                        urqmdCounted.getHistogram();

tests/modules/testObservationPlane.cpp

@@ -17,6 +17,8 @@
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
+#include <corsika/output/NoOutput.hpp>
+
 #include <SetupTestEnvironment.hpp>
 #include <SetupTestStack.hpp>
 #include <SetupTestTrajectory.hpp>
@@ -54,7 +56,7 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
   SECTION("horizontal plane") {
 
     Plane const obsPlane(Point(cs, {10_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
-    ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
+    ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
 
     LengthType const length = obs.getMaxStepLength(particle, no_used_track);
     ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
@@ -65,10 +67,10 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
 
   SECTION("transparent plane") {
     Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
-    ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 0., 1.}));
+    ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 0., 1.}));
 
     LengthType const length = obs.getMaxStepLength(particle, no_used_track);
-    ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
+    ProcessReturn const ret = obs.doContinuous(particle, no_used_track, false);
 
     CHECK(length / 1_m == Approx(1).margin(1e-4));
     CHECK(ret == ProcessReturn::Ok);
@@ -83,7 +85,7 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
 
     Plane const obsPlane(Point(cs, {10_m, 5_m, 5_m}),
                          DirectionVector(cs, {1, 0.1, -0.05}));
-    ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
+    ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
 
     LengthType const length = obs.getMaxStepLength(particle, no_used_track);
     ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);