smash::CrossSections Class Reference

#include <crosssections.h>

The CrossSections class assembles everything that is needed to calculate cross sections and gathers a list of all possible reactions for the incoming particles at the given energy with the calculated cross sections.

Definition at line 31 of file crosssections.h.

Collaboration diagram for smash::CrossSections:

Public Member Functions
	CrossSections (const ParticleList &incoming_particles, double sqrt_s, const std::pair< FourVector, FourVector > potentials)
	Construct CrossSections instance. More...
CollisionBranchList	generate_collision_list (const ScatterActionsFinderParameters &finder_parameters, StringProcess *string_process) const
	Generate a list of all possible collisions between the incoming particles with the given c.m. More...
double	parametrized_total (const ScatterActionsFinderParameters &finder_parameters) const
	Select the parametrization for the total cross section, given the types of incoming particles. More...
CollisionBranchPtr	elastic (const ScatterActionsFinderParameters &finder_parameters) const
	Determine the elastic cross section for this collision. More...
CollisionBranchList	two_to_one (CharmRescattering charm_rescattering) const
	Find all resonances that can be produced in a 2->1 collision of the two input particles and the production cross sections of these resonances. More...
double	formation (const ParticleType &type_resonance, double cm_momentum_sqr) const
	Calculates the 2-to-1 resonance production cross section for a given resonance using the Breit-Wigner distribution as probability amplitude. More...
CollisionBranchList	rare_two_to_two () const
	Find all 2->2 processes which are suppressed at high energies when strings are turned on with probabilites, but important for the production of rare species such as strange particles. More...
CollisionBranchList	two_to_two (const ReactionsBitSet &included_2to2, double KN_offset, CharmRescattering charm_rescattering) const
	Find all inelastic 2->2 processes for the given scattering. More...
CollisionBranchList	two_to_three () const
	Find all 2->3 processes for the given scattering. More...
CollisionBranchList	two_to_four () const
	Find all 2->4 processes for the given scattering. More...
CollisionBranchList	string_excitation (double total_string_xs, StringProcess *string_process, const ScatterActionsFinderParameters &finder_parameters) const
	Determine the cross section for string excitations, which is given by the difference between the parametrized total cross section and all the explicitly implemented channels at low energy (elastic, resonance excitation, etc). More...
CollisionBranchPtr	NNbar_annihilation (double current_xs, double scale_xs) const
	Determine the cross section for NNbar annihilation, which is given by the difference between the parametrized total cross section and all the explicitly implemented channels at low energy (in this case only elastic). More...
CollisionBranchList	NNbar_creation () const
	Determine the cross section for NNbar creation, which is given by detailed balance from the reverse reaction. More...
CollisionBranchPtr	NNbar_to_5pi (double scale_xs) const
	Create collision branch for NNbar annihilation going directly into 5 pions. More...
double	high_energy (const ScatterActionsFinderParameters &finder_parameters) const
	Determine the parametrized total cross section at high energies for the given collision, which is non-zero for Baryon-Baryon and Nucleon-Pion scatterings currently. More...
double	string_probability (const ScatterActionsFinderParameters &finder_parameters) const
double	transition_probability_at_sqrts (double region_lower, double region_upper) const
	Computes a smooth transition probability as a function of sqrt(s). More...

Static Public Member Functions
static double	sum_xs_of (const CollisionBranchList &list)
	Helper function: Sum all cross sections of the given process list. More...
static double	two_to_three_xs (const ParticleType &type_in1, const ParticleType &type_in2, double sqrts)
	Determine 2->3 cross section for the scattering of the given particle types. More...
static double	two_to_four_xs (const ParticleType &type_in1, const ParticleType &type_in2, double sqrts)
	Determine 2->4 cross section for the scattering of the given particle types. More...

Private Member Functions
double	elastic_parametrization (const ScatterActionsFinderParameters &finder_parameters) const
	Choose the appropriate parametrizations for given incoming particles and return the (parametrized) elastic cross section. More...
double	nn_el () const
	Determine the (parametrized) elastic cross section for a nucleon-nucleon (NN) collision. More...
double	npi_el () const
	Determine the elastic cross section for a nucleon-pion (Npi) collision. More...
double	nk_el () const
	Determine the elastic cross section for a nucleon-kaon (NK) collision. More...
std::optional< double >	Dpi_and_Dstarpi_elastic () const
	Determine the elastic cross section for a D meson-pion (Dpi) or a D*-pion (D*pi) collision. More...
double	Dpi_and_Dstarpi_inelastic () const
	Determine the inelastic cross section for a D meson-pion (Dpi) or a D*-pion (D*pi) collision. More...
std::optional< double >	Deta_and_Dstareta_elastic () const
	Determine the elastic cross section for a D meson-eta (Deta) or a D*-eta (D*eta) collision. More...
std::optional< double >	DK_and_DstarK_elastic () const
	Determine the elastic cross section for a D meson-kaon (DK) or a D*-kaon (D*K) collision. More...
double	DK_and_DstarK_inelastic () const
	Determine the inelastic cross section for a D meson-kaon (DK) or a D*-kaon (D*K) collision. More...
std::optional< double >	DN_elastic () const
	Determine the elastic cross section for a D meson-nucleon (DN) collision, If the center-of-mass energy for the collision is below the lower bound of the energy range of the underlying cross section data, the return value is zero. More...
double	DN_inelastic () const
	Determine the inelastic cross section for a D meson-nucleon (DN) collision. More...
std::optional< double >	DDelta_elastic () const
	Determine the elastic cross section for a D meson-Delta (DΔ) collision, If the center-of-mass energy for the collision is below the lower bound of the energy range of the underlying cross section data, the return value is zero. More...
double	DDelta_inelastic () const
	Determine the inelastic cross section for a D meson-Delta (DΔ) collision. More...
CollisionBranchList	npi_yk () const
	Find all processes for Nucleon-Pion to Hyperon-Kaon Scattering. More...
CollisionBranchList	bb_xx_except_nn (const ReactionsBitSet &included_2to2) const
	Find all inelastic 2->2 processes for Baryon-Baryon (BB) Scattering except the more specific Nucleon-Nucleon Scattering. More...
CollisionBranchList	nn_xx (const ReactionsBitSet &included_2to2) const
	Find all inelastic 2->2 processes for Nucelon-Nucelon Scattering. More...
CollisionBranchList	nk_xx (const ReactionsBitSet &included_2to2, double KN_offset) const
	Find all inelastic 2->2 background processes for Nucleon-Kaon (NK) Scattering. More...
CollisionBranchList	Dpi_and_Dstarpi_xx (const ReactionsBitSet &included_2to2, CharmRescattering charm_rescattering) const
	Find all inelastic 2->2 processes for D meson-pion (Dpi) and D*-pion (D*pi) scattering. More...
CollisionBranchList	DK_and_DstarK_xx (const ReactionsBitSet &included_2to2, CharmRescattering charm_rescattering) const
	Find all inelastic 2->2 processes for D meson-kaon (DK) and D*-kaon (D*K) scattering. More...
CollisionBranchList	DN_xx (const ReactionsBitSet &included_2to2, CharmRescattering charm_rescattering) const
	Find all inelastic 2->2 processes for D meson-nucleon (DN) scatterings. More...
CollisionBranchList	DDelta_xx (const ReactionsBitSet &included_2to2, CharmRescattering charm_rescattering) const
	Find all inelastic 2->2 processes for D meson-Delta (DΔ) scatterings. More...
CollisionBranchList	deltak_xx (const ReactionsBitSet &included_2to2) const
	Find all inelastic 2->2 processes for Delta-Kaon (DeltaK) Scattering. More...
CollisionBranchList	ypi_xx (const ReactionsBitSet &included_2to2) const
	Find all inelastic 2->2 processes for Hyperon-Pion (Ypi) Scattering. More...
CollisionBranchList	dpi_xx (const ReactionsBitSet &included_2to2) const
	Find all inelastic 2->2 processes involving Pion and (anti-) Deuteron (dpi), specifically dπ→ NN, d̅π→ N̅N̅; πd→ πd' (mockup for πd→ πnp), πd̅→ πd̅' and reverse. More...
