Merge branch '351-multiple-showers' into 'master'

Resolve "Multiple showers"

Closes #351

See merge request AirShowerPhysics/corsika!326

examples/CMakeLists.txt

@@ -7,7 +7,7 @@ set (CMAKE_VERBOSE_MAKEFILE OFF) # this can be changed with `make VERBOSE=1`
 find_package (corsika CONFIG REQUIRED)
 
 # this is only for CORSIKA_REGISTER_EXAMPLE
-include ("${CMAKE_CURRENT_SOURCE_DIR}/CMakeHelper.cmake") 
+include ("${CMAKE_CURRENT_SOURCE_DIR}/CMakeHelper.cmake")
 
 add_executable (helix_example helix_example.cpp)
 target_link_libraries (helix_example CORSIKA8::CORSIKA8)
@@ -35,7 +35,7 @@ CORSIKA_REGISTER_EXAMPLE (cascade_proton_example)
 
 add_executable (vertical_EAS vertical_EAS.cpp)
 target_link_libraries (vertical_EAS CORSIKA8::CORSIKA8)
-CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
+CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000. 1)
 
 add_executable (stopping_power stopping_power.cpp)
 target_link_libraries (stopping_power CORSIKA8::CORSIKA8)

examples/vertical_EAS.cpp

@@ -91,6 +91,11 @@ void registerRandomStreams(const int seed) {
 template <typename T>
 using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
 
+// argv : 1.number of nucleons, 2.number of protons,
+//        3.total energy in GeV, 4.number of showers,
+//        5.output directory
+//        6.seed (0 by default to generate random values for all)
+
 int main(int argc, char** argv) {
 
   corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
@@ -98,7 +103,7 @@ int main(int argc, char** argv) {
 
   CORSIKA_LOG_INFO("vertical_EAS");
 
-  if (argc < 4) {
+  if (argc < 5) {
     std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed] \n"
                  "       if A=0, Z is interpreted as PDG code \n"
                  "       if no seed is given, a random seed is chosen \n"
@@ -107,8 +112,16 @@ int main(int argc, char** argv) {
   }
   feenableexcept(FE_INVALID);
 
+  string output_dir = "";
   int seed = 0;
-  if (argc > 4) seed = std::stoi(std::string(argv[4]));
+  int number_showers = std::stoi(std::string(argv[4]));
+  if (argc > 5) {
+    output_dir = argv[5];
+    if (output_dir.back() != '/' && output_dir != "") output_dir += '/';
+  }
+
+  if (argc > 6) { seed = std::stoi(std::string(argv[6])); }
+
   // initialize random number sequence(s)
   registerRandomStreams(seed);
 
@@ -150,173 +163,190 @@ int main(int argc, char** argv) {
       fmt::ptr(env.getUniverse()->getContainingNode(
           Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m}))));
 
-  // 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 = -observationHeight * cos(thetaRad) +
-                 sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
-                      static_pow<2>(injectionHeight));
-  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) {
-        stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
-      } else if (Z == 0) {
-        stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
-      } else {
-        std::cerr << "illegal parameters" << std::endl;
-        return EXIT_FAILURE;
-      }
+  for (int i_shower = 1; i_shower < number_showers + 1; i_shower++) {
+
+    // directory for outputs
+    string labHist_dir =
+        output_dir + "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
+    string cMSHist_dir =
+        output_dir + "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
+    string longprof_dir =
+        output_dir + "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
+    string tracks_dir = output_dir + "tracks_" + to_string(i_shower) + ".dat";
+    string particles_dir = output_dir + "particles_" + to_string(i_shower) + ".dat";
+
+    std::cout << std::endl;
+    std::cout << "Shower " << i_shower << "/" << number_showers << std::endl;
+
+    // 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 {
-      stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
+      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 = -observationHeight * cos(thetaRad) +
+                   sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
+                        static_pow<2>(injectionHeight));
+    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));
 
