ProposalProcessBase.hpp

/*
 * (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

#include <PROPOSAL/PROPOSAL.h>

#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/random/RNGManager.hpp>

#include <array>

namespace corsika::proposal {

  static constexpr std::array<Code, 7> tracked{
      Code::Photon, Code::Electron, Code::Positron, Code::MuMinus,
      Code::MuPlus, Code::TauPlus,  Code::TauMinus,
  };

  static constexpr double v_cut = 0.01;

  static std::map<Code, PROPOSAL::ParticleDef> particle = {
      {Code::Photon, PROPOSAL::GammaDef()},   {Code::Electron, PROPOSAL::EMinusDef()},
      {Code::Positron, PROPOSAL::EPlusDef()}, {Code::MuMinus, PROPOSAL::MuMinusDef()},
      {Code::MuPlus, PROPOSAL::MuPlusDef()},  {Code::TauMinus, PROPOSAL::TauMinusDef()},
      {Code::TauPlus, PROPOSAL::TauPlusDef()}};


  // Cross sections for muons and taus
  template <typename T>
  auto cross_builder = [](PROPOSAL::Medium& medium,
                          corsika::units::si::HEPEnergyType
                              emCut) { 
    auto particle_def = T();
    CORSIKA_LOG_DEBUG("PROPOSAL: Generate general cross sections for particle type {} ",
                      particle_def.name);
    auto interpolate = true;

    auto p_cut =
        std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, v_cut, false);

    // do not use v_cut for ionization due to numerical problems
    auto p_cut_no_vcut =
        std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, false);

    auto cross_vec = std::vector<std::shared_ptr<PROPOSAL::CrossSectionBase>>();
    cross_vec.push_back(PROPOSAL::make_crosssection(
        PROPOSAL::crosssection::BremsKelnerKokoulinPetrukhin{false}, particle_def, medium,
        p_cut, interpolate));
    cross_vec.push_back(PROPOSAL::make_crosssection(
        PROPOSAL::crosssection::EpairKelnerKokoulinPetrukhin{false}, particle_def, medium,
        p_cut, interpolate));
    cross_vec.push_back(PROPOSAL::make_crosssection(
        PROPOSAL::crosssection::IonizBetheBlochRossi{*p_cut_no_vcut}, particle_def,
        medium, p_cut_no_vcut, interpolate));
    cross_vec.push_back(PROPOSAL::make_crosssection(
        PROPOSAL::crosssection::PhotoAbramowiczLevinLevyMaor97{
            std::make_unique<PROPOSAL::crosssection::ShadowButkevichMikheyev>()},
        particle_def, medium, p_cut, interpolate));

    return cross_vec;
  };

  // cross sections for photons
  template <>
  auto cross_builder<PROPOSAL::GammaDef> =
      [](PROPOSAL::Medium& medium,
         corsika::units::si::HEPEnergyType) { 
        auto particle_def = PROPOSAL::GammaDef();
        CORSIKA_LOG_DEBUG(
            "PROPOSAL: Generate photon cross sections for particle type {} ",
            particle_def.name);

        auto interpolate = true;
        auto cross_vec = std::vector<std::shared_ptr<PROPOSAL::CrossSectionBase>>();
        auto photopair =
            PROPOSAL::make_crosssection(PROPOSAL::crosssection::PhotoPairKochMotz{false},
                                        particle_def, medium, nullptr, interpolate);
        cross_vec.push_back(std::move(photopair));
        auto compton =
            PROPOSAL::make_crosssection(PROPOSAL::crosssection::ComptonKleinNishina{},
                                        particle_def, medium, nullptr, interpolate);
        cross_vec.push_back(std::move(compton));
        auto photoproduction = PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::PhotoproductionHeckC7Shadowing{}, particle_def,
            medium, nullptr, interpolate);
        cross_vec.push_back(std::move(photoproduction));
        auto photoeffect =
            PROPOSAL::make_crosssection(PROPOSAL::crosssection::PhotoeffectSauter{},
                                        particle_def, medium, nullptr, interpolate);
        cross_vec.push_back(std::move(photoeffect));
        return cross_vec;
      };

  // cross sections for electrons
  template <>
  auto cross_builder<PROPOSAL::EMinusDef> =
      [](PROPOSAL::Medium& medium,
         corsika::units::si::HEPEnergyType
             emCut) { 
        auto particle_def = PROPOSAL::EMinusDef();
        CORSIKA_LOG_DEBUG(
            "PROPOSAL: Generate electron cross sections for particle type {} ",
            particle_def.name);
        auto interpolate = true;

        auto p_cut = std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV,
                                                                         v_cut, false);

        // do not use v_cut for ionization due to numerical problems
        auto p_cut_no_vcut =
            std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, false);

        auto cross_vec = std::vector<std::shared_ptr<PROPOSAL::CrossSectionBase>>();
        cross_vec.push_back(PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::BremsElectronScreening{false}, particle_def, medium,
            p_cut, interpolate));
        cross_vec.push_back(PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::EpairForElectronPositron{false}, particle_def, medium,
            p_cut, interpolate));
        cross_vec.push_back(PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::IonizBergerSeltzerMoller{*p_cut_no_vcut},
            particle_def, medium, p_cut_no_vcut, interpolate));
        cross_vec.push_back(PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::PhotoAbramowiczLevinLevyMaor97{
                std::make_unique<PROPOSAL::crosssection::ShadowButkevichMikheyev>()},
            particle_def, medium, p_cut, interpolate));

        return cross_vec;
      };

  // cross sections for positrons
  template <>
  auto cross_builder<PROPOSAL::EPlusDef> =
      [](PROPOSAL::Medium& medium,
         corsika::units::si::HEPEnergyType
             emCut) { 
        auto particle_def = PROPOSAL::EPlusDef();
        CORSIKA_LOG_DEBUG(
            "PROPOSAL: Generate positron cross sections for particle type {} ",
            particle_def.name);
        auto interpolate = true;

        auto p_cut = std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV,
                                                                         v_cut, false);

        // do not use v_cut for ionization due to numerical problems
        auto p_cut_no_vcut =
            std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, false);

        auto cross_vec = std::vector<std::shared_ptr<PROPOSAL::CrossSectionBase>>();
        cross_vec.push_back(PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::BremsElectronScreening{false}, particle_def, medium,
            p_cut, interpolate));
        cross_vec.push_back(PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::EpairForElectronPositron{false}, particle_def, medium,
            p_cut, interpolate));
        cross_vec.push_back(PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::IonizBergerSeltzerBhabha{*p_cut_no_vcut},
            particle_def, medium, p_cut_no_vcut, interpolate));
        cross_vec.push_back(PROPOSAL::make_crosssection(
            PROPOSAL::crosssection::PhotoAbramowiczLevinLevyMaor97{
                std::make_unique<PROPOSAL::crosssection::ShadowButkevichMikheyev>()},
            particle_def, medium, p_cut, interpolate));
        cross_vec.push_back(
            PROPOSAL::make_crosssection(PROPOSAL::crosssection::AnnihilationHeitler{},
                                        particle_def, medium, nullptr, interpolate));

        return cross_vec;
      };

  static std::map<Code, std::function<PROPOSAL::crosssection_list_t(
                            PROPOSAL::Medium&, corsika::units::si::HEPEnergyType)>>
      cross = {{Code::Photon, cross_builder<PROPOSAL::GammaDef>},
               {Code::Electron, cross_builder<PROPOSAL::EMinusDef>},
               {Code::Positron, cross_builder<PROPOSAL::EPlusDef>},
               {Code::MuMinus, cross_builder<PROPOSAL::MuMinusDef>},
               {Code::MuPlus, cross_builder<PROPOSAL::MuPlusDef>},
               {Code::TauMinus, cross_builder<PROPOSAL::TauMinusDef>},
               {Code::TauPlus, cross_builder<PROPOSAL::TauPlusDef>}};

  class ProposalProcessBase {
  protected:
    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("proposal");

    std::unordered_map<std::size_t, PROPOSAL::Medium>
        media; 

    template <typename TEnvironment>
    ProposalProcessBase(TEnvironment const& _env);

    bool canInteract(Code pcode) const;

    using calc_key_t = std::pair<std::size_t, Code>;

    struct hash {
      size_t operator()(const calc_key_t& p) const noexcept;
    };

    virtual void buildCalculator(Code, NuclearComposition const&) = 0;

    template <typename Particle, typename Calculators>
    auto getCalculator(Particle& vP, Calculators& calc) {
      const auto& comp = vP.getNode()->getModelProperties().getNuclearComposition();
      auto calc_it = calc.find(std::make_pair(comp.getHash(), vP.getPID()));
      if (calc_it != calc.end()) return calc_it;
      buildCalculator(vP.getPID(), comp);
      return getCalculator(vP, calc);
    }
  };
} // namespace corsika::proposal

#include <corsika/detail/modules/proposal/ProposalProcessBase.inl>