CollisionBranchList	dn_xx (const ReactionsBitSet &included_2to2) const
	Find all inelastic 2->2 processes involving Nucleon and (anti-) Deuteron (dN), specifically Nd → Nd', N̅d → N̅d', N̅d̅→ N̅d̅', Nd̅→ Nd̅' and reverse (e.g. More...
double	string_hard_cross_section () const
	Determine the (parametrized) hard non-diffractive string cross section for this collision. More...
CollisionBranchList	bar_bar_to_nuc_nuc (bool is_anti_particles) const
	Calculate cross sections for 2 → 2 resonance absorption (i.e. More...
template<class IntegrationMethod >
CollisionBranchList	find_nn_xsection_from_type (const ParticleTypePtrList &type_res_1, const ParticleTypePtrList &type_res_2, const IntegrationMethod integrator) const
	Utility function to avoid code replication in nn_xx(). More...
double	cm_momentum () const
	Determine the momenta of the incoming particles in the center-of-mass system. More...
template<typename F >
void	add_channel (CollisionBranchList &process_list, F &&get_xsection, double sqrts, const ParticleType &type_a, const ParticleType &type_b) const
	Helper function: Add a 2-to-2 channel to a collision branch list given a cross section. More...

Static Private Member Functions
static double	xs_dpi_dprimepi (double sqrts, double cm_mom, ParticleTypePtr produced_nucleus, const ParticleType &type_pi)
	Parametrized cross section for πd→ πd' (mockup for πd→ πnp), πd̅→ πd̅' and reverse, see Oliinychenko:2018ugs [49] for details. More...
static double	xs_dn_dprimen (double sqrts, double cm_mom, ParticleTypePtr produced_nucleus, const ParticleType &type_nucleus, const ParticleType &type_N)
	Parametrized cross section for Nd → Nd', N̅d → N̅d', N̅d̅→ N̅d̅', Nd̅→ Nd̅' and reverse (e.g. More...
static double	nn_to_resonance_matrix_element (double sqrts, const ParticleType &type_a, const ParticleType &type_b, int twoI)
	Scattering matrix amplitude squared (divided by 16π) for resonance production processes like NN → NR and NN → ΔR, where R is a baryon resonance (Δ, N*, Δ*). More...

Private Attributes
const ParticleList	incoming_particles_
	List with data of scattering particles. More...
const double	sqrt_s_
	Total energy in the center-of-mass frame. More...
const std::pair< FourVector, FourVector >	potentials_
	Potentials at the interacting point. More...
const bool	is_BBbar_pair_
	Whether incoming particles are a pair of a baryon and an antibaryon (could be different baryon types) More...
const bool	is_NNbar_pair_
	Whether incoming particles are a nulecon-antinucleon pair (same isospin) More...

Constructor & Destructor Documentation

CrossSections()

smash::CrossSections::CrossSections	(	const ParticleList &	incoming_particles,
		double	sqrt_s,
		const std::pair< FourVector, FourVector >	potentials
	)

Construct CrossSections instance.

Parameters

[in]	incoming_particles	Particles that are interacting.
[in]	sqrt_s	Center-of-mass energy of the reaction.
[in]	potentials	Potentials at the interacting point. They are used to calculate the corrections on the thresholds.

Definition at line 219 of file crosssections.cc.

    : incoming_particles_(incoming_particles),
      sqrt_s_(sqrt_s),
      potentials_(potentials),
      is_BBbar_pair_(incoming_particles_[0].type().is_baryon() &&
                     incoming_particles_[1].type().is_baryon() &&
                     incoming_particles_[0].type().antiparticle_sign() ==
                         -incoming_particles_[1].type().antiparticle_sign()),
      is_NNbar_pair_(
          incoming_particles_[0].type().is_nucleon() &&
          incoming_particles_[1].pdgcode() ==
              incoming_particles_[0].type().get_antiparticle()->pdgcode()) {}

Member Function Documentation

generate_collision_list()

CollisionBranchList smash::CrossSections::generate_collision_list	(	const ScatterActionsFinderParameters &	finder_parameters,
		StringProcess *	string_process
	)		const

Generate a list of all possible collisions between the incoming particles with the given c.m.

energy and the calculated cross sections. The string processes are not added at this step if it is not triggerd according to the probability. It will then be added in add_all_scatterings in scatteraction.cc

Parameters

[in]	finder_parameters	parameters for collision finding.
[in]	string_process	a pointer to the StringProcess object, which is used for string excitation and fragmentation.

Returns

List of all possible collisions.

Definition at line 234 of file crosssections.cc.

                                         {
  CollisionBranchList process_list;
  const ParticleType& t1 = incoming_particles_[0].type();
  const ParticleType& t2 = incoming_particles_[1].type();
 
  double p_pythia = 0.;
  if (finder_parameters.strings_with_probability) {
    p_pythia = string_probability(finder_parameters);
  }
 
  /* Elastic collisions between two nucleons with sqrt_s below
   * low_snn_cut can not happen. */
  const bool reject_by_nucleon_elastic_cutoff =
      t1.is_nucleon() && t2.is_nucleon() &&
      t1.antiparticle_sign() == t2.antiparticle_sign() &&
      sqrt_s_ < finder_parameters.low_snn_cut;
  bool incl_elastic =
      finder_parameters.included_2to2[IncludedReactions::Elastic];
  if (incl_elastic && !reject_by_nucleon_elastic_cutoff) {
    process_list.emplace_back(elastic(finder_parameters));
  }
  if (incoming_particles_[0].is_core() != incoming_particles_[1].is_core()) {
    return process_list;
  }
  if (p_pythia > 0.) {
    /* String-excitation cross section =
     * Parametrized total cross - the contributions
     * from all other present channels. */
    const double sig_current = sum_xs_of(process_list);
    const double sig_string = std::max(
        0., finder_parameters.scale_xs * high_energy(finder_parameters) -
                sig_current);
    append_list(
        process_list,
        string_excitation(sig_string, string_process, finder_parameters),
        p_pythia);
    append_list(process_list, rare_two_to_two(),
                p_pythia * finder_parameters.scale_xs);
  }
  if (p_pythia < 1.) {
    if (finder_parameters.two_to_one) {
      // resonance formation (2->1)
      append_list(process_list,
                  two_to_one(finder_parameters.charm_rescattering),
                  (1. - p_pythia) * finder_parameters.scale_xs);
    }
    if (finder_parameters.included_2to2.any()) {
      // 2->2 (inelastic)
      append_list(process_list,
                  two_to_two(finder_parameters.included_2to2,
                             finder_parameters.transition_high_energy.KN_offset,
                             finder_parameters.charm_rescattering),
                  (1. - p_pythia) * finder_parameters.scale_xs);
    }
    if (finder_parameters
            .included_multi[IncludedMultiParticleReactions::Deuteron_3to2] ==
        1) {
      // 2->3 (deuterons only 2-to-3 reaction at the moment)
      append_list(process_list, two_to_three(),
                  (1. - p_pythia) * finder_parameters.scale_xs);
    }
    if (finder_parameters
            .included_multi[IncludedMultiParticleReactions::A3_Nuclei_4to2] ==
        1) {
      // 2->4
      append_list(process_list, two_to_four(),
                  (1. - p_pythia) * finder_parameters.scale_xs);
    }
  }
  if (finder_parameters.nnbar_treatment == NNbarTreatment::TwoToFive &&
      is_NNbar_pair_) {
    // NNbar directly to 5 pions (2-to-5)
    process_list.emplace_back(NNbar_to_5pi(finder_parameters.scale_xs));
  }
 
  /* NNbar annihilation thru NNbar → ρh₁(1170); combined with the decays
   * ρ → ππ and h₁(1170) → πρ, this gives a final state of 5 pions.
   * Only use in cases when detailed balance MUST happen, i.e. in a box! */
  if (finder_parameters.nnbar_treatment == NNbarTreatment::Resonances) {
    if (is_NNbar_pair_) {
      /* Has to be called after the other processes are already determined,
       * so that the sum of the cross sections includes all other processes. */
      process_list.emplace_back(NNbar_annihilation(sum_xs_of(process_list),
                                                   finder_parameters.scale_xs));
    } else {
      append_list(process_list, NNbar_creation(), finder_parameters.scale_xs);
    }
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

parametrized_total()

double smash::CrossSections::parametrized_total

(

const ScatterActionsFinderParameters &

finder_parameters

)

const

Select the parametrization for the total cross section, given the types of incoming particles.

Parameters

[in]

finder_parameters

Parameters for collision finding, containing cut for low energy NN interactions.

Returns

The appropriate total cross section value.

Definition at line 327 of file crosssections.cc.

                                                                   {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  double total_xs = 0.;
  if (pdg_a.is_baryon() && pdg_b.is_baryon() &&
      sqrt_s_ > finder_parameters.low_snn_cut) {
    if (pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
      // NN
      total_xs = (pdg_a == pdg_b) ? pp_total(sqrt_s_ * sqrt_s_)
                                  : np_total(sqrt_s_ * sqrt_s_);
    } else {
      // NNbar
      total_xs = ppbar_total(sqrt_s_ * sqrt_s_);
    }
    total_xs *= finder_parameters.AQM_scaling_factor(pdg_a) *
                finder_parameters.AQM_scaling_factor(pdg_b);
  } else if ((pdg_a.is_baryon() && pdg_b.is_meson()) ||
             (pdg_a.is_meson() && pdg_b.is_baryon())) {
    const PdgCode& meson = pdg_a.is_meson() ? pdg_a : pdg_b;
    const PdgCode& baryon = pdg_a.is_meson() ? pdg_b : pdg_a;
    if (meson.is_kaon() && baryon.is_nucleon()) {
      if ((meson.code() == pdg::K_p && baryon.code() == pdg::p) ||
          (meson.code() == pdg::K_z && baryon.code() == pdg::n) ||
          (meson.code() == pdg::K_m && baryon.code() == -pdg::p) ||
          (meson.code() == pdg::Kbar_z && baryon.code() == -pdg::n)) {
        // K⁺p, K⁰n, and anti-processes
        total_xs = kplusp_total(sqrt_s_ * sqrt_s_);
      } else if ((meson.code() == pdg::K_p && baryon.code() == -pdg::p) ||
                 (meson.code() == pdg::K_z && baryon.code() == -pdg::n) ||
                 (meson.code() == pdg::K_m && baryon.code() == pdg::p) ||
                 (meson.code() == pdg::Kbar_z && baryon.code() == pdg::n)) {
        // K⁻p, K̅⁰n, and anti-processes
        total_xs = kminusp_total(sqrt_s_ * sqrt_s_);
      } else if ((meson.code() == pdg::K_p && baryon.code() == pdg::n) ||
                 (meson.code() == pdg::K_z && baryon.code() == pdg::p) ||
                 (meson.code() == pdg::K_m && baryon.code() == -pdg::n) ||
                 (meson.code() == pdg::Kbar_z && baryon.code() == -pdg::p)) {
        // K⁺n, K⁰p, and anti-processes
        total_xs = kplusn_total(sqrt_s_ * sqrt_s_);
      } else if ((meson.code() == pdg::K_p && baryon.code() == -pdg::n) ||
                 (meson.code() == pdg::K_z && baryon.code() == -pdg::p) ||
                 (meson.code() == pdg::K_m && baryon.code() == pdg::n) ||
                 (meson.code() == pdg::Kbar_z && baryon.code() == pdg::p)) {
        // K⁻n, K̅⁰p and anti-processes
        total_xs = kminusn_total(sqrt_s_ * sqrt_s_);
      }
    } else if (meson.is_pion() && baryon.is_nucleon()) {
      // π⁺(p,nbar), π⁻(n,pbar)
      if ((meson.code() == pdg::pi_p &&
           (baryon.code() == pdg::p || baryon.code() == -pdg::n)) ||
          (meson.code() == pdg::pi_m &&
           (baryon.code() == pdg::n || baryon.code() == -pdg::p))) {
        total_xs = piplusp_total(sqrt_s_);
      } else if (meson.code() == pdg::pi_z) {
        // π⁰N
        total_xs = 0.5 * (piplusp_total(sqrt_s_) + piminusp_total(sqrt_s_));
      } else {
        // π⁻(p,nbar), π⁺(n,pbar)
        total_xs = piminusp_total(sqrt_s_);
      }
    } else if (meson.is_Dmeson()) {
      const CharmRescattering& charm_rescattering =
          finder_parameters.charm_rescattering;
      std::optional<double> elastic_xs = std::nullopt;
      double inelastic_xs = 0.;
      if (baryon.is_nucleon()) {
        elastic_xs = DN_elastic();
        inelastic_xs = DN_inelastic();
      } else if (baryon.is_Delta()) {
        elastic_xs = DDelta_elastic();
        inelastic_xs = DDelta_inelastic();
      }
      if ((charm_rescattering == CharmRescattering::T_Matrix) &&
          elastic_xs.has_value()) {
        total_xs = elastic_xs.value() + inelastic_xs;
      } else {
        /* use AQM if charm_rescattering == CharmRescattering::Resonances or if
         * tmp_elastic_xs has no value, which happens either when sqrts is above
         * the upper bound of the energy range of the underlying cross section
         * data or there is no underlying data for the two colliding particles
         */
        total_xs = AQM_based_on_piminusp_high_energy(
            sqrt_s_, pdg_a, pdg_b, finder_parameters.AQM_scaling_factor(pdg_a),
            finder_parameters.AQM_scaling_factor(pdg_b));
      }
    } else {
      // M*+B* goes to AQM high energy π⁻p
      total_xs = AQM_based_on_piminusp_high_energy(
          sqrt_s_, pdg_a, pdg_b, finder_parameters.AQM_scaling_factor(pdg_a),
          finder_parameters.AQM_scaling_factor(pdg_b));
    }
  } else if (pdg_a.is_meson() && pdg_b.is_meson()) {
    if (pdg_a.is_pion() && pdg_b.is_pion()) {
      switch (pdg_a.isospin3() * pdg_b.isospin3() / 4) {
        // π⁺π⁻
        case -1:
          total_xs = pipluspiminus_total(sqrt_s_);
          break;
        case 0:
          // π⁰π⁰
          if (pdg_a.isospin3() + pdg_b.isospin3() == 0) {
            total_xs = pizeropizero_total(sqrt_s_);
          } else {
            // π⁺π⁰: similar to π⁺π⁻
            total_xs = pipluspiminus_total(sqrt_s_);
          }
          break;
        // π⁺π⁺ goes to π⁻p AQM
        case 1:
          total_xs = (2. / 3.) * piminusp_high_energy(sqrt_s_ * sqrt_s_);
          break;
        default:
          throw std::runtime_error("wrong isospin in ππ scattering");
      }
    } else if ((pdg_a.is_Dmeson() || pdg_b.is_Dmeson()) ||
               (pdg_a.is_Dstar2007() || pdg_b.is_Dstar2007())) {
      const CharmRescattering& charm_rescattering =
          finder_parameters.charm_rescattering;
      std::optional<double> elastic_xs = std::nullopt;
      double inelastic_xs = 0.;
      if (pdg_a.is_pion() || pdg_b.is_pion()) {
        elastic_xs = Dpi_and_Dstarpi_elastic();
        inelastic_xs = Dpi_and_Dstarpi_inelastic();
      } else if (pdg_a.is_eta() || pdg_b.is_eta()) {
        elastic_xs = Deta_and_Dstareta_elastic();
        // no inelastic scattering for Deta or D*eta
      } else if (pdg_a.is_kaon() || pdg_b.is_kaon()) {
        elastic_xs = DK_and_DstarK_elastic();
        inelastic_xs = DK_and_DstarK_inelastic();
      }
      if ((charm_rescattering == CharmRescattering::T_Matrix) &&
          elastic_xs.has_value()) {
        total_xs = elastic_xs.value() + inelastic_xs;
      } else {
        /* use AQM if charm_rescattering == CharmRescattering::Resonances or if
         * tmp_elastic_xs has no value, which happens either when sqrts is above
         * the upper bound of the energy range of the underlying cross section
         * data or there is no underlying data for the two colliding particles
         */
        total_xs = AQM_based_on_piminusp_high_energy(
            sqrt_s_, pdg_a, pdg_b, finder_parameters.AQM_scaling_factor(pdg_a),
            finder_parameters.AQM_scaling_factor(pdg_b));
      }
    } else {
      // M*+M* goes to AQM high energy π⁻p
      total_xs = AQM_based_on_piminusp_high_energy(
          sqrt_s_, pdg_a, pdg_b, finder_parameters.AQM_scaling_factor(pdg_a),
          finder_parameters.AQM_scaling_factor(pdg_b));
    }
  }
  return (total_xs + finder_parameters.additional_el_xs) *
         finder_parameters.scale_xs;
}

Here is the call graph for this function:

Here is the caller graph for this function:

sum_xs_of()

static double smash::CrossSections::sum_xs_of ( const CollisionBranchList &  list )

inlinestatic

Helper function: Sum all cross sections of the given process list.

Definition at line 77 of file crosssections.h.

                                                           {
    double xs_sum = 0.0;
    for (auto& proc : list) {
      xs_sum += proc->weight();
    }
    return xs_sum;
  }

Here is the caller graph for this function:

elastic()

CollisionBranchPtr smash::CrossSections::elastic

(

const ScatterActionsFinderParameters &

finder_parameters

)

const

Determine the elastic cross section for this collision.

If elastic_par is given (and positive), we just use a constant cross section of that size, otherwise a parametrization of the elastic cross section is used (if available). Optional a constant additional elastic cross section is added

Parameters

[in]

finder_parameters

parameters for collision finding, including cross section modifications from config file.

Note

The additional constant elastic cross section contribution is added after the scaling of the cross section.

Returns

A ProcessBranch object containing the cross section and final-state IDs.

Definition at line 482 of file crosssections.cc.

                                                                   {
  double elastic_xs = 0.;
 
  if (finder_parameters.elastic_parameter >= 0.) {
    // use constant elastic cross section from config file
    elastic_xs = finder_parameters.elastic_parameter;
  } else {
    // use parametrization
    elastic_xs = elastic_parametrization(finder_parameters);
  }
  /* when using a factor to scale the cross section and an additional
   * contribution to the elastic cross section, the contribution is added first
   * and then everything is scaled */
  return std::make_unique<CollisionBranch>(
      incoming_particles_[0].type(), incoming_particles_[1].type(),
      (elastic_xs + finder_parameters.additional_el_xs) *
          finder_parameters.scale_xs,
      ProcessType::Elastic);
}

Here is the call graph for this function:

Here is the caller graph for this function:

two_to_one()

CollisionBranchList smash::CrossSections::two_to_one

(

CharmRescattering

charm_rescattering

)

const

Find all resonances that can be produced in a 2->1 collision of the two input particles and the production cross sections of these resonances.

Given the data and type information of two colliding particles, create a list of possible resonance production processes and their cross sections.

If Charm_Rescattering_Method is not set to "Resonances" then interactions of charmed hadrons will not be considered in two to one processes.

Parameters

[in]

charm_rescattering

Type of charm rescattering

Returns

A list of processes with resonance in the final state. Each element in the list contains the type of the final-state particle and the cross section for that particular process.

Definition at line 1476 of file crosssections.cc.

                                                      {
  CollisionBranchList resonance_process_list;
  const ParticleType& type_particle_a = incoming_particles_[0].type();
  const ParticleType& type_particle_b = incoming_particles_[1].type();
 
  if (charm_rescattering == CharmRescattering::T_Matrix) {
    const PdgCode& pdg_a = type_particle_a.pdgcode();
    const PdgCode& pdg_b = type_particle_b.pdgcode();
    const bool Dmeson_present = pdg_a.is_Dmeson() || pdg_b.is_Dmeson();
    const bool Dstar_present = pdg_a.is_Dstar2007() || pdg_b.is_Dstar2007();
    const bool light_meson_present = pdg_a.is_pion() || pdg_b.is_pion() ||
                                     pdg_a.is_eta() || pdg_b.is_eta() ||
                                     pdg_a.is_kaon() || pdg_b.is_kaon();
    const bool nucleon_or_Delta_present = pdg_a.is_nucleon() ||
                                          pdg_b.is_nucleon() ||
                                          pdg_a.is_Delta() || pdg_b.is_Delta();
    if ((Dmeson_present && (light_meson_present || nucleon_or_Delta_present)) ||
        (Dstar_present && light_meson_present)) {
      return resonance_process_list;
    }
  }
 
  const double m1 = incoming_particles_[0].effective_mass();
  const double m2 = incoming_particles_[1].effective_mass();
  const double p_cm_sqr = pCM_sqr(sqrt_s_, m1, m2);
 
  ParticleTypePtrList possible_resonances =
      list_possible_resonances(&type_particle_a, &type_particle_b);
 
  // Find all the possible resonances
  for (const ParticleTypePtr type_resonance : possible_resonances) {
    double resonance_xsection = formation(*type_resonance, p_cm_sqr);
 
    // If cross section is non-negligible, add resonance to the list
    if (resonance_xsection > really_small) {
      resonance_process_list.push_back(std::make_unique<CollisionBranch>(
          *type_resonance, resonance_xsection, ProcessType::TwoToOne));
      logg[LCrossSections].debug("Found resonance: ", *type_resonance);
      logg[LCrossSections].debug(type_particle_a.name(), type_particle_b.name(),
                                 "->", type_resonance->name(),
                                 " at sqrt(s)[GeV] = ", sqrt_s_,
                                 " with xs[mb] = ", resonance_xsection);
    }
  }
  return resonance_process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

formation()

double smash::CrossSections::formation	(	const ParticleType &	type_resonance,
		double	cm_momentum_sqr
	)		const

Calculates the 2-to-1 resonance production cross section for a given resonance using the Breit-Wigner distribution as probability amplitude.

See eq. (176) in Buss:2011mx [16].

Parameters

[in]	type_resonance	Type information for the resonance to be produced.
[in]	cm_momentum_sqr	Square of the center-of-mass momentum of the two initial particles.

Returns

The cross section for the process [initial particle a] + [initial particle b] -> resonance.

Definition at line 1524 of file crosssections.cc.

                                                              {
  const ParticleType& type_particle_a = incoming_particles_[0].type();
  const ParticleType& type_particle_b = incoming_particles_[1].type();
 
  // Calculate partial in-width.
  const double partial_width = type_resonance.get_partial_in_width(
      sqrt_s_, incoming_particles_[0], incoming_particles_[1]);
  if (partial_width <= 0.) {
    return 0.;
  }
 
  assert(type_resonance.charge() ==
         type_particle_a.charge() + type_particle_b.charge());
  assert(type_resonance.baryon_number() ==
         type_particle_a.baryon_number() + type_particle_b.baryon_number());
 
  const double spinfactor =
      static_cast<double>(type_resonance.spin() + 1) /
      ((type_particle_a.spin() + 1) * (type_particle_b.spin() + 1));
  const int sym_factor =
      (type_particle_a.pdgcode() == type_particle_b.pdgcode()) ? 2 : 1;
  return spinfactor * sym_factor * 2. * M_PI * M_PI / cm_momentum_sqr *
         type_resonance.full_spectral_function(sqrt_s_) * partial_width *
         hbarc * hbarc / fm2_mb;
}

Here is the call graph for this function:

Here is the caller graph for this function:

rare_two_to_two()

CollisionBranchList smash::CrossSections::rare_two_to_two

(

)

const

Find all 2->2 processes which are suppressed at high energies when strings are turned on with probabilites, but important for the production of rare species such as strange particles.

This function should call the different, more specific functions for the different scatterings. But so far, only Nucleon-Pion to Hyperon- Kaon scattering is implemented.

Returns

List of all possible rare 2->2 processes.

Definition at line 503 of file crosssections.cc.

                                                         {
  CollisionBranchList process_list;
  const ParticleData& data_a = incoming_particles_[0];
  const ParticleData& data_b = incoming_particles_[1];
  const auto& pdg_a = data_a.pdgcode();
  const auto& pdg_b = data_b.pdgcode();
  if ((pdg_a.is_nucleon() && pdg_b.is_pion()) ||
      (pdg_b.is_nucleon() && pdg_a.is_pion())) {
    process_list = npi_yk();
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

two_to_two()

CollisionBranchList smash::CrossSections::two_to_two	(	const ReactionsBitSet &	included_2to2,
		double	KN_offset,
		CharmRescattering	charm_rescattering
	)		const

Find all inelastic 2->2 processes for the given scattering.

This function calls the different, more specific functions for the different scatterings.

Parameters

[in]	included_2to2	Which 2->2 reactions are enabled?
[in]	KN_offset	Offset to the minimum energy for string production in KN scatterings
[in]	charm_rescattering	Type of charm rescattering

Returns

List of all possible inelastic 2->2 processes.

Definition at line 1551 of file crosssections.cc.

                                                      {
  CollisionBranchList process_list;
  const ParticleData& data_a = incoming_particles_[0];
  const ParticleData& data_b = incoming_particles_[1];
  const ParticleType& type_a = data_a.type();
  const ParticleType& type_b = data_b.type();
  const auto& pdg_a = data_a.pdgcode();
  const auto& pdg_b = data_b.pdgcode();
 
  if (data_a.is_baryon() && data_b.is_baryon()) {
    if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
        pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
      // Nucleon Nucleon Scattering
      process_list = nn_xx(included_2to2);
    } else {
      // Baryon Baryon Scattering
      process_list = bb_xx_except_nn(included_2to2);
    }
  } else if ((type_a.is_baryon() && type_b.is_meson()) ||
             (type_a.is_meson() && type_b.is_baryon())) {
    // Baryon Meson Scattering
    if ((pdg_a.is_nucleon() && pdg_b.is_kaon()) ||
        (pdg_b.is_nucleon() && pdg_a.is_kaon())) {
      // Nucleon Kaon Scattering
      process_list = nk_xx(included_2to2, KN_offset);
    } else if ((pdg_a.is_hyperon() && pdg_b.is_pion()) ||
               (pdg_b.is_hyperon() && pdg_a.is_pion())) {
      // Hyperon Pion Scattering
      process_list = ypi_xx(included_2to2);
    } else if ((pdg_a.is_Delta() && pdg_b.is_kaon()) ||
               (pdg_b.is_Delta() && pdg_a.is_kaon())) {
      // Delta Kaon Scattering
      process_list = deltak_xx(included_2to2);
    } else if ((pdg_a.is_nucleon() && pdg_b.is_Dmeson()) ||
               (pdg_b.is_nucleon() && pdg_a.is_Dmeson())) {
      // Nucleon D meson Scattering
      process_list = DN_xx(included_2to2, charm_rescattering);
    } else if ((pdg_a.is_Delta() && pdg_b.is_Dmeson()) ||
               (pdg_b.is_Delta() && pdg_a.is_Dmeson())) {
      // Delta D meson Scattering
      process_list = DDelta_xx(included_2to2, charm_rescattering);
    }
  } else if (type_a.is_meson() && type_b.is_meson()) {
    if ((pdg_a.is_Dmeson() || pdg_b.is_Dmeson()) ||
        (pdg_a.is_Dstar2007() || pdg_b.is_Dstar2007())) {
      if (pdg_a.is_pion() || pdg_b.is_pion()) {
        // D or D* - Pion Scattering
        process_list = Dpi_and_Dstarpi_xx(included_2to2, charm_rescattering);
      } else if (pdg_a.is_eta() || pdg_b.is_eta()) {
        // D or D* - Eta Scattering (inelastic) not existent in T-matrix method
        return process_list;
      } else if (pdg_a.is_kaon() || pdg_b.is_kaon()) {
        // D or D* - Kaon Scattering
        process_list = DK_and_DstarK_xx(included_2to2, charm_rescattering);
      }
    }
  } else if (type_a.is_nucleus() || type_b.is_nucleus()) {
    if ((type_a.is_nucleon() && type_b.is_nucleus()) ||
        (type_b.is_nucleon() && type_a.is_nucleus())) {
      // Nucleon Deuteron and Nucleon d' Scattering
      process_list = dn_xx(included_2to2);
    } else if (((type_a.is_deuteron() || type_a.is_dprime()) &&
                pdg_b.is_pion()) ||
               ((type_b.is_deuteron() || type_b.is_dprime()) &&
                pdg_a.is_pion())) {
      // Pion Deuteron and Pion d' Scattering
      process_list = dpi_xx(included_2to2);
    }
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

two_to_three()

CollisionBranchList smash::CrossSections::two_to_three

(

)

const

Find all 2->3 processes for the given scattering.

This function calls the different, more specific functions for the different scatterings.

Returns

List of all possible 2->3 processes.

Definition at line 1625 of file crosssections.cc.

                                                      {
  CollisionBranchList process_list;
  const ParticleType& type_a = incoming_particles_[0].type();
  const ParticleType& type_b = incoming_particles_[1].type();
 
  if ((type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
      (type_b.is_deuteron() && type_a.pdgcode().is_pion())) {
    const ParticleType& type_pi = type_a.pdgcode().is_pion() ? type_a : type_b;
    const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
 
    if (type_nucleus.baryon_number() > 0) {
      // πd → πpn
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_n = ParticleType::find(pdg::n);
 
      process_list.push_back(std::make_unique<CollisionBranch>(
          type_pi, type_p, type_n, two_to_three_xs(type_a, type_b, sqrt_s_),
          ProcessType::TwoToThree));
    } else {
      // πd̅ → πp̅n̅
      const auto& type_anti_p = ParticleType::find(-pdg::p);
      const auto& type_anti_n = ParticleType::find(-pdg::n);
 
      process_list.push_back(std::make_unique<CollisionBranch>(
          type_pi, type_anti_p, type_anti_n,
          two_to_three_xs(type_a, type_b, sqrt_s_), ProcessType::TwoToThree));
    }
  }
 
  if ((type_a.is_nucleon() && type_b.is_deuteron()) ||
      (type_b.is_nucleon() && type_a.is_deuteron())) {
    const ParticleType& type_N = type_a.is_nucleon() ? type_a : type_b;
    const ParticleType& type_nucleus = type_a.is_deuteron() ? type_a : type_b;
 
    if (type_nucleus.baryon_number() > 0) {
      // Nd → Nnp, N̅d → N̅np
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_n = ParticleType::find(pdg::n);
 
      process_list.push_back(std::make_unique<CollisionBranch>(
          type_N, type_p, type_n, two_to_three_xs(type_a, type_b, sqrt_s_),
          ProcessType::TwoToThree));
    } else {
      // Nd̅ → Np̅n̅, N̅d̅ → N̅p̅n̅
      const auto& type_anti_p = ParticleType::find(-pdg::p);
      const auto& type_anti_n = ParticleType::find(-pdg::n);
 
      process_list.push_back(std::make_unique<CollisionBranch>(
          type_N, type_anti_p, type_anti_n,
          two_to_three_xs(type_a, type_b, sqrt_s_), ProcessType::TwoToThree));
    }
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

two_to_four()

CollisionBranchList smash::CrossSections::two_to_four

(

)

const

Find all 2->4 processes for the given scattering.

This function calls the different, more specific functions for the different scatterings.

Returns

List of all possible 2->4 processes.

Definition at line 1680 of file crosssections.cc.

                                                     {
  CollisionBranchList process_list;
  ParticleTypePtr type_nucleus = &(incoming_particles_[0].type());
  ParticleTypePtr type_catalyzer = &(incoming_particles_[1].type());
  if (!type_nucleus->is_nucleus()) {
    type_nucleus = &(incoming_particles_[1].type());
    type_catalyzer = &(incoming_particles_[0].type());
  }
 
  if (type_nucleus->is_nucleus() &&
      std::abs(type_nucleus->baryon_number()) == 3 &&
      (type_catalyzer->is_pion() || type_catalyzer->is_nucleon())) {
    const ParticleTypePtr type_p = ParticleType::try_find(pdg::p);
    const ParticleTypePtr type_n = ParticleType::try_find(pdg::n);
    const ParticleTypePtr type_anti_p = ParticleType::try_find(-pdg::p);
    const ParticleTypePtr type_anti_n = ParticleType::try_find(-pdg::n);
    const ParticleTypePtr type_la = ParticleType::try_find(pdg::Lambda);
    const ParticleTypePtr type_anti_la = ParticleType::try_find(-pdg::Lambda);
 
    // Find nucleus components
    ParticleTypePtrList components;
    components.reserve(3);
    const PdgCode nucleus_pdg = type_nucleus->pdgcode();
    for (int i = 0; i < nucleus_pdg.nucleus_p(); i++) {
      components.push_back(type_p);
    }
    for (int i = 0; i < nucleus_pdg.nucleus_n(); i++) {
      components.push_back(type_n);
    }
    for (int i = 0; i < nucleus_pdg.nucleus_ap(); i++) {
      components.push_back(type_anti_p);
    }
    for (int i = 0; i < nucleus_pdg.nucleus_an(); i++) {
      components.push_back(type_anti_n);
    }
    for (int i = 0; i < nucleus_pdg.nucleus_La(); i++) {
      components.push_back(type_la);
    }
    for (int i = 0; i < nucleus_pdg.nucleus_aLa(); i++) {
      components.push_back(type_anti_la);
    }
    if (sqrt_s_ > type_catalyzer->mass() + components[0]->mass() +
                      components[1]->mass() + components[2]->mass()) {
      process_list.push_back(std::make_unique<CollisionBranch>(
          *type_catalyzer, *(components[0]), *(components[1]), *(components[2]),
          two_to_four_xs(*type_nucleus, *type_catalyzer, sqrt_s_),
          ProcessType::TwoToFour));
    }
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

string_excitation()

CollisionBranchList smash::CrossSections::string_excitation	(	double	total_string_xs,
		StringProcess *	string_process,
		const ScatterActionsFinderParameters &	finder_parameters
	)		const

Determine the cross section for string excitations, which is given by the difference between the parametrized total cross section and all the explicitly implemented channels at low energy (elastic, resonance excitation, etc).

Parameters

[in]	total_string_xs	Total cross section for the string process [mb]
[in]	string_process	a pointer to the StringProcess object, which is used for string excitation and fragmentation
[in]	finder_parameters	parameters for collision finding and cross sections

Returns

List of subprocesses (single-diffractive, double-diffractive and non-diffractive) with their cross sections

Exceptions

std::runtime_error

if string_process is a null pointer

This method has to be called after all other processes have been determined.

Todo:

Same assumption made by NNbar_annihilation. Resolve.

Definition at line 4135 of file crosssections.cc.

                                                                   {
  if (!string_process) {
    throw std::runtime_error("string_process should be initialized.");
  }
 
  CollisionBranchList channel_list;
  if (total_string_xs <= 0.) {
    return channel_list;
  }
 
  double mandelstam_s = sqrt_s_ * sqrt_s_;
  /* Get mapped PDG id for evaluation of the parametrized cross sections for
   * diffractive processes. This must be rescaled according to additive quark
   * model in the case of exotic hadrons. Also calculate the multiplicative
   * factor for AQM based on the quark contents. */
  std::array<int, 2> pdgid;
  double AQM_scaling = 1.;
  for (int i = 0; i < 2; i++) {
    PdgCode pdg = incoming_particles_[i].type().pdgcode();
    pdgid[i] = StringProcess::pdg_map_for_pythia(pdg);
    AQM_scaling *= finder_parameters.AQM_scaling_factor(pdg);
  }
 
  /* Determine if the initial state is a baryon-antibaryon pair,
   * which can annihilate. */
  bool can_annihilate = false;
  if (is_BBbar_pair_) {
    int n_q_types = 5;  // u, d, s, c, b
    for (int iq = 1; iq <= n_q_types; iq++) {
      std::array<int, 2> nquark;
      for (int i = 0; i < 2; i++) {
        nquark[i] =
            incoming_particles_[i].type().pdgcode().net_quark_number(iq);
      }
      if (nquark[0] != 0 && nquark[1] != 0) {
        can_annihilate = true;
        break;
      }
    }
  }
 
  /* The case for baryon/anti-baryon annihilation is treated separately,
   * as in this case we use only one way to break up the particles, namely
   * into 2 mesonic strings of equal mass after annihilating one quark-
   * anti-quark pair. See StringProcess::next_BBbarAnn() */
  double sig_annihilation = 0.0;
  if (can_annihilate) {
    /* In the case of baryon-antibaryon pair, the parametrized cross section for
     * annihilation will be added. See xs_ppbar_annihilation(). */
    mandelstam_s = effective_AQM_s(
        mandelstam_s, incoming_particles_[0].effective_mass(),
        incoming_particles_[1].effective_mass(), nucleon_mass, nucleon_mass);
    double xs_param = xs_ppbar_annihilation(mandelstam_s);
    if (finder_parameters.use_AQM) {
      xs_param *= AQM_scaling;
    }
    sig_annihilation = std::min(total_string_xs, xs_param);
  }
 
  /* Total parametrized cross-section (I) and pythia-produced total
   * cross-section (II) do not necessarily coincide. If I > II then
   * non-diffractive cross-section is reinforced to get I == II.
   * If I < II then partial cross-sections are drained one-by-one
   * to reduce II until I == II:
   * first non-diffractive, then double-diffractive, then
   * single-diffractive AB->AX and AB->XB in equal proportion.
   * The way it is done here is not unique. I (ryu) think that at high energy
   * collision this is not an issue, but at sqrt_s < 10 GeV it may matter. */
  std::array<double, 3> xs_diffractive =
      string_process->cross_sections_diffractive(pdgid[0], pdgid[1],
                                                 std::sqrt(mandelstam_s));
 
  if (finder_parameters.use_AQM) {
    for (double& x : xs_diffractive) {
      x *= AQM_scaling;
    }
  }
 
  double single_diffr_AX = xs_diffractive[0];
  double single_diffr_XB = xs_diffractive[1];
  double double_diffr = xs_diffractive[2];
 
  double single_diffr = single_diffr_AX + single_diffr_XB;
  double diffractive = single_diffr + double_diffr;
 
  const double nondiffractive =
      std::max(0., total_string_xs - sig_annihilation - diffractive);
 
  diffractive = total_string_xs - sig_annihilation - nondiffractive;
  double_diffr = std::max(0., diffractive - single_diffr);
 
  const double a =
      single_diffr > 0.0 ? (diffractive - double_diffr) / single_diffr : 0.0;
 
  single_diffr_AX *= a;
  single_diffr_XB *= a;
 
  assert(std::abs(single_diffr_AX + single_diffr_XB + double_diffr +
                  sig_annihilation + nondiffractive - total_string_xs) < 1.e-6);
  enum class Proc { ND, SD_AX, SD_XB, DD, N };
  enum class Comp { Soft, Hard, N };
 
  constexpr std::size_t n_proc = static_cast<std::size_t>(Proc::N);
  constexpr std::size_t n_comp = static_cast<std::size_t>(Comp::N);
 
  std::array<std::array<double, n_comp>, n_proc> split_xs{};
 
  auto proc_idx = [](Proc p) { return static_cast<std::size_t>(p); };
 
  auto comp_idx = [](Comp c) { return static_cast<std::size_t>(c); };
 
  auto soft = [&](Proc p) -> double& {
    return split_xs[proc_idx(p)][comp_idx(Comp::Soft)];
  };
 
  auto hard = [&](Proc p) -> double& {
    return split_xs[proc_idx(p)][comp_idx(Comp::Hard)];
  };
 
  auto set_all_soft = [&](Proc p, double total) {
    soft(p) = total;
    hard(p) = 0.0;
  };
 
  auto split = [&](Proc p, double total, double weight_hard) {
    hard(p) = total * weight_hard;
    soft(p) = total - hard(p);
  };
  auto split_all_string_processes = [&](double weight_hard) {
    split(Proc::ND, nondiffractive, weight_hard);
    split(Proc::SD_AX, single_diffr_AX, weight_hard);
    split(Proc::SD_XB, single_diffr_XB, weight_hard);
    split(Proc::DD, double_diffr, weight_hard);
  };
 
  if (finder_parameters.hard_string_transition_mode ==
      HardStringTransitionMode::Custom_Range) {
    const auto& [hard_transition_start, hard_transition_end] =
        finder_parameters.hard_string_transition_energy_range;
 
    const double weight_hard = transition_probability_at_sqrts(
        hard_transition_start, hard_transition_end);
    split_all_string_processes(weight_hard);
  } else if (nondiffractive > 0.0) {
    const double hard_xsec = AQM_scaling * string_hard_cross_section();
 
    /* Use the non-diffractive exponential transition probability for all
 
    * string-excitation channels, so that soft strings are suppressed at high
    * energies also for diffractive processes. */
    const double weight_soft = std::exp(-hard_xsec / nondiffractive);
    const double weight_hard = std::clamp(1.0 - weight_soft, 0.0, 1.0);
    split_all_string_processes(weight_hard);
  } else {
    set_all_soft(Proc::ND, nondiffractive);
    set_all_soft(Proc::SD_AX, single_diffr_AX);
    set_all_soft(Proc::SD_XB, single_diffr_XB);
    set_all_soft(Proc::DD, double_diffr);
  }
  logg[LCrossSections].debug("Soft string cross sections [mb] are");
  logg[LCrossSections].debug("Soft single-diffractive AB->AX: ",
                             soft(Proc::SD_AX));
  logg[LCrossSections].debug("Soft single-diffractive AB->XB: ",
                             soft(Proc::SD_XB));
  logg[LCrossSections].debug("Soft double-diffractive AB->XX: ",
                             soft(Proc::DD));
  logg[LCrossSections].debug("Soft non-diffractive: ", soft(Proc::ND));
  logg[LCrossSections].debug("B-Bbar annihilation: ", sig_annihilation);
 
  logg[LCrossSections].debug("Hard string cross sections [mb] are");
  logg[LCrossSections].debug("Hard single-diffractive AB->AX: ",
                             hard(Proc::SD_AX));
  logg[LCrossSections].debug("Hard single-diffractive AB->XB: ",
                             hard(Proc::SD_XB));
  logg[LCrossSections].debug("Hard double-diffractive AB->XX: ",
                             hard(Proc::DD));
  logg[LCrossSections].debug("Hard non-diffractive: ", hard(Proc::ND));
 
  // cross section of soft string excitation including annihilation
  const double sig_string_soft = soft(Proc::SD_AX) + soft(Proc::SD_XB) +
                                 soft(Proc::DD) + soft(Proc::ND) +
                                 sig_annihilation;
 
  // fill the list of process channels
  if (sig_string_soft > 0.) {
    channel_list.push_back(std::make_unique<CollisionBranch>(
        soft(Proc::SD_AX), ProcessType::StringSoftSingleDiffractiveAX));
    channel_list.push_back(std::make_unique<CollisionBranch>(
        soft(Proc::SD_XB), ProcessType::StringSoftSingleDiffractiveXB));
    channel_list.push_back(std::make_unique<CollisionBranch>(
        soft(Proc::DD), ProcessType::StringSoftDoubleDiffractive));
    channel_list.push_back(std::make_unique<CollisionBranch>(
        soft(Proc::ND), ProcessType::StringSoftNonDiffractive));
 
    if (can_annihilate) {
      channel_list.push_back(std::make_unique<CollisionBranch>(
          sig_annihilation, ProcessType::StringSoftAnnihilation));
    }
  }
 
  if (hard(Proc::SD_AX) > 0.) {
    channel_list.push_back(std::make_unique<CollisionBranch>(
        hard(Proc::SD_AX), ProcessType::StringHardSingleDiffractiveAX));
  }
  if (hard(Proc::SD_XB) > 0.) {
    channel_list.push_back(std::make_unique<CollisionBranch>(
        hard(Proc::SD_XB), ProcessType::StringHardSingleDiffractiveXB));
  }
  if (hard(Proc::DD) > 0.) {
    channel_list.push_back(std::make_unique<CollisionBranch>(
        hard(Proc::DD), ProcessType::StringHardDoubleDiffractive));
  }
  if (hard(Proc::ND) > 0.) {
    channel_list.push_back(std::make_unique<CollisionBranch>(
        hard(Proc::ND), ProcessType::StringHardNonDiffractive));
  }
 
  return channel_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

NNbar_annihilation()

CollisionBranchPtr smash::CrossSections::NNbar_annihilation	(	double	current_xs,
		double	scale_xs
	)		const

Determine the cross section for NNbar annihilation, which is given by the difference between the parametrized total cross section and all the explicitly implemented channels at low energy (in this case only elastic).

Parameters

[in]	current_xs	Sum of all cross sections of already determined processes
[in]	scale_xs	Factor by which all (partial) cross sections are scaled

Returns

Collision Branch with NNbar annihilation process and its cross section

This method has to be called after all other processes have been determined.

Todo:

Same assumption made by string_excitation. Resolve.

Definition at line 4468 of file crosssections.cc.

                                                          {
  /* Calculate NNbar cross section:
   * Parametrized total minus all other present channels.*/
  const double s = sqrt_s_ * sqrt_s_;
  double nnbar_xsec = std::max(0., ppbar_total(s) * scale_xs - current_xs);
  logg[LCrossSections].debug("NNbar cross section is: ", nnbar_xsec);
  // Make collision channel NNbar -> ρh₁(1170); eventually decays into 5π
  return std::make_unique<CollisionBranch>(ParticleType::find(pdg::h1),
                                           ParticleType::find(pdg::rho_z),
                                           nnbar_xsec, ProcessType::TwoToTwo);
}

Here is the call graph for this function:

Here is the caller graph for this function:

NNbar_creation()

CollisionBranchList smash::CrossSections::NNbar_creation

(

)

const

Determine the cross section for NNbar creation, which is given by detailed balance from the reverse reaction.

See NNbar_annihilation.

Returns

Collision Branch with NNbar creation process and its cross section

Definition at line 4481 of file crosssections.cc.

                                                        {
  CollisionBranchList channel_list;
  const ParticleType& type_a = incoming_particles_[0].type();
  const ParticleType& type_b = incoming_particles_[1].type();
  if ((type_a.pdgcode() == pdg::rho_z && type_b.pdgcode() == pdg::h1) ||
      (type_a.pdgcode() == pdg::h1 && type_b.pdgcode() == pdg::rho_z)) {
    /* Calculate NNbar reverse cross section:
     * from reverse reaction (see NNbar_annihilation_cross_section).*/
    const double s = sqrt_s_ * sqrt_s_;
    const double pcm = cm_momentum();
 
    const auto& type_N = ParticleType::find(pdg::p);
    const auto& type_Nbar = ParticleType::find(-pdg::p);
 
    // Check available energy
    if (sqrt_s_ - 2 * type_N.mass() < 0) {
      return channel_list;
    }
 
    double xsection = detailed_balance_factor_RR(sqrt_s_, pcm, type_a, type_b,
                                                 type_N, type_Nbar) *
                      std::max(0., ppbar_total(s) - ppbar_elastic(s));
    logg[LCrossSections].debug("NNbar reverse cross section is: ", xsection);
    channel_list.push_back(std::make_unique<CollisionBranch>(
        type_N, type_Nbar, xsection, ProcessType::TwoToTwo));
    channel_list.push_back(std::make_unique<CollisionBranch>(
        ParticleType::find(pdg::n), ParticleType::find(-pdg::n), xsection,
        ProcessType::TwoToTwo));
  }
  return channel_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

NNbar_to_5pi()

CollisionBranchPtr smash::CrossSections::NNbar_to_5pi

(

double

scale_xs

)

const

Create collision branch for NNbar annihilation going directly into 5 pions.

The cross section is given by the parametrized ppbar cross section, which is also used for the reverse 5-to-2 process.

Parameters

[in]

scale_xs

Factor by which all (partial) cross sections are scaled

Returns

Collision Branch with NNbar annihilation process

Definition at line 4451 of file crosssections.cc.

                                                                          {
  const double s = sqrt_s_ * sqrt_s_;
  /* Use difference between total and elastic in order to conserve detailed
   * balance for all inelastoc NNbar processes. */
  const double nnbar_xsec = std::max(0., ppbar_total(s) - ppbar_elastic(s));
  logg[LCrossSections].debug("NNbar cross section for 2-to-5 is: ", nnbar_xsec);
 
  /* Make collision channel NNbar -> 5π (with same final state as resonance
   * approach). */
  const auto& type_piz = ParticleType::find(pdg::pi_z);
  const auto& type_pip = ParticleType::find(pdg::pi_p);
  const auto& type_pim = ParticleType::find(pdg::pi_m);
  return std::make_unique<CollisionBranch>(
      type_pip, type_pim, type_pip, type_pim, type_piz, nnbar_xsec * scale_xs,
      ProcessType::TwoToFive);
}

Here is the call graph for this function:

Here is the caller graph for this function:

two_to_three_xs()

double smash::CrossSections::two_to_three_xs ( const ParticleType &  type_in1, const ParticleType &  type_in2, double  sqrts  )

static

Determine 2->3 cross section for the scattering of the given particle types.

That the function only depends on the types of particles (plus sqrt(s)) and not on the specific particles, is an assumption needed in order to treat the 3->2 back-reaction with the stochastic criterion, where this function also needs to be called for 3-to-2 collision probability with only types and sqrt(s) known at this point. Therefore the function is also made static.

Parameters

[in]	type_in1	first scatterning particle type
[in]	type_in2	second scatterning particle type
[in]	sqrts	center-of-mass energy of scattering

Returns

cross section for 2->3 process

Definition at line 1732 of file crosssections.cc.

                                                    {
  double xs = 0.0;
  ParticleTypePtr type_nucleus = &type_a, type_catalyzer = &type_b;
  if (!type_nucleus->is_nucleus()) {
    type_nucleus = &type_b;
    type_catalyzer = &type_a;
  }
 
  bool nonzero_xs = type_nucleus->is_nucleus() &&
                    (type_catalyzer->is_pion() || type_catalyzer->is_nucleon());
  if (!nonzero_xs) {
    return 0.0;
  }
 
  const double md = type_nucleus->mass(), mcat = type_catalyzer->mass();
  const double Tkin = (sqrts * sqrts - (md + mcat) * (md + mcat)) / (2.0 * md);
 
  // Should normally never happen, but may be a useful safeguard
  if (Tkin <= 0.0) {
    return 0.0;
  }
 
  if (type_catalyzer->is_pion()) {
    xs = deuteron_pion_inelastic(Tkin);
  } else if (type_catalyzer->is_nucleon()) {
    if (type_nucleus->pdgcode().antiparticle_sign() ==
        type_catalyzer->pdgcode().antiparticle_sign()) {
      // Nd and N̅d̅
      xs = deuteron_nucleon_inelastic(Tkin);
    } else {
      // N̅d and Nd̅
      xs = deuteron_antinucleon_inelastic(Tkin);
    }
  }
  return xs;
}

Here is the call graph for this function:

Here is the caller graph for this function:

two_to_four_xs()

double smash::CrossSections::two_to_four_xs ( const ParticleType &  type_in1, const ParticleType &  type_in2, double  sqrts  )

static

Determine 2->4 cross section for the scattering of the given particle types.

Same assumptions as for 2->3 cross section, see respective documentation.

Parameters

[in]	type_in1	first scatterning particle type
[in]	type_in2	second scatterning particle type
[in]	sqrts	center-of-mass energy of scattering

Returns

cross section for 2->4 process

Definition at line 1771 of file crosssections.cc.

                                                                               {
  double xs = 0.0;
  ParticleTypePtr type_nucleus = &type_a, type_catalyzer = &type_b;
  if (!type_nucleus->is_nucleus()) {
    type_nucleus = &type_b;
    type_catalyzer = &type_a;
  }
  bool nonzero_xs = type_nucleus->is_nucleus() &&
                    (type_catalyzer->is_pion() || type_catalyzer->is_nucleon());
  if (!nonzero_xs) {
    return 0.0;
  }
 
  const double mA = type_nucleus->mass(), mcat = type_catalyzer->mass();
  const double Tkin = (sqrts * sqrts - (mA + mcat) * (mA + mcat)) / (2.0 * mA);
  const int A = type_nucleus->pdgcode().nucleus_A();
  // Should normally never happen, but may be a useful safeguard
  if (A != 3 || Tkin <= 0.0) {
    return 0.0;
  }
 
  if (type_catalyzer->is_pion()) {
    xs = A / 2. * deuteron_pion_inelastic(Tkin);
  } else if (type_catalyzer->is_nucleon()) {
    if (type_nucleus->pdgcode().antiparticle_sign() ==
        type_catalyzer->pdgcode().antiparticle_sign()) {
      // N + A, anti-N + anti-A
      xs = A / 2. * deuteron_nucleon_inelastic(Tkin);
    } else {
      // N̅ + A and N + anti-A
      xs = A / 2. * deuteron_antinucleon_inelastic(Tkin);
    }
  }
  return xs;
}

Here is the call graph for this function:

Here is the caller graph for this function:

high_energy()

double smash::CrossSections::high_energy

(

const ScatterActionsFinderParameters &

finder_parameters

)

const

Determine the parametrized total cross section at high energies for the given collision, which is non-zero for Baryon-Baryon and Nucleon-Pion scatterings currently.

This is rescaled by AQM factors.

Parameters

[in]

finder_parameters

parameters for collision finding and cross sections.

Definition at line 4356 of file crosssections.cc.

                                                                   {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
 
  const double s = sqrt_s_ * sqrt_s_;
  double xs = 0.;
 
  // Currently all BB collisions use the nucleon-nucleon parametrizations.
  if (pdg_a.is_baryon() && pdg_b.is_baryon()) {
    const double eff_s = effective_AQM_s(
        s, incoming_particles_[0].effective_mass(),
        incoming_particles_[1].effective_mass(), nucleon_mass, nucleon_mass);
    if (pdg_a == pdg_b) {
      xs = pp_high_energy(eff_s);  // pp, nn
    } else if (pdg_a.antiparticle_sign() * pdg_b.antiparticle_sign() == 1) {
      xs = np_high_energy(eff_s);  // np, nbarpbar
    } else if (pdg_a.antiparticle_sign() * pdg_b.antiparticle_sign() == -1) {
      /* In the case of baryon-antibaryon interactions,
       * the low-energy cross section must be involved
       * due to annihilation processes (via strings). */
      double xs_l = ppbar_total(eff_s);
      double xs_h = 0.;
      if (pdg_a.is_antiparticle_of(pdg_b)) {
        xs_h = ppbar_high_energy(eff_s);  // ppbar, nnbar
      } else {
        xs_h = npbar_high_energy(eff_s);  // npbar, nbarp
      }
      /* Transition between low and high energy is set to be consistent with
       * that defined in string_probability(). */
      auto [region_lower, region_upper] =
          finder_parameters.transition_high_energy.sqrts_range_NN;
      double prob_high =
          transition_probability_at_sqrts(region_lower, region_upper);
      xs = xs_l * (1. - prob_high) + xs_h * prob_high;
    }
  }
 
  // Pion nucleon interaction / baryon-meson
  if ((pdg_a == pdg::pi_p && pdg_b == pdg::p) ||
      (pdg_b == pdg::pi_p && pdg_a == pdg::p) ||
      (pdg_a == pdg::pi_m && pdg_b == pdg::n) ||
      (pdg_b == pdg::pi_m && pdg_a == pdg::n)) {
    xs = piplusp_high_energy(s);  // pi+ p, pi- n
  } else if ((pdg_a == pdg::pi_m && pdg_b == pdg::p) ||
             (pdg_b == pdg::pi_m && pdg_a == pdg::p) ||
             (pdg_a == pdg::pi_p && pdg_b == pdg::n) ||
             (pdg_b == pdg::pi_p && pdg_a == pdg::n)) {
    xs = piminusp_high_energy(s);  // pi- p, pi+ n
  } else if ((pdg_a.is_meson() && pdg_b.is_baryon()) ||
             (pdg_b.is_meson() && pdg_a.is_baryon())) {
    xs = piminusp_high_energy(s);  // default for baryon-meson
  }
 
  /* Meson-meson interaction goes through AQM from pi+p,
   * see user guide "Use_AQM" */
  if (pdg_a.is_meson() && pdg_b.is_meson()) {
    /* 2/3 factor since difference of 1 meson between meson-meson
     * and baryon-meson */
    xs = 2. / 3. * piplusp_high_energy(s);
  }
 
  // AQM scaling for cross-sections
  xs *= finder_parameters.AQM_scaling_factor(pdg_a) *
        finder_parameters.AQM_scaling_factor(pdg_b);
 
  return xs;
}

Here is the call graph for this function:

Here is the caller graph for this function:

string_probability()

double smash::CrossSections::string_probability

(

const ScatterActionsFinderParameters &

finder_parameters

)

const

Returns

the probability whether the scattering between the incoming particles is via string fragmentation or not.

If use_transition_probability is true: The string fragmentation is implemented in the same way in GiBUU (Physics Reports 512(2012), 1-124, pg. 33). If the center of mass energy is low, two particles scatter through the resonance channels. If high, the outgoing particles are generated by string fragmentation. If in between, the out- going particles are generated either through the resonance channels or string fragmentation by chance. In detail, the low energy region is from the threshold to (mix_scatter_type_energy - mix_scatter_type_window_width), while the high energy region is from (mix_scatter_type_energy + mix_scatter_type_window_width) to infinity. In between, the probability for string fragmentation increases smoothly from 0 to 1 as the c.m. energy.

If use_transition_probability is false: The string fragmentation is implemented similarly to what is in UrQMD (Bass:1998ca [8]). If sqrts is lower than some cutoff value, there are no strings. If higher, strings are allowed, with the cross-section being the difference between some parametrized total cross-section and the sum of all other channels, if this parametrization is larger than the sum of the channels. If not, strings are not allowed (this cross-section check is performed directly after the function is called, for technical reasons).

Both of these methods are initially implemented for NN and Npi cross- sections, and extended using the AQM to all BB, BM and MM interactions.

Baryon-antibaryon annihilation also uses this function to decide whether to produce strings or not. Since there are no other contributions for this process, there are no cutoffs or gradual increase in the probability of this process happening or not, it just requires the proper combination of incoming particles and config parameters.

Parameters

[in]

finder_parameters

parameters for collision finding and cross sections.

Definition at line 4742 of file crosssections.cc.

                                                                   {
  /* string fragmentation is enabled when strings_switch is on and the process
   * is included in pythia. */
  if (!finder_parameters.strings_switch) {
    return 0.;
  }
 
  const ParticleType& t1 = incoming_particles_[0].type();
  const ParticleType& t2 = incoming_particles_[1].type();
  const bool treat_BBbar_with_strings =
      (finder_parameters.nnbar_treatment == NNbarTreatment::Strings);
  const bool is_NN_scattering =
      t1.is_nucleon() && t2.is_nucleon() &&
      t1.antiparticle_sign() == t2.antiparticle_sign();
  const bool is_BBbar_scattering =
      (treat_BBbar_with_strings && is_BBbar_pair_ &&
       finder_parameters.use_AQM) ||
      (t1.is_nucleon() && t2.is_nucleon() &&
       t1.antiparticle_sign() != t2.antiparticle_sign());
  const bool is_Npi_scattering = (t1.pdgcode().is_pion() && t2.is_nucleon()) ||
                                 (t1.is_nucleon() && t2.pdgcode().is_pion());
  /* True for baryon-baryon, anti-baryon-anti-baryon, baryon-meson,
   * anti-baryon-meson and meson-meson*/
  const bool is_AQM_scattering =
      finder_parameters.use_AQM &&
      ((t1.is_baryon() && t2.is_baryon() &&
        t1.antiparticle_sign() == t2.antiparticle_sign()) ||
       ((t1.is_baryon() && t2.is_meson()) ||
        (t2.is_baryon() && t1.is_meson())) ||
       (t1.is_meson() && t2.is_meson()));
  const double mass_sum =
      incoming_particles_[0].pole_mass() + incoming_particles_[1].pole_mass();
 
  if (!is_NN_scattering && !is_BBbar_scattering && !is_Npi_scattering &&
      !is_AQM_scattering) {
    return 0.;
  } else if (is_NNbar_pair_ && !treat_BBbar_with_strings) {
    return 0.;
  } else if (is_BBbar_scattering) {
    // BBbar only goes through strings, so there are no "window" considerations
    return 1.;
  } else {
    /* true for K+ p and K0 p (+ antiparticles), which have special treatment
     * to fit data */
    const PdgCode pdg1 = t1.pdgcode(), pdg2 = t2.pdgcode();
    const bool is_KplusP =
        ((pdg1 == pdg::K_p || pdg1 == pdg::K_z) && (pdg2 == pdg::p)) ||
        ((pdg2 == pdg::K_p || pdg2 == pdg::K_z) && (pdg1 == pdg::p)) ||
        ((pdg1 == -pdg::K_p || pdg1 == -pdg::K_z) && (pdg2 == -pdg::p)) ||
        ((pdg2 == -pdg::K_p || pdg2 == -pdg::K_z) && (pdg1 == -pdg::p));
    // where to start the AQM strings above mass sum
    double aqm_offset =
        finder_parameters.transition_high_energy.sqrts_add_lower;
    if (is_KplusP) {
      /* for this specific case we have data. This corresponds to the point
       * where the AQM parametrization is smaller than the current 2to2
       * parametrization, which starts growing and diverges from exp. data */
      aqm_offset = finder_parameters.transition_high_energy.KN_offset;
    } else if (pdg1.is_pion() && pdg2.is_pion()) {
      aqm_offset = finder_parameters.transition_high_energy.pipi_offset;
    }
    /* if we do not use the probability transition algorithm, this is always a
     * string contribution if the energy is large enough */
    if (!finder_parameters.strings_with_probability) {
      return static_cast<double>(sqrt_s_ > mass_sum + aqm_offset);
    }
    /* No strings at low energy, only strings at high energy and
     * a transition region in the middle. Determine transition region: */
    double region_lower, region_upper;
    if (is_Npi_scattering) {
      std::tie(region_lower, region_upper) =
          finder_parameters.transition_high_energy.sqrts_range_Npi;
    } else if (is_NN_scattering) {
      std::tie(region_lower, region_upper) =
          finder_parameters.transition_high_energy.sqrts_range_NN;
    } else {  // AQM - Additive Quark Model
      /* Transition region around 0.9 larger than the sum of pole masses;
       * highly arbitrary, feel free to improve */
      region_lower = mass_sum + aqm_offset;
      region_upper = mass_sum + aqm_offset +
                     finder_parameters.transition_high_energy.sqrts_range_width;
    }
 
    return transition_probability_at_sqrts(region_lower, region_upper);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

transition_probability_at_sqrts()

double smash::CrossSections::transition_probability_at_sqrts	(	double	region_lower,
		double	region_upper
	)		const

Computes a smooth transition probability as a function of sqrt(s).

The probability is:

• 0 for sqrt(s) < region_lower
• 1 for sqrt(s) > region_upper
• smoothly varying between 0 and 1 inside the transition region according to a sinusoidal profile

This probability can be used to determine the relative contribution of two processes or regimes (e.g. soft ↔ hard string excitation).

Parameters

[in]	region_lower	Lower bound of the transition region in sqrt(s) [GeV]
[in]	region_upper	Upper bound of the transition region in sqrt(s) [GeV]

Returns

Transition probability in the interval [0,1]

Definition at line 4830 of file crosssections.cc.

                                                    {
  if (sqrt_s_ < region_lower) {
    return 0.;
  } else if (sqrt_s_ > region_upper) {
    return 1.;
  }
 
  /* Map sqrt_s_ from [region_lower, region_upper] to [-0.5, 0.5] that
   * sin(pi * x) goes from -1 to 1 leading to a probability within 0 and 1. */
  const double x = (sqrt_s_ - 0.5 * (region_lower + region_upper)) /
                   (region_upper - region_lower);
  assert(x >= -0.5 && x <= 0.5);
  double prob = 0.5 * (std::sin(M_PI * x) + 1.0);
  assert(prob >= 0. && prob <= 1.);
 
  return prob;
}

Here is the caller graph for this function:

elastic_parametrization()

double smash::CrossSections::elastic_parametrization ( const ScatterActionsFinderParameters &  finder_parameters ) const

private

Choose the appropriate parametrizations for given incoming particles and return the (parametrized) elastic cross section.

Parameters

[in]

finder_parameters

parameters for collision finding and cross sections, containing whether to extend string cross-sections with AQM and the offset to the minimum energy for string production in \(\pi\pi \) scatterings

Returns

Elastic cross section

Exceptions

std::runtime_error

if elastic cross section is negative.

Definition at line 516 of file crosssections.cc.

                                                                   {
  const bool use_AQM = finder_parameters.use_AQM;
  const double pipi_offset =
      finder_parameters.transition_high_energy.pipi_offset;
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  double elastic_xs = 0.0;
  if ((pdg_a.is_nucleon() && pdg_b.is_pion()) ||
      (pdg_b.is_nucleon() && pdg_a.is_pion())) {
    // Elastic Nucleon Pion Scattering
    elastic_xs = npi_el();
  } else if ((pdg_a.is_nucleon() && pdg_b.is_kaon()) ||
             (pdg_b.is_nucleon() && pdg_a.is_kaon())) {
    // Elastic Nucleon Kaon Scattering
    elastic_xs = nk_el();
  } else if (pdg_a.is_Dmeson() || pdg_b.is_Dmeson()) {
    const CharmRescattering& charm_rescattering =
        finder_parameters.charm_rescattering;
    std::optional<double> tmp_elastic_xs = std::nullopt;
    if (pdg_a.is_nucleon() || pdg_b.is_nucleon()) {
      tmp_elastic_xs = DN_elastic();
    } else if (pdg_a.is_Delta() || pdg_b.is_Delta()) {
      tmp_elastic_xs = DDelta_elastic();
    }
    if ((charm_rescattering == CharmRescattering::T_Matrix) &&
        tmp_elastic_xs.has_value()) {
      elastic_xs = tmp_elastic_xs.value();
    } else if (use_AQM) {
      /* use AQM if charm_rescattering == CharmRescattering::Resonances or if
       * tmp_elastic_xs has no value, which happens either when sqrts is above
       * the upper bound of the energy range of the underlying cross section
       * data or there is no underlying data for the two colliding particles */
      const double m1 = incoming_particles_[0].effective_mass();
      const double m2 = incoming_particles_[1].effective_mass();
      const double s = sqrt_s_ * sqrt_s_;
      elastic_xs = 2. / 3. * piplusp_elastic_AQM(s, m1, m2) *
                   finder_parameters.AQM_scaling_factor(pdg_a) *
                   finder_parameters.AQM_scaling_factor(pdg_b);
    } else {
      const ParticleType& a = incoming_particles_[0].type();
      const ParticleType& b = incoming_particles_[1].type();
      warn_if_charm_rescattering_enabled_and_AQM_disabled(sqrt_s_, a, b,
                                                          charm_rescattering);
      return 0.;
    }
  } else if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
             pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
    // Elastic Nucleon Nucleon Scattering
    elastic_xs = nn_el();
  } else if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
             pdg_a.antiparticle_sign() == -pdg_b.antiparticle_sign()) {
    // Elastic Nucleon anti-Nucleon Scattering
    elastic_xs = ppbar_elastic(sqrt_s_ * sqrt_s_);
  } else if (pdg_a.is_nucleus() || pdg_b.is_nucleus()) {
    const PdgCode& pdg_nucleus = pdg_a.is_nucleus() ? pdg_a : pdg_b;
    const PdgCode& pdg_other = pdg_a.is_nucleus() ? pdg_b : pdg_a;
    const bool is_deuteron = pdg_nucleus.is_deuteron();  // d or anti-d
    if (is_deuteron && pdg_other.is_pion()) {
      // Elastic (Anti-)deuteron Pion Scattering
      elastic_xs = deuteron_pion_elastic(sqrt_s_ * sqrt_s_);
    } else if (is_deuteron && pdg_other.is_nucleon()) {
      // Elastic (Anti-)deuteron (Anti-)Nucleon Scattering
      elastic_xs = deuteron_nucleon_elastic(sqrt_s_ * sqrt_s_);
    }
  } else if ((pdg_a.is_Dmeson() || pdg_b.is_Dmeson()) ||
             (pdg_a.is_Dstar2007() || pdg_b.is_Dstar2007())) {
    const CharmRescattering& charm_rescattering =
        finder_parameters.charm_rescattering;
    std::optional<double> tmp_elastic_xs = std::nullopt;
    if (pdg_a.is_pion() || pdg_b.is_pion()) {
      tmp_elastic_xs = Dpi_and_Dstarpi_elastic();
    } else if (pdg_a.is_eta() || pdg_b.is_eta()) {
      tmp_elastic_xs = Deta_and_Dstareta_elastic();
    } else if (pdg_a.is_kaon() || pdg_b.is_kaon()) {
      tmp_elastic_xs = DK_and_DstarK_elastic();
    }
    if ((charm_rescattering == CharmRescattering::T_Matrix) &&
        tmp_elastic_xs.has_value()) {
      elastic_xs = tmp_elastic_xs.value();
    } else if (use_AQM) {
      /* use AQM if charm_rescattering == CharmRescattering::Resonances or if
       * tmp_elastic_xs has no value, which happens either when sqrts is above
       * the upper bound of the energy range of the underlying cross section
       * data or there is no underlying data for the two colliding particles */
      const double m1 = incoming_particles_[0].effective_mass();
      const double m2 = incoming_particles_[1].effective_mass();
      const double s = sqrt_s_ * sqrt_s_;
      elastic_xs = 2. / 3. * piplusp_elastic_AQM(s, m1, m2) *
                   finder_parameters.AQM_scaling_factor(pdg_a) *
                   finder_parameters.AQM_scaling_factor(pdg_b);
    } else {
      const ParticleType& a = incoming_particles_[0].type();
      const ParticleType& b = incoming_particles_[1].type();
      warn_if_charm_rescattering_enabled_and_AQM_disabled(sqrt_s_, a, b,
                                                          charm_rescattering);
      return 0.;
    }
  } else if (use_AQM) {
    const double m1 = incoming_particles_[0].effective_mass();
    const double m2 = incoming_particles_[1].effective_mass();
    const double s = sqrt_s_ * sqrt_s_;
    if (pdg_a.is_baryon() && pdg_b.is_baryon()) {
      elastic_xs = nn_el();  // valid also for annihilation
    } else if ((pdg_a.is_meson() && pdg_b.is_baryon()) ||
               (pdg_b.is_meson() && pdg_a.is_baryon())) {
      elastic_xs = piplusp_elastic_AQM(s, m1, m2);
    } else if (pdg_a.is_meson() && pdg_b.is_meson()) {
      /* Special case: the pi+pi- elastic cross-section goes through resonances
       * at low sqrt_s, so we turn it off for this region so as not to destroy
       * the agreement with experimental data; this does not
       * apply to other pi pi cross-sections, which do not have any data */
      if (((pdg_a == pdg::pi_p && pdg_b == pdg::pi_m) ||
           (pdg_a == pdg::pi_m && pdg_b == pdg::pi_p)) &&
          (m1 + m2 + pipi_offset) > sqrt_s_) {
        elastic_xs = 0.0;
      } else {
        // meson-meson goes through scaling from π+p parametrization
        elastic_xs = 2. / 3. * piplusp_elastic_AQM(s, m1, m2);
      }
    }
    elastic_xs *= finder_parameters.AQM_scaling_factor(pdg_a) *
                  finder_parameters.AQM_scaling_factor(pdg_b);
  }
  if (elastic_xs < 0.) {
    throw_xsec_is_negative(sqrt_s_, elastic_xs, incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return elastic_xs;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

nn_el()

double smash::CrossSections::nn_el ( ) const

private

Determine the (parametrized) elastic cross section for a nucleon-nucleon (NN) collision.

Returns

Elastic cross section for NN

Exceptions

std::runtime_error

if positive cross section cannot be specified.

Definition at line 648 of file crosssections.cc.

                                  {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
 
  // Use parametrized cross sections.
  double sig_el = -1.;
  const double s = sqrt_s_ * sqrt_s_;
  const bool is_NN_pair = pdg_a.is_nucleon() && pdg_b.is_nucleon();
  if (is_NN_pair) {
    if (is_BBbar_pair_) {
      // npbar and ppbar
      sig_el = ppbar_elastic(s);
    } else {
      sig_el = (pdg_a == pdg_b) ? pp_elastic(s) : np_elastic(s);
    }
  } else {
    // AQM - Additive Quark Model
    const double m1 = incoming_particles_[0].effective_mass();
    const double m2 = incoming_particles_[1].effective_mass();
    if (is_BBbar_pair_) {
      sig_el =
          ppbar_elastic(effective_AQM_s(s, m1, m2, nucleon_mass, nucleon_mass));
    } else {
      sig_el = pp_elastic_high_energy(s, m1, m2);
    }
  }
 
  if (sig_el > 0.) {
    return sig_el;
  } else {
    throw_xsec_is_negative(sqrt_s_, sig_el, incoming_particles_[0],
                           incoming_particles_[1], __func__);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

npi_el()

double smash::CrossSections::npi_el ( ) const

private

Determine the elastic cross section for a nucleon-pion (Npi) collision.

It is given by a parametrization of experimental data.

Returns

Elastic cross section for Npi

Exceptions

std::runtime_error	if incoming particles are not nucleon+pion.
std::runtime_error	if positive cross section cannot be specified.

Definition at line 683 of file crosssections.cc.

                                   {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
 
  const PdgCode& nucleon = pdg_a.is_nucleon() ? pdg_a : pdg_b;
  const PdgCode& pion = pdg_a.is_nucleon() ? pdg_b : pdg_a;
  assert(pion != nucleon);
 
  const double s = sqrt_s_ * sqrt_s_;
 
  double sig_el = 0.;
  switch (nucleon.code()) {
    case pdg::p:
      switch (pion.code()) {
        case pdg::pi_p:
          sig_el = piplusp_elastic(s);
          break;
        case pdg::pi_m:
          sig_el = piminusp_elastic(s);
          break;
        case pdg::pi_z:
          sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
          break;
      }
      break;
    case pdg::n:
      switch (pion.code()) {
        case pdg::pi_p:
          sig_el = piminusp_elastic(s);
          break;
        case pdg::pi_m:
          sig_el = piplusp_elastic(s);
          break;
        case pdg::pi_z:
          sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
          break;
      }
      break;
    case -pdg::p:
      switch (pion.code()) {
        case pdg::pi_p:
          sig_el = piminusp_elastic(s);
          break;
        case pdg::pi_m:
          sig_el = piplusp_elastic(s);
          break;
        case pdg::pi_z:
          sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
          break;
      }
      break;
    case -pdg::n:
      switch (pion.code()) {
        case pdg::pi_p:
          sig_el = piplusp_elastic(s);
          break;
        case pdg::pi_m:
          sig_el = piminusp_elastic(s);
          break;
        case pdg::pi_z:
          sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
          break;
      }
      break;
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_el > 0) {
    return sig_el;
  } else {
    throw_xsec_is_negative(sqrt_s_, sig_el, incoming_particles_[0],
                           incoming_particles_[1], __func__);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

nk_el()

double smash::CrossSections::nk_el ( ) const

private

Determine the elastic cross section for a nucleon-kaon (NK) collision.

It is given by a parametrization of experimental data.

Returns

Elastic cross section for NK

Exceptions

std::runtime_error	if incoming particles are not nucleon+kaon.
std::runtime_error	if positive cross section cannot be specified.

Definition at line 990 of file crosssections.cc.

                                  {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
 
  const PdgCode& nucleon = pdg_a.is_nucleon() ? pdg_a : pdg_b;
  const PdgCode& kaon = pdg_a.is_nucleon() ? pdg_b : pdg_a;
  assert(kaon != nucleon);
 
  const double s = sqrt_s_ * sqrt_s_;
 
  double sig_el = 0.;
  switch (nucleon.code()) {
    case pdg::p:
      switch (kaon.code()) {
        case pdg::K_p:
          sig_el = kplusp_elastic_background(s);
          break;
        case pdg::K_m:
          sig_el = kminusp_elastic_background(s);
          break;
        case pdg::K_z:
          sig_el = k0p_elastic_background(s);
          break;
        case pdg::Kbar_z:
          sig_el = kbar0p_elastic_background(s);
          break;
      }
      break;
    case pdg::n:
      switch (kaon.code()) {
        case pdg::K_p:
          sig_el = kplusn_elastic_background(s);
          break;
        case pdg::K_m:
          sig_el = kminusn_elastic_background(s);
          break;
        case pdg::K_z:
          sig_el = k0n_elastic_background(s);
          break;
        case pdg::Kbar_z:
          sig_el = kbar0n_elastic_background(s);
          break;
      }
      break;
    case -pdg::p:
      switch (kaon.code()) {
        case pdg::K_p:
          sig_el = kminusp_elastic_background(s);
          break;
        case pdg::K_m:
          sig_el = kplusp_elastic_background(s);
          break;
        case pdg::K_z:
          sig_el = kbar0p_elastic_background(s);
          break;
        case pdg::Kbar_z:
          sig_el = k0p_elastic_background(s);
          break;
      }
      break;
    case -pdg::n:
      switch (kaon.code()) {
        case pdg::K_p:
          sig_el = kminusn_elastic_background(s);
          break;
        case pdg::K_m:
          sig_el = kplusn_elastic_background(s);
          break;
        case pdg::K_z:
          sig_el = kbar0n_elastic_background(s);
          break;
        case pdg::Kbar_z:
          sig_el = k0n_elastic_background(s);
          break;
      }
      break;
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_el > 0) {
    return sig_el;
  } else {
    throw_xsec_is_negative(sqrt_s_, sig_el, incoming_particles_[0],
                           incoming_particles_[1], __func__);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

Dpi_and_Dstarpi_elastic()

std::optional< double > smash::CrossSections::Dpi_and_Dstarpi_elastic ( ) const

private

Determine the elastic cross section for a D meson-pion (Dpi) or a D*-pion (D*pi) collision.

If the center-of-mass energy for the collision is below the lower bound of the energy range of the underlying cross section data, the return value is zero. If it is above the upper bound, the return value will be std::nullopt.

Returns

Elastic cross section for Dpi or D*pi.

Exceptions

std::runtime_error	if incoming particles are not Dpi or D*pi.
std::runtime_error	if cross section is negative.

Definition at line 1079 of file crosssections.cc.

                                                                 {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
  const auto pdg_pion =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
 
  std::optional<double> sig_el = std::nullopt;
  switch (pack(pdg_D, pdg_pion)) {
    // Checks for D mesons scatterings
    case pack(pdg::D_z, pdg::pi_p):
    case pack(pdg::Dbar_z, pdg::pi_m): {  // Same xsec for charge conjugation.
      sig_el = Dzeropiplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::pi_m):
    case pack(pdg::Dbar_z, pdg::pi_p): {  // Same xsec for charge conjugation.
      sig_el = Dzeropiminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::pi_z):
    case pack(pdg::Dbar_z, pdg::pi_z): {  // Same xsec for charge conjugation.
      sig_el = Dzeropizero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::pi_p):
    case pack(pdg::D_m, pdg::pi_m): {  // Same xsec for charge conjugation.
      sig_el = Dpluspiplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::pi_m):
    case pack(pdg::D_m, pdg::pi_p): {  // Same xsec for charge conjugation.
      sig_el = Dpluspiminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::pi_z):
    case pack(pdg::D_m, pdg::pi_z): {  // Same xsec for charge conjugation.
      sig_el = Dpluspizero_elastic(sqrt_s_);
      break;
    }
    // Checks for D* mesons scatterings
    case pack(pdg::Dstar_z, pdg::pi_p):
    case pack(pdg::Dstarbar_z, pdg::pi_m): {  // Same xs for charge conjugation.
      sig_el = Dstarzeropiplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::pi_m):
    case pack(pdg::Dstarbar_z, pdg::pi_p): {  // Same xs for charge conjugation.
      sig_el = Dstarzeropiminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::pi_z):
    case pack(pdg::Dstarbar_z, pdg::pi_z): {  // Same xs for charge conjugation.
      sig_el = Dstarzeropizero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::pi_p):
    case pack(pdg::Dstar_m, pdg::pi_m): {  // Same xsec for charge conjugation.
      sig_el = Dstarpluspiplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::pi_m):
    case pack(pdg::Dstar_m, pdg::pi_p): {  // Same xsec for charge conjugation.
      sig_el = Dstarpluspiminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::pi_z):
    case pack(pdg::Dstar_m, pdg::pi_z): {  // Same xsec for charge conjugation.
      sig_el = Dstarpluspizero_elastic(sqrt_s_);
      break;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_el.has_value() && sig_el.value() < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_el.value(), incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_el;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

Dpi_and_Dstarpi_inelastic()

double smash::CrossSections::Dpi_and_Dstarpi_inelastic ( ) const

private

Determine the inelastic cross section for a D meson-pion (Dpi) or a D*-pion (D*pi) collision.

Returns

Inlastic cross section for Dpi or D*pi.

Exceptions

std::runtime_error	if incoming particles are not Dpi or D*pi.
std::runtime_error	if cross section is negative.

Definition at line 3116 of file crosssections.cc.

                                                      {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
  const auto pdg_pion =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
 
  double sig_inel = -1.;
  switch (pack(pdg_D, pdg_pion)) {
    // Checks for D mesons scatterings
    case pack(pdg::D_z, pdg::pi_p):
    case pack(pdg::Dbar_z, pdg::pi_m): {  // Same xsec for charge conjugation.
      sig_inel = Dzeropiplus_Dpluspizero(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::pi_z):
    case pack(pdg::Dbar_z, pdg::pi_z): {  // Same xsec for charge conjugation.
      sig_inel = Dzeropizero_Dpluspiminus(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::pi_m):
    case pack(pdg::D_m, pdg::pi_p): {  // Same xsec for charge conjugation.
      sig_inel = Dpluspiminus_Dzeropizero(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::pi_z):
    case pack(pdg::D_m, pdg::pi_z): {  // Same xsec for charge conjugation.
      sig_inel = Dpluspizero_Dzeropiplus(sqrt_s_);
      break;
    }
    // Checks for D* mesons scatterings
    case pack(pdg::Dstar_z, pdg::pi_p):
    case pack(pdg::Dstarbar_z, pdg::pi_m): {  // Same xs for charge conjugation.
      sig_inel = Dstarzeropiplus_Dstarpluspizero(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::pi_z):
    case pack(pdg::Dstarbar_z, pdg::pi_z): {  // Same xs for charge conjugation.
      sig_inel = Dstarzeropizero_Dstarpluspiminus(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::pi_m):
    case pack(pdg::Dstar_m, pdg::pi_p): {  // Same xsec for charge conjugation.
      sig_inel = Dstarpluspiminus_Dstarzeropizero(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::pi_z):
    case pack(pdg::Dstar_m, pdg::pi_z): {  // Same xsec for charge conjugation.
      sig_inel = Dstarpluspizero_Dstarzeropiplus(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::pi_m):
    case pack(pdg::Dbar_z, pdg::pi_p):
    case pack(pdg::D_p, pdg::pi_p):
    case pack(pdg::D_m, pdg::pi_m):
    case pack(pdg::Dstar_z, pdg::pi_m):
    case pack(pdg::Dstarbar_z, pdg::pi_p):
    case pack(pdg::Dstar_p, pdg::pi_p):
    case pack(pdg::Dstar_m, pdg::pi_m): {
      // These combinations can only scatter elastically.
      return 0.;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_inel < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_inel, incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_inel;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

Deta_and_Dstareta_elastic()

std::optional< double > smash::CrossSections::Deta_and_Dstareta_elastic ( ) const

private

Determine the elastic cross section for a D meson-eta (Deta) or a D*-eta (D*eta) collision.

If the center-of-mass energy for the collision is below the lower bound of the energy range of the underlying cross section data, the return value is zero. If it is above the upper bound, the return value will be std::nullopt.

Returns

Elastic cross section for Deta or D*eta.

Exceptions

std::runtime_error	if incoming particles are not Deta or D*eta.
std::runtime_error	if cross section is negative.

Definition at line 1164 of file crosssections.cc.

                                                                   {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
  const auto pdg_eta =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
 
  std::optional<double> sig_el = std::nullopt;
  switch (pack(pdg_D, pdg_eta)) {
    // Checks for D mesons scatterings
    case pack(pdg::D_p, pdg::eta):
    case pack(pdg::D_m, pdg::eta): {  // Same xsec for charge conjugation.
      sig_el = Dpluseta_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::eta):
    case pack(pdg::Dbar_z, pdg::eta): {  // Same xsec for charge conjugation.
      sig_el = Dzeroeta_elastic(sqrt_s_);
      break;
    }
    // Checks for D* mesons scatterings
    case pack(pdg::Dstar_p, pdg::eta):
    case pack(pdg::Dstar_m, pdg::eta): {  // Same xsec for charge conjugation.
      sig_el = Dstarpluseta_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::eta):
    case pack(pdg::Dstarbar_z, pdg::eta): {  // Same xs for charge conjugation.
      sig_el = Dstarzeroeta_elastic(sqrt_s_);
      break;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_el.has_value() && sig_el.value() < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_el.value(), incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_el;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

DK_and_DstarK_elastic()

std::optional< double > smash::CrossSections::DK_and_DstarK_elastic ( ) const

private

Determine the elastic cross section for a D meson-kaon (DK) or a D*-kaon (D*K) collision.

If the center-of-mass energy for the collision is below the lower bound of the energy range of the underlying cross section data, the return value is zero. If it is above the upper bound, the return value will be std::nullopt.

Returns

Elastic cross section for DK or D*K.

Exceptions

std::runtime_error	if incoming particles are not DK or D*K.
std::runtime_error	if cross section is negative.

Definition at line 1209 of file crosssections.cc.

                                                               {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
  const auto pdg_kaon =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
 
  std::optional<double> sig_el = std::nullopt;
  switch (pack(pdg_D, pdg_kaon)) {
    // Checks for D mesons scatterings
    case pack(pdg::D_p, pdg::K_p):
    case pack(pdg::D_m, pdg::K_m): {  // Same xsec for charge conjugation.
      sig_el = DplusKplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::K_z):
    case pack(pdg::D_m, pdg::Kbar_z): {  // Same xsec for charge conjugation.
      sig_el = DplusKzero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::K_p):
    case pack(pdg::Dbar_z, pdg::K_m): {  // Same xsec for charge conjugation.
      sig_el = DzeroKplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::K_z):
    case pack(pdg::Dbar_z, pdg::Kbar_z): {  // Same xsec for charge conjugation.
      sig_el = DzeroKzero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::Kbar_z):
    case pack(pdg::D_m, pdg::K_z): {  // Same xsec for charge conjugation.
      sig_el = DplusKbarzero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::K_m):
    case pack(pdg::D_m, pdg::K_p): {  // Same xsec for charge conjugation.
      sig_el = DplusKminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Kbar_z):
    case pack(pdg::Dbar_z, pdg::K_z): {  // Same xsec for charge conjugation.
      sig_el = DzeroKbarzero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::K_m):
    case pack(pdg::Dbar_z, pdg::K_p): {  // Same xsec for charge conjugation.
      sig_el = DzeroKminus_elastic(sqrt_s_);
      break;
    }
    // Checks for D* mesons scatterings
    case pack(pdg::Dstar_p, pdg::K_p):
    case pack(pdg::Dstar_m, pdg::K_m): {  // Same xsec for charge conjugation.
      sig_el = DstarplusKplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::K_z):
    case pack(pdg::Dstar_m, pdg::Kbar_z): {  // Same xs for charge conjugation.
      sig_el = DstarplusKzero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::K_p):
    case pack(pdg::Dstarbar_z, pdg::K_m): {  // Same xs for charge conjugation.
      sig_el = DstarzeroKplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::K_z):
    case pack(pdg::Dstarbar_z, pdg::Kbar_z): {  // Same xs for charge conjugat.
      sig_el = DstarzeroKzero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::Kbar_z):
    case pack(pdg::Dstar_m, pdg::K_z): {  // Same xsec for charge conjugation.
      sig_el = DstarplusKbarzero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::K_m):
    case pack(pdg::Dstar_m, pdg::K_p): {  // Same xsec for charge conjugation.
      sig_el = DstarplusKminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::Kbar_z):
    case pack(pdg::Dstarbar_z, pdg::K_z): {  // Same xs for charge conjugation.
      sig_el = DstarzeroKbarzero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::K_m):
    case pack(pdg::Dstarbar_z, pdg::K_p): {  // Same xs for charge conjugation.
      sig_el = DstarzeroKminus_elastic(sqrt_s_);
      break;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_el.has_value() && sig_el.value() < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_el.value(), incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_el;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

DK_and_DstarK_inelastic()

double smash::CrossSections::DK_and_DstarK_inelastic ( ) const

private

Determine the inelastic cross section for a D meson-kaon (DK) or a D*-kaon (D*K) collision.

Returns

Inlastic cross section for DK or D*K.

Exceptions

std::runtime_error	if incoming particles are not DK or D*K.
std::runtime_error	if cross section is negative.

Definition at line 3192 of file crosssections.cc.

                                                    {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
  const auto pdg_kaon =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
 
  double sig_inel = -1.;
  switch (pack(pdg_D, pdg_kaon)) {
    // Checks for D mesons scatterings
    case pack(pdg::D_p, pdg::K_z):
    case pack(pdg::D_m, pdg::Kbar_z): {  // Same xsec for charge conjugation.
      sig_inel = DplusKzero_DzeroKplus(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::K_p):
    case pack(pdg::Dbar_z, pdg::K_m): {  // Same xsec for charge conjugation.
      sig_inel = DzeroKplus_DplusKzero(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::K_m):
    case pack(pdg::D_m, pdg::K_p): {  // Same xsec for charge conjugation.
      sig_inel = DplusKminus_DzeroKbarzero(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Kbar_z):
    case pack(pdg::Dbar_z, pdg::K_z): {  // Same xsec for charge conjugation.
      sig_inel = DzeroKbarzero_DplusKminus(sqrt_s_);
      break;
    }
    // Checks for D* mesons scatterings
    case pack(pdg::Dstar_p, pdg::K_z):
    case pack(pdg::Dstar_m, pdg::Kbar_z): {  // Same xs for charge conjugation.
      sig_inel = DstarplusKzero_DstarzeroKplus(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::K_p):
    case pack(pdg::Dstarbar_z, pdg::K_m): {  // Same xs for charge conjugation.
      sig_inel = DstarzeroKplus_DstarplusKzero(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_p, pdg::K_m):
    case pack(pdg::Dstar_m, pdg::K_p): {  // Same xsec for charge conjugation.
      sig_inel = DstarplusKminus_DstarzeroKbarzero(sqrt_s_);
      break;
    }
    case pack(pdg::Dstar_z, pdg::Kbar_z):
    case pack(pdg::Dstarbar_z, pdg::K_z): {  // Same xs for charge conjugation.
      sig_inel = DstarzeroKbarzero_DstarplusKminus(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::K_p):
    case pack(pdg::D_p, pdg::Kbar_z):
    case pack(pdg::D_z, pdg::K_z):
    case pack(pdg::D_z, pdg::K_m):
    case pack(pdg::D_m, pdg::K_z):
    case pack(pdg::D_m, pdg::K_m):
    case pack(pdg::Dbar_z, pdg::K_p):
    case pack(pdg::Dbar_z, pdg::Kbar_z):
    case pack(pdg::Dstar_p, pdg::K_p):
    case pack(pdg::Dstar_p, pdg::Kbar_z):
    case pack(pdg::Dstar_z, pdg::K_z):
    case pack(pdg::Dstar_z, pdg::K_m):
    case pack(pdg::Dstar_m, pdg::K_z):
    case pack(pdg::Dstar_m, pdg::K_m):
    case pack(pdg::Dstarbar_z, pdg::K_p):
    case pack(pdg::Dstarbar_z, pdg::Kbar_z): {
      // These combinations can only scatter elastically.
      return 0.;
      break;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_inel < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_inel, incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_inel;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

DN_elastic()

std::optional< double > smash::CrossSections::DN_elastic ( ) const

private

Determine the elastic cross section for a D meson-nucleon (DN) collision, If the center-of-mass energy for the collision is below the lower bound of the energy range of the underlying cross section data, the return value is zero.

If it is above the upper bound, the return value will be std::nullopt.

Returns

Elastic cross section for DN.

Exceptions

std::runtime_error	if incoming particles are not DN.
std::runtime_error	if cross section is negative.

Definition at line 1314 of file crosssections.cc.

                                                    {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
  const auto pdg_nucleon = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
 
  std::optional<double> sig_el = std::nullopt;
  switch (pack(pdg_D, pdg_nucleon)) {
    case pack(pdg::D_p, pdg::n):
    case pack(pdg::D_m, -pdg::n): {  // Same xsec for charge conjugation.
      sig_el = Dplusn_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::p):
    case pack(pdg::D_m, -pdg::p): {  // Same xsec for charge conjugation.
      sig_el = Dplusp_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::n):
    case pack(pdg::Dbar_z, -pdg::n): {  // Same xsec for charge conjugation.
      sig_el = Dzeron_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::p):
    case pack(pdg::Dbar_z, -pdg::p): {  // Same xsec for charge conjugation.
      sig_el = Dzerop_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::n):
    case pack(pdg::D_p, -pdg::n): {  // Same xsec for charge conjugation.
      sig_el = Dminusn_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::p):
    case pack(pdg::D_p, -pdg::p): {  // Same xsec for charge conjugation.
      sig_el = Dminusp_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::n):
    case pack(pdg::D_z, -pdg::n): {  // Same xsec for charge conjugation.
      sig_el = Dbarzeron_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::p):
    case pack(pdg::D_z, -pdg::p): {  // Same xsec for charge conjugation.
      sig_el = Dbarzerop_elastic(sqrt_s_);
      break;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_el.has_value() && sig_el.value() < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_el.value(), incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_el;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

DN_inelastic()

double smash::CrossSections::DN_inelastic ( ) const

private

Determine the inelastic cross section for a D meson-nucleon (DN) collision.

Returns

Inlastic cross section for DN.

Exceptions

std::runtime_error	if incoming particles are not DN.
std::runtime_error	if cross section is negative.

Definition at line 3277 of file crosssections.cc.

                                         {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
  const auto pdg_nucleon = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
 
  double sig_inel = -1.;
  switch (pack(pdg_D, pdg_nucleon)) {
    case pack(pdg::D_p, pdg::n):
    case pack(pdg::D_m, -pdg::n): {  // Same xsec for charge conjugation.
      sig_inel = Dplusn_Dzerop(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::p):
    case pack(pdg::Dbar_z, -pdg::p): {  // Same xsec for charge conjugation.
      sig_inel = Dzerop_Dplusn(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::p):
    case pack(pdg::D_p, -pdg::p): {  // Same xsec for charge conjugation.
      sig_inel = Dminusp_Dbarzeron(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::n):
    case pack(pdg::D_z, -pdg::n): {  // Same xsec for charge conjugation.
      sig_inel = Dbarzeron_Dminusp(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, -pdg::n):
    case pack(pdg::D_p, pdg::p):
    case pack(pdg::D_z, -pdg::p):
    case pack(pdg::D_z, pdg::n):
    case pack(pdg::D_m, pdg::n):
    case pack(pdg::D_m, -pdg::p):
    case pack(pdg::Dbar_z, pdg::p):
    case pack(pdg::Dbar_z, -pdg::n): {
      // These combinations can only scatter elastically.
      return 0.;
      break;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_inel < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_inel, incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_inel;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

DDelta_elastic()

std::optional< double > smash::CrossSections::DDelta_elastic ( ) const

private

Determine the elastic cross section for a D meson-Delta (DΔ) collision, If the center-of-mass energy for the collision is below the lower bound of the energy range of the underlying cross section data, the return value is zero.

If it is above the upper bound, the return value will be std::nullopt.

Returns

Elastic cross section for DΔ.

Exceptions

std::runtime_error	if incoming particles are not DΔ.
std::runtime_error	if cross section is negative.

Definition at line 1375 of file crosssections.cc.

                                                        {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
  const auto pdg_Delta = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
 
  std::optional<double> sig_el = std::nullopt;
  switch (pack(pdg_D, pdg_Delta)) {
    case pack(pdg::D_p, pdg::Delta_p):
    case pack(pdg::D_m, -pdg::Delta_p): {  // Same xsec for charge conjugation.
      sig_el = DplusDeltaplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::Delta_pp):
    case pack(pdg::D_m, -pdg::Delta_pp): {  // Same xsec for charge conjugation.
      sig_el = DplusDeltaplusplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::Delta_m):
    case pack(pdg::D_m, -pdg::Delta_m): {  // Same xsec for charge conjugation.
      sig_el = DplusDeltaminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::Delta_z):
    case pack(pdg::D_m, -pdg::Delta_z): {  // Same xsec for charge conjugation.
      sig_el = DplusDeltazero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_p):
    case pack(pdg::Dbar_z, -pdg::Delta_p): {  // Same xs for charge conjugation.
      sig_el = DzeroDeltaplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_pp):
    case pack(pdg::Dbar_z, -pdg::Delta_pp): {  // Same xs for charge conjugate.
      sig_el = DzeroDeltaplusplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_m):
    case pack(pdg::Dbar_z, -pdg::Delta_m): {  // Same xs for charge conjugation.
      sig_el = DzeroDeltaminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_z):
    case pack(pdg::Dbar_z, -pdg::Delta_z): {  // Same xs for charge conjugation.
      sig_el = DzeroDeltazero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_p):
    case pack(pdg::D_p, -pdg::Delta_p): {  // Same xsec for charge conjugation.
      sig_el = DminusDeltaplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_pp):
    case pack(pdg::D_p, -pdg::Delta_pp): {  // Same xsec for charge conjugation.
      sig_el = DminusDeltaplusplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_m):
    case pack(pdg::D_p, -pdg::Delta_m): {  // Same xsec for charge conjugation.
      sig_el = DminusDeltaminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_z):
    case pack(pdg::D_p, -pdg::Delta_z): {  // Same xsec for charge conjugation.
      sig_el = DminusDeltazero_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_p):
    case pack(pdg::D_z, -pdg::Delta_p): {  // Same xsec for charge conjugation.
      sig_el = DbarzeroDeltaplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_pp):
    case pack(pdg::D_z, -pdg::Delta_pp): {  // Same xsec for charge conjugate.
      sig_el = DbarzeroDeltaplusplus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_m):
    case pack(pdg::D_z, -pdg::Delta_m): {  // Same xsec for charge conjugation.
      sig_el = DbarzeroDeltaminus_elastic(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_z):
    case pack(pdg::D_z, -pdg::Delta_z): {  // Same xsec for charge conjugation.
      sig_el = DbarzeroDeltazero_elastic(sqrt_s_);
      break;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_el.has_value() && sig_el.value() < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_el.value(), incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_el;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

DDelta_inelastic()

double smash::CrossSections::DDelta_inelastic ( ) const

private

Determine the inelastic cross section for a D meson-Delta (DΔ) collision.

Returns

Inlastic cross section for DΔ.

Exceptions

std::runtime_error	if incoming particles are not DΔ.
std::runtime_error	if cross section is negative.

Definition at line 3330 of file crosssections.cc.

                                             {
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
  const auto pdg_Delta = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
 
  double sig_inel = -1.;
  switch (pack(pdg_D, pdg_Delta)) {
    case pack(pdg::D_p, pdg::Delta_p):
    case pack(pdg::D_m, -pdg::Delta_p): {  // Same xsec for charge conjugation.
      sig_inel = DplusDeltaplus_DzeroDeltaplusplus(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::Delta_m):
    case pack(pdg::D_m, -pdg::Delta_m): {  // Same xsec for charge conjugation.
      sig_inel = DplusDeltaminus_DzeroDeltazero(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::Delta_z):
    case pack(pdg::D_m, -pdg::Delta_z): {  // Same xsec for charge conjugation.
      sig_inel = DplusDeltazero_DzeroDeltaplus(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_p):
    case pack(pdg::Dbar_z, -pdg::Delta_p): {  // Same xs for charge conjugation.
      sig_inel = DzeroDeltaplus_DplusDeltazero(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_pp):
    case pack(pdg::Dbar_z, -pdg::Delta_pp): {  // Same xs for charge conjugate.
      sig_inel = DzeroDeltaplusplus_DplusDeltaplus(sqrt_s_);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_z):
    case pack(pdg::Dbar_z, -pdg::Delta_z): {  // Same xs for charge conjugation.
      sig_inel = DzeroDeltazero_DplusDeltaminus(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_p):
    case pack(pdg::D_p, -pdg::Delta_p): {  // Same xsec for charge conjugation.
      sig_inel = DminusDeltaplus_DbarzeroDeltazero(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_pp):
    case pack(pdg::D_p, -pdg::Delta_pp): {  // Same xsec for charge conjugation.
      sig_inel = DminusDeltaplusplus_DbarzeroDeltaplus(sqrt_s_);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_z):
    case pack(pdg::D_p, -pdg::Delta_z): {  // Same xsec for charge conjugation.
      sig_inel = DminusDeltazero_DbarzeroDeltaminus(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_p):
    case pack(pdg::D_z, -pdg::Delta_p): {  // Same xsec for charge conjugation.
      sig_inel = DbarzeroDeltaplus_DminusDeltaplusplus(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_m):
    case pack(pdg::D_z, -pdg::Delta_m): {  // Same xsec for charge conjugation.
      sig_inel = DbarzeroDeltaminus_DminusDeltazero(sqrt_s_);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_z):
    case pack(pdg::D_z, -pdg::Delta_z): {  // Same xsec for charge conjugation.
      sig_inel = DbarzeroDeltazero_DminusDeltaplus(sqrt_s_);
      break;
    }
    case pack(pdg::D_p, pdg::Delta_pp):
    case pack(pdg::D_p, -pdg::Delta_m):
    case pack(pdg::D_z, -pdg::Delta_pp):
    case pack(pdg::D_z, pdg::Delta_m):
    case pack(pdg::D_m, -pdg::Delta_pp):
    case pack(pdg::D_m, pdg::Delta_m):
    case pack(pdg::Dbar_z, pdg::Delta_pp):
    case pack(pdg::Dbar_z, -pdg::Delta_m): {
      // These combinations can only scatter elastically.
      return 0.;
      break;
    }
    default:
      throw_xsec_is_not_implemented(incoming_particles_[0],
                                    incoming_particles_[1], __func__);
  }
 
  if (sig_inel < 0.) {
    throw_xsec_is_negative(sqrt_s_, sig_inel, incoming_particles_[0],
                           incoming_particles_[1], __func__);
  } else {
    return sig_inel;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

npi_yk()

CollisionBranchList smash::CrossSections::npi_yk ( ) const

private

Find all processes for Nucleon-Pion to Hyperon-Kaon Scattering.

These scatterings are suppressed at high energies when strings are turned on with probabilities, so they need to be added back manually.

Returns

List of all possible Npi -> YK reactions with their cross sections

Definition at line 760 of file crosssections.cc.

                                                {
  const ParticleType& a = incoming_particles_[0].type();
  const ParticleType& b = incoming_particles_[1].type();
  const ParticleType& type_nucleon = a.pdgcode().is_nucleon() ? a : b;
  const ParticleType& type_pion = a.pdgcode().is_nucleon() ? b : a;
 
  const auto pdg_nucleon = type_nucleon.pdgcode().code();
  const auto pdg_pion = type_pion.pdgcode().code();
 
  const double s = sqrt_s_ * sqrt_s_;
 
  /* The cross sections are paramectrized for four isospin channels. The
   * cross sections of the rest isospin channels are obtained using
   * Clebsch-Gordan coefficients */
 
  CollisionBranchList process_list;
  switch (pdg_nucleon) {
    case pdg::p: {
      switch (pdg_pion) {
        case pdg::pi_p: {
          const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
          const auto& type_K_p = ParticleType::find(pdg::K_p);
          add_channel(
              process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
              sqrt_s_, type_K_p, type_Sigma_p);
          break;
        }
        case pdg::pi_m: {
          const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
          const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
          const auto& type_Lambda = ParticleType::find(pdg::Lambda);
          const auto& type_K_p = ParticleType::find(pdg::K_p);
          const auto& type_K_z = ParticleType::find(pdg::K_z);
          add_channel(
              process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
              sqrt_s_, type_K_p, type_Sigma_m);
          add_channel(
              process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
              type_K_z, type_Sigma_z);
          add_channel(
              process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
              type_K_z, type_Lambda);
          break;
        }
        case pdg::pi_z: {
          const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
          const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
          const auto& type_Lambda = ParticleType::find(pdg::Lambda);
          const auto& type_K_p = ParticleType::find(pdg::K_p);
          const auto& type_K_z = ParticleType::find(pdg::K_z);
          add_channel(
              process_list,
              [&] {
                return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
                              piminusp_sigma0k0_res(s) +
                              piminusp_sigmaminuskplus_pdg(s));
              },
              sqrt_s_, type_K_p, type_Sigma_z);
          add_channel(
              process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
              type_K_z, type_Sigma_p);
          add_channel(
              process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
              sqrt_s_, type_K_p, type_Lambda);
          break;
        }
      }
      break;
    }
    case pdg::n: {
      switch (pdg_pion) {
        case pdg::pi_p: {
          const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
          const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
          const auto& type_Lambda = ParticleType::find(pdg::Lambda);
          const auto& type_K_p = ParticleType::find(pdg::K_p);
          const auto& type_K_z = ParticleType::find(pdg::K_z);
          add_channel(
              process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
              sqrt_s_, type_K_z, type_Sigma_p);
          add_channel(
              process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
              type_K_p, type_Sigma_z);
          add_channel(
              process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
              type_K_p, type_Lambda);
          break;
        }
        case pdg::pi_m: {
          const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
          const auto& type_K_z = ParticleType::find(pdg::K_z);
          add_channel(
              process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
              sqrt_s_, type_K_z, type_Sigma_m);
          break;
        }
        case pdg::pi_z: {
          const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
          const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
          const auto& type_Lambda = ParticleType::find(pdg::Lambda);
          const auto& type_K_p = ParticleType::find(pdg::K_p);
          const auto& type_K_z = ParticleType::find(pdg::K_z);
          add_channel(
              process_list,
              [&] {
                return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
                              piminusp_sigma0k0_res(s) +
                              piminusp_sigmaminuskplus_pdg(s));
              },
              sqrt_s_, type_K_z, type_Sigma_z);
          add_channel(
              process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
              type_K_p, type_Sigma_m);
          add_channel(
              process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
              sqrt_s_, type_K_z, type_Lambda);
          break;
        }
      }
      break;
    }
    case -pdg::p: {
      switch (pdg_pion) {
        case pdg::pi_p: {
          const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
          const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
          const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
          const auto& type_K_m = ParticleType::find(-pdg::K_p);
          const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
          add_channel(
              process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
              sqrt_s_, type_K_m, type_Sigma_m_bar);
          add_channel(
              process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
              type_Kbar_z, type_Sigma_z_bar);
          add_channel(
              process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
              type_Kbar_z, type_Lambda_bar);
          break;
        }
        case pdg::pi_m: {
          const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
          const auto& type_K_m = ParticleType::find(-pdg::K_p);
          add_channel(
              process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
              sqrt_s_, type_K_m, type_Sigma_p_bar);
          break;
        }
        case pdg::pi_z: {
          const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
          const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
          const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
          const auto& type_K_m = ParticleType::find(-pdg::K_p);
          const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
          add_channel(
              process_list,
              [&] {
                return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
                              piminusp_sigma0k0_res(s) +
                              piminusp_sigmaminuskplus_pdg(s));
              },
              sqrt_s_, type_K_m, type_Sigma_z_bar);
          add_channel(
              process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
              type_Kbar_z, type_Sigma_p_bar);
          add_channel(
              process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
              sqrt_s_, type_K_m, type_Lambda_bar);
          break;
        }
      }
      break;
    }
    case -pdg::n: {
      switch (pdg_pion) {
        case pdg::pi_p: {
          const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
          const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
          add_channel(
              process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
              sqrt_s_, type_Kbar_z, type_Sigma_m_bar);
          break;
        }
        case pdg::pi_m: {
          const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
          const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
          const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
          const auto& type_K_m = ParticleType::find(-pdg::K_p);
          const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
          add_channel(
              process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
              sqrt_s_, type_Kbar_z, type_Sigma_p_bar);
          add_channel(
              process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
              type_K_m, type_Sigma_z_bar);
          add_channel(
              process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
              type_K_m, type_Lambda_bar);
          break;
        }
        case pdg::pi_z: {
          const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
          const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
          const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
          const auto& type_K_m = ParticleType::find(-pdg::K_p);
          const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
          add_channel(
              process_list,
              [&] {
                return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
                              piminusp_sigma0k0_res(s) +
                              piminusp_sigmaminuskplus_pdg(s));
              },
              sqrt_s_, type_Kbar_z, type_Sigma_z_bar);
          add_channel(
              process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
              type_K_m, type_Sigma_m_bar);
          add_channel(
              process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
              sqrt_s_, type_Kbar_z, type_Lambda_bar);
          break;
        }
      }
      break;
    }
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

bb_xx_except_nn()

CollisionBranchList smash::CrossSections::bb_xx_except_nn ( const ReactionsBitSet &  included_2to2 ) const

private

Find all inelastic 2->2 processes for Baryon-Baryon (BB) Scattering except the more specific Nucleon-Nucleon Scattering.

Parameters

[in]

included_2to2

Which 2->2 reactions are enabled?

Returns

List of all possible BB reactions with their cross sections

Definition at line 1808 of file crosssections.cc.

                                                {
  CollisionBranchList process_list;
  const ParticleType& type_a = incoming_particles_[0].type();
  const ParticleType& type_b = incoming_particles_[1].type();
 
  bool same_sign = type_a.antiparticle_sign() == type_b.antiparticle_sign();
  bool any_nucleus = type_a.is_nucleus() || type_b.is_nucleus();
  if (!same_sign && !any_nucleus) {
    return process_list;
  }
  bool anti_particles = type_a.antiparticle_sign() == -1;
  if (type_a.is_nucleon() || type_b.is_nucleon()) {
    // N R → N N, N̅ R → N̅ N̅
    if (included_2to2[IncludedReactions::NN_to_NR] == 1) {
      process_list = bar_bar_to_nuc_nuc(anti_particles);
    }
  } else if (type_a.is_Delta() || type_b.is_Delta()) {
    // Δ R → N N, Δ̅ R → N̅ N̅
    if (included_2to2[IncludedReactions::NN_to_DR] == 1) {
      process_list = bar_bar_to_nuc_nuc(anti_particles);
    }
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

nn_xx()

CollisionBranchList smash::CrossSections::nn_xx ( const ReactionsBitSet &  included_2to2 ) const

private

Find all inelastic 2->2 processes for Nucelon-Nucelon Scattering.

Calculate cross sections for resonance production from nucleon-nucleon collisions (i.e. N N -> N R, N N -> Delta R).

Checks are processed in the following order:

1. Charge conservation
2. Isospin factors (Clebsch-Gordan)
3. Enough energy for all decay channels to be available for the resonance

Parameters

[in]

included_2to2

Which 2->2 reactions are enabled?

Returns

List of resonance production processes possible in the collision of the two nucleons. Each element in the list contains the type(s) of the final state particle(s) and the cross section for that particular process.

Definition at line 1835 of file crosssections.cc.

                                                {
  CollisionBranchList process_list, channel_list;
 
  const double sqrts = sqrt_s_;
 
  /* Find whether colliding particles are nucleons or anti-nucleons;
   * adjust lists of produced particles. */
  bool both_antinucleons =
      (incoming_particles_[0].type().antiparticle_sign() == -1) &&
      (incoming_particles_[1].type().antiparticle_sign() == -1);
  const ParticleTypePtrList& nuc_or_anti_nuc =
      both_antinucleons ? ParticleType::list_anti_nucleons()
                        : ParticleType::list_nucleons();
  const ParticleTypePtrList& delta_or_anti_delta =
      both_antinucleons ? ParticleType::list_anti_Deltas()
                        : ParticleType::list_Deltas();
  // Find N N → N R channels.
  if (included_2to2[IncludedReactions::NN_to_NR] == 1) {
    channel_list = find_nn_xsection_from_type(
        ParticleType::list_baryon_resonances(), nuc_or_anti_nuc,
        [&sqrts](const ParticleType& type_res_1, const ParticleType&) {
          return type_res_1.iso_multiplet()->get_integral_NR(sqrts);
        });
    process_list.reserve(process_list.size() + channel_list.size());
    std::move(channel_list.begin(), channel_list.end(),
              std::inserter(process_list, process_list.end()));
    channel_list.clear();
  }
 
  // Find N N → Δ R channels.
  if (included_2to2[IncludedReactions::NN_to_DR] == 1) {
    channel_list = find_nn_xsection_from_type(
        ParticleType::list_baryon_resonances(), delta_or_anti_delta,
        [&sqrts](const ParticleType& type_res_1,
                 const ParticleType& type_res_2) {
          return type_res_1.iso_multiplet()->get_integral_RR(
              type_res_2.iso_multiplet(), sqrts);
        });
    process_list.reserve(process_list.size() + channel_list.size());
    std::move(channel_list.begin(), channel_list.end(),
              std::inserter(process_list, process_list.end()));
    channel_list.clear();
  }
 
  // Find N N → dπ and N̅ N̅→ d̅π channels.
  ParticleTypePtr deuteron = ParticleType::try_find(pdg::deuteron);
  ParticleTypePtr antideutron = ParticleType::try_find(pdg::antideuteron);
  ParticleTypePtr pim = ParticleType::try_find(pdg::pi_m);
  ParticleTypePtr pi0 = ParticleType::try_find(pdg::pi_z);
  ParticleTypePtr pip = ParticleType::try_find(pdg::pi_p);
  // Make sure all the necessary particle types are found
  if (deuteron && antideutron && pim && pi0 && pip &&
      included_2to2[IncludedReactions::PiDeuteron_to_NN] == 1) {
    const ParticleTypePtrList deutron_list = {deuteron};
    const ParticleTypePtrList antideutron_list = {antideutron};
    const ParticleTypePtrList pion_list = {pim, pi0, pip};
    channel_list = find_nn_xsection_from_type(
        (both_antinucleons ? antideutron_list : deutron_list), pion_list,
        [&sqrts](const ParticleType& type_res_1,
                 const ParticleType& type_res_2) {
          return pCM(sqrts, type_res_1.mass(), type_res_2.mass());
        });
    process_list.reserve(process_list.size() + channel_list.size());
    std::move(channel_list.begin(), channel_list.end(),
              std::inserter(process_list, process_list.end()));
    channel_list.clear();
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

nk_xx()

CollisionBranchList smash::CrossSections::nk_xx ( const ReactionsBitSet &  included_2to2, double  KN_offset  ) const

private

Find all inelastic 2->2 background processes for Nucleon-Kaon (NK) Scattering.

Parameters

[in]	included_2to2	Which 2->2 reactions are enabled?
[in]	KN_offset	Offset to the minimum energy for string production in KN scatterings

Returns

List of all possible NK reactions with their cross sections

Definition at line 1907 of file crosssections.cc.

                                                                       {
  const ParticleType& a = incoming_particles_[0].type();
  const ParticleType& b = incoming_particles_[1].type();
  const ParticleType& type_nucleon = a.pdgcode().is_nucleon() ? a : b;
  const ParticleType& type_kaon = a.pdgcode().is_nucleon() ? b : a;
 
  const auto pdg_nucleon = type_nucleon.pdgcode().code();
  const auto pdg_kaon = type_kaon.pdgcode().code();
 
  const double s = sqrt_s_ * sqrt_s_;
 
  // Some variable declarations for frequently used quantities
  const auto sigma_kplusp = kplusp_inelastic_background(s);
  const auto sigma_kplusn = kplusn_inelastic_background(s);
 
  /* At high energy, the parametrization we use diverges from experimental
   * data. This cutoff represents the point where the AQM cross section
   * becomes smaller than this parametrization, so we cut it here, and fully
   * switch to AQM beyond this point. */
  const double KN_to_KDelta_cutoff = KN_offset +
                                     incoming_particles_[0].pole_mass() +
                                     incoming_particles_[1].pole_mass();
 
  bool incl_KN_to_KN = included_2to2[IncludedReactions::KN_to_KN] == 1;
  bool incl_KN_to_KDelta =
      included_2to2[IncludedReactions::KN_to_KDelta] == 1 &&
      sqrt_s_ < KN_to_KDelta_cutoff;
  bool incl_Strangeness_exchange =
      included_2to2[IncludedReactions::Strangeness_exchange] == 1;
 
  CollisionBranchList process_list;
  switch (pdg_kaon) {
    case pdg::K_m: {
      /* All inelastic K- N channels here are strangeness exchange, plus one
       * charge exchange. */
      switch (pdg_nucleon) {
        case pdg::p: {
          if (incl_Strangeness_exchange) {
            const auto& type_pi_z = ParticleType::find(pdg::pi_z);
            const auto& type_pi_m = ParticleType::find(pdg::pi_m);
            const auto& type_pi_p = ParticleType::find(pdg::pi_p);
            const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
            const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
            const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
            const auto& type_Lambda = ParticleType::find(pdg::Lambda);
            add_channel(
                process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
                sqrt_s_, type_pi_m, type_Sigma_p);
            add_channel(
                process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
                sqrt_s_, type_pi_p, type_Sigma_m);
            add_channel(
                process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Sigma_z);
            add_channel(
                process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Lambda);
          }
          if (incl_KN_to_KN) {
            const auto& type_n = ParticleType::find(pdg::n);
            const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
            add_channel(
                process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
                type_Kbar_z, type_n);
          }
          break;
        }
        case pdg::n: {
          if (incl_Strangeness_exchange) {
            const auto& type_pi_z = ParticleType::find(pdg::pi_z);
            const auto& type_pi_m = ParticleType::find(pdg::pi_m);
            const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
            const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
            const auto& type_Lambda = ParticleType::find(pdg::Lambda);
            add_channel(
                process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
                sqrt_s_, type_pi_m, type_Sigma_z);
            add_channel(
                process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Sigma_m);
            add_channel(
                process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
                sqrt_s_, type_pi_m, type_Lambda);
          }
          break;
        }
        case -pdg::p: {
          if (incl_KN_to_KDelta) {
            const auto& type_K_m = ParticleType::find(pdg::K_m);
            const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
            const auto& type_Delta_pp_bar = ParticleType::find(-pdg::Delta_pp);
            const auto& type_Delta_p_bar = ParticleType::find(-pdg::Delta_p);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
                                            type_nucleon, type_kaon,
                                            type_Kbar_z, type_Delta_pp_bar);
                },
                sqrt_s_, type_Kbar_z, type_Delta_pp_bar);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
                                            type_nucleon, type_kaon, type_K_m,
                                            type_Delta_p_bar);
                },
                sqrt_s_, type_K_m, type_Delta_p_bar);
          }
          break;
        }
        case -pdg::n: {
          if (incl_KN_to_KDelta) {
            const auto& type_K_m = ParticleType::find(pdg::K_m);
            const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
            const auto& type_Delta_p_bar = ParticleType::find(-pdg::Delta_p);
            const auto& type_Delta_z_bar = ParticleType::find(-pdg::Delta_z);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
                                            type_nucleon, type_kaon,
                                            type_Kbar_z, type_Delta_p_bar);
                },
                sqrt_s_, type_Kbar_z, type_Delta_p_bar);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
                                            type_nucleon, type_kaon, type_K_m,
                                            type_Delta_z_bar);
                },
                sqrt_s_, type_K_m, type_Delta_z_bar);
          }
          if (incl_KN_to_KN) {
            const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
            const auto& type_p_bar = ParticleType::find(-pdg::p);
            add_channel(
                process_list, [&] { return kplusn_k0p(s); }, sqrt_s_,
                type_Kbar_z, type_p_bar);
          }
          break;
        }
      }
      break;
    }
    case pdg::K_p: {
      /* All inelastic channels are K+ N -> K Delta -> K pi N, with identical
       * cross section, weighted by the isospin factor. */
      switch (pdg_nucleon) {
        case pdg::p: {
          if (incl_KN_to_KDelta) {
            const auto& type_K_p = ParticleType::find(pdg::K_p);
            const auto& type_K_z = ParticleType::find(pdg::K_z);
            const auto& type_Delta_pp = ParticleType::find(pdg::Delta_pp);
            const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusp *
                         kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
                                                       type_K_z, type_Delta_pp);
                },
                sqrt_s_, type_K_z, type_Delta_pp);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusp *
                         kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
                                                       type_K_p, type_Delta_p);
                },
                sqrt_s_, type_K_p, type_Delta_p);
          }
          break;
        }
        case pdg::n: {
          if (incl_KN_to_KDelta) {
            const auto& type_K_p = ParticleType::find(pdg::K_p);
            const auto& type_K_z = ParticleType::find(pdg::K_z);
            const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
            const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusn *
                         kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
                                                       type_K_z, type_Delta_p);
                },
                sqrt_s_, type_K_z, type_Delta_p);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusn *
                         kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
                                                       type_K_p, type_Delta_z);
                },
                sqrt_s_, type_K_p, type_Delta_z);
          }
          if (incl_KN_to_KN) {
            const auto& type_K_z = ParticleType::find(pdg::K_z);
            const auto& type_p = ParticleType::find(pdg::p);
            add_channel(
                process_list, [&] { return kplusn_k0p(s); }, sqrt_s_, type_K_z,
                type_p);
          }
          break;
        }
        case -pdg::p: {
          if (incl_Strangeness_exchange) {
            const auto& type_pi_z = ParticleType::find(pdg::pi_z);
            const auto& type_pi_m = ParticleType::find(pdg::pi_m);
            const auto& type_pi_p = ParticleType::find(pdg::pi_p);
            const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
            const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
            const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
            const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
            add_channel(
                process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
                sqrt_s_, type_pi_p, type_Sigma_p_bar);
            add_channel(
                process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
                sqrt_s_, type_pi_m, type_Sigma_m_bar);
            add_channel(
                process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Sigma_z_bar);
            add_channel(
                process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Lambda_bar);
          }
          if (incl_KN_to_KN) {
            const auto& type_n_bar = ParticleType::find(-pdg::n);
            const auto& type_K_z = ParticleType::find(pdg::K_z);
            add_channel(
                process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
                type_K_z, type_n_bar);
          }
          break;
        }
        case -pdg::n: {
          if (incl_Strangeness_exchange) {
            const auto& type_pi_z = ParticleType::find(pdg::pi_z);
            const auto& type_pi_p = ParticleType::find(pdg::pi_p);
            const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
            const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
            const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
            add_channel(
                process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
                sqrt_s_, type_pi_p, type_Sigma_z_bar);
            add_channel(
                process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Sigma_m_bar);
            add_channel(
                process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
                sqrt_s_, type_pi_p, type_Lambda_bar);
          }
          break;
        }
      }
      break;
    }
    case pdg::K_z: {
      /* K+ and K0 have the same mass and spin, so their cross sections are
       * assumed to only differ in isospin factors. For the initial state, we
       * assume that K0 p is equivalent to K+ n and K0 n equivalent to K+ p,
       * like for the elastic background. */
      switch (pdg_nucleon) {
        case pdg::p: {
          if (incl_KN_to_KDelta) {
            const auto& type_K_p = ParticleType::find(pdg::K_p);
            const auto& type_K_z = ParticleType::find(pdg::K_z);
            const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
            const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusn *
                         kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
                                                       type_K_z, type_Delta_p);
                },
                sqrt_s_, type_K_z, type_Delta_p);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusn *
                         kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
                                                       type_K_p, type_Delta_z);
                },
                sqrt_s_, type_K_p, type_Delta_z);
          }
          if (incl_KN_to_KN) {
            const auto& type_K_p = ParticleType::find(pdg::K_p);
            const auto& type_n = ParticleType::find(pdg::n);
            add_channel(
                process_list,
                [&] {
                  // The isospin factor is 1, see the parametrizations tests.
                  return kplusn_k0p(s);
                },
                sqrt_s_, type_K_p, type_n);
          }
          break;
        }
        case pdg::n: {
          if (incl_KN_to_KDelta) {
            const auto& type_K_p = ParticleType::find(pdg::K_p);
            const auto& type_K_z = ParticleType::find(pdg::K_z);
            const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
            const auto& type_Delta_m = ParticleType::find(pdg::Delta_m);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusp *
                         kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
                                                       type_K_z, type_Delta_z);
                },
                sqrt_s_, type_K_z, type_Delta_z);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusp *
                         kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
                                                       type_K_p, type_Delta_m);
                },
                sqrt_s_, type_K_p, type_Delta_m);
          }
          break;
        }
        case -pdg::p: {
          if (incl_Strangeness_exchange) {
            const auto& type_pi_z = ParticleType::find(pdg::pi_z);
            const auto& type_pi_m = ParticleType::find(pdg::pi_m);
            const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
            const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
            const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
            add_channel(
                process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
                sqrt_s_, type_pi_m, type_Sigma_z_bar);
            add_channel(
                process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Sigma_p_bar);
            add_channel(
                process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
                sqrt_s_, type_pi_m, type_Lambda_bar);
          }
          break;
        }
        case -pdg::n: {
          if (incl_Strangeness_exchange) {
            const auto& type_pi_z = ParticleType::find(pdg::pi_z);
            const auto& type_pi_m = ParticleType::find(pdg::pi_m);
            const auto& type_pi_p = ParticleType::find(pdg::pi_p);
            const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
            const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
            const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
            const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
            add_channel(
                process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
                sqrt_s_, type_pi_m, type_Sigma_m_bar);
            add_channel(
                process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
                sqrt_s_, type_pi_p, type_Sigma_p_bar);
            add_channel(
                process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Sigma_z_bar);
            add_channel(
                process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Lambda_bar);
          }
          if (incl_KN_to_KN) {
            const auto& type_K_p = ParticleType::find(pdg::K_p);
            const auto& type_p_bar = ParticleType::find(-pdg::p);
            add_channel(
                process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
                type_K_p, type_p_bar);
          }
          break;
        }
      }
      break;
    }
    case pdg::Kbar_z:
      switch (pdg_nucleon) {
        case pdg::p: {
          if (incl_Strangeness_exchange) {
            const auto& type_pi_z = ParticleType::find(pdg::pi_z);
            const auto& type_pi_p = ParticleType::find(pdg::pi_p);
            const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
            const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
            const auto& type_Lambda = ParticleType::find(pdg::Lambda);
            add_channel(
                process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Sigma_p);
            add_channel(
                process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
                sqrt_s_, type_pi_p, type_Sigma_z);
            add_channel(
                process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
                sqrt_s_, type_pi_p, type_Lambda);
          }
          break;
        }
        case pdg::n: {
          if (incl_Strangeness_exchange) {
            const auto& type_pi_z = ParticleType::find(pdg::pi_z);
            const auto& type_pi_m = ParticleType::find(pdg::pi_m);
            const auto& type_pi_p = ParticleType::find(pdg::pi_p);
            const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
            const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
            const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
            const auto& type_Lambda = ParticleType::find(pdg::Lambda);
            add_channel(
                process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
                sqrt_s_, type_pi_p, type_Sigma_m);
            add_channel(
                process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
                sqrt_s_, type_pi_m, type_Sigma_p);
            add_channel(
                process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Sigma_z);
            add_channel(
                process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
                sqrt_s_, type_pi_z, type_Lambda);
          }
          if (incl_KN_to_KN) {
            const auto& type_p = ParticleType::find(pdg::p);
            const auto& type_K_m = ParticleType::find(pdg::K_m);
            add_channel(
                process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
                type_K_m, type_p);
          }
          break;
        }
        case -pdg::p: {
          if (incl_KN_to_KDelta) {
            const auto& type_K_m = ParticleType::find(pdg::K_m);
            const auto& type_Kbar_z = type_kaon;
            const auto& type_Delta_bar_m = ParticleType::find(-pdg::Delta_p);
            const auto& type_Delta_bar_z = ParticleType::find(-pdg::Delta_z);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
                                            type_nucleon, type_kaon,
                                            type_Kbar_z, type_Delta_bar_m);
                },
                sqrt_s_, type_Kbar_z, type_Delta_bar_m);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
                                            type_nucleon, type_kaon, type_K_m,
                                            type_Delta_bar_z);
                },
                sqrt_s_, type_K_m, type_Delta_bar_z);
          }
          if (incl_KN_to_KN) {
            const auto& type_K_m = ParticleType::find(pdg::K_m);
            const auto& type_n_bar = ParticleType::find(-pdg::n);
            add_channel(
                process_list,
                [&] {
                  // The isospin factor is 1, see the parametrizations tests.
                  return kplusn_k0p(s);
                },
                sqrt_s_, type_K_m, type_n_bar);
          }
          break;
        }
        case -pdg::n: {
          if (incl_KN_to_KDelta) {
            const auto& type_K_m = ParticleType::find(pdg::K_m);
            const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
            const auto& type_Delta_z_bar = ParticleType::find(-pdg::Delta_z);
            const auto& type_Delta_m_bar = ParticleType::find(-pdg::Delta_m);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
                                            type_nucleon, type_kaon,
                                            type_Kbar_z, type_Delta_z_bar);
                },
                sqrt_s_, type_Kbar_z, type_Delta_z_bar);
            add_channel(
                process_list,
                [&] {
                  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
                                            type_nucleon, type_kaon, type_K_m,
                                            type_Delta_m_bar);
                },
                sqrt_s_, type_K_m, type_Delta_m_bar);
          }
          break;
        }
      }
      break;
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Dpi_and_Dstarpi_xx()

CollisionBranchList smash::CrossSections::Dpi_and_Dstarpi_xx ( const ReactionsBitSet &  included_2to2, CharmRescattering  charm_rescattering  ) const

private

Find all inelastic 2->2 processes for D meson-pion (Dpi) and D*-pion (D*pi) scattering.

Parameters

[in]	included_2to2	Which 2->2 reactions are enabled?
[in]	charm_rescattering	Type of charm rescattering

Returns

List of all possible Dpi or D*pi reactions with their cross sections

Definition at line 3423 of file crosssections.cc.

                                                      {
  CollisionBranchList process_list;
  if ((included_2to2[IncludedReactions::Charm_T_matrix] == 0) ||
      !(charm_rescattering == CharmRescattering::T_Matrix)) {
    return process_list;
  }
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
  const auto pdg_pion =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
 
  /* Adding the following channels, a check for detailed balance is not needed
   * because the parametrization is explicit. */
  switch (pack(pdg_D, pdg_pion)) {
    // Channels for Dpi scatterings
    case pack(pdg::D_z, pdg::pi_p): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_pi_z = ParticleType::find(pdg::pi_z);
      // D0pi+ -> D+pi0
      add_channel(
          process_list, [&] { return Dzeropiplus_Dpluspizero(sqrt_s_); },
          sqrt_s_, type_D_p, type_pi_z);
      break;
    }
    case pack(pdg::Dbar_z, pdg::pi_m): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_pi_z = ParticleType::find(pdg::pi_z);
      // D0barpi- -> D-pi0 (charge conjugation of D0pi+ -> D+pi0)
      add_channel(
          process_list, [&] { return Dzeropiplus_Dpluspizero(sqrt_s_); },
          sqrt_s_, type_D_m, type_pi_z);
      break;
    }
    case pack(pdg::D_z, pdg::pi_z): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_pi_m = ParticleType::find(pdg::pi_m);
      // D0pi0 -> D+pi-
      add_channel(
          process_list, [&] { return Dzeropizero_Dpluspiminus(sqrt_s_); },
          sqrt_s_, type_D_p, type_pi_m);
      break;
    }
    case pack(pdg::Dbar_z, pdg::pi_z): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_pi_p = ParticleType::find(pdg::pi_p);
      // Dbar0pi0 -> D-pi+ (charge conjugation of D0pi0 -> D+pi-)
      add_channel(
          process_list, [&] { return Dzeropizero_Dpluspiminus(sqrt_s_); },
          sqrt_s_, type_D_m, type_pi_p);
      break;
    }
    case pack(pdg::D_p, pdg::pi_m): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_pi_z = ParticleType::find(pdg::pi_z);
      // D+pi- -> D0pi0
      add_channel(
          process_list, [&] { return Dpluspiminus_Dzeropizero(sqrt_s_); },
          sqrt_s_, type_D_z, type_pi_z);
      break;
    }
    case pack(pdg::D_m, pdg::pi_p): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_pi_z = ParticleType::find(pdg::pi_z);
      // D-pi+ -> Dbar0pi0 (charge conjugation of D+pi- -> D0pi0)
      add_channel(
          process_list, [&] { return Dpluspiminus_Dzeropizero(sqrt_s_); },
          sqrt_s_, type_Dbar_z, type_pi_z);
      break;
    }
    case pack(pdg::D_p, pdg::pi_z): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_pi_p = ParticleType::find(pdg::pi_p);
      // D+pi0 -> D0pi+
      add_channel(
          process_list, [&] { return Dpluspizero_Dzeropiplus(sqrt_s_); },
          sqrt_s_, type_D_z, type_pi_p);
      break;
    }
    case pack(pdg::D_m, pdg::pi_z): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_pi_m = ParticleType::find(pdg::pi_m);
      // D-pi0 -> Dbar0pi- (charge conjugation of D+pi0 -> D0pi+)
      add_channel(
          process_list, [&] { return Dpluspizero_Dzeropiplus(sqrt_s_); },
          sqrt_s_, type_Dbar_z, type_pi_m);
      break;
    }
    // Channels for D*pi scatterings
    case pack(pdg::Dstar_z, pdg::pi_p): {
      const auto& type_Dstar_p = ParticleType::find(pdg::Dstar_p);
      const auto& type_pi_z = ParticleType::find(pdg::pi_z);
      // D*(2007)0pi+ -> D*(2010)+pi0
      add_channel(
          process_list,
          [&] { return Dstarzeropiplus_Dstarpluspizero(sqrt_s_); }, sqrt_s_,
          type_Dstar_p, type_pi_z);
      break;
    }
    case pack(pdg::Dstarbar_z, pdg::pi_m): {
      const auto& type_Dstar_m = ParticleType::find(pdg::Dstar_m);
      const auto& type_pi_z = ParticleType::find(pdg::pi_z);
      /* D*(2007)bar0pi- -> D*(2010)-pi0
       * (charge conjugation of D*(2007)0pi+ -> D*(2010)+pi0) */
      add_channel(
          process_list,
          [&] { return Dstarzeropiplus_Dstarpluspizero(sqrt_s_); }, sqrt_s_,
          type_Dstar_m, type_pi_z);
      break;
    }
    case pack(pdg::Dstar_z, pdg::pi_z): {
      const auto& type_Dstar_p = ParticleType::find(pdg::Dstar_p);
      const auto& type_pi_m = ParticleType::find(pdg::pi_m);
      // D*(2007)0pi0 -> D*(2010)+pi-
      add_channel(
          process_list,
          [&] { return Dstarzeropizero_Dstarpluspiminus(sqrt_s_); }, sqrt_s_,
          type_Dstar_p, type_pi_m);
      break;
    }
    case pack(pdg::Dstarbar_z, pdg::pi_z): {
      const auto& type_Dstar_m = ParticleType::find(pdg::Dstar_m);
      const auto& type_pi_p = ParticleType::find(pdg::pi_p);
      /* D*(2007)bar0pi0 -> D*(2010)-pi+
       * (charge conjugation of D*(2007)0pi0 -> D*(2010)+pi-) */
      add_channel(
          process_list,
          [&] { return Dstarzeropizero_Dstarpluspiminus(sqrt_s_); }, sqrt_s_,
          type_Dstar_m, type_pi_p);
      break;
    }
    case pack(pdg::Dstar_p, pdg::pi_m): {
      const auto& type_Dstar_z = ParticleType::find(pdg::Dstar_z);
      const auto& type_pi_z = ParticleType::find(pdg::pi_z);
      // D*(2010)+pi- -> D*(2007)0pi0
      add_channel(
          process_list,
          [&] { return Dstarpluspiminus_Dstarzeropizero(sqrt_s_); }, sqrt_s_,
          type_Dstar_z, type_pi_z);
      break;
    }
    case pack(pdg::Dstar_m, pdg::pi_p): {
      const auto& type_Dstarbar_z = ParticleType::find(pdg::Dstarbar_z);
      const auto& type_pi_z = ParticleType::find(pdg::pi_z);
      /* D*(2010)-pi+ -> D*(2007)bar0pi0
       * (charge conjugation of D*(2010)+pi- -> D*(2007)0pi0) */
      add_channel(
          process_list,
          [&] { return Dstarpluspiminus_Dstarzeropizero(sqrt_s_); }, sqrt_s_,
          type_Dstarbar_z, type_pi_z);
      break;
    }
    case pack(pdg::Dstar_p, pdg::pi_z): {
      const auto& type_Dstar_z = ParticleType::find(pdg::Dstar_z);
      const auto& type_pi_p = ParticleType::find(pdg::pi_p);
      // D*(2010)+pi0 -> D*(2007)0pi+
      add_channel(
          process_list,
          [&] { return Dstarpluspizero_Dstarzeropiplus(sqrt_s_); }, sqrt_s_,
          type_Dstar_z, type_pi_p);
      break;
    }
    case pack(pdg::Dstar_m, pdg::pi_z): {
      const auto& type_Dstarbar_z = ParticleType::find(pdg::Dstarbar_z);
      const auto& type_pi_m = ParticleType::find(pdg::pi_m);
      /* D*(2010)-pi0 -> D*(2007)bar0pi-
       * (charge conjugation of D*(2010)+pi0 -> D*(2007)0pi+) */
      add_channel(
          process_list,
          [&] { return Dstarpluspizero_Dstarzeropiplus(sqrt_s_); }, sqrt_s_,
          type_Dstarbar_z, type_pi_m);
      break;
    }
    default:
      break;
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

DK_and_DstarK_xx()

CollisionBranchList smash::CrossSections::DK_and_DstarK_xx ( const ReactionsBitSet &  included_2to2, CharmRescattering  charm_rescattering  ) const

private

Find all inelastic 2->2 processes for D meson-kaon (DK) and D*-kaon (D*K) scattering.

Parameters

[in]	included_2to2	Which 2->2 reactions are enabled?
[in]	charm_rescattering	Type of charm rescattering

Returns

List of all possible DK or D*K reactions with their cross sections

Definition at line 3606 of file crosssections.cc.

                                                      {
  CollisionBranchList process_list;
  if ((included_2to2[IncludedReactions::Charm_T_matrix] == 0) ||
      !(charm_rescattering == CharmRescattering::T_Matrix)) {
    return process_list;
  }
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
  const auto pdg_kaon =
      (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
 
  /* Adding the following channels, a check for detailed balance is not needed
   * because the parametrization is explicit. */
  switch (pack(pdg_D, pdg_kaon)) {
    // Channels for DK scatterings
    case pack(pdg::D_p, pdg::K_z): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      // D+K0 -> D0K+
      add_channel(
          process_list, [&] { return DplusKzero_DzeroKplus(sqrt_s_); }, sqrt_s_,
          type_D_z, type_K_p);
      break;
    }
    case pack(pdg::D_m, pdg::Kbar_z): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      // D-Kbar0 -> Dbar0K- (charge conjugation of D+K0 -> D0K+)
      add_channel(
          process_list, [&] { return DplusKzero_DzeroKplus(sqrt_s_); }, sqrt_s_,
          type_Dbar_z, type_K_m);
      break;
    }
    case pack(pdg::D_z, pdg::K_p): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      // D0K+ -> D+K0
      add_channel(
          process_list, [&] { return DzeroKplus_DplusKzero(sqrt_s_); }, sqrt_s_,
          type_D_p, type_K_z);
      break;
    }
    case pack(pdg::Dbar_z, pdg::K_m): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      // Dbar0K- -> D-Kbar0 (charge conjugation of D0K+ -> D+K0)
      add_channel(
          process_list, [&] { return DzeroKplus_DplusKzero(sqrt_s_); }, sqrt_s_,
          type_D_m, type_Kbar_z);
      break;
    }
    case pack(pdg::D_p, pdg::K_m): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      // D+K- -> D0Kbar0
      add_channel(
          process_list, [&] { return DplusKminus_DzeroKbarzero(sqrt_s_); },
          sqrt_s_, type_D_z, type_Kbar_z);
      break;
    }
    case pack(pdg::D_m, pdg::K_p): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      // D-K+ -> Dbar0K0 (charge conjugation of D+K- -> D0Kbar0)
      add_channel(
          process_list, [&] { return DplusKminus_DzeroKbarzero(sqrt_s_); },
          sqrt_s_, type_Dbar_z, type_K_z);
      break;
    }
    case pack(pdg::D_z, pdg::Kbar_z): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      // D0Kbar0 -> D+K-
      add_channel(
          process_list, [&] { return DzeroKbarzero_DplusKminus(sqrt_s_); },
          sqrt_s_, type_D_p, type_K_m);
      break;
    }
    case pack(pdg::Dbar_z, pdg::K_z): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      // Dbar0K0 -> D-K+ (charge conjugation of D0Kbar0 -> D+K-)
      add_channel(
          process_list, [&] { return DzeroKbarzero_DplusKminus(sqrt_s_); },
          sqrt_s_, type_D_m, type_K_p);
      break;
    }
    // Channels for D*K scatterings
    case pack(pdg::Dstar_p, pdg::K_z): {
      const auto& type_Dstar_z = ParticleType::find(pdg::Dstar_z);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      // D*(2010)+K0 -> D*(2007)0K+
      add_channel(
          process_list, [&] { return DstarplusKzero_DstarzeroKplus(sqrt_s_); },
          sqrt_s_, type_Dstar_z, type_K_p);
      break;
    }
    case pack(pdg::Dstar_m, pdg::Kbar_z): {
      const auto& type_Dstarbar_z = ParticleType::find(pdg::Dstarbar_z);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      /* D*(2010)-Kbar0 -> D*(2007)bar0K-
       * (charge conjugation of D*(2010)+K0 -> D*(2007)0K+) */
      add_channel(
          process_list, [&] { return DstarplusKzero_DstarzeroKplus(sqrt_s_); },
          sqrt_s_, type_Dstarbar_z, type_K_m);
      break;
    }
    case pack(pdg::Dstar_z, pdg::K_p): {
      const auto& type_Dstar_p = ParticleType::find(pdg::Dstar_p);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      // D*(2007)0K+ -> D*(2010)+K0
      add_channel(
          process_list, [&] { return DstarzeroKplus_DstarplusKzero(sqrt_s_); },
          sqrt_s_, type_Dstar_p, type_K_z);
      break;
    }
    case pack(pdg::Dstarbar_z, pdg::K_m): {
      const auto& type_Dstar_m = ParticleType::find(pdg::Dstar_m);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      /* D*(2007)bar0K- -> D*(2010)-Kbar0
       * (charge conjugation of D*(2007)0K+ -> D*(2010)+K0) */
      add_channel(
          process_list, [&] { return DstarzeroKplus_DstarplusKzero(sqrt_s_); },
          sqrt_s_, type_Dstar_m, type_Kbar_z);
      break;
    }
    case pack(pdg::Dstar_p, pdg::K_m): {
      const auto& type_Dstar_z = ParticleType::find(pdg::Dstar_z);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      // D*(2010)+K- -> D*(2007)0Kbar0
      add_channel(
          process_list,
          [&] { return DstarplusKminus_DstarzeroKbarzero(sqrt_s_); }, sqrt_s_,
          type_Dstar_z, type_Kbar_z);
      break;
    }
    case pack(pdg::Dstar_m, pdg::K_p): {
      const auto& type_Dstarbar_z = ParticleType::find(pdg::Dstarbar_z);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      /* D*(2010)-K+ -> D*(2007)bar0K0
       * (charge conjugation of D*(2010)+K- -> D*(2007)0Kbar0) */
      add_channel(
          process_list,
          [&] { return DstarplusKminus_DstarzeroKbarzero(sqrt_s_); }, sqrt_s_,
          type_Dstarbar_z, type_K_z);
      break;
    }
    case pack(pdg::Dstar_z, pdg::Kbar_z): {
      const auto& type_Dstar_p = ParticleType::find(pdg::Dstar_p);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      // D*(2007)0Kbar0 -> D*(2010)+K-
      add_channel(
          process_list,
          [&] { return DstarzeroKbarzero_DstarplusKminus(sqrt_s_); }, sqrt_s_,
          type_Dstar_p, type_K_m);
      break;
    }
    case pack(pdg::Dstarbar_z, pdg::K_z): {
      const auto& type_Dstar_m = ParticleType::find(pdg::Dstar_m);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      /* D*(2007)bar0K0 -> D*(2010)-K+
       * (charge conjugation of D*(2007)0Kbar0 -> D*(2010)+K-) */
      add_channel(
          process_list,
          [&] { return DstarzeroKbarzero_DstarplusKminus(sqrt_s_); }, sqrt_s_,
          type_Dstar_m, type_K_p);
      break;
    }
    default:
      break;
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

DN_xx()

CollisionBranchList smash::CrossSections::DN_xx ( const ReactionsBitSet &  included_2to2, CharmRescattering  charm_rescattering  ) const

private

Find all inelastic 2->2 processes for D meson-nucleon (DN) scatterings.

Parameters

[in]	included_2to2	Which 2->2 reactions are enabled?
[in]	charm_rescattering	Type of charm rescattering

Returns

List of all possible DN reactions with their cross sections

Definition at line 3785 of file crosssections.cc.

                                                      {
  CollisionBranchList process_list;
  if ((included_2to2[IncludedReactions::Charm_T_matrix] == 0) ||
      !(charm_rescattering == CharmRescattering::T_Matrix)) {
    return process_list;
  }
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
  const auto pdg_nucleon = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
 
  /* Adding the following channels, a check for detailed balance is not needed
   * because the parametrization is explicit. */
  switch (pack(pdg_D, pdg_nucleon)) {
    case pack(pdg::D_p, pdg::n): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_p = ParticleType::find(pdg::p);
      // D+n -> D0p
      add_channel(
          process_list, [&] { return Dplusn_Dzerop(sqrt_s_); }, sqrt_s_,
          type_D_z, type_p);
      break;
    }
    case pack(pdg::D_m, -pdg::n): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_pbar = ParticleType::find(-pdg::p);
      // D-nbar -> Dbar0pbar (charge conjugation of D+n -> D0p)
      add_channel(
          process_list, [&] { return Dplusn_Dzerop(sqrt_s_); }, sqrt_s_,
          type_Dbar_z, type_pbar);
      break;
    }
    case pack(pdg::D_z, pdg::p): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_n = ParticleType::find(pdg::n);
      // D0p -> D+n
      add_channel(
          process_list, [&] { return Dzerop_Dplusn(sqrt_s_); }, sqrt_s_,
          type_D_p, type_n);
      break;
    }
    case pack(pdg::Dbar_z, -pdg::p): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_nbar = ParticleType::find(-pdg::n);
      // Dbar0pbar -> D-nbar (charge conjugation of D0p -> D+n)
      add_channel(
          process_list, [&] { return Dzerop_Dplusn(sqrt_s_); }, sqrt_s_,
          type_D_m, type_nbar);
      break;
    }
    case pack(pdg::D_m, pdg::p): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_n = ParticleType::find(pdg::n);
      // D-p -> Dbar0n
      add_channel(
          process_list, [&] { return Dminusp_Dbarzeron(sqrt_s_); }, sqrt_s_,
          type_Dbar_z, type_n);
      break;
    }
    case pack(pdg::D_p, -pdg::p): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_nbar = ParticleType::find(-pdg::n);
      // D+pbar -> D0nbar (charge conjugation of D-p -> Dbar0n)
      add_channel(
          process_list, [&] { return Dminusp_Dbarzeron(sqrt_s_); }, sqrt_s_,
          type_D_z, type_nbar);
      break;
    }
    case pack(pdg::Dbar_z, pdg::n): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_p = ParticleType::find(pdg::p);
      // Dbar0n -> D-p
      add_channel(
          process_list, [&] { return Dbarzeron_Dminusp(sqrt_s_); }, sqrt_s_,
          type_D_m, type_p);
      break;
    }
    case pack(pdg::D_z, -pdg::n): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_pbar = ParticleType::find(-pdg::p);
      // D0nbar -> D+pbar (charge conjugation of Dbar0n -> D-p)
      add_channel(
          process_list, [&] { return Dbarzeron_Dminusp(sqrt_s_); }, sqrt_s_,
          type_D_p, type_pbar);
      break;
    }
    default:
      break;
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

DDelta_xx()

CollisionBranchList smash::CrossSections::DDelta_xx ( const ReactionsBitSet &  included_2to2, CharmRescattering  charm_rescattering  ) const

private

Find all inelastic 2->2 processes for D meson-Delta (DΔ) scatterings.

Parameters

[in]	included_2to2	Which 2->2 reactions are enabled?
[in]	charm_rescattering	Type of charm rescattering

Returns

List of all possible DΔ reactions with their cross sections

Definition at line 3880 of file crosssections.cc.

                                                      {
  CollisionBranchList process_list;
  if ((included_2to2[IncludedReactions::Charm_T_matrix] == 0) ||
      !(charm_rescattering == CharmRescattering::T_Matrix)) {
    return process_list;
  }
  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
  const auto pdg_Delta = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
 
  /* Adding the following channels, a check for detailed balance is not needed
   * because the parametrization is explicit. */
  switch (pack(pdg_D, pdg_Delta)) {
    case pack(pdg::D_p, pdg::Delta_p): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_Delta_pp = ParticleType::find(pdg::Delta_pp);
      // D+Δ+ -> D0Δ++
      add_channel(
          process_list,
          [&] { return DplusDeltaplus_DzeroDeltaplusplus(sqrt_s_); }, sqrt_s_,
          type_D_z, type_Delta_pp);
      break;
    }
    case pack(pdg::D_m, -pdg::Delta_p): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_Deltabar_mm = ParticleType::find(-pdg::Delta_pp);
      // D-Δbar- -> Dbar0Δbar-- (charge conjugation of D+Δ+ -> D0Δ++)
      add_channel(
          process_list,
          [&] { return DplusDeltaplus_DzeroDeltaplusplus(sqrt_s_); }, sqrt_s_,
          type_Dbar_z, type_Deltabar_mm);
      break;
    }
    case pack(pdg::D_p, pdg::Delta_m): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
      // D+Δ- -> D0Δ0
      add_channel(
          process_list, [&] { return DplusDeltaminus_DzeroDeltazero(sqrt_s_); },
          sqrt_s_, type_D_z, type_Delta_z);
      break;
    }
    case pack(pdg::D_m, -pdg::Delta_m): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_Deltabar_z = ParticleType::find(-pdg::Delta_z);
      // D-Δbar- -> Dbar0Δbar0 (charge conjugation of D+Δ- -> D0Δ0)
      add_channel(
          process_list, [&] { return DplusDeltaminus_DzeroDeltazero(sqrt_s_); },
          sqrt_s_, type_Dbar_z, type_Deltabar_z);
      break;
    }
    case pack(pdg::D_p, pdg::Delta_z): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
      // D+Δ0 -> D0Δ+
      add_channel(
          process_list, [&] { return DplusDeltazero_DzeroDeltaplus(sqrt_s_); },
          sqrt_s_, type_D_z, type_Delta_p);
      break;
    }
    case pack(pdg::D_m, -pdg::Delta_z): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_Deltabar_m = ParticleType::find(-pdg::Delta_p);
      // D-Δbar0 -> Dbar0Δbar- (charge conjugation of D+Δ0 -> D0Δ+)
      add_channel(
          process_list, [&] { return DplusDeltazero_DzeroDeltaplus(sqrt_s_); },
          sqrt_s_, type_Dbar_z, type_Deltabar_m);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_p): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
      // D0Δ+ -> D+Δ0
      add_channel(
          process_list, [&] { return DzeroDeltaplus_DplusDeltazero(sqrt_s_); },
          sqrt_s_, type_D_p, type_Delta_z);
      break;
    }
    case pack(pdg::Dbar_z, -pdg::Delta_p): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_Deltabar_z = ParticleType::find(-pdg::Delta_z);
      // Dbar0Δbar- -> D-Δbar0 (charge conjugation of D0Δ+ -> D+Δ0)
      add_channel(
          process_list, [&] { return DzeroDeltaplus_DplusDeltazero(sqrt_s_); },
          sqrt_s_, type_D_m, type_Deltabar_z);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_pp): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
      // D0Δ++ -> D+Δ+
      add_channel(
          process_list,
          [&] { return DzeroDeltaplusplus_DplusDeltaplus(sqrt_s_); }, sqrt_s_,
          type_D_p, type_Delta_p);
      break;
    }
    case pack(pdg::Dbar_z, -pdg::Delta_pp): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_Deltabar_m = ParticleType::find(-pdg::Delta_p);
      // Dbar0Δbar-- -> D-Δbar- (charge conjugation of D0Δ++ -> D+Δ+)
      add_channel(
          process_list,
          [&] { return DzeroDeltaplusplus_DplusDeltaplus(sqrt_s_); }, sqrt_s_,
          type_D_m, type_Deltabar_m);
      break;
    }
    case pack(pdg::D_z, pdg::Delta_z): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_Delta_m = ParticleType::find(pdg::Delta_m);
      // D0Δ0 -> D+Δ-
      add_channel(
          process_list, [&] { return DzeroDeltazero_DplusDeltaminus(sqrt_s_); },
          sqrt_s_, type_D_p, type_Delta_m);
      break;
    }
    case pack(pdg::Dbar_z, -pdg::Delta_z): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_Deltabar_p = ParticleType::find(-pdg::Delta_m);
      // Dbar0Δbar0 -> D-Δbar+ (charge conjugation of D0Δ0 -> D+Δ-)
      add_channel(
          process_list, [&] { return DzeroDeltazero_DplusDeltaminus(sqrt_s_); },
          sqrt_s_, type_D_m, type_Deltabar_p);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_p): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
      // D-Δ+ -> Dbar0Δ0
      add_channel(
          process_list,
          [&] { return DminusDeltaplus_DbarzeroDeltazero(sqrt_s_); }, sqrt_s_,
          type_Dbar_z, type_Delta_z);
      break;
    }
    case pack(pdg::D_p, -pdg::Delta_p): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_Deltabar_z = ParticleType::find(-pdg::Delta_z);
      // D+Δbar- -> D0Δbar0 (charge conjugation of D-Δ+ -> Dbar0Δ0)
      add_channel(
          process_list,
          [&] { return DminusDeltaplus_DbarzeroDeltazero(sqrt_s_); }, sqrt_s_,
          type_D_z, type_Deltabar_z);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_pp): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
      // D-Δ++ -> Dbar0Δ+
      add_channel(
          process_list,
          [&] { return DminusDeltaplusplus_DbarzeroDeltaplus(sqrt_s_); },
          sqrt_s_, type_Dbar_z, type_Delta_p);
      break;
    }
    case pack(pdg::D_p, -pdg::Delta_pp): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_Deltabar_m = ParticleType::find(-pdg::Delta_p);
      // D+Δbar-- -> D0Δbar- (charge conjugation of D-Δ++ -> Dbar0Δ+)
      add_channel(
          process_list,
          [&] { return DminusDeltaplusplus_DbarzeroDeltaplus(sqrt_s_); },
          sqrt_s_, type_D_z, type_Deltabar_m);
      break;
    }
    case pack(pdg::D_m, pdg::Delta_z): {
      const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
      const auto& type_Delta_m = ParticleType::find(pdg::Delta_m);
      // D-Δ0 -> Dbar0Δ-
      add_channel(
          process_list,
          [&] { return DminusDeltazero_DbarzeroDeltaminus(sqrt_s_); }, sqrt_s_,
          type_Dbar_z, type_Delta_m);
      break;
    }
    case pack(pdg::D_p, -pdg::Delta_z): {
      const auto& type_D_z = ParticleType::find(pdg::D_z);
      const auto& type_Deltabar_p = ParticleType::find(-pdg::Delta_m);
      // D+Δbar0 -> D0Δbar+ (charge conjugation of D-Δ0 -> Dbar0Δ-)
      add_channel(
          process_list,
          [&] { return DminusDeltazero_DbarzeroDeltaminus(sqrt_s_); }, sqrt_s_,
          type_D_z, type_Deltabar_p);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_p): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_Delta_pp = ParticleType::find(pdg::Delta_pp);
      // Dbar0Δ+ -> D-Δ++
      add_channel(
          process_list,
          [&] { return DbarzeroDeltaplus_DminusDeltaplusplus(sqrt_s_); },
          sqrt_s_, type_D_m, type_Delta_pp);
      break;
    }
    case pack(pdg::D_z, -pdg::Delta_p): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_Deltabar_mm = ParticleType::find(-pdg::Delta_pp);
      // D0Δbar- -> D+Δbar-- (charge conjugation of Dbar0Δ+ -> D-Δ++)
      add_channel(
          process_list,
          [&] { return DbarzeroDeltaplus_DminusDeltaplusplus(sqrt_s_); },
          sqrt_s_, type_D_p, type_Deltabar_mm);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_m): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
      // Dbar0Δ- -> D-Δ0
      add_channel(
          process_list,
          [&] { return DbarzeroDeltaminus_DminusDeltazero(sqrt_s_); }, sqrt_s_,
          type_D_m, type_Delta_z);
      break;
    }
    case pack(pdg::D_z, -pdg::Delta_m): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_Deltabar_z = ParticleType::find(-pdg::Delta_z);
      // D0Δbar+ -> D+Δbar0 (charge conjugation of Dbar0Δ- -> D-Δ0)
      add_channel(
          process_list,
          [&] { return DbarzeroDeltaminus_DminusDeltazero(sqrt_s_); }, sqrt_s_,
          type_D_p, type_Deltabar_z);
      break;
    }
    case pack(pdg::Dbar_z, pdg::Delta_z): {
      const auto& type_D_m = ParticleType::find(pdg::D_m);
      const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
      // Dbar0Δ0 -> D-Δ+
      add_channel(
          process_list,
          [&] { return DbarzeroDeltazero_DminusDeltaplus(sqrt_s_); }, sqrt_s_,
          type_D_m, type_Delta_p);
      break;
    }
    case pack(pdg::D_z, -pdg::Delta_z): {
      const auto& type_D_p = ParticleType::find(pdg::D_p);
      const auto& type_Deltabar_m = ParticleType::find(-pdg::Delta_p);
      // D0Δbar0 -> D+Δbar- (charge conjugation of Dbar0Δ0 -> D-Δ+)
      add_channel(
          process_list,
          [&] { return DbarzeroDeltazero_DminusDeltaplus(sqrt_s_); }, sqrt_s_,
          type_D_p, type_Deltabar_m);
      break;
    }
    default:
      break;
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

