Interaction.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 <tuple>
#include <boost/filesystem/path.hpp>

#include <corsika/framework/utility/CorsikaData.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/process/InteractionProcess.hpp>
#include <corsika/modules/pythia8/Pythia8.hpp>

namespace corsika::pythia8 {
  class Interaction : public InteractionProcess<Interaction> {

  public:
    Interaction(
        boost::filesystem::path const& mpiInitFile = corsika_data("Pythia/main184.mpi"),
        bool const print_listing = false);
    ~Interaction();

    bool canInteract(Code const) const;

    bool isValid(Code const projectileId, Code const targetId,
                 HEPEnergyType const sqrtS) const;
    std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
        Code const projectile, Code const target, FourMomentum const& projectileP4,
        FourMomentum const& targetP4) const;

    CrossSectionType getCrossSection(
        Code const projectile, Code const target, FourMomentum const& projectileP4,
        FourMomentum const& targetP4) const; // non-const for now

    template <typename TView>
    void doInteraction(TView& output, Code const projectileId, Code const targetId,
                       FourMomentum const& projectileP4, FourMomentum const& targetP4);

    double getAverageSubcollisions(Code targetId, CrossSectionType sigTot) const;

    static std::array constexpr validTargets_{Code::Oxygen, Code::Nitrogen,
                                              Code::Argon,  Code::Hydrogen,
                                              Code::Proton, Code::Neutron};

  private:
    default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("pythia");
    //~ bool const internalDecays_ = true;
    int count_ = 0;
    bool const print_listing_ = false;
    Pythia8::Pythia pythiaMain_;
    Pythia8::Pythia mutable pythiaColl_; // Pythia's getSigma...() are not marked const...
    double const probSD_ =
        0.3; // Fraction of single diffractive events beyond first collision in nucleus.
    HEPEnergyType const eMaxLab_ = 1e18_eV;
    HEPEnergyType const eKinMinLab_ = 0.2_GeV;
  };

} // namespace corsika::pythia8

#include <corsika/detail/modules/pythia8/Interaction.inl>