crosssections.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2015-2022
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_CROSSSECTIONS_H_
#define SRC_INCLUDE_SMASH_CROSSSECTIONS_H_
 
#include <memory>
#include <utility>
 
#include "forwarddeclarations.h"
#include "isoparticletype.h"
#include "particles.h"
#include "potential_globals.h"
#include "stringprocess.h"
 
namespace smash {
 
namespace transit_high_energy {
const std::array<double, 2> sqrts_range_Npi = {1.9, 2.2};
const std::array<double, 2> sqrts_range_NN = {3.5, 4.5};
const double sqrts_add_lower = 0.9;
const double sqrts_range = 1.0;
 
const double pipi_offset = 1.12;
 
const double KN_offset = 15.15;
}  // namespace transit_high_energy
 
class CrossSections {
 public:
  CrossSections(const ParticleList& incoming_particles, const double sqrt_s,
                const std::pair<FourVector, FourVector> potentials);
 
  CollisionBranchList generate_collision_list(
      double elastic_parameter, bool two_to_one_switch,
      ReactionsBitSet included_2to2,
      MultiParticleReactionsBitSet included_multi, double low_snn_cut,
      bool strings_switch, bool use_AQM, bool strings_with_probability,
      NNbarTreatment nnbar_treatment, StringProcess* string_process,
      double scale_xs, double additional_el_xs) const;
 
  static double sum_xs_of(const CollisionBranchList& list) {
    double xs_sum = 0.0;
    for (auto& proc : list) {
      xs_sum += proc->weight();
    }
    return xs_sum;
  }
 
  CollisionBranchPtr elastic(double elast_par, bool use_AQM, double add_el_xs,
                             double scale_xs) const;
 
  CollisionBranchList two_to_one() const;
 
  double formation(const ParticleType& type_resonance,
                   double cm_momentum_sqr) const;
 
  CollisionBranchList rare_two_to_two() const;
 
  CollisionBranchList two_to_two(ReactionsBitSet included_2to2) const;
 
  CollisionBranchList two_to_three() const;
 
  CollisionBranchList two_to_four() const;
 
  CollisionBranchList string_excitation(double total_string_xs,
                                        StringProcess* string_process,
                                        bool use_AQM) const;
 
  CollisionBranchPtr NNbar_annihilation(const double current_xs,
                                        const double scale_xs) const;
 
  CollisionBranchList NNbar_creation() const;
 
  CollisionBranchPtr NNbar_to_5pi(const double scale_xs) const;
 
  static double d_pi_inelastic_xs(double pion_kinetic_energy);
 
  static double d_N_inelastic_xs(double N_kinetic_energy);
 
  static double d_aN_inelastic_xs(double aN_kinetic_energy);
 
  static double two_to_three_xs(const ParticleType& type_in1,
                                const ParticleType& type_in2, double sqrts);
 
  static double two_to_four_xs(const ParticleType& type_in1,
                               const ParticleType& type_in2, double sqrts);
 
  double high_energy() const;
 
  double string_probability(bool strings_switch,
                            bool use_transition_probability, bool use_AQM,
                            bool treat_nnbar_with_strings) const;
 
  double probability_transit_high(const double region_lower,
                                  const double region_upper) const;
 
 private:
  double elastic_parametrization(bool use_AQM) const;
 
  double nn_el() const;
 
  double npi_el() const;
 
  double nk_el() const;
 
  CollisionBranchList npi_yk() const;
 
  CollisionBranchList bb_xx_except_nn(ReactionsBitSet included_2to2) const;
 
  CollisionBranchList nn_xx(ReactionsBitSet included_2to2) const;
 
  CollisionBranchList nk_xx(ReactionsBitSet included_2to2) const;
 
  CollisionBranchList deltak_xx(ReactionsBitSet included_2to2) const;
 
  CollisionBranchList ypi_xx(ReactionsBitSet included_2to2) const;
 
  CollisionBranchList dpi_xx(ReactionsBitSet included_2to2) const;
 
  CollisionBranchList dn_xx(ReactionsBitSet included_2to2) const;
 
  static double xs_dpi_dprimepi(const double sqrts, const double cm_mom,
                                ParticleTypePtr produced_nucleus,
                                const ParticleType& type_pi);
 
  static double xs_dn_dprimen(const double sqrts, const double cm_mom,
                              ParticleTypePtr produced_nucleus,
                              const ParticleType& type_nucleus,
                              const ParticleType& type_N);
 
  double string_hard_cross_section() const;
 
  CollisionBranchList bar_bar_to_nuc_nuc(const bool is_anti_particles) const;
 
  static double nn_to_resonance_matrix_element(double sqrts,
                                               const ParticleType& type_a,
                                               const ParticleType& type_b,
                                               const int twoI);
 
  template <class IntegrationMethod>
  CollisionBranchList find_nn_xsection_from_type(
      const ParticleTypePtrList& type_res_1,
      const ParticleTypePtrList& type_res_2,
      const IntegrationMethod integrator) const;
 
  double cm_momentum() const {
    const double m1 = incoming_particles_[0].effective_mass();
    const double m2 = incoming_particles_[1].effective_mass();
    return pCM(sqrt_s_, m1, m2);
  }
 
  const ParticleList incoming_particles_;
 
  const double sqrt_s_;
 
  const std::pair<FourVector, FourVector> potentials_;
 
  const bool is_BBbar_pair_;
 
  const bool is_NNbar_pair_;
 
  template <typename F>
  void add_channel(CollisionBranchList& process_list, F&& get_xsection,
                   double sqrts, const ParticleType& type_a,
                   const ParticleType& type_b) const {
    const double sqrt_s_min =
        type_a.min_mass_spectral() + type_b.min_mass_spectral();
    /* Determine wether the process is below the threshold. */
    double scale_B = 0.0;
    double scale_I3 = 0.0;
    bool is_below_threshold;
    FourVector incoming_momentum = FourVector();
    if (pot_pointer != nullptr) {
      for (const auto& p : incoming_particles_) {
        incoming_momentum += p.momentum();
        scale_B += pot_pointer->force_scale(p.type()).first;
        scale_I3 +=
            pot_pointer->force_scale(p.type()).second * p.type().isospin3_rel();
      }
      scale_B -= pot_pointer->force_scale(type_a).first;
      scale_I3 -=
          pot_pointer->force_scale(type_a).second * type_a.isospin3_rel();
      scale_B -= pot_pointer->force_scale(type_b).first;
      scale_I3 -=
          pot_pointer->force_scale(type_b).second * type_b.isospin3_rel();
      is_below_threshold = (incoming_momentum + potentials_.first * scale_B +
                            potentials_.second * scale_I3)
                               .abs() <= sqrt_s_min;
    } else {
      is_below_threshold = (sqrts <= sqrt_s_min);
    }
    if (is_below_threshold) {
      return;
    }
    const auto xsection = get_xsection();
    if (xsection > really_small) {
      process_list.push_back(make_unique<CollisionBranch>(
          type_a, type_b, xsection, ProcessType::TwoToTwo));
    }
  }
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_CROSSSECTIONS_H_