deltak_xx()

CollisionBranchList smash::CrossSections::deltak_xx ( const ReactionsBitSet &  included_2to2 ) const

private

Find all inelastic 2->2 processes for Delta-Kaon (DeltaK) Scattering.

Parameters

[in]

included_2to2

Which 2->2 reactions are enabled?

Returns

List of all possible DeltaK reactions with their cross sections

Definition at line 2408 of file crosssections.cc.

                                                {
  CollisionBranchList process_list;
  if (included_2to2[IncludedReactions::KN_to_KDelta] == 0) {
    return process_list;
  }
  const ParticleType& a = incoming_particles_[0].type();
  const ParticleType& b = incoming_particles_[1].type();
  const ParticleType& type_delta = a.pdgcode().is_Delta() ? a : b;
  const ParticleType& type_kaon = a.pdgcode().is_Delta() ? b : a;
 
  const auto pdg_delta = type_delta.pdgcode().code();
  const auto pdg_kaon = type_kaon.pdgcode().code();
 
  const double s = sqrt_s_ * sqrt_s_;
  const double pcm = cm_momentum();
  /* The cross sections are determined from the backward reactions via detailed
   * balance. The same isospin factors as for the backward reaction are used. */
  switch (pack(pdg_delta, pdg_kaon)) {
    case pack(pdg::Delta_pp, pdg::K_z):
    case pack(pdg::Delta_p, pdg::K_p): {
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
                                              type_kaon, type_p, type_K_p) *
                   kaon_nucleon_ratios.get_ratio(type_p, type_K_p, type_kaon,
                                                 type_delta) *
                   kplusp_inelastic_background(s);
          },
          sqrt_s_, type_p, type_K_p);
      break;
    }
    case pack(-pdg::Delta_pp, pdg::Kbar_z):
    case pack(-pdg::Delta_p, pdg::K_m): {
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
                                              type_kaon, type_p_bar, type_K_m) *
                   kaon_nucleon_ratios.get_ratio(type_p_bar, type_K_m,
                                                 type_kaon, type_delta) *
                   kplusp_inelastic_background(s);
          },
          sqrt_s_, type_p_bar, type_K_m);
      break;
    }
    case pack(pdg::Delta_p, pdg::K_z):
    case pack(pdg::Delta_z, pdg::K_p): {
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
                                              type_kaon, type_n, type_K_p) *
                   kaon_nucleon_ratios.get_ratio(type_n, type_K_p, type_kaon,
                                                 type_delta) *
                   kplusn_inelastic_background(s);
          },
          sqrt_s_, type_n, type_K_p);
 
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
                                              type_kaon, type_p, type_K_z) *
                   kaon_nucleon_ratios.get_ratio(type_p, type_K_z, type_kaon,
                                                 type_delta) *
                   kplusn_inelastic_background(s);
          },
          sqrt_s_, type_p, type_K_z);
      break;
    }
    case pack(-pdg::Delta_p, pdg::Kbar_z):
    case pack(-pdg::Delta_z, pdg::K_m): {
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
                                              type_kaon, type_n_bar, type_K_m) *
                   kaon_nucleon_ratios.get_ratio(type_n_bar, type_K_m,
                                                 type_kaon, type_delta) *
                   kplusn_inelastic_background(s);
          },
          sqrt_s_, type_n_bar, type_K_m);
 
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
                                              type_kaon, type_p_bar,
                                              type_Kbar_z) *
                   kaon_nucleon_ratios.get_ratio(type_p_bar, type_Kbar_z,
                                                 type_kaon, type_delta) *
                   kplusn_inelastic_background(s);
          },
          sqrt_s_, type_p_bar, type_Kbar_z);
      break;
    }
    case pack(pdg::Delta_z, pdg::K_z):
    case pack(pdg::Delta_m, pdg::K_p): {
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
                                              type_kaon, type_n, type_K_z) *
                   kaon_nucleon_ratios.get_ratio(type_n, type_K_z, type_kaon,
                                                 type_delta) *
                   kplusp_inelastic_background(s);
          },
          sqrt_s_, type_n, type_K_z);
      break;
    }
    case pack(-pdg::Delta_z, pdg::Kbar_z):
    case pack(-pdg::Delta_m, pdg::K_m): {
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
                                              type_kaon, type_n_bar,
                                              type_Kbar_z) *
                   kaon_nucleon_ratios.get_ratio(type_n_bar, type_Kbar_z,
                                                 type_kaon, type_delta) *
                   kplusp_inelastic_background(s);
          },
          sqrt_s_, type_n_bar, type_Kbar_z);
      break;
    }
    default:
      break;
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ypi_xx()