-  // we make the axis much longer than the inj-core distance since the
-  // profile will go beyond the core, depending on zenith angle
-  std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
-            << std::endl;
-
-  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
-
-  // setup processes, decays and interactions
-
-  // corsika::qgsjetII::Interaction qgsjet;
-  corsika::sibyll::Interaction sibyll;
-  InteractionCounter sibyllCounted(sibyll);
-
-  corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
-  InteractionCounter sibyllNucCounted(sibyllNuc);
-
-  corsika::pythia8::Decay decayPythia;
-
-  // use sibyll decay routine for decays of particles unknown to pythia
-  corsika::sibyll::Decay decaySibyll{{
-      Code::N1440Plus,
-      Code::N1440MinusBar,
-      Code::N1440_0,
-      Code::N1440_0Bar,
-      Code::N1710Plus,
-      Code::N1710MinusBar,
-      Code::N1710_0,
-      Code::N1710_0Bar,
-
-      Code::Pi1300Plus,
-      Code::Pi1300Minus,
-      Code::Pi1300_0,
-
-      Code::KStar0_1430_0,
-      Code::KStar0_1430_0Bar,
-      Code::KStar0_1430_Plus,
-      Code::KStar0_1430_MinusBar,
-  }};
-
-  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);
-
-  OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
-  TrackWriter trackWriter("tracks.dat");
-
-  LongitudinalProfile longprof{showerAxis};
-
-  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
-  ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
-                                    "particles.dat");
-
-  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;
+    } else {
+      if (A == 1) {
+        if (Z == 1) {
+          stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
+        } else if (Z == 0) {
+          stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
+        } else {
+          std::cerr << "illegal parameters" << std::endl;
+          return EXIT_FAILURE;
+        }
+      } else {
+        stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns));
+      }
     }
-  };
-  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);
-
-  // define air shower object, run simulation
-  setup::Tracking tracking;
-  Cascade EAS(env, tracking, sequence, stack);
-
-  // to fix the point of first interaction, uncomment the following two lines:
-  //  EAS.forceInteraction();
-
-  EAS.run();
-
-  cut.showResults();
-  emContinuous.showResults();
-  observationLevel.showResults();
-  const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
-                               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();
-
-  auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
-                     urqmdCounted.getHistogram();
-
-  save_hist(hists.labHist(), "inthist_lab_verticalEAS.npz", true);
-  save_hist(hists.CMSHist(), "inthist_cms_verticalEAS.npz", true);
-  longprof.save("longprof_verticalEAS.txt");
+
+    // we make the axis much longer than the inj-core distance since the
+    // profile will go beyond the core, depending on zenith angle
+    std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
+              << std::endl;
+
+    ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
+
+    // setup processes, decays and interactions
+
+    // corsika::qgsjetII::Interaction qgsjet;
+    corsika::sibyll::Interaction sibyll;
+    InteractionCounter sibyllCounted(sibyll);
+
+    corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+    InteractionCounter sibyllNucCounted(sibyllNuc);
+
+    corsika::pythia8::Decay decayPythia;
+
+    // use sibyll decay routine for decays of particles unknown to pythia
+    corsika::sibyll::Decay decaySibyll{{
+        Code::N1440Plus,
+        Code::N1440MinusBar,
+        Code::N1440_0,
+        Code::N1440_0Bar,
+        Code::N1710Plus,
+        Code::N1710MinusBar,
+        Code::N1710_0,
+        Code::N1710_0Bar,
+
+        Code::Pi1300Plus,
+        Code::Pi1300Minus,
+        Code::Pi1300_0,
+
+        Code::KStar0_1430_0,
+        Code::KStar0_1430_0Bar,
+        Code::KStar0_1430_Plus,
+        Code::KStar0_1430_MinusBar,
+    }};
+
+    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);
+
+    OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
+    TrackWriter trackWriter(tracks_dir);
+
+    LongitudinalProfile longprof{showerAxis};
+
+    Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+    ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
+                                      particles_dir);
+
+    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);
+
+    // define air shower object, run simulation
+    setup::Tracking tracking;
+    Cascade EAS(env, tracking, sequence, stack);
+
+    // to fix the point of first interaction, uncomment the following two lines:
+    //  EAS.forceInteraction();
+
+    EAS.run();
+
+    cut.showResults();
+    emContinuous.showResults();
+    observationLevel.showResults();
+    const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
+                                 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();
+
+    auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
+                       urqmdCounted.getHistogram();
+
+    save_hist(hists.labHist(), labHist_dir, true);
+    save_hist(hists.CMSHist(), cMSHist_dir, true);
+    longprof.save(longprof_dir);
+  }
 }