CONEX in vertial_EAS

Documentation/Examples/CMakeLists.txt

@@ -93,6 +93,7 @@ if (Pythia8_FOUND)
     ProcessUrQMD
     ProcessSwitch
     CORSIKAcascade
+    ProcessCONEXSourceCut
     ProcessEnergyLoss
     ProcessObservationPlane
     ProcessInteractionCounter

Documentation/Examples/vertical_EAS.cc

@@ -21,6 +21,7 @@
 #include <corsika/process/energy_loss/EnergyLoss.h>
 #include <corsika/process/interaction_counter/InteractionCounter.h>
 #include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
 #include <corsika/process/observation_plane/ObservationPlane.h>
 #include <corsika/process/on_shell_check/OnShellCheck.h>
 #include <corsika/process/particle_cut/ParticleCut.h>
@@ -118,7 +119,7 @@ int main(int argc, char** argv) {
   cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
        << endl;
 
-  auto const observationHeight = 1.4_km + builder.getEarthRadius();
+  auto const observationHeight = 0_km + builder.getEarthRadius();
   auto const injectionHeight = 112.75_km + builder.getEarthRadius();
   auto const t = -observationHeight * cos(thetaRad) +
                  sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
@@ -192,11 +193,18 @@ int main(int argc, char** argv) {
   Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
   process::observation_plane::ObservationPlane observationLevel(obsPlane,
                                                                 "particles.dat");
-
-  // assemble all processes into an ordered process list
-
+  
   process::UrQMD::UrQMD urqmd;
   process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
+  
+  sibyllNuc.Init();
+  sibyll.Init();
+                                                                
+  process::conex_source_cut::CONEXSourceCut conexSource(
+      center, showerAxis, t, injectionHeight, E0,
+      particles::GetPDG(particles::Code::Proton));
+
+  // assemble all processes into an ordered process list
 
   auto sibyllSequence = sibyllNucCounted << sibyllCounted;
   process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
@@ -218,6 +226,7 @@ int main(int argc, char** argv) {
   EAS.Run();
 
   eLoss.PrintProfile(); // print longitudinal profile
+  conexSource.SolveCE();
 
   cut.ShowResults();
   const HEPEnergyType Efinal =