CollisionBranchList smash::CrossSections::ypi_xx ( const ReactionsBitSet &  included_2to2 ) const

private

Find all inelastic 2->2 processes for Hyperon-Pion (Ypi) Scattering.

Parameters

[in]

included_2to2

Which 2->2 reactions are enabled?

Returns

List of all possible Ypi reactions with their cross sections

Definition at line 2558 of file crosssections.cc.

                                                {
  CollisionBranchList process_list;
  if (included_2to2[IncludedReactions::Strangeness_exchange] == 0) {
    return process_list;
  }
  const ParticleType& a = incoming_particles_[0].type();
  const ParticleType& b = incoming_particles_[1].type();
  const ParticleType& type_hyperon = a.pdgcode().is_hyperon() ? a : b;
  const ParticleType& type_pion = a.pdgcode().is_hyperon() ? b : a;
 
  const auto pdg_hyperon = type_hyperon.pdgcode().code();
  const auto pdg_pion = type_pion.pdgcode().code();
 
  const double s = sqrt_s_ * sqrt_s_;
 
  switch (pack(pdg_hyperon, pdg_pion)) {
    case pack(pdg::Sigma_z, pdg::pi_m): {
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n, type_K_m) *
                   kminusn_piminussigma0(sqrt_s_);
          },
          sqrt_s_, type_n, type_K_m);
      break;
    }
    case pack(pdg::Sigma_z, pdg::pi_p): {
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p, type_Kbar_z) *
                   kminusn_piminussigma0(sqrt_s_);
          },
          sqrt_s_, type_p, type_Kbar_z);
      break;
    }
    case pack(-pdg::Sigma_z, pdg::pi_p): {
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n_bar, type_K_p) *
                   kminusn_piminussigma0(sqrt_s_);
          },
          sqrt_s_, type_n_bar, type_K_p);
      break;
    }
    case pack(-pdg::Sigma_z, pdg::pi_m): {
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p_bar, type_K_z) *
                   kminusn_piminussigma0(sqrt_s_);
          },
          sqrt_s_, type_p_bar, type_K_z);
      break;
    }
    case pack(pdg::Sigma_m, pdg::pi_z): {
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n, type_K_m) *
                   kminusn_piminussigma0(sqrt_s_);
          },
          sqrt_s_, type_n, type_K_m);
      break;
    }
    case pack(pdg::Sigma_p, pdg::pi_z): {
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p, type_Kbar_z) *
                   kminusn_piminussigma0(sqrt_s_);
          },
          sqrt_s_, type_p, type_Kbar_z);
      break;
    }
    case pack(-pdg::Sigma_m, pdg::pi_z): {
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n_bar, type_K_p) *
                   kminusn_piminussigma0(sqrt_s_);
          },
          sqrt_s_, type_n_bar, type_K_p);
      break;
    }
    case pack(-pdg::Sigma_p, pdg::pi_z): {
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p_bar, type_K_z) *
                   kminusn_piminussigma0(sqrt_s_);
          },
          sqrt_s_, type_p_bar, type_K_z);
      break;
    }
    case pack(pdg::Lambda, pdg::pi_m): {
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n, type_K_m) *
                   kminusn_piminuslambda(sqrt_s_);
          },
          sqrt_s_, type_n, type_K_m);
      break;
    }
    case pack(pdg::Lambda, pdg::pi_p): {
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p, type_Kbar_z) *
                   kminusn_piminuslambda(sqrt_s_);
          },
          sqrt_s_, type_p, type_Kbar_z);
      break;
    }
    case pack(-pdg::Lambda, pdg::pi_p): {
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n_bar, type_K_p) *
                   kminusn_piminuslambda(sqrt_s_);
          },
          sqrt_s_, type_n_bar, type_K_p);
      break;
    }
    case pack(-pdg::Lambda, pdg::pi_m): {
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p_bar, type_K_z) *
                   kminusn_piminuslambda(sqrt_s_);
          },
          sqrt_s_, type_p_bar, type_K_z);
      break;
    }
    case pack(pdg::Sigma_z, pdg::pi_z): {
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p, type_K_m) *
                   kminusp_pi0sigma0(sqrt_s_);
          },
          sqrt_s_, type_p, type_K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n, type_Kbar_z) *
                   kminusp_pi0sigma0(sqrt_s_);
          },
          sqrt_s_, type_n, type_Kbar_z);
      break;
    }
    case pack(-pdg::Sigma_z, pdg::pi_z): {
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p_bar, type_K_p) *
                   kminusp_pi0sigma0(sqrt_s_);
          },
          sqrt_s_, type_p_bar, type_K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n_bar, type_K_z) *
                   kminusp_pi0sigma0(sqrt_s_);
          },
          sqrt_s_, type_n_bar, type_K_z);
      break;
    }
    case pack(pdg::Sigma_m, pdg::pi_p): {
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p, type_K_m) *
                   kminusp_piplussigmaminus(sqrt_s_);
          },
          sqrt_s_, type_p, type_K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n, type_Kbar_z) *
                   kminusp_piminussigmaplus(sqrt_s_);
          },
          sqrt_s_, type_n, type_Kbar_z);
      break;
    }
    case pack(-pdg::Sigma_m, pdg::pi_m): {
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p_bar, type_K_p) *
                   kminusp_piplussigmaminus(sqrt_s_);
          },
          sqrt_s_, type_p_bar, type_K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n_bar, type_K_z) *
                   kminusp_piminussigmaplus(sqrt_s_);
          },
          sqrt_s_, type_n_bar, type_K_z);
      break;
    }
    case pack(pdg::Lambda, pdg::pi_z): {
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p, type_K_m) *
                   kminusp_pi0lambda(sqrt_s_);
          },
          sqrt_s_, type_p, type_K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n, type_Kbar_z) *
                   kminusp_pi0lambda(sqrt_s_);
          },
          sqrt_s_, type_n, type_Kbar_z);
      break;
    }
    case pack(-pdg::Lambda, pdg::pi_z): {
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p_bar, type_K_p) *
                   kminusp_pi0lambda(sqrt_s_);
          },
          sqrt_s_, type_p_bar, type_K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n_bar, type_K_z) *
                   kminusp_pi0lambda(sqrt_s_);
          },
          sqrt_s_, type_n_bar, type_K_z);
      break;
    }
    case pack(pdg::Sigma_p, pdg::pi_m): {
      const auto& type_p = ParticleType::find(pdg::p);
      const auto& type_n = ParticleType::find(pdg::n);
      const auto& type_K_m = ParticleType::find(pdg::K_m);
      const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p, type_K_m) *
                   kminusp_piminussigmaplus(sqrt_s_);
          },
          sqrt_s_, type_p, type_K_m);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n, type_Kbar_z) *
                   kminusp_piplussigmaminus(sqrt_s_);
          },
          sqrt_s_, type_n, type_Kbar_z);
      break;
    }
    case pack(-pdg::Sigma_p, pdg::pi_p): {
      const auto& type_p_bar = ParticleType::find(-pdg::p);
      const auto& type_n_bar = ParticleType::find(-pdg::n);
      const auto& type_K_p = ParticleType::find(pdg::K_p);
      const auto& type_K_z = ParticleType::find(pdg::K_z);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_p_bar, type_K_p) *
                   kminusp_piminussigmaplus(sqrt_s_);
          },
          sqrt_s_, type_p_bar, type_K_p);
      add_channel(
          process_list,
          [&] {
            return detailed_balance_factor_stable(s, type_hyperon, type_pion,
                                                  type_n_bar, type_K_z) *
                   kminusp_piplussigmaminus(sqrt_s_);
          },
          sqrt_s_, type_n_bar, type_K_z);
      break;
    }
    default:
      break;
  }
 
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

