smash::CrossSections Class Reference

#include <crosssections.h>

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

Definition at line 62 of file crosssections.h.

Collaboration diagram for smash::CrossSections:

Public Member Functions
	CrossSections (const ParticleList &incoming_particles, const double sqrt_s, const std::pair< FourVector, FourVector > potentials)
	Construct CrossSections instance. More...
CollisionBranchList	generate_collision_list (double elastic_parameter, bool two_to_one_switch, ReactionsBitSet included_2to2, double low_snn_cut, bool strings_switch, bool use_AQM, bool strings_with_probability, NNbarTreatment nnbar_treatment, StringProcess *string_process) const
	Generate a list of all possible collisions between the incoming particles with the given c.m. More...
CollisionBranchPtr	elastic (double elast_par, bool use_AQM) const
	Determine the elastic cross section for this collision. More...
CollisionBranchList	two_to_one () 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
	Return the 2-to-1 resonance production cross section for a given resonance. 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 (ReactionsBitSet included_2to2) const
	Find all inelastic 2->2 processes for the given scattering. More...
CollisionBranchList	string_excitation (double total_string_xs, StringProcess *string_process, bool use_AQM) 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 (const double current_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...
double	high_energy () 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 (bool strings_switch, bool use_transition_probability, bool use_AQM, bool treat_nnbar_with_strings) const
double	probability_transit_high (const double region_lower, const double region_upper) const

Private Member Functions
double	elastic_parametrization (bool use_AQM) 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...
CollisionBranchList	npi_yk () const
	Find all processes for Nucleon-Pion to Hyperon-Kaon Scattering. More...
CollisionBranchList	bb_xx_except_nn (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 (ReactionsBitSet included_2to2) const
	Find all inelastic 2->2 processes for Nucelon-Nucelon Scattering. More...
CollisionBranchList	nk_xx (ReactionsBitSet included_2to2) const
	Find all inelastic 2->2 background processes for Nucleon-Kaon (NK) Scattering. More...
CollisionBranchList	deltak_xx (ReactionsBitSet included_2to2) const
	Find all inelastic 2->2 processes for Delta-Kaon (DeltaK) Scattering. More...
CollisionBranchList	ypi_xx (ReactionsBitSet included_2to2) const
	Find all inelastic 2->2 processes for Hyperon-Pion (Ypi) Scattering. More...
CollisionBranchList	dpi_xx (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 (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 (const bool is_anti_particles) const
	Calculate cross sections for 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	nn_to_resonance_matrix_element (double sqrts, const ParticleType &type_a, const ParticleType &type_b, const 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 baryon-antibaryon pair. More...

Constructor & Destructor Documentation

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

Construct CrossSections instance.

Parameters

[in]	incoming_particles	Particles that are reacting.
[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 112 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()) {}

Member Function Documentation

CollisionBranchList smash::CrossSections::generate_collision_list	(	double	elastic_parameter,
		bool	two_to_one_switch,
		ReactionsBitSet	included_2to2,
		double	low_snn_cut,
		bool	strings_switch,
		bool	use_AQM,
		bool	strings_with_probability,
		NNbarTreatment	nnbar_treatment,
		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's not triggerd according to the probability. It will then be added in add_all_scatterings in scatteraction.cc

Parameters

[in]	elastic_parameter	Value of the constant global elastic cross section, if it is non-zero. The parametrized elastic cross section is used otherwise.
[in]	two_to_one_switch	2->1 reactions enabled?
[in]	included_2to2	Which 2->2 ractions are enabled?
[in]	low_snn_cut	Elastic collisions with CME below are forbidden.
[in]	strings_switch	Are string processes enabled?
[in]	use_AQM	Is the Additive Quark Model enabled?
[in]	strings_with_probability	Are string processes triggered according to a probability?
[in]	nnbar_treatment	NNbar treatment through resonance, strings or none
[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 123 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 (strings_with_probability) {
    p_pythia =
        string_probability(strings_switch, strings_with_probability, use_AQM,
                           nnbar_treatment == NNbarTreatment::Strings);
  }

  /* 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_ < low_snn_cut;
  bool incl_elastic = included_2to2[IncludedReactions::Elastic];
  if (incl_elastic && !reject_by_nucleon_elastic_cutoff) {
    process_list.emplace_back(elastic(elastic_parameter, use_AQM));
  }
  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., high_energy() - sig_current);
    append_list(process_list,
                string_excitation(sig_string, string_process, use_AQM),
                p_pythia);
    append_list(process_list, rare_two_to_two(), p_pythia);
  }
  if (p_pythia < 1.) {
    if (two_to_one_switch) {
      // resonance formation (2->1)
      append_list(process_list, two_to_one(), 1. - p_pythia);
    }
    if (included_2to2.any()) {
      // 2->2 (inelastic)
      append_list(process_list, two_to_two(included_2to2), 1. - p_pythia);
    }
  }
  /* 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 (nnbar_treatment == NNbarTreatment::Resonances) {
    if (included_2to2[IncludedReactions::NNbar] != 1) {
      throw std::runtime_error(
          "'NNbar' has to be in the list of allowed 2 to 2 processes "
          "to enable annihilation to go through resonances");
    }
    if (t1.is_nucleon() && t2.pdgcode() == t1.get_antiparticle()->pdgcode()) {
      /* 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)));
    }
    if ((t1.pdgcode() == pdg::rho_z && t2.pdgcode() == pdg::h1) ||
        (t1.pdgcode() == pdg::h1 && t2.pdgcode() == pdg::rho_z)) {
      append_list(process_list, NNbar_creation());
    }
  }
  return process_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CollisionBranchPtr smash::CrossSections::elastic	(	double	elast_par,
		bool	use_AQM
	)		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).

Parameters

[in]	elast_par	Elastic cross section parameter from the input file.
[in]	use_AQM	Whether to extend elastic cross-sections with AQM.

Returns

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

Definition at line 191 of file crosssections.cc.

                                                              {
  double elastic_xs = 0.;
  if (elast_par >= 0.) {
    // use constant elastic cross section from config file
    elastic_xs = elast_par;
  } else {
    // use parametrization
    elastic_xs = elastic_parametrization(use_AQM);
  }
  return make_unique<CollisionBranch>(incoming_particles_[0].type(),
                                      incoming_particles_[1].type(), elastic_xs,
                                      ProcessType::Elastic);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CollisionBranchList smash::CrossSections::two_to_one

(

)

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.

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 703 of file crosssections.cc.

                                                    {
  const auto& log = logger<LogArea::CrossSections>();
  CollisionBranchList resonance_process_list;
  const ParticleType& type_particle_a = incoming_particles_[0].type();
  const ParticleType& type_particle_b = incoming_particles_[1].type();

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

  // Find all the possible resonances
  for (const ParticleType& type_resonance : ParticleType::list_all()) {
    /* Not a resonance, go to next type of particle */
    if (type_resonance.is_stable()) {
      continue;
    }

    // Same resonance as in the beginning, ignore
    if ((!type_particle_a.is_stable() &&
         type_resonance.pdgcode() == type_particle_a.pdgcode()) ||
        (!type_particle_b.is_stable() &&
         type_resonance.pdgcode() == type_particle_b.pdgcode())) {
      continue;
    }

    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(make_unique<CollisionBranch>(
          type_resonance, resonance_xsection, ProcessType::TwoToOne));
      log.debug("Found resonance: ", type_resonance);
      log.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:

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

Return the 2-to-1 resonance production cross section for a given resonance.

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.

Calculate resonance production cross section using the Breit-Wigner distribution as probability amplitude. See Eq. (176) in Buss:2011mx.

Definition at line 743 of file crosssections.cc.

                                                              {
  const ParticleType& type_particle_a = incoming_particles_[0].type();
  const ParticleType& type_particle_b = incoming_particles_[1].type();
  // Check for charge conservation.
  if (type_resonance.charge() !=
      type_particle_a.charge() + type_particle_b.charge()) {
    return 0.;
  }

  // Check for baryon-number conservation.
  if (type_resonance.baryon_number() !=
      type_particle_a.baryon_number() + type_particle_b.baryon_number()) {
    return 0.;
  }

  // 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.;
  }

  // Calculate spin factor
  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.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:

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 possibe rare 2->2 processes.

Definition at line 206 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:

CollisionBranchList smash::CrossSections::two_to_two

(

ReactionsBitSet

included_2to2

)

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?

Returns

List of all possibe inelastic 2->2 processes.

Definition at line 780 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);
    } 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 (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:

CollisionBranchList smash::CrossSections::string_excitation	(	double	total_string_xs,
		StringProcess *	string_process,
		bool	use_AQM
	)		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]	use_AQM	whether to extend string cross-sections with AQM

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 2091 of file crosssections.cc.

                                                                               {
  const auto& log = logger<LogArea::CrossSections>();

  if (!string_process) {
    throw std::runtime_error("string_process should be initialized.");
  }

  CollisionBranchList channel_list;
  if (total_string_xs <= 0.) {
    return channel_list;
  }

  /* 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_factor = 1.;
  for (int i = 0; i < 2; i++) {
    PdgCode pdg = incoming_particles_[i].type().pdgcode();
    pdgid[i] = StringProcess::pdg_map_for_pythia(pdg);
    AQM_factor *= (1. - 0.4 * pdg.frac_strange());
  }

  /* 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;
      }
    }
  }

  /* 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 =
      string_process->cross_sections_diffractive(pdgid[0], pdgid[1], sqrt_s_);
  if (use_AQM) {
    for (int ip = 0; ip < 3; ip++) {
      xs[ip] *= AQM_factor;
    }
  }
  double single_diffr_AX = xs[0], single_diffr_XB = xs[1], double_diffr = xs[2];
  double single_diffr = single_diffr_AX + single_diffr_XB;
  double diffractive = single_diffr + double_diffr;

  /* 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(). */
    double xs_param = xs_ppbar_annihilation(sqrt_s_ * sqrt_s_);
    if (use_AQM) {
      xs_param *= AQM_factor;
    }
    sig_annihilation = std::min(total_string_xs, xs_param);
  }

  const double nondiffractive_all =
      std::max(0., total_string_xs - sig_annihilation - diffractive);
  diffractive = total_string_xs - sig_annihilation - nondiffractive_all;
  double_diffr = std::max(0., diffractive - single_diffr);
  const double a = (diffractive - double_diffr) / single_diffr;
  single_diffr_AX *= a;
  single_diffr_XB *= a;
  assert(std::abs(single_diffr_AX + single_diffr_XB + double_diffr +
                  sig_annihilation + nondiffractive_all - total_string_xs) <
         1.e-6);

  double nondiffractive_soft = 0.;
  double nondiffractive_hard = 0.;
  if (nondiffractive_all > 0.) {
    /* Hard string process is added by hard cross section
     * in conjunction with multipartion interaction picture
     * \iref{Sjostrand:1987su}. */
    const double hard_xsec = AQM_factor * string_hard_cross_section();
    nondiffractive_soft =
        nondiffractive_all * std::exp(-hard_xsec / nondiffractive_all);
    nondiffractive_hard = nondiffractive_all - nondiffractive_soft;
  }
  log.debug("String cross sections [mb] are");
  log.debug("Single-diffractive AB->AX: ", single_diffr_AX);
  log.debug("Single-diffractive AB->XB: ", single_diffr_XB);
  log.debug("Double-diffractive AB->XX: ", double_diffr);
  log.debug("Soft non-diffractive: ", nondiffractive_soft);
  log.debug("Hard non-diffractive: ", nondiffractive_hard);
  log.debug("B-Bbar annihilation: ", sig_annihilation);

  /* cross section of soft string excitation including annihilation */
  const double sig_string_soft = total_string_xs - nondiffractive_hard;

  /* fill the list of process channels */
  if (sig_string_soft > 0.) {
    channel_list.push_back(make_unique<CollisionBranch>(
        single_diffr_AX, ProcessType::StringSoftSingleDiffractiveAX));
    channel_list.push_back(make_unique<CollisionBranch>(
        single_diffr_XB, ProcessType::StringSoftSingleDiffractiveXB));
    channel_list.push_back(make_unique<CollisionBranch>(
        double_diffr, ProcessType::StringSoftDoubleDiffractive));
    channel_list.push_back(make_unique<CollisionBranch>(
        nondiffractive_soft, ProcessType::StringSoftNonDiffractive));
    if (can_annihilate) {
      channel_list.push_back(make_unique<CollisionBranch>(
          sig_annihilation, ProcessType::StringSoftAnnihilation));
    }
  }
  if (nondiffractive_hard > 0.) {
    channel_list.push_back(make_unique<CollisionBranch>(
        nondiffractive_hard, ProcessType::StringHard));
  }
  return channel_list;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CollisionBranchPtr smash::CrossSections::NNbar_annihilation

(

const double

current_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

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 2315 of file crosssections.cc.

                                   {
  const auto& log = logger<LogArea::CrossSections>();
  /* 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) - current_xs);
  log.debug("NNbar cross section is: ", nnbar_xsec);
  // Make collision channel NNbar -> ρh₁(1170); eventually decays into 5π
  return 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:

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_cross_section

Returns

Collision Branch with NNbar creation process and its cross section

Definition at line 2329 of file crosssections.cc.

                                                        {
  const auto& log = logger<LogArea::CrossSections>();
  CollisionBranchList channel_list;
  /* 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, incoming_particles_[0].type(),
                        incoming_particles_[1].type(), type_N, type_Nbar) *
                    std::max(0., ppbar_total(s) - ppbar_elastic(s));
  log.debug("NNbar reverse cross section is: ", xsection);
  channel_list.push_back(make_unique<CollisionBranch>(
      type_N, type_Nbar, xsection, ProcessType::TwoToTwo));
  channel_list.push_back(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:

double smash::CrossSections::high_energy

(

)

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.

Definition at line 2228 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()) {
    if (pdg_a == pdg_b) {
      xs = pp_high_energy(s);  // pp, nn
    } else if (pdg_a.antiparticle_sign() * pdg_b.antiparticle_sign() == 1) {
      xs = np_high_energy(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(s);
      double xs_h = 0.;
      if (pdg_a.is_antiparticle_of(pdg_b)) {
        xs_h = ppbar_high_energy(s);  // ppbar, nnbar
      } else {
        xs_h = npbar_high_energy(s);  // npbar, nbarp
      }
      /* Transition between low and high energy is set to be consistent with
       * that defined in string_probability(). */
      double region_lower = transit_high_energy::sqrts_range_NN[0];
      double region_upper = transit_high_energy::sqrts_range_NN[1];
      double prob_high = probability_transit_high(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 *= (1. - 0.4 * pdg_a.frac_strange()) * (1. - 0.4 * pdg_b.frac_strange());

  return xs;
}

Here is the call graph for this function:

Here is the caller graph for this function:

double smash::CrossSections::string_probability	(	bool	strings_switch,
		bool	use_transition_probability,
		bool	use_AQM,
		bool	treat_nnbar_with_strings
	)		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). 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]	strings_switch	Is string fragmentation enabled?
[in]	use_transition_probability	which algorithm to use for string treatment (see Switch_on_String_with_Probability)
[in]	use_AQM	whether AQM is activated
[in]	treat_nnbar_with_strings	use strings for nnbar treatment?

Definition at line 2587 of file crosssections.cc.

                                                                              {
  /* string fragmentation is enabled when strings_switch is on and the process
   * is included in pythia. */
  if (!strings_switch) {
    return 0.;
  }

  const ParticleType& t1 = incoming_particles_[0].type();
  const ParticleType& t2 = incoming_particles_[1].type();

  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_ && 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 =
      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_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 = transit_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 = transit_high_energy::KN_offset;
    } else if (pdg1.is_pion() && pdg2.is_pion()) {
      aqm_offset = transit_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 (!use_transition_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) {
      region_lower = transit_high_energy::sqrts_range_Npi[0];
      region_upper = transit_high_energy::sqrts_range_Npi[1];
    } else if (is_NN_scattering) {
      region_lower = transit_high_energy::sqrts_range_NN[0];
      region_upper = transit_high_energy::sqrts_range_NN[1];
    } 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 + transit_high_energy::sqrts_range;
    }

    if (sqrt_s_ > region_upper) {
      return 1.;
    } else if (sqrt_s_ < region_lower) {
      return 0.;
    } else {
      // Rescale transition region to [-1, 1]
      return probability_transit_high(region_lower, region_upper);
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

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

Parameters

[in]	region_lower	the lowest sqrts in the transition region [GeV]
[in]	region_upper	the highest sqrts in the transition region [GeV]

Returns

probability to have the high energy interaction (via string)

Definition at line 2677 of file crosssections.cc.

                                                                {
  if (sqrt_s_ < region_lower) {
    return 0.;
  }

  if (sqrt_s_ > region_upper) {
    return 1.;
  }

  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:

double smash::CrossSections::elastic_parametrization ( bool  use_AQM ) const

private

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

Parameters

[in]

use_AQM

whether AQM is activated

Returns

Elastic cross section

Definition at line 219 of file crosssections.cc.

                                                                {
  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_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 =
        std::abs(pdg_nucleus.get_decimal()) == pdg::decimal_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 (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_high_energy(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 + transit_high_energy::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 *=
        (1. - 0.4 * pdg_a.frac_strange()) * (1. - 0.4 * pdg_b.frac_strange());
  }
  return elastic_xs;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 280 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_;

  // Use parametrized cross sections.
  double sig_el = -1.;
  if (pdg_a.antiparticle_sign() == -pdg_b.antiparticle_sign()) {
    sig_el = ppbar_elastic(s);
  } else if (pdg_a.is_nucleon() && pdg_b.is_nucleon()) {
    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();
    sig_el = pp_elastic_high_energy(s, m1, m2);
  }
  if (sig_el > 0.) {
    return sig_el;
  } else {
    std::stringstream ss;
    const auto name_a = incoming_particles_[0].type().name();
    const auto name_b = incoming_particles_[1].type().name();
    ss << "problem in CrossSections::elastic: a=" << name_a << " b=" << name_b
       << " j_a=" << pdg_a.spin() << " j_b=" << pdg_b.spin()
       << " sigma=" << sig_el << " s=" << s;
    throw std::runtime_error(ss.str());
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 311 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 std::runtime_error(
          "only the elastic cross section for proton-pion "
          "is implemented");
  }

  if (sig_el > 0) {
    return sig_el;
  } else {
    std::stringstream ss;
    const auto name_a = incoming_particles_[0].type().name();
    const auto name_b = incoming_particles_[1].type().name();
    ss << "problem in CrossSections::elastic: a=" << name_a << " b=" << name_b
       << " j_a=" << pdg_a.spin() << " j_b=" << pdg_b.spin()
       << " sigma=" << sig_el << " s=" << s;
    throw std::runtime_error(ss.str());
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 608 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 std::runtime_error(
          "elastic cross section for antinucleon-kaon "
          "not implemented");
  }

  if (sig_el > 0) {
    return sig_el;
  } else {
    std::stringstream ss;
    const auto name_a = incoming_particles_[0].type().name();
    const auto name_b = incoming_particles_[1].type().name();
    ss << "problem in CrossSections::elastic: a=" << name_a << " b=" << name_b
       << " j_a=" << pdg_a.spin() << " j_b=" << pdg_b.spin()
       << " sigma=" << sig_el << " s=" << s;
    throw std::runtime_error(ss.str());
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 394 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:

CollisionBranchList smash::CrossSections::bb_xx_except_nn ( 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 830 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:

CollisionBranchList smash::CrossSections::nn_xx ( 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 857 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, 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 deutron =
      ParticleType::try_find(PdgCode::from_decimal(pdg::decimal_d));
  ParticleTypePtr antideutron =
      ParticleType::try_find(PdgCode::from_decimal(pdg::decimal_antid));
  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 (deutron && antideutron && pim && pi0 && pip &&
      included_2to2[IncludedReactions::PiDeuteron_to_NN] == 1) {
    const ParticleTypePtrList deutron_list = {deutron};
    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:

CollisionBranchList smash::CrossSections::nk_xx ( ReactionsBitSet  included_2to2 ) const

private

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

Parameters

[in]

included_2to2

Which 2->2 reactions are enabled?

Returns

List of all possible NK reactions with their cross sections

Definition at line 929 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 = transit_high_energy::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:

CollisionBranchList smash::CrossSections::deltak_xx ( 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 1412 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:

CollisionBranchList smash::CrossSections::ypi_xx ( 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 1561 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:

CollisionBranchList smash::CrossSections::dpi_xx ( 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 1917 of file crosssections.cc.

                                                                             {
  const auto& log = logger<LogArea::ScatterAction>();
  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(make_unique<CollisionBranch>(
                *nuc_a, *nuc_b, xsection, ProcessType::TwoToTwo));
            log.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();
    const double s = sqrt_s_ * sqrt_s_;
    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;
      }
      // 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_momentum());
      if (produced_nucleus->is_stable()) {
        assert(!type_nucleus.is_stable());
        xsection *= pCM_from_s(s, type_pi.mass(), produced_nucleus->mass());
      } else {
        assert(type_nucleus.is_stable());
        const double resonance_integral =
            produced_nucleus->iso_multiplet()->get_integral_piR(sqrts);
        xsection *= resonance_integral;
        log.debug("Resonance integral ", resonance_integral,
                  ", matrix element: ", matrix_element,
                  ", cm_momentum: ", cm_momentum());
      }
      process_list.push_back(make_unique<CollisionBranch>(
          type_pi, *produced_nucleus, xsection, ProcessType::TwoToTwo));
      log.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:

CollisionBranchList smash::CrossSections::dn_xx ( 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 2031 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();
  const double s = sqrt_s_ * sqrt_s_;
  const double sqrts = sqrt_s_;

  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;
    }
    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
       * [\iref{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_momentum());
    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;
    }
    process_list.push_back(make_unique<CollisionBranch>(
        type_N, *produced_nucleus, xsection, ProcessType::TwoToTwo));
    const auto& log = logger<LogArea::ScatterAction>();
    log.debug(type_N.name(), type_nucleus.name(), "→ ", type_N.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:

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 2291 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.
    cross_sec = NN_string_hard(sqrt_s_ * sqrt_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:

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

private

Calculate cross sections for resonance absorption (i.e.

NR->NN and ΔR->NN).

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.

Cross section for 2->2 resonance absorption, obtained via detailed balance from the inverse reaction. See eqs. (B.6), (B.9) and (181) in Buss:2011mx. There are factors for spin, isospin and symmetry involved.

Definition at line 2358 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;
        }

        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(make_unique<CollisionBranch>(
              *nuc_a, *nuc_b, xsection, ProcessType::TwoToTwo));
          const auto& log = logger<LogArea::CrossSections>();
          log.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:

double smash::CrossSections::nn_to_resonance_matrix_element ( double  sqrts, const ParticleType &  type_a, const ParticleType &  type_b, const 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]

All other processes use a constant matrix element, similar to Bass:1998ca, equ. (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]. 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], eq. 29.

Definition at line 2425 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 of 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.;
  }
  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);
  } 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;
      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:

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 and Eq. (3.29) in Bass:1998ca

Definition at line 2513 of file crosssections.cc.

                                              {
  const ParticleType& type_particle_a = incoming_particles_[0].type();
  const ParticleType& type_particle_b = incoming_particles_[1].type();

  const auto& log = logger<LogArea::CrossSections>();
  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);

        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(make_unique<CollisionBranch>(
              *type_res_1, *type_res_2, xsection, ProcessType::TwoToTwo));
          log.debug("Found 2->2 creation process for resonance ", type_res_1,
                    ", ", type_res_2);
          log.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:

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 458 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:

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 487 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(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

const ParticleList smash::CrossSections::incoming_particles_

private

List with data of scattering particles.

Definition at line 465 of file crosssections.h.

const double smash::CrossSections::sqrt_s_

private

Total energy in the center-of-mass frame.

Definition at line 468 of file crosssections.h.

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 474 of file crosssections.h.

const bool smash::CrossSections::is_BBbar_pair_

private

Whether incoming particles are a baryon-antibaryon pair.

Definition at line 477 of file crosssections.h.

---

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

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