Compare revisions

Felix Riehn · Marvin Gottowik · Marvin Gottowik · 6a0185b6 · 6a0185b6 · 6a0185b6
--- a/applications/c8_air_shower.cpp
+++ b/applications/c8_air_shower.cpp
@@ -195,6 +195,10 @@ int main(int argc, char** argv) {
      ->default_val(0.3)
      ->check(CLI::Range(0.000001, 1.e13))
      ->group("Config");
+  app.add_option("--taucut", "Min. kin. energy of tau leptons in tracking (GeV)")
+      ->default_val(0.3)
+      ->check(CLI::Range(0.000001, 1.e13))
+      ->group("Config");
  bool track_neutrinos = false;
  app.add_flag("--track-neutrinos", track_neutrinos, "switch on tracking of neutrinos")
      ->group("Config");
@@ -459,17 +463,19 @@ int main(int argc, char** argv) {
  HEPEnergyType const emcut = 1_GeV * app["--emcut"]->as<double>();
  HEPEnergyType const hadcut = 1_GeV * app["--hadcut"]->as<double>();
  HEPEnergyType const mucut = 1_GeV * app["--mucut"]->as<double>();
-  ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut,
+  HEPEnergyType const taucut = 1_GeV * app["--taucut"]->as<double>();
+  ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut, taucut,
                                             !track_neutrinos, dEdX);

  // tell proposal that we are interested in all energy losses above the particle cut
-  set_energy_production_threshold(Code::Electron, std::min({emcut, hadcut, mucut}));
-  set_energy_production_threshold(Code::Positron, std::min({emcut, hadcut, mucut}));
-  set_energy_production_threshold(Code::Photon, std::min({emcut, hadcut, mucut}));
-  set_energy_production_threshold(Code::MuMinus, std::min({emcut, hadcut, mucut}));
-  set_energy_production_threshold(Code::MuPlus, std::min({emcut, hadcut, mucut}));
-  set_energy_production_threshold(Code::TauMinus, std::min({emcut, hadcut, mucut}));
-  set_energy_production_threshold(Code::TauPlus, std::min({emcut, hadcut, mucut}));
+  auto const prod_threshold = std::min({emcut, hadcut, mucut, taucut});
+  set_energy_production_threshold(Code::Electron, prod_threshold);
+  set_energy_production_threshold(Code::Positron, prod_threshold);
+  set_energy_production_threshold(Code::Photon, prod_threshold);
+  set_energy_production_threshold(Code::MuMinus, prod_threshold);
+  set_energy_production_threshold(Code::MuPlus, prod_threshold);
+  set_energy_production_threshold(Code::TauMinus, prod_threshold);
+  set_energy_production_threshold(Code::TauPlus, prod_threshold);

  // energy threshold for high energy hadronic model. Affects LE/HE switch for
  // hadron interactions and the hadronic photon model in proposal

--- a/corsika/detail/modules/ParticleCut.inl
+++ b/corsika/detail/modules/ParticleCut.inl
@@ -17,19 +17,23 @@ namespace corsika {
  inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eEleCut,
                                           HEPEnergyType const ePhoCut,
                                           HEPEnergyType const eHadCut,
-                                           HEPEnergyType const eMuCut, bool const inv,
+                                           HEPEnergyType const eMuCut,
+                                           HEPEnergyType const eTauCut, bool const inv,
                                           TArgs&&... outputArgs)
      : TOutput(std::forward<TArgs>(outputArgs)...)
      , cut_electrons_(eEleCut)
      , cut_photons_(ePhoCut)
-      , cut_muons_(eMuCut)
      , cut_hadrons_(eHadCut)
+      , cut_muons_(eMuCut)
+      , cut_tau_(eTauCut)
      , doCutInv_(inv) {
    for (auto p : get_all_particles()) {
      if (is_hadron(p)) // nuclei are also hadrons
        set_kinetic_energy_propagation_threshold(p, eHadCut);
      else if (is_muon(p))
        set_kinetic_energy_propagation_threshold(p, eMuCut);
+      else if (p == Code::TauMinus || p == Code::TauPlus)
+        set_kinetic_energy_propagation_threshold(p, eTauCut);
      else if (p == Code::Electron || p == Code::Positron)
        set_kinetic_energy_propagation_threshold(p, eEleCut);
      else if (p == Code::Photon)
@@ -40,7 +44,8 @@ namespace corsika {
        "setting kinetic energy thresholds: electrons = {} GeV, photons = {} GeV, "
        "hadrons = {} GeV, "
        "muons = {} GeV",
-        eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV);
+        "tau = {} GeV", eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV,
+        eTauCut / 1_GeV);
  }

  template <typename TOutput>
@@ -152,6 +157,7 @@ namespace corsika {
    CORSIKA_LOG_DEBUG("kinetic energy threshold for photons is {} GeV",
                      cut_photons_ / 1_GeV);
    CORSIKA_LOG_DEBUG("kinetic energy threshold for muons is {} GeV", cut_muons_ / 1_GeV);
+    CORSIKA_LOG_DEBUG("kinetic energy threshold for tau is {} GeV", cut_tau_ / 1_GeV);
    CORSIKA_LOG_DEBUG("kinetic energy threshold for hadrons is {} GeV",
                      cut_hadrons_ / 1_GeV);

@@ -171,6 +177,7 @@ namespace corsika {
    node["cut_photons"] = cut_photons_ / 1_GeV;
    node["cut_muons"] = cut_muons_ / 1_GeV;
    node["cut_hadrons"] = cut_hadrons_ / 1_GeV;
+    node["cut_tau"] = cut_tau_ / 1_GeV;
    node["cut_invisibles"] = doCutInv_;
    for (auto const& cut : cuts_) {
      node[fmt::format("cut_{}", cut.first)] = cut.second / 1_GeV;

--- a/corsika/modules/ParticleCut.hpp
+++ b/corsika/modules/ParticleCut.hpp
@@ -44,8 +44,8 @@ namespace corsika {
     */
    template <typename... TArgs>
    ParticleCut(HEPEnergyType const eEleCut, HEPEnergyType const ePhoCut,
-                HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv,
-                TArgs&&... args);
+                HEPEnergyType const eHadCut, HEPEnergyType const eMuCut,
+                HEPEnergyType const eTauCut, bool const inv, TArgs&&... args);

    /**
     * particle cut with kinetic energy thresholds for all particles.
@@ -103,6 +103,7 @@ namespace corsika {
    HEPEnergyType getElectronKineticECut() const { return cut_electrons_; }
    HEPEnergyType getPhotonKineticECut() const { return cut_photons_; }
    HEPEnergyType getMuonKineticECut() const { return cut_muons_; }
+    HEPEnergyType getTauKineticECut() const { return cut_tau_; }
    HEPEnergyType getHadronKineticECut() const { return cut_hadrons_; }

    //! get configuration of this node, for output
@@ -116,8 +117,9 @@ namespace corsika {
  private:
    HEPEnergyType cut_electrons_;
    HEPEnergyType cut_photons_;
-    HEPEnergyType cut_muons_;
    HEPEnergyType cut_hadrons_;
+    HEPEnergyType cut_muons_;
+    HEPEnergyType cut_tau_;
    bool doCutInv_;
    std::unordered_map<Code const, HEPEnergyType const> cuts_;
  }; // namespace corsika

--- a/examples/cascade_examples/em_shower.cpp
+++ b/examples/cascade_examples/em_shower.cpp
@@ -146,8 +146,8 @@ int main(int argc, char** argv) {

  // setup processes, decays and interactions
  EnergyLossWriter energyloss{showerAxis, dX};
-  ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
-                                                   energyloss);
+  ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV,
+                                                   100_GeV, true, energyloss);

  corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
                                      corsika::setup::C7trackedParticles);

--- a/examples/cascade_examples/mars.cpp
+++ b/examples/cascade_examples/mars.cpp
@@ -313,7 +313,8 @@ int main(int argc, char** argv) {

  HEPEnergyType const emcut = 1_GeV;
  HEPEnergyType const hadcut = 1_GeV;
-  ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, dEdX);
+  ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, hadcut, true,
+                                             dEdX);

  // tell proposal that we are interested in all energy losses above the particle cut
  set_energy_production_threshold(Code::Electron, std::min({emcut, hadcut}));

--- a/examples/cascade_examples/radio_em_shower.cpp
+++ b/examples/cascade_examples/radio_em_shower.cpp
@@ -231,8 +231,8 @@ int main(int argc, char** argv) {

  // setup processes, decays and interactions
  EnergyLossWriter energyloss{showerAxis, dX};
-  ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
-                                                   energyloss);
+  ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV,
+                                                   100_GeV, true, energyloss);

  corsika::sibyll::Interaction sibyll(corsika::get_all_elements_in_universe(env),
                                      corsika::setup::C7trackedParticles);

--- a/examples/cascade_examples/water.cpp
+++ b/examples/cascade_examples/water.cpp
@@ -211,7 +211,8 @@ int main(int argc, char** argv) {
  // particle production threshold
  HEPEnergyType const emCut = eCut;
  HEPEnergyType const hadCut = eCut;
-  ParticleCut<SubWriter<decltype(dEdX)>> cut(emCut, emCut, hadCut, hadCut, true, dEdX);
+  ParticleCut<SubWriter<decltype(dEdX)>> cut(emCut, emCut, hadCut, hadCut, hadCut, true,
+                                             dEdX);

  // tell proposal that we are interested in all energy losses above the particle cut
  set_energy_production_threshold(Code::Electron, std::min({emCut, hadCut}));

--- a/tests/modules/testParticleCut.cpp
+++ b/tests/modules/testParticleCut.cpp
@@ -56,7 +56,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
  SECTION("cut on particle type: inv") {

    // particle cut with 20GeV threshold for all, also cut invisible
-    ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+    ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
    CHECK(cut.getHadronKineticECut() == 20_GeV);

    // add primary particle to stack
@@ -83,7 +83,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {

  SECTION("cut on particle type: em") {

-    ParticleCut cut(1_EeV, 1_EeV, 1_GeV, 1_GeV, false);
+    ParticleCut cut(1_EeV, 1_EeV, 1_GeV, 1_GeV, 1_GeV, false);

    // add primary particle to stack
    auto particle = stack.addParticle(std::make_tuple(
@@ -105,7 +105,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
  }

  SECTION("cut low energy") {
-    ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+    ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, 20_GeV, true);

    // add primary particle to stack
    auto particle = stack.addParticle(std::make_tuple(
@@ -133,8 +133,8 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
    CHECK(view.getSize() == 13);
  }

-  SECTION("cut low energy: electrons, photons, hadrons and muons") {
-    ParticleCut cut(5_MeV, 5_MeV, 5_GeV, 5_GeV, true);
+  SECTION("cut low energy: electrons, photons, hadrons, muons and tau") {
+    ParticleCut cut(5_MeV, 5_MeV, 5_GeV, 5_GeV, 5_MeV, true);

    // add primary particle to stack
    auto particle = stack.addParticle(std::make_tuple(Code::Proton, Eabove - Proton::mass,
@@ -177,7 +177,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
  }

  SECTION("cut on time") {
-    ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, false);
+    ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, 20_GeV, false);
    const TimeType too_late = 1_s;

    // add primary particle to stack
@@ -205,7 +205,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {

  SECTION("cut on doContinous, just invisibles") {

-    ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true);
+    ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, 20_GeV, true);

    // add particles, all with energies above the threshold
    // only cut is by species
No results found