dpi_xx()

CollisionBranchList smash::CrossSections::dpi_xx ( const ReactionsBitSet &  included_2to2 ) const

private

Find all inelastic 2->2 processes involving Pion and (anti-) Deuteron (dpi), specifically dπ→ NN, d̅π→ N̅N̅; πd→ πd' (mockup for πd→ πnp), πd̅→ πd̅' and reverse.

Parameters

[in]

included_2to2

Which 2->2 reactions are enabled?

Returns

List of all possible dpi reactions with their cross sections

Definition at line 2955 of file crosssections.cc.

                                                {
  CollisionBranchList process_list;
  const double sqrts = sqrt_s_;
  const ParticleType& type_a = incoming_particles_[0].type();
  const ParticleType& type_b = incoming_particles_[1].type();
 
  // pi d -> N N
  bool is_pid = (type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
                (type_b.is_deuteron() && type_a.pdgcode().is_pion());
  if (is_pid && included_2to2[IncludedReactions::PiDeuteron_to_NN] == 1) {
    const int baryon_number = type_a.baryon_number() + type_b.baryon_number();
    ParticleTypePtrList nuc = (baryon_number > 0)
                                  ? ParticleType::list_nucleons()
                                  : ParticleType::list_anti_nucleons();
    const double s = sqrt_s_ * sqrt_s_;
    for (ParticleTypePtr nuc_a : nuc) {
      for (ParticleTypePtr nuc_b : nuc) {
        if (type_a.charge() + type_b.charge() !=
            nuc_a->charge() + nuc_b->charge()) {
          continue;
        }
        // loop over total isospin
        for (const int twoI : I_tot_range(*nuc_a, *nuc_b)) {
          const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
              type_a, type_b, *nuc_a, *nuc_b, twoI);
          // If Clebsch-Gordan coefficient = 0, don't bother with the rest.
          if (std::abs(isospin_factor) < really_small) {
            continue;
          }
 
          // Calculate matrix element for inverse process.
          const double matrix_element =
              nn_to_resonance_matrix_element(sqrts, type_a, type_b, twoI);
          if (matrix_element <= 0.) {
            continue;
          }
 
          const double spin_factor = (nuc_a->spin() + 1) * (nuc_b->spin() + 1);
          const int sym_fac_in =
              (type_a.iso_multiplet() == type_b.iso_multiplet()) ? 2 : 1;
          const int sym_fac_out =
              (nuc_a->iso_multiplet() == nuc_b->iso_multiplet()) ? 2 : 1;
          double p_cm_final = pCM_from_s(s, nuc_a->mass(), nuc_b->mass());
          const double xsection = isospin_factor * spin_factor * sym_fac_in /
                                  sym_fac_out * p_cm_final * matrix_element /
                                  (s * cm_momentum());
 
          if (xsection > really_small) {
            process_list.push_back(std::make_unique<CollisionBranch>(
                *nuc_a, *nuc_b, xsection, ProcessType::TwoToTwo));
            logg[LScatterAction].debug(type_a.name(), type_b.name(), "->",
                                       nuc_a->name(), nuc_b->name(),
                                       " at sqrts [GeV] = ", sqrts,
                                       " with cs[mb] = ", xsection);
          }
        }
      }
    }
  }
 
  // pi d -> pi d' (effectively pi d -> pi p n)  AND reverse, pi d' -> pi d
  bool is_pid_or_pidprime = ((type_a.is_deuteron() || type_a.is_dprime()) &&
                             type_b.pdgcode().is_pion()) ||
                            ((type_b.is_deuteron() || type_b.is_dprime()) &&
                             type_a.pdgcode().is_pion());
  if (is_pid_or_pidprime &&
      included_2to2[IncludedReactions::PiDeuteron_to_pidprime] == 1) {
    const ParticleType& type_pi = type_a.pdgcode().is_pion() ? type_a : type_b;
    const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
    ParticleTypePtrList nuclei = ParticleType::list_light_nuclei();
    for (ParticleTypePtr produced_nucleus : nuclei) {
      // Elastic collisions are treated in a different function
      if (produced_nucleus == &type_nucleus ||
          produced_nucleus->charge() != type_nucleus.charge() ||
          produced_nucleus->baryon_number() != type_nucleus.baryon_number()) {
        continue;
      }
      const double xsection =
          xs_dpi_dprimepi(sqrts, cm_momentum(), produced_nucleus, type_pi);
      process_list.push_back(std::make_unique<CollisionBranch>(
          type_pi, *produced_nucleus, xsection, ProcessType::TwoToTwo));
      logg[LScatterAction].debug(type_pi.name(), type_nucleus.name(), "→ ",
                                 type_pi.name(), produced_nucleus->name(),
                                 " at ", sqrts, " GeV, xs[mb] = ", xsection);
    }
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

dn_xx()

CollisionBranchList smash::CrossSections::dn_xx ( const ReactionsBitSet &  included_2to2 ) const

private

Find all inelastic 2->2 processes involving Nucleon and (anti-) Deuteron (dN), specifically Nd → Nd', N̅d → N̅d', N̅d̅→ N̅d̅', Nd̅→ Nd̅' and reverse (e.g.

Nd'→ Nd).

Parameters

[in]

included_2to2

Which 2->2 reactions are enabled?

Returns

List of all possible dN reactions with their cross sections

Definition at line 3086 of file crosssections.cc.

                                                {
  const ParticleType& type_a = incoming_particles_[0].type();
  const ParticleType& type_b = incoming_particles_[1].type();
  const ParticleType& type_N = type_a.is_nucleon() ? type_a : type_b;
  const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
  CollisionBranchList process_list;
  if (included_2to2[IncludedReactions::NDeuteron_to_Ndprime] == 0) {
    return process_list;
  }
  ParticleTypePtrList nuclei = ParticleType::list_light_nuclei();
 
  for (ParticleTypePtr produced_nucleus : nuclei) {
    // No elastic collisions for now, respect conservation laws
    if (produced_nucleus == &type_nucleus ||
        produced_nucleus->charge() != type_nucleus.charge() ||
        produced_nucleus->baryon_number() != type_nucleus.baryon_number()) {
      continue;
    }
    const double xsection = xs_dn_dprimen(
        sqrt_s_, cm_momentum(), produced_nucleus, type_nucleus, type_N);
    process_list.push_back(std::make_unique<CollisionBranch>(
        type_N, *produced_nucleus, xsection, ProcessType::TwoToTwo));
    logg[LScatterAction].debug(type_N.name(), type_nucleus.name(), "→ ",
                               type_N.name(), produced_nucleus->name(), " at ",
                               sqrt_s_, " GeV, xs[mb] = ", xsection);
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

xs_dpi_dprimepi()

double smash::CrossSections::xs_dpi_dprimepi ( double  sqrts, double  cm_mom, ParticleTypePtr  produced_nucleus, const ParticleType &  type_pi  )

staticprivate

Parametrized cross section for πd→ πd' (mockup for πd→ πnp), πd̅→ πd̅' and reverse, see Oliinychenko:2018ugs [49] for details.

Parameters

[in]	sqrts	square-root of mandelstam s
[in]	cm_mom	center of mass momentum of incoming particles
[in]	produced_nucleus	type of outgoing deuteron or d-prime
[in]	type_pi	type of scattering pion

Returns

cross section for given scattering

Matrix element is fit to match the inelastic pi+ d -> pi+ n p cross-section from the Fig. 5 of Arndt:1994bs [7].

Definition at line 2922 of file crosssections.cc.

                                                                   {
  const double s = sqrts * sqrts;
  // same matrix element for πd and πd̅
  const double tmp = sqrts - pion_mass - deuteron_mass;
  /**
   * Matrix element is fit to match the inelastic pi+ d -> pi+ n p cross-section
   * from the Fig. 5 of \iref{Arndt:1994bs}.
   */
  const double matrix_element =
      295.5 + 2.862 / (0.00283735 + pow_int(sqrts - 2.181, 2)) +
      0.0672 / pow_int(tmp, 2) - 6.61753 / tmp;
 
  const double spin_factor =
      (produced_nucleus->spin() + 1) * (type_pi.spin() + 1);
  /* Isospin factor is always the same, so it is included into the matrix
   * element. Symmetry factor is always 1 here. The (hbarc)^2/16 pi factor is
   * absorbed into matrix element. */
  double xsection = matrix_element * spin_factor / (s * cm_mom);
  if (produced_nucleus->is_stable()) {
    xsection *= pCM_from_s(s, type_pi.mass(), produced_nucleus->mass());
  } else {
    const double resonance_integral =
        produced_nucleus->iso_multiplet()->get_integral_piR(sqrts);
    xsection *= resonance_integral;
    logg[LScatterAction].debug("Resonance integral ", resonance_integral,
                               ", matrix element: ", matrix_element,
                               ", cm_momentum: ", cm_mom);
  }
  return xsection;
}

Here is the call graph for this function:

Here is the caller graph for this function:

xs_dn_dprimen()

double smash::CrossSections::xs_dn_dprimen ( double  sqrts, double  cm_mom, ParticleTypePtr  produced_nucleus, const ParticleType &  type_nucleus, const ParticleType &  type_N  )

staticprivate

Parametrized cross section for Nd → Nd', N̅d → N̅d', N̅d̅→ N̅d̅', Nd̅→ Nd̅' and reverse (e.g.

Nd'→ Nd), see Oliinychenko:2018ugs [49] for details.

Parameters

[in]	sqrts	square-root of mandelstam s
[in]	cm_mom	center of mass momentum of incoming particles
[in]	produced_nucleus	type of outgoing deuteron or d-prime
[in]	type_nucleus	type of scattering (incoming) deuteron or d-prime
[in]	type_N	type of scattering nucleon

Returns

cross section for given scattering

Nd → Nd', N̅d̅→ N̅d̅' and reverse: Fit to match experimental cross-section Nd -> Nnp from [19].

N̅d → N̅d', Nd̅→ Nd̅' and reverse: Fit to roughly match experimental cross-section N̅d -> N̅ np from Bizzarri:1973sp [12].

Definition at line 3045 of file crosssections.cc.

                                                                {
  const double s = sqrts * sqrts;
  double matrix_element = 0.0;
  double tmp = sqrts - nucleon_mass - deuteron_mass;
  assert(tmp >= 0.0);
  if (std::signbit(type_N.baryon_number()) ==
      std::signbit(type_nucleus.baryon_number())) {
    /**
     * Nd → Nd', N̅d̅→ N̅d̅' and reverse:
     * Fit to match experimental cross-section Nd -> Nnp from \cite Carlson1973.
     */
    matrix_element = 79.0474 / std::pow(tmp, 0.7897) + 654.596 * tmp;
  } else {
    /**
     * N̅d →  N̅d', Nd̅→ Nd̅' and reverse:
     * Fit to roughly match experimental cross-section N̅d -> N̅ np from
     * \iref{Bizzarri:1973sp}.
     */
    matrix_element = 342.572 / std::pow(tmp, 0.6);
  }
  const double spin_factor =
      (produced_nucleus->spin() + 1) * (type_N.spin() + 1);
  /* Isospin factor is always the same, so it is included into matrix element
   * Symmetry factor is always 1 here. Absorb (hbarc)^2/16 pi factor into matrix
   * element. */
  double xsection = matrix_element * spin_factor / (s * cm_mom);
  if (produced_nucleus->is_stable()) {
    assert(!type_nucleus.is_stable());
    xsection *= pCM_from_s(s, type_N.mass(), produced_nucleus->mass());
  } else {
    assert(type_nucleus.is_stable());
    const double resonance_integral =
        produced_nucleus->iso_multiplet()->get_integral_NR(sqrts);
    xsection *= resonance_integral;
  }
  return xsection;
}

Here is the call graph for this function:

Here is the caller graph for this function:

string_hard_cross_section()

double smash::CrossSections::string_hard_cross_section ( ) const

private

Determine the (parametrized) hard non-diffractive string cross section for this collision.

Returns

Parametrized cross section (without AQM scaling).

Definition at line 4425 of file crosssections.cc.

                                                      {
  double cross_sec = 0.;
  // Hard strings can only be excited if the lower cutoff by Pythia is fulfilled
  if (sqrt_s_ <= minimum_sqrts_pythia_can_handle) {
    return cross_sec;
  }
  const ParticleData& data_a = incoming_particles_[0];
  const ParticleData& data_b = incoming_particles_[1];
 
  if (data_a.is_baryon() && data_b.is_baryon()) {
    // Nucleon-nucleon cross section is used for all baryon-baryon cases.
    const double eff_s =
        effective_AQM_s(sqrt_s_ * sqrt_s_, data_a.effective_mass(),
                        data_b.effective_mass(), nucleon_mass, nucleon_mass);
    cross_sec = NN_string_hard(eff_s);
  } else if (data_a.is_baryon() || data_b.is_baryon()) {
    // Nucleon-pion cross section is used for all baryon-meson cases.
    cross_sec = Npi_string_hard(sqrt_s_ * sqrt_s_);
  } else {
    // Pion-pion cross section is used for all meson-meson cases.
    cross_sec = pipi_string_hard(sqrt_s_ * sqrt_s_);
  }
 
  return cross_sec;
}

Here is the call graph for this function:

Here is the caller graph for this function:

bar_bar_to_nuc_nuc()

CollisionBranchList smash::CrossSections::bar_bar_to_nuc_nuc ( bool  is_anti_particles ) const

private

Calculate cross sections for 2 → 2 resonance absorption (i.e.

NR → NN and ΔR → NN). See eqs. (B.6), (B.9) and (181) in Buss:2011mx [16].

Parameters

[in]

is_anti_particles

Whether the colliding particles are antiparticles

Returns

List of possible resonance absorption processes. Each element of the list contains the types of the final-state particles and the cross section for that particular process.

Definition at line 4513 of file crosssections.cc.

                                        {
  const ParticleType& type_a = incoming_particles_[0].type();
  const ParticleType& type_b = incoming_particles_[1].type();
  CollisionBranchList process_list;
 
  const double s = sqrt_s_ * sqrt_s_;
  // CM momentum in final state
  double p_cm_final = std::sqrt(s - 4. * nucleon_mass * nucleon_mass) / 2.;
 
  ParticleTypePtrList nuc_or_anti_nuc;
  if (is_anti_particles) {
    nuc_or_anti_nuc = ParticleType::list_anti_nucleons();
  } else {
    nuc_or_anti_nuc = ParticleType::list_nucleons();
  }
 
  // Loop over all nucleon or anti-nucleon charge states.
  for (ParticleTypePtr nuc_a : nuc_or_anti_nuc) {
    for (ParticleTypePtr nuc_b : nuc_or_anti_nuc) {
      /* Check for charge conservation. */
      if (type_a.charge() + type_b.charge() !=
          nuc_a->charge() + nuc_b->charge()) {
        continue;
      }
      // loop over total isospin
      for (const int twoI : I_tot_range(*nuc_a, *nuc_b)) {
        const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
            type_a, type_b, *nuc_a, *nuc_b, twoI);
        // If Clebsch-Gordan coefficient is zero, don't bother with the rest
        if (std::abs(isospin_factor) < really_small) {
          continue;
        }
 
        // Calculate matrix element for inverse process.
        const double matrix_element =
            nn_to_resonance_matrix_element(sqrt_s_, type_a, type_b, twoI);
        if (matrix_element <= 0.) {
          continue;
        }
 
        /* Cross section for 2->2 resonance absorption, obtained via detailed
         * balance from the inverse reaction. */
        const double spin_factor = (nuc_a->spin() + 1) * (nuc_b->spin() + 1);
        const int sym_fac_in =
            (type_a.iso_multiplet() == type_b.iso_multiplet()) ? 2 : 1;
        const int sym_fac_out =
            (nuc_a->iso_multiplet() == nuc_b->iso_multiplet()) ? 2 : 1;
        const double xsection = isospin_factor * spin_factor * sym_fac_in /
                                sym_fac_out * p_cm_final * matrix_element /
                                (s * cm_momentum());
 
        if (xsection > really_small) {
          process_list.push_back(std::make_unique<CollisionBranch>(
              *nuc_a, *nuc_b, xsection, ProcessType::TwoToTwo));
          logg[LCrossSections].debug(
              "2->2 absorption with original particles: ", type_a, type_b);
        }
      }
    }
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

nn_to_resonance_matrix_element()

double smash::CrossSections::nn_to_resonance_matrix_element ( double  sqrts, const ParticleType &  type_a, const ParticleType &  type_b, int  twoI  )

staticprivate

Scattering matrix amplitude squared (divided by 16π) for resonance production processes like NN → NR and NN → ΔR, where R is a baryon resonance (Δ, N*, Δ*).

Includes no spin or isospin factors.

Parameters

[in]	sqrts	sqrt(Mandelstam-s), i.e. collision CMS energy.
[in]	type_a	Type information for the first final-state particle.
[in]	type_b	Type information for the second final-state particle.
[in]	twoI	Twice the total isospin of the involved state.

Returns

Matrix amplitude squared \( |\mathcal{M}(\sqrt{s})|^2/16\pi \).

NN → NΔ: fit sqrt(s)-dependence to OBE model [Dmitriev:1986st [23]]

All other processes use a constant matrix element, similar to Bass:1998ca [8], eq. (3.35).

pn → pnη cross section is known to be larger than the corresponding pp → ppη cross section by a factor of 6.5 [Calen:1998vh [17]]. Since the eta is mainly produced by an intermediate N*(1535) we introduce an explicit isospin asymmetry for the production of N*(1535) produced in pn vs. pp similar to [Teis:1996kx [67]], eq. (29).

Definition at line 4577 of file crosssections.cc.

                                                                     {
  const double m_a = type_a.mass();
  const double m_b = type_b.mass();
  const double msqr = 2. * (m_a * m_a + m_b * m_b);
  /* If the c.m. energy is larger than the sum of the pole masses of the
   * outgoing particles plus three times the sum of the widths plus 3 GeV, the
   * collision will be neglected. This can be problematic for some final-state
   * cross sections, but at energies that high strings are used anyway. */
  const double w_a = type_a.width_at_pole();
  const double w_b = type_b.width_at_pole();
  const double uplmt = m_a + m_b + 3.0 * (w_a + w_b) + 3.0;
  if (sqrts > uplmt) {
    return 0.;
  }
  /// NN → NΔ: fit sqrt(s)-dependence to OBE model [\iref{Dmitriev:1986st}]
  if (((type_a.is_Delta() && type_b.is_nucleon()) ||
       (type_b.is_Delta() && type_a.is_nucleon())) &&
      (type_a.antiparticle_sign() == type_b.antiparticle_sign())) {
    return 68. / std::pow(sqrts - 1.104, 1.951);
    /**
     * All other processes use a constant matrix element, similar to
     * \iref{Bass:1998ca}, eq. (3.35).
     */
  } else if (((type_a.is_Nstar() && type_b.is_nucleon()) ||
              (type_b.is_Nstar() && type_a.is_nucleon())) &&
             type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
    // NN → NN*
    if (twoI == 2) {
      return 4.5 / msqr;
    } else if (twoI == 0) {
      const double parametrization = 14. / msqr;
      /**
       * pn → pnη cross section is known to be larger than the corresponding
       * pp → ppη cross section by a factor of 6.5 [\iref{Calen:1998vh}].
       * Since the eta is mainly produced by an intermediate N*(1535) we
       * introduce an explicit isospin asymmetry for the production of N*(1535)
       * produced in pn vs. pp similar to [\iref{Teis:1996kx}], eq. (29).
       */
      if (type_a.is_Nstar1535() || type_b.is_Nstar1535()) {
        return 6.5 * parametrization;
      } else {
        return parametrization;
      }
    }
  } else if (((type_a.is_Deltastar() && type_b.is_nucleon()) ||
              (type_b.is_Deltastar() && type_a.is_nucleon())) &&
             type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
    // NN → NΔ*
    return 15. / msqr;
  } else if ((type_a.is_Delta() && type_b.is_Delta()) &&
             (type_a.antiparticle_sign() == type_b.antiparticle_sign())) {
    // NN → ΔΔ
    if (twoI == 2) {
      return 45. / msqr;
    } else if (twoI == 0) {
      return 120. / msqr;
    }
  } else if (((type_a.is_Nstar() && type_b.is_Delta()) ||
              (type_b.is_Nstar() && type_a.is_Delta())) &&
             type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
    // NN → ΔN*
    return 7. / msqr;
  } else if (((type_a.is_Deltastar() && type_b.is_Delta()) ||
              (type_b.is_Deltastar() && type_a.is_Delta())) &&
             type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
    // NN → ΔΔ*
    if (twoI == 2) {
      return 15. / msqr;
    } else if (twoI == 0) {
      return 25. / msqr;
    }
  } else if ((type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
             (type_b.is_deuteron() && type_a.pdgcode().is_pion())) {
    /* This parametrization is the result of fitting d+pi->NN cross-section.
     * Already Breit-Wigner-like part provides a good fit, exponential fixes
     * behaviour around the treshold. The d+pi experimental cross-section
     * was taken from Fig. 2 of [\iref{Tanabe:1987vg}]. */
    return 0.055 / (pow_int(sqrts - 2.145, 2) + pow_int(0.065, 2)) *
           (1.0 - std::exp(-(sqrts - 2.0) * 20.0));
  }
 
  // all cases not listed: zero!
  return 0.;
}

Here is the call graph for this function:

Here is the caller graph for this function:

find_nn_xsection_from_type()

template<class IntegrationMethod >

CollisionBranchList smash::CrossSections::find_nn_xsection_from_type ( const ParticleTypePtrList &  type_res_1, const ParticleTypePtrList &  type_res_2, const IntegrationMethod  integrator  ) const

private

Utility function to avoid code replication in nn_xx().

Parameters

[in]	type_res_1	List of possible first final resonance types
[in]	type_res_2	List of possible second final resonance types
[in]	integrator	Used to integrate over the kinematically allowed mass range of the Breit-Wigner distribution

Returns

List of all possible NN reactions with their cross sections with different final states

Cross section for 2->2 process with 1/2 resonance(s) in final state. Based on eq. (46) in Weil:2013mya [72] and eq. (3.29) in Bass:1998ca [8]

Definition at line 4666 of file crosssections.cc.

                                              {
  const ParticleType& type_particle_a = incoming_particles_[0].type();
  const ParticleType& type_particle_b = incoming_particles_[1].type();
 
  CollisionBranchList channel_list;
  const double s = sqrt_s_ * sqrt_s_;
 
  // Loop over specified first resonance list
  for (ParticleTypePtr type_res_1 : list_res_1) {
    // Loop over specified second resonance list
    for (ParticleTypePtr type_res_2 : list_res_2) {
      // Check for charge conservation.
      if (type_res_1->charge() + type_res_2->charge() !=
          type_particle_a.charge() + type_particle_b.charge()) {
        continue;
      }
 
      // loop over total isospin
      for (const int twoI : I_tot_range(type_particle_a, type_particle_b)) {
        const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
            type_particle_a, type_particle_b, *type_res_1, *type_res_2, twoI);
        // If Clebsch-Gordan coefficient is zero, don't bother with the rest.
        if (std::abs(isospin_factor) < really_small) {
          continue;
        }
 
        // Integration limits.
        const double lower_limit = type_res_1->min_mass_kinematic();
        const double upper_limit = sqrt_s_ - type_res_2->mass();
        /* Check the available energy (requiring it to be a little above the
         * threshold, because the integration will not work if it's too close).
         */
        if (upper_limit - lower_limit < 1E-3) {
          continue;
        }
 
        // Calculate matrix element.
        const double matrix_element = nn_to_resonance_matrix_element(
            sqrt_s_, *type_res_1, *type_res_2, twoI);
        if (matrix_element <= 0.) {
          continue;
        }
 
        /* Calculate resonance production cross section
         * using the Breit-Wigner distribution as probability amplitude.
         * Integrate over the allowed resonance mass range. */
        const double resonance_integral = integrator(*type_res_1, *type_res_2);
 
        /**
         * Cross section for 2->2 process with 1/2 resonance(s) in final state.
         * Based on eq. (46) in \iref{Weil:2013mya} and eq. (3.29) in
         * \iref{Bass:1998ca}
         */
        const double spin_factor =
            (type_res_1->spin() + 1) * (type_res_2->spin() + 1);
        const double xsection = isospin_factor * spin_factor * matrix_element *
                                resonance_integral / (s * cm_momentum());
 
        if (xsection > really_small) {
          channel_list.push_back(std::make_unique<CollisionBranch>(
              *type_res_1, *type_res_2, xsection, ProcessType::TwoToTwo));
          logg[LCrossSections].debug(
              "Found 2->2 creation process for resonance ", type_res_1, ", ",
              type_res_2);
          logg[LCrossSections].debug("2->2 with original particles: ",
                                     type_particle_a, type_particle_b);
        }
      }
    }
  }
  return channel_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

cm_momentum()

double smash::CrossSections::cm_momentum ( ) const

inlineprivate

Determine the momenta of the incoming particles in the center-of-mass system.

Returns

Center-of-mass momentum

Definition at line 735 of file crosssections.h.

                             {
    const double m1 = incoming_particles_[0].effective_mass();
    const double m2 = incoming_particles_[1].effective_mass();
    return pCM(sqrt_s_, m1, m2);
  }

Here is the call graph for this function:

Here is the caller graph for this function:

add_channel()

template<typename F >

void smash::CrossSections::add_channel ( CollisionBranchList &  process_list, F &&  get_xsection, double  sqrts, const ParticleType &  type_a, const ParticleType &  type_b  ) const

inlineprivate

Helper function: Add a 2-to-2 channel to a collision branch list given a cross section.

The cross section is only calculated if there is enough energy for the process. If the cross section is small, the branch is not added.

Definition at line 770 of file crosssections.h.

                                                     {
    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));
    }
  }

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

incoming_particles_

const ParticleList smash::CrossSections::incoming_particles_

private

List with data of scattering particles.

Definition at line 742 of file crosssections.h.

sqrt_s_

const double smash::CrossSections::sqrt_s_

private

Total energy in the center-of-mass frame.

Definition at line 745 of file crosssections.h.

potentials_

const std::pair<FourVector, FourVector> smash::CrossSections::potentials_

private

Potentials at the interacting point.

They are used to calculate the corrections on the threshold energies.

Definition at line 751 of file crosssections.h.

is_BBbar_pair_

const bool smash::CrossSections::is_BBbar_pair_

private

Whether incoming particles are a pair of a baryon and an antibaryon (could be different baryon types)

Definition at line 757 of file crosssections.h.

is_NNbar_pair_

const bool smash::CrossSections::is_NNbar_pair_

private

Whether incoming particles are a nulecon-antinucleon pair (same isospin)

Definition at line 760 of file crosssections.h.

---

The documentation for this class was generated from the following files:

• src/include/smash/crosssections.h
• src/crosssections.cc