Merge branch 'proposal_tables_v2' into 'master'

Restructuring of PROPOSAL table management

See merge request !557

corsika/detail/modules/proposal/ContinuousProcess.inl

@@ -28,11 +28,9 @@ namespace corsika::proposal {
     if (code == Code::Photon) return; // no continuous builders needed for photons
 
     // interpolate the crosssection for given media and energy cut. These may
-    // take some minutes if you have to build the tables and cannot read the
+    // take some minutes if you have to build the tables and cannot read the tables
     // from disk
-    auto const emCut = get_energy_production_threshold(
-        code); //! energy resolutions globally defined for individual particles
-    auto c = p_cross->second(media.at(comp.getHash()), emCut);
+    auto c = p_cross->second(media.at(comp.getHash()), proposal_energycutsettings[code]);
 
     // choose multiple scattering model
     static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::MoliereInterpol;

corsika/detail/modules/proposal/InteractionModel.inl

@@ -37,12 +37,9 @@ namespace corsika::proposal {
       throw std::runtime_error("PROPOSAL could not find corresponding builder");
 
     // interpolate the crosssection for given media and energy cut. These may
-    // take some minutes if you have to build the tables and cannot read the
+    // take some minutes if you have to build the tables and cannot read the tables
     // from disk
-    auto const emCut = get_energy_production_threshold(
-        code); //! energy resolutions globally defined for individual particles
-
-    auto c = p_cross->second(media.at(comp.getHash()), emCut);
+    auto c = p_cross->second(media.at(comp.getHash()), proposal_energycutsettings[code]);
 
     // Look which interactions take place and build the corresponding
     // interaction and secondary builder. The interaction integral will

corsika/detail/modules/proposal/ProposalProcessBase.inl

@@ -27,6 +27,27 @@ namespace corsika::proposal {
     return false;
   }
 
+  inline HEPEnergyType ProposalProcessBase::getOptimizedEmCut(Code code) const {
+    // get energy above which energy losses need to be considered
+    auto const production_threshold = get_energy_production_threshold(code);
+
+    HEPEnergyType lowest_table_value = 0_GeV;
+
+    // find tables for EnergyCuts closest (but still smaller than) production_threshold
+    for (auto table_energy : energycut_table_values) {
+      if (table_energy <= production_threshold && table_energy > lowest_table_value) {
+        lowest_table_value = table_energy;
+      }
+    }
+
+    if (lowest_table_value == 0_GeV) {
+      // no appropriate table available
+      return production_threshold;
+    }
+
+    return lowest_table_value;
+  };
+
   template <typename TEnvironment>
   inline ProposalProcessBase::ProposalProcessBase(TEnvironment const& _env) {
     _env.getUniverse()->walk([&](auto& vtn) {
@@ -54,6 +75,34 @@ namespace corsika::proposal {
     //! path, otherwise interpolation tables would only stored in main memory if
     //! no explicit intrpolation def is specified.
     PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str();
+
+    //! Initialize EnergyCutSettings
+    for (auto particle_code : tracked) {
+      if (particle_code == Code::Photon) {
+        // no EnergyCut for photon, only-stochastic propagation
+        continue;
+      }
+      // use optimized emcut for which tables should be available
+      auto optimized_ecut = getOptimizedEmCut(particle_code);
+      CORSIKA_LOG_INFO("PROPOSAL: Use tables with EnergyCut of {} for particle {}.",
+                       optimized_ecut, particle_code);
+      proposal_energycutsettings[particle_code] = optimized_ecut;
+    }
+
+    //! Initialize PROPOSAL tables for all media and all particles
+    for (auto medium : media) {
+      for (auto particle_code : tracked) {
+        buildTables(medium.second, particle_code,
+                    proposal_energycutsettings[particle_code]);
+      }
+    }
+  }
+
+  void ProposalProcessBase::buildTables(PROPOSAL::Medium medium, Code particle_code,
+                                        HEPEnergyType emcut) {
+    CORSIKA_LOG_DEBUG("PROPOSAL: Initialize tables for particle {} in {}", particle_code,
+                      medium.GetName());
+    cross[particle_code](medium, emcut);
   }
 
   inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const

corsika/modules/proposal/ProposalProcessBase.hpp

@@ -32,6 +32,14 @@ namespace corsika::proposal {
   //!
   static constexpr double v_cut = 0.01;
 
+  //!
+  //! List of EnergyCut values for which CORSIKA will, by default, create and provide
+  //! PROPOSAL tables.
+  //!
+  static constexpr std::array<HEPEnergyType, 10> energycut_table_values{
+      1000_MeV, 100_MeV, 20_MeV,   10_MeV,   3_MeV,
+      1_MeV,    0.4_MeV, 0.25_MeV, 0.15_MeV, 0.05_MeV};
+
   //!
   //! Internal map from particle codes to particle properties required for
   //! crosssections, decay and scattering algorithms. In the future the
@@ -223,6 +231,10 @@ namespace corsika::proposal {
         media; //!< maps nuclear composition from univers to media to produce
                //!< crosssections, which requires further ionization constants.
 
+    //!< save emcut for tracked particles
+    std::unordered_map<Code, corsika::units::si::HEPEnergyType>
+        proposal_energycutsettings;
+
     //!
     //! Store cut and  nuclear composition of the whole universe in media which are
     //! required for creating crosssections by proposal.
@@ -235,6 +247,12 @@ namespace corsika::proposal {
     //!
     bool canInteract(Code pcode) const;
 
+    //!
+    //! Finds the optimal EnergyCut for which PROPOSAL tables should (by default) be
+    //! available.
+    //!
+    HEPEnergyType getOptimizedEmCut(Code code) const;
+
     using calc_key_t = std::pair<std::size_t, Code>;
 
     //!
@@ -250,6 +268,11 @@ namespace corsika::proposal {
     //!
     virtual void buildCalculator(Code, NuclearComposition const&) = 0;
 
+    //!
+    //! Initialize PROPOSAL tables for given medium, code, and energy cut
+    //!
+    void buildTables(PROPOSAL::Medium, Code, HEPEnergyType);
+
     //!
     //! Searches the particle dependet calculator dependent of actuall medium composition
     //! and particle type. If no calculator is found, the corresponding new calculator is

examples/corsika.cpp

@@ -398,12 +398,14 @@ int main(int argc, char** argv) {
   HEPEnergyType const mucut = 1_GeV * app["--mucut"]->as<double>();
   ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut, true, dEdX);
 
-  // tell proposal that we are interested in all energy losses above the emcut
-  set_energy_production_threshold(Code::Electron, emcut);
-  set_energy_production_threshold(Code::Positron, emcut);
-  set_energy_production_threshold(Code::Photon, emcut);
-  set_energy_production_threshold(Code::MuMinus, mucut);
-  set_energy_production_threshold(Code::MuPlus, mucut);
+  // 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}));
 
   // energy threshold for high energy hadronic model. Affects LE/HE switch for
   // hadron interactions and the hadronic photon model in proposal