crosssections.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2013-2014,2018-2023
 *      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 "scatteractionsfinderparameters.h"
#include "stringprocess.h"
 
namespace smash {
 
class CrossSections {
 public:
  CrossSections(const ParticleList& incoming_particles, double sqrt_s,
                const std::pair<FourVector, FourVector> potentials);
 
  CollisionBranchList generate_collision_list(
      const ScatterActionsFinderParameters& finder_parameters,
      StringProcess* string_process) const;
 
  double parametrized_total(
      const ScatterActionsFinderParameters& finder_parameters) 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(
      const ScatterActionsFinderParameters& finder_parameters) 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(const ReactionsBitSet& included_2to2,
                                 double KN_offset) 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(double current_xs,
                                        double scale_xs) const;
 
  CollisionBranchList NNbar_creation() const;
 
  CollisionBranchPtr NNbar_to_5pi(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 StringTransitionParameters& transition_high_energy) const;
 
  double string_probability(
      const ScatterActionsFinderParameters& finder_parameters) const;
 
  double probability_transit_high(double region_lower,
                                  double region_upper) const;
 
 private:
  double elastic_parametrization(bool use_AQM, double pipi_offset) const;
 
  double nn_el() const;
 
  double npi_el() const;
 
  double nk_el() const;
 
  CollisionBranchList npi_yk() const;
 
  CollisionBranchList bb_xx_except_nn(
      const ReactionsBitSet& included_2to2) const;
 
  CollisionBranchList nn_xx(const ReactionsBitSet& included_2to2) const;
 
  CollisionBranchList nk_xx(const ReactionsBitSet& included_2to2,
                            double KN_offset) const;
 
  CollisionBranchList deltak_xx(const ReactionsBitSet& included_2to2) const;
 
  CollisionBranchList ypi_xx(const ReactionsBitSet& included_2to2) const;
 
  CollisionBranchList dpi_xx(const ReactionsBitSet& included_2to2) const;
 
  CollisionBranchList dn_xx(const ReactionsBitSet& included_2to2) const;
 
  static double xs_dpi_dprimepi(double sqrts, double cm_mom,
                                ParticleTypePtr produced_nucleus,
                                const ParticleType& type_pi);
 
  static double xs_dn_dprimen(double sqrts, 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(std::make_unique<CollisionBranch>(
          type_a, type_b, xsection, ProcessType::TwoToTwo));
    }
  }
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_CROSSSECTIONS_H_