From 3f6d2b517d918afd1103b6a9b8282aa6ddbe0c04 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 7 Oct 2020 22:10:02 +0200
Subject: [PATCH] better interface for CONEX, better output

---
 Documentation/Examples/hybrid_MC.cc  | 23 ++++++++++-------------
 Environment/MediumTypes.h            | 22 ++++++++++++++++++++++
 Processes/NullModel/testNullModel.cc | 12 ++++++------
 3 files changed, 38 insertions(+), 19 deletions(-)
 create mode 100644 Environment/MediumTypes.h

diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc
index 93531a2a0..036de755a 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/Documentation/Examples/hybrid_MC.cc
@@ -23,18 +23,17 @@
 #include <corsika/logging/Logging.h>
 #include <corsika/process/ProcessSequence.h>
 #include <corsika/process/StackProcess.h>
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
 #include <corsika/process/energy_loss/EnergyLoss.h>
 #include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
 #include <corsika/process/observation_plane/ObservationPlane.h>
 #include <corsika/process/on_shell_check/OnShellCheck.h>
-#include <corsika/process/energy_loss/EnergyLoss.h>
 #include <corsika/process/particle_cut/ParticleCut.h>
 #include <corsika/process/pythia/Decay.h>
 #include <corsika/process/sibyll/Decay.h>
 #include <corsika/process/sibyll/Interaction.h>
 #include <corsika/process/sibyll/NuclearInteraction.h>
 #include <corsika/process/switch_process/SwitchProcess.h>
-#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
 #include <corsika/process/tracking_line/TrackingLine.h>
 #include <corsika/process/urqmd/UrQMD.h>
 #include <corsika/random/RNGManager.h>
@@ -76,7 +75,6 @@ void registerRandomStreams(const int seed) {
 template <typename T>
 using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
 
-
 int main(int argc, char** argv) {
 
   logging::SetLevel(logging::level::info);
@@ -100,7 +98,8 @@ int main(int argc, char** argv) {
   EnvType env;
   const CoordinateSystem& rootCS = env.GetCoordinateSystem();
   Point const center{rootCS, 0_m, 0_m, 0_m};
-  environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{center};
+  environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{
+      center};
   builder.setNuclearComposition(
       {{particles::Code::Nitrogen, particles::Code::Oxygen},
        {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
@@ -216,8 +215,7 @@ int main(int argc, char** argv) {
   process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
 
   corsika::process::conex_source_cut::CONEXSourceCut conex(
-      center, showerAxis, t, injectionHeight, E0,
-      particles::GetPDG(particles::Code::Proton));
+      center, showerAxis, t, injectionHeight, E0, particles::Code::Proton);
 
   process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
 
@@ -237,8 +235,8 @@ int main(int argc, char** argv) {
                                                        55_GeV);
   auto decaySequence = decayPythia << decaySibyll;
 
-  auto sequence = switchProcess << reset_particle_mass << decaySequence 
-                                <<  eLoss << cut << conex << longprof << observationLevel;
+  auto sequence = switchProcess << reset_particle_mass << decaySequence << eLoss << cut
+                                << conex << longprof << observationLevel;
 
   // define air shower object, run simulation
   tracking_line::TrackingLine tracking;
@@ -251,22 +249,21 @@ int main(int argc, char** argv) {
   EAS.Run();
 
   cut.ShowResults();
-  //eLoss.ShowResults();
+  eLoss.showResults();
   observationLevel.ShowResults();
   const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
-    cut.GetEmEnergy() + //eLoss.GetEnergyLost() +
+                               cut.GetEmEnergy() + eLoss.energyLost() +
                                observationLevel.GetEnergyGround();
   cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
        << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
   observationLevel.Reset();
   cut.Reset();
-  //eLoss.Reset();
+  eLoss.reset();
 
-  //conex.addParticle(0, E0, 0_eV, emPosition, momentum.normalized(), 0_s);
   conex.SolveCE();
 
   auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
-    urqmdCounted.GetHistogram();
+                     urqmdCounted.GetHistogram();
 
   hists.saveLab("inthist_lab.txt");
   hists.saveCMS("inthist_cms.txt");
diff --git a/Environment/MediumTypes.h b/Environment/MediumTypes.h
new file mode 100644
index 000000000..bc44560fd
--- /dev/null
+++ b/Environment/MediumTypes.h
@@ -0,0 +1,22 @@
+/*
+ * (c) Copyright 2020 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::environment {
+
+  /**
+   * Medium types are useful most importantly for effective models
+   * like energy losses. a particular medium (mixture of components)
+   * may have specif properties not reflected by its mixture of
+   * components.
+   */
+
+  enum class EMediumType { eUnknown, eAir, eWater, eIce, eRock };
+
+} // namespace corsika::environment
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
index 0440bad4b..cc142637f 100644
--- a/Processes/NullModel/testNullModel.cc
+++ b/Processes/NullModel/testNullModel.cc
@@ -30,14 +30,14 @@ TEST_CASE("NullModel", "[processes]") {
   [[maybe_unused]] auto const& env_dummy = env;
   [[maybe_unused]] auto const& node_dummy = nodePtr;
 
-
-  auto [stack, view] = setup::testing::setupStack(particles::Code::Electron, 0,0, 100_GeV, nodePtr, cs);
-  [[maybe_unused]] const auto& dummyView = view; 
+  auto [stack, view] =
+      setup::testing::setupStack(particles::Code::Electron, 0, 0, 100_GeV, nodePtr, cs);
+  [[maybe_unused]] const auto& dummyView = view;
   auto particle = stack->first();
-  
+
   geometry::Point const origin(cs, {0_m, 0_m, 0_m});
-  geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second,
-                                                           0_m / second, 1_m / second);
+  geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
+                                                           1_m / second);
   geometry::Line line(origin, v);
   geometry::Trajectory<geometry::Line> track(line, 10_s);
 
-- 
GitLab