Version: SMASH-3.1
crosssections.cc
Go to the documentation of this file.
1 /*
2  *
3  * Copyright (c) 2018-2024
4  * SMASH Team
5  *
6  * GNU General Public License (GPLv3 or later)
7  *
8  */
9 
10 #include "smash/crosssections.h"
11 
12 #include "smash/clebschgordan.h"
13 #include "smash/constants.h"
14 #include "smash/logging.h"
15 #include "smash/parametrizations.h"
16 #include "smash/pow.h"
17 
18 namespace smash {
19 static constexpr int LCrossSections = LogArea::CrossSections::id;
20 static constexpr int LScatterAction = LogArea::ScatterAction::id;
21 
28 static double detailed_balance_factor_stable(double s, const ParticleType& a,
29  const ParticleType& b,
30  const ParticleType& c,
31  const ParticleType& d) {
32  double spin_factor = (c.spin() + 1) * (d.spin() + 1);
33  spin_factor /= (a.spin() + 1) * (b.spin() + 1);
34  double symmetry_factor = (1 + (a == b));
35  symmetry_factor /= (1 + (c == d));
36  const double momentum_factor = pCM_sqr_from_s(s, c.mass(), d.mass()) /
37  pCM_sqr_from_s(s, a.mass(), b.mass());
38  return spin_factor * symmetry_factor * momentum_factor;
39 }
40 
47 static double detailed_balance_factor_RK(double sqrts, double pcm,
48  const ParticleType& a,
49  const ParticleType& b,
50  const ParticleType& c,
51  const ParticleType& d) {
52  assert(!a.is_stable());
53  assert(b.pdgcode().is_kaon());
54  double spin_factor = (c.spin() + 1) * (d.spin() + 1);
55  spin_factor /= (a.spin() + 1) * (b.spin() + 1);
56  double symmetry_factor = (1 + (a == b));
57  symmetry_factor /= (1 + (c == d));
58  const double momentum_factor =
59  pCM_sqr(sqrts, c.mass(), d.mass()) /
60  (pcm * a.iso_multiplet()->get_integral_RK(sqrts));
61  return spin_factor * symmetry_factor * momentum_factor;
62 }
63 
70 static double detailed_balance_factor_RR(double sqrts, double pcm,
71  const ParticleType& a,
72  const ParticleType& b,
73  const ParticleType& c,
74  const ParticleType& d) {
75  assert(!a.is_stable());
76  assert(!b.is_stable());
77  double spin_factor = (c.spin() + 1) * (d.spin() + 1);
78  spin_factor /= (a.spin() + 1) * (b.spin() + 1);
79  double symmetry_factor = (1 + (a == b));
80  symmetry_factor /= (1 + (c == d));
81  const double momentum_factor =
82  pCM_sqr(sqrts, c.mass(), d.mass()) /
83  (pcm * a.iso_multiplet()->get_integral_RR(b.iso_multiplet(), sqrts));
84  return spin_factor * symmetry_factor * momentum_factor;
85 }
86 
91 static void append_list(CollisionBranchList& main_list,
92  CollisionBranchList in_list, double weight = 1.) {
93  main_list.reserve(main_list.size() + in_list.size());
94  for (auto& proc : in_list) {
95  proc->set_weight(proc->weight() * weight);
96  main_list.emplace_back(std::move(proc));
97  }
98 }
99 
100 CrossSections::CrossSections(const ParticleList& incoming_particles,
101  const double sqrt_s,
102  const std::pair<FourVector, FourVector> potentials)
103  : incoming_particles_(incoming_particles),
104  sqrt_s_(sqrt_s),
105  potentials_(potentials),
106  is_BBbar_pair_(incoming_particles_[0].type().is_baryon() &&
107  incoming_particles_[1].type().is_baryon() &&
108  incoming_particles_[0].type().antiparticle_sign() ==
109  -incoming_particles_[1].type().antiparticle_sign()),
110  is_NNbar_pair_(
111  incoming_particles_[0].type().is_nucleon() &&
112  incoming_particles_[1].pdgcode() ==
113  incoming_particles_[0].type().get_antiparticle()->pdgcode()) {}
114 
115 CollisionBranchList CrossSections::generate_collision_list(
116  const ScatterActionsFinderParameters& finder_parameters,
117  StringProcess* string_process) const {
118  CollisionBranchList process_list;
119  const ParticleType& t1 = incoming_particles_[0].type();
120  const ParticleType& t2 = incoming_particles_[1].type();
121 
122  double p_pythia = 0.;
123  if (finder_parameters.strings_with_probability) {
124  p_pythia = string_probability(finder_parameters);
125  }
126 
127  /* Elastic collisions between two nucleons with sqrt_s below
128  * low_snn_cut can not happen. */
129  const bool reject_by_nucleon_elastic_cutoff =
130  t1.is_nucleon() && t2.is_nucleon() &&
131  t1.antiparticle_sign() == t2.antiparticle_sign() &&
132  sqrt_s_ < finder_parameters.low_snn_cut;
133  bool incl_elastic =
134  finder_parameters.included_2to2[IncludedReactions::Elastic];
135  if (incl_elastic && !reject_by_nucleon_elastic_cutoff) {
136  process_list.emplace_back(elastic(finder_parameters));
137  }
138  if (p_pythia > 0.) {
139  /* String-excitation cross section =
140  * Parametrized total cross - the contributions
141  * from all other present channels. */
142  const double sig_current = sum_xs_of(process_list);
143  const double sig_string = std::max(
144  0., finder_parameters.scale_xs *
145  high_energy(finder_parameters.transition_high_energy) -
146  sig_current);
147  append_list(process_list,
148  string_excitation(sig_string, string_process,
149  finder_parameters.use_AQM),
150  p_pythia);
151  append_list(process_list, rare_two_to_two(),
152  p_pythia * finder_parameters.scale_xs);
153  }
154  if (p_pythia < 1.) {
155  if (finder_parameters.two_to_one) {
156  // resonance formation (2->1)
157  append_list(process_list, two_to_one(),
158  (1. - p_pythia) * finder_parameters.scale_xs);
159  }
160  if (finder_parameters.included_2to2.any()) {
161  // 2->2 (inelastic)
162  append_list(
163  process_list,
164  two_to_two(finder_parameters.included_2to2,
165  finder_parameters.transition_high_energy.KN_offset),
166  (1. - p_pythia) * finder_parameters.scale_xs);
167  }
168  if (finder_parameters
169  .included_multi[IncludedMultiParticleReactions::Deuteron_3to2] ==
170  1) {
171  // 2->3 (deuterons only 2-to-3 reaction at the moment)
172  append_list(process_list, two_to_three(),
173  (1. - p_pythia) * finder_parameters.scale_xs);
174  }
175  if (finder_parameters
176  .included_multi[IncludedMultiParticleReactions::A3_Nuclei_4to2] ==
177  1) {
178  // 2->4
179  append_list(process_list, two_to_four(),
180  (1. - p_pythia) * finder_parameters.scale_xs);
181  }
182  }
183  if (finder_parameters.nnbar_treatment == NNbarTreatment::TwoToFive &&
184  is_NNbar_pair_) {
185  // NNbar directly to 5 pions (2-to-5)
186  process_list.emplace_back(NNbar_to_5pi(finder_parameters.scale_xs));
187  }
188 
189  /* NNbar annihilation thru NNbar → ρh₁(1170); combined with the decays
190  * ρ → ππ and h₁(1170) → πρ, this gives a final state of 5 pions.
191  * Only use in cases when detailed balance MUST happen, i.e. in a box! */
192  if (finder_parameters.nnbar_treatment == NNbarTreatment::Resonances) {
193  if (is_NNbar_pair_) {
194  /* Has to be called after the other processes are already determined,
195  * so that the sum of the cross sections includes all other processes. */
196  process_list.emplace_back(NNbar_annihilation(sum_xs_of(process_list),
197  finder_parameters.scale_xs));
198  } else {
199  append_list(process_list, NNbar_creation(), finder_parameters.scale_xs);
200  }
201  }
202  return process_list;
203 }
204 
205 double CrossSections::parametrized_total(
206  const ScatterActionsFinderParameters& finder_parameters) const {
207  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
208  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
209  double total_xs = 0.;
210  if (pdg_a.is_baryon() && pdg_b.is_baryon() &&
211  sqrt_s_ > finder_parameters.low_snn_cut) {
212  if (pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
213  // NN
214  total_xs = (pdg_a == pdg_b) ? pp_total(sqrt_s_ * sqrt_s_)
215  : np_total(sqrt_s_ * sqrt_s_);
216  } else {
217  // NNbar
218  total_xs = ppbar_total(sqrt_s_ * sqrt_s_);
219  }
220  total_xs *=
221  (1 - 0.4 * pdg_a.frac_strange()) * (1 - 0.4 * pdg_b.frac_strange());
222  } else if ((pdg_a.is_baryon() && pdg_b.is_meson()) ||
223  (pdg_a.is_meson() && pdg_b.is_baryon())) {
224  const PdgCode& meson = pdg_a.is_meson() ? pdg_a : pdg_b;
225  const PdgCode& baryon = pdg_a.is_meson() ? pdg_b : pdg_a;
226  if (meson.is_kaon() && baryon.is_nucleon()) {
227  if ((meson.code() == pdg::K_p && baryon.code() == pdg::p) ||
228  (meson.code() == pdg::K_z && baryon.code() == pdg::n) ||
229  (meson.code() == pdg::K_m && baryon.code() == -pdg::p) ||
230  (meson.code() == pdg::Kbar_z && baryon.code() == -pdg::n)) {
231  // K⁺p, K⁰n, and anti-processes
232  total_xs = kplusp_total(sqrt_s_ * sqrt_s_);
233  } else if ((meson.code() == pdg::K_p && baryon.code() == -pdg::p) ||
234  (meson.code() == pdg::K_z && baryon.code() == -pdg::n) ||
235  (meson.code() == pdg::K_m && baryon.code() == pdg::p) ||
236  (meson.code() == pdg::Kbar_z && baryon.code() == pdg::n)) {
237  // K⁻p, K̅⁰n, and anti-processes
238  total_xs = kminusp_total(sqrt_s_ * sqrt_s_);
239  } else if ((meson.code() == pdg::K_p && baryon.code() == pdg::n) ||
240  (meson.code() == pdg::K_z && baryon.code() == pdg::p) ||
241  (meson.code() == pdg::K_m && baryon.code() == -pdg::n) ||
242  (meson.code() == pdg::Kbar_z && baryon.code() == -pdg::p)) {
243  // K⁺n, K⁰p, and anti-processes
244  total_xs = kplusn_total(sqrt_s_ * sqrt_s_);
245  } else if ((meson.code() == pdg::K_p && baryon.code() == -pdg::n) ||
246  (meson.code() == pdg::K_z && baryon.code() == -pdg::p) ||
247  (meson.code() == pdg::K_m && baryon.code() == pdg::n) ||
248  (meson.code() == pdg::Kbar_z && baryon.code() == pdg::p)) {
249  // K⁻n, K̅⁰p and anti-processes
250  total_xs = kminusn_total(sqrt_s_ * sqrt_s_);
251  }
252  } else if (meson.is_pion() && baryon.is_nucleon()) {
253  // π⁺(p,nbar), π⁻(n,pbar)
254  if ((meson.code() == pdg::pi_p &&
255  (baryon.code() == pdg::p || baryon.code() == -pdg::n)) ||
256  (meson.code() == pdg::pi_m &&
257  (baryon.code() == pdg::n || baryon.code() == -pdg::p))) {
258  total_xs = piplusp_total(sqrt_s_);
259  } else if (meson.code() == pdg::pi_z) {
260  // π⁰N
261  total_xs = 0.5 * (piplusp_total(sqrt_s_) + piminusp_total(sqrt_s_));
262  } else {
263  // π⁻(p,nbar), π⁺(n,pbar)
264  total_xs = piminusp_total(sqrt_s_);
265  }
266  } else {
267  // M*+B* goes to AQM high energy π⁻p
268  total_xs = piminusp_high_energy(sqrt_s_ * sqrt_s_) *
269  (1 - 0.4 * pdg_a.frac_strange()) *
270  (1 - 0.4 * pdg_b.frac_strange());
271  }
272  } else if (pdg_a.is_meson() && pdg_b.is_meson()) {
273  if (pdg_a.is_pion() && pdg_b.is_pion()) {
274  switch (pdg_a.isospin3() * pdg_b.isospin3() / 4) {
275  // π⁺π⁻
276  case -1:
277  total_xs = pipluspiminus_total(sqrt_s_);
278  break;
279  case 0:
280  // π⁰π⁰
281  if (pdg_a.isospin3() + pdg_b.isospin3() == 0) {
282  total_xs = pizeropizero_total(sqrt_s_);
283  } else {
284  // π⁺π⁰: similar to π⁺π⁻
285  total_xs = pipluspiminus_total(sqrt_s_);
286  }
287  break;
288  // π⁺π⁺ goes to π⁻p AQM
289  case 1:
290  total_xs = (2. / 3.) * piminusp_high_energy(sqrt_s_ * sqrt_s_);
291  break;
292  default:
293  throw std::runtime_error("wrong isospin in ππ scattering");
294  }
295  } else {
296  // M*+M* goes to AQM high energy π⁻p
297  total_xs = (2. / 3.) * piminusp_high_energy(sqrt_s_ * sqrt_s_) *
298  (1 - 0.4 * pdg_a.frac_strange()) *
299  (1 - 0.4 * pdg_b.frac_strange());
300  }
301  }
302  return (total_xs + finder_parameters.additional_el_xs) *
303  finder_parameters.scale_xs;
304 }
305 
306 CollisionBranchPtr CrossSections::elastic(
307  const ScatterActionsFinderParameters& finder_parameters) const {
308  double elastic_xs = 0.;
309  if (finder_parameters.elastic_parameter >= 0.) {
310  // use constant elastic cross section from config file
311  elastic_xs = finder_parameters.elastic_parameter;
312  } else {
313  // use parametrization
314  elastic_xs = elastic_parametrization(
315  finder_parameters.use_AQM,
316  finder_parameters.transition_high_energy.pipi_offset);
317  }
318  /* when using a factor to scale the cross section and an additional
319  * contribution to the elastic cross section, the contribution is added first
320  * and then everything is scaled */
321  return std::make_unique<CollisionBranch>(
322  incoming_particles_[0].type(), incoming_particles_[1].type(),
323  (elastic_xs + finder_parameters.additional_el_xs) *
324  finder_parameters.scale_xs,
325  ProcessType::Elastic);
326 }
327 
328 CollisionBranchList CrossSections::rare_two_to_two() const {
329  CollisionBranchList process_list;
330  const ParticleData& data_a = incoming_particles_[0];
331  const ParticleData& data_b = incoming_particles_[1];
332  const auto& pdg_a = data_a.pdgcode();
333  const auto& pdg_b = data_b.pdgcode();
334  if ((pdg_a.is_nucleon() && pdg_b.is_pion()) ||
335  (pdg_b.is_nucleon() && pdg_a.is_pion())) {
336  process_list = npi_yk();
337  }
338  return process_list;
339 }
340 
341 double CrossSections::elastic_parametrization(const bool use_AQM,
342  const double pipi_offset) const {
343  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
344  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
345  double elastic_xs = 0.0;
346  if ((pdg_a.is_nucleon() && pdg_b.is_pion()) ||
347  (pdg_b.is_nucleon() && pdg_a.is_pion())) {
348  // Elastic Nucleon Pion Scattering
349  elastic_xs = npi_el();
350  } else if ((pdg_a.is_nucleon() && pdg_b.is_kaon()) ||
351  (pdg_b.is_nucleon() && pdg_a.is_kaon())) {
352  // Elastic Nucleon Kaon Scattering
353  elastic_xs = nk_el();
354  } else if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
355  pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
356  // Elastic Nucleon Nucleon Scattering
357  elastic_xs = nn_el();
358  } else if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
359  pdg_a.antiparticle_sign() == -pdg_b.antiparticle_sign()) {
360  // Elastic Nucleon anti-Nucleon Scattering
361  elastic_xs = ppbar_elastic(sqrt_s_ * sqrt_s_);
362  } else if (pdg_a.is_nucleus() || pdg_b.is_nucleus()) {
363  const PdgCode& pdg_nucleus = pdg_a.is_nucleus() ? pdg_a : pdg_b;
364  const PdgCode& pdg_other = pdg_a.is_nucleus() ? pdg_b : pdg_a;
365  const bool is_deuteron = pdg_nucleus.is_deuteron(); // d or anti-d
366  if (is_deuteron && pdg_other.is_pion()) {
367  // Elastic (Anti-)deuteron Pion Scattering
368  elastic_xs = deuteron_pion_elastic(sqrt_s_ * sqrt_s_);
369  } else if (is_deuteron && pdg_other.is_nucleon()) {
370  // Elastic (Anti-)deuteron (Anti-)Nucleon Scattering
371  elastic_xs = deuteron_nucleon_elastic(sqrt_s_ * sqrt_s_);
372  }
373  } else if (use_AQM) {
374  const double m1 = incoming_particles_[0].effective_mass();
375  const double m2 = incoming_particles_[1].effective_mass();
376  const double s = sqrt_s_ * sqrt_s_;
377  if (pdg_a.is_baryon() && pdg_b.is_baryon()) {
378  elastic_xs = nn_el(); // valid also for annihilation
379  } else if ((pdg_a.is_meson() && pdg_b.is_baryon()) ||
380  (pdg_b.is_meson() && pdg_a.is_baryon())) {
381  elastic_xs = piplusp_elastic_high_energy(s, m1, m2);
382  } else if (pdg_a.is_meson() && pdg_b.is_meson()) {
383  /* Special case: the pi+pi- elastic cross-section goes through resonances
384  * at low sqrt_s, so we turn it off for this region so as not to destroy
385  * the agreement with experimental data; this does not
386  * apply to other pi pi cross-sections, which do not have any data */
387  if (((pdg_a == pdg::pi_p && pdg_b == pdg::pi_m) ||
388  (pdg_a == pdg::pi_m && pdg_b == pdg::pi_p)) &&
389  (m1 + m2 + pipi_offset) > sqrt_s_) {
390  elastic_xs = 0.0;
391  } else {
392  // meson-meson goes through scaling from π+p parametrization
393  elastic_xs = 2. / 3. * piplusp_elastic_AQM(s, m1, m2);
394  }
395  }
396  elastic_xs *=
397  (1. - 0.4 * pdg_a.frac_strange()) * (1. - 0.4 * pdg_b.frac_strange());
398  }
399  return elastic_xs;
400 }
401 
402 double CrossSections::nn_el() const {
403  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
404  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
405 
406  const double s = sqrt_s_ * sqrt_s_;
407 
408  // Use parametrized cross sections.
409  double sig_el = -1.;
410  if (pdg_a.antiparticle_sign() == -pdg_b.antiparticle_sign()) {
411  sig_el = ppbar_elastic(s);
412  } else if (pdg_a.is_nucleon() && pdg_b.is_nucleon()) {
413  sig_el = (pdg_a == pdg_b) ? pp_elastic(s) : np_elastic(s);
414  } else {
415  // AQM - Additive Quark Model
416  const double m1 = incoming_particles_[0].effective_mass();
417  const double m2 = incoming_particles_[1].effective_mass();
418  sig_el = pp_elastic_high_energy(s, m1, m2);
419  }
420  if (sig_el > 0.) {
421  return sig_el;
422  } else {
423  std::stringstream ss;
424  const auto name_a = incoming_particles_[0].type().name();
425  const auto name_b = incoming_particles_[1].type().name();
426  ss << "problem in CrossSections::elastic: a=" << name_a << " b=" << name_b
427  << " j_a=" << pdg_a.spin() << " j_b=" << pdg_b.spin()
428  << " sigma=" << sig_el << " s=" << s;
429  throw std::runtime_error(ss.str());
430  }
431 }
432 
433 double CrossSections::npi_el() const {
434  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
435  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
436 
437  const PdgCode& nucleon = pdg_a.is_nucleon() ? pdg_a : pdg_b;
438  const PdgCode& pion = pdg_a.is_nucleon() ? pdg_b : pdg_a;
439  assert(pion != nucleon);
440 
441  const double s = sqrt_s_ * sqrt_s_;
442 
443  double sig_el = 0.;
444  switch (nucleon.code()) {
445  case pdg::p:
446  switch (pion.code()) {
447  case pdg::pi_p:
448  sig_el = piplusp_elastic(s);
449  break;
450  case pdg::pi_m:
451  sig_el = piminusp_elastic(s);
452  break;
453  case pdg::pi_z:
454  sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
455  break;
456  }
457  break;
458  case pdg::n:
459  switch (pion.code()) {
460  case pdg::pi_p:
461  sig_el = piminusp_elastic(s);
462  break;
463  case pdg::pi_m:
464  sig_el = piplusp_elastic(s);
465  break;
466  case pdg::pi_z:
467  sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
468  break;
469  }
470  break;
471  case -pdg::p:
472  switch (pion.code()) {
473  case pdg::pi_p:
474  sig_el = piminusp_elastic(s);
475  break;
476  case pdg::pi_m:
477  sig_el = piplusp_elastic(s);
478  break;
479  case pdg::pi_z:
480  sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
481  break;
482  }
483  break;
484  case -pdg::n:
485  switch (pion.code()) {
486  case pdg::pi_p:
487  sig_el = piplusp_elastic(s);
488  break;
489  case pdg::pi_m:
490  sig_el = piminusp_elastic(s);
491  break;
492  case pdg::pi_z:
493  sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
494  break;
495  }
496  break;
497  default:
498  throw std::runtime_error(
499  "only the elastic cross section for proton-pion "
500  "is implemented");
501  }
502 
503  if (sig_el > 0) {
504  return sig_el;
505  } else {
506  std::stringstream ss;
507  const auto name_a = incoming_particles_[0].type().name();
508  const auto name_b = incoming_particles_[1].type().name();
509  ss << "problem in CrossSections::elastic: a=" << name_a << " b=" << name_b
510  << " j_a=" << pdg_a.spin() << " j_b=" << pdg_b.spin()
511  << " sigma=" << sig_el << " s=" << s;
512  throw std::runtime_error(ss.str());
513  }
514 }
515 
516 CollisionBranchList CrossSections::npi_yk() const {
517  const ParticleType& a = incoming_particles_[0].type();
518  const ParticleType& b = incoming_particles_[1].type();
519  const ParticleType& type_nucleon = a.pdgcode().is_nucleon() ? a : b;
520  const ParticleType& type_pion = a.pdgcode().is_nucleon() ? b : a;
521 
522  const auto pdg_nucleon = type_nucleon.pdgcode().code();
523  const auto pdg_pion = type_pion.pdgcode().code();
524 
525  const double s = sqrt_s_ * sqrt_s_;
526 
527  /* The cross sections are paramectrized for four isospin channels. The
528  * cross sections of the rest isospin channels are obtained using
529  * Clebsch-Gordan coefficients */
530 
531  CollisionBranchList process_list;
532  switch (pdg_nucleon) {
533  case pdg::p: {
534  switch (pdg_pion) {
535  case pdg::pi_p: {
536  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
537  const auto& type_K_p = ParticleType::find(pdg::K_p);
538  add_channel(
539  process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
540  sqrt_s_, type_K_p, type_Sigma_p);
541  break;
542  }
543  case pdg::pi_m: {
544  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
545  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
546  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
547  const auto& type_K_p = ParticleType::find(pdg::K_p);
548  const auto& type_K_z = ParticleType::find(pdg::K_z);
549  add_channel(
550  process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
551  sqrt_s_, type_K_p, type_Sigma_m);
552  add_channel(
553  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
554  type_K_z, type_Sigma_z);
555  add_channel(
556  process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
557  type_K_z, type_Lambda);
558  break;
559  }
560  case pdg::pi_z: {
561  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
562  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
563  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
564  const auto& type_K_p = ParticleType::find(pdg::K_p);
565  const auto& type_K_z = ParticleType::find(pdg::K_z);
566  add_channel(
567  process_list,
568  [&] {
569  return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
570  piminusp_sigma0k0_res(s) +
571  piminusp_sigmaminuskplus_pdg(s));
572  },
573  sqrt_s_, type_K_p, type_Sigma_z);
574  add_channel(
575  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
576  type_K_z, type_Sigma_p);
577  add_channel(
578  process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
579  sqrt_s_, type_K_p, type_Lambda);
580  break;
581  }
582  }
583  break;
584  }
585  case pdg::n: {
586  switch (pdg_pion) {
587  case pdg::pi_p: {
588  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
589  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
590  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
591  const auto& type_K_p = ParticleType::find(pdg::K_p);
592  const auto& type_K_z = ParticleType::find(pdg::K_z);
593  add_channel(
594  process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
595  sqrt_s_, type_K_z, type_Sigma_p);
596  add_channel(
597  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
598  type_K_p, type_Sigma_z);
599  add_channel(
600  process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
601  type_K_p, type_Lambda);
602  break;
603  }
604  case pdg::pi_m: {
605  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
606  const auto& type_K_z = ParticleType::find(pdg::K_z);
607  add_channel(
608  process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
609  sqrt_s_, type_K_z, type_Sigma_m);
610  break;
611  }
612  case pdg::pi_z: {
613  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
614  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
615  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
616  const auto& type_K_p = ParticleType::find(pdg::K_p);
617  const auto& type_K_z = ParticleType::find(pdg::K_z);
618  add_channel(
619  process_list,
620  [&] {
621  return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
622  piminusp_sigma0k0_res(s) +
623  piminusp_sigmaminuskplus_pdg(s));
624  },
625  sqrt_s_, type_K_z, type_Sigma_z);
626  add_channel(
627  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
628  type_K_p, type_Sigma_m);
629  add_channel(
630  process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
631  sqrt_s_, type_K_z, type_Lambda);
632  break;
633  }
634  }
635  break;
636  }
637  case -pdg::p: {
638  switch (pdg_pion) {
639  case pdg::pi_p: {
640  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
641  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
642  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
643  const auto& type_K_m = ParticleType::find(-pdg::K_p);
644  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
645  add_channel(
646  process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
647  sqrt_s_, type_K_m, type_Sigma_m_bar);
648  add_channel(
649  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
650  type_Kbar_z, type_Sigma_z_bar);
651  add_channel(
652  process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
653  type_Kbar_z, type_Lambda_bar);
654  break;
655  }
656  case pdg::pi_m: {
657  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
658  const auto& type_K_m = ParticleType::find(-pdg::K_p);
659  add_channel(
660  process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
661  sqrt_s_, type_K_m, type_Sigma_p_bar);
662  break;
663  }
664  case pdg::pi_z: {
665  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
666  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
667  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
668  const auto& type_K_m = ParticleType::find(-pdg::K_p);
669  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
670  add_channel(
671  process_list,
672  [&] {
673  return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
674  piminusp_sigma0k0_res(s) +
675  piminusp_sigmaminuskplus_pdg(s));
676  },
677  sqrt_s_, type_K_m, type_Sigma_z_bar);
678  add_channel(
679  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
680  type_Kbar_z, type_Sigma_p_bar);
681  add_channel(
682  process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
683  sqrt_s_, type_K_m, type_Lambda_bar);
684  break;
685  }
686  }
687  break;
688  }
689  case -pdg::n: {
690  switch (pdg_pion) {
691  case pdg::pi_p: {
692  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
693  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
694  add_channel(
695  process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
696  sqrt_s_, type_Kbar_z, type_Sigma_m_bar);
697  break;
698  }
699  case pdg::pi_m: {
700  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
701  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
702  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
703  const auto& type_K_m = ParticleType::find(-pdg::K_p);
704  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
705  add_channel(
706  process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
707  sqrt_s_, type_Kbar_z, type_Sigma_p_bar);
708  add_channel(
709  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
710  type_K_m, type_Sigma_z_bar);
711  add_channel(
712  process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
713  type_K_m, type_Lambda_bar);
714  break;
715  }
716  case pdg::pi_z: {
717  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
718  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
719  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
720  const auto& type_K_m = ParticleType::find(-pdg::K_p);
721  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
722  add_channel(
723  process_list,
724  [&] {
725  return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
726  piminusp_sigma0k0_res(s) +
727  piminusp_sigmaminuskplus_pdg(s));
728  },
729  sqrt_s_, type_Kbar_z, type_Sigma_z_bar);
730  add_channel(
731  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
732  type_K_m, type_Sigma_m_bar);
733  add_channel(
734  process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
735  sqrt_s_, type_Kbar_z, type_Lambda_bar);
736  break;
737  }
738  }
739  break;
740  }
741  }
742 
743  return process_list;
744 }
745 
746 double CrossSections::nk_el() const {
747  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
748  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
749 
750  const PdgCode& nucleon = pdg_a.is_nucleon() ? pdg_a : pdg_b;
751  const PdgCode& kaon = pdg_a.is_nucleon() ? pdg_b : pdg_a;
752  assert(kaon != nucleon);
753 
754  const double s = sqrt_s_ * sqrt_s_;
755 
756  double sig_el = 0.;
757  switch (nucleon.code()) {
758  case pdg::p:
759  switch (kaon.code()) {
760  case pdg::K_p:
761  sig_el = kplusp_elastic_background(s);
762  break;
763  case pdg::K_m:
764  sig_el = kminusp_elastic_background(s);
765  break;
766  case pdg::K_z:
767  sig_el = k0p_elastic_background(s);
768  break;
769  case pdg::Kbar_z:
770  sig_el = kbar0p_elastic_background(s);
771  break;
772  }
773  break;
774  case pdg::n:
775  switch (kaon.code()) {
776  case pdg::K_p:
777  sig_el = kplusn_elastic_background(s);
778  break;
779  case pdg::K_m:
780  sig_el = kminusn_elastic_background(s);
781  break;
782  case pdg::K_z:
783  sig_el = k0n_elastic_background(s);
784  break;
785  case pdg::Kbar_z:
786  sig_el = kbar0n_elastic_background(s);
787  break;
788  }
789  break;
790  case -pdg::p:
791  switch (kaon.code()) {
792  case pdg::K_p:
793  sig_el = kminusp_elastic_background(s);
794  break;
795  case pdg::K_m:
796  sig_el = kplusp_elastic_background(s);
797  break;
798  case pdg::K_z:
799  sig_el = kbar0p_elastic_background(s);
800  break;
801  case pdg::Kbar_z:
802  sig_el = k0p_elastic_background(s);
803  break;
804  }
805  break;
806  case -pdg::n:
807  switch (kaon.code()) {
808  case pdg::K_p:
809  sig_el = kminusn_elastic_background(s);
810  break;
811  case pdg::K_m:
812  sig_el = kplusn_elastic_background(s);
813  break;
814  case pdg::K_z:
815  sig_el = kbar0n_elastic_background(s);
816  break;
817  case pdg::Kbar_z:
818  sig_el = k0n_elastic_background(s);
819  break;
820  }
821  break;
822  default:
823  throw std::runtime_error(
824  "elastic cross section for antinucleon-kaon "
825  "not implemented");
826  }
827 
828  if (sig_el > 0) {
829  return sig_el;
830  } else {
831  std::stringstream ss;
832  const auto name_a = incoming_particles_[0].type().name();
833  const auto name_b = incoming_particles_[1].type().name();
834  ss << "problem in CrossSections::elastic: a=" << name_a << " b=" << name_b
835  << " j_a=" << pdg_a.spin() << " j_b=" << pdg_b.spin()
836  << " sigma=" << sig_el << " s=" << s;
837  throw std::runtime_error(ss.str());
838  }
839 }
840 
841 CollisionBranchList CrossSections::two_to_one() const {
842  CollisionBranchList resonance_process_list;
843  const ParticleType& type_particle_a = incoming_particles_[0].type();
844  const ParticleType& type_particle_b = incoming_particles_[1].type();
845 
846  const double m1 = incoming_particles_[0].effective_mass();
847  const double m2 = incoming_particles_[1].effective_mass();
848  const double p_cm_sqr = pCM_sqr(sqrt_s_, m1, m2);
849 
850  ParticleTypePtrList possible_resonances =
851  list_possible_resonances(&type_particle_a, &type_particle_b);
852 
853  // Find all the possible resonances
854  for (const ParticleTypePtr type_resonance : possible_resonances) {
855  double resonance_xsection = formation(*type_resonance, p_cm_sqr);
856 
857  // If cross section is non-negligible, add resonance to the list
858  if (resonance_xsection > really_small) {
859  resonance_process_list.push_back(std::make_unique<CollisionBranch>(
860  *type_resonance, resonance_xsection, ProcessType::TwoToOne));
861  logg[LCrossSections].debug("Found resonance: ", *type_resonance);
862  logg[LCrossSections].debug(type_particle_a.name(), type_particle_b.name(),
863  "->", type_resonance->name(),
864  " at sqrt(s)[GeV] = ", sqrt_s_,
865  " with xs[mb] = ", resonance_xsection);
866  }
867  }
868  return resonance_process_list;
869 }
870 
871 double CrossSections::formation(const ParticleType& type_resonance,
872  double cm_momentum_sqr) const {
873  const ParticleType& type_particle_a = incoming_particles_[0].type();
874  const ParticleType& type_particle_b = incoming_particles_[1].type();
875 
876  // Calculate partial in-width.
877  const double partial_width = type_resonance.get_partial_in_width(
878  sqrt_s_, incoming_particles_[0], incoming_particles_[1]);
879  if (partial_width <= 0.) {
880  return 0.;
881  }
882 
883  assert(type_resonance.charge() ==
884  type_particle_a.charge() + type_particle_b.charge());
885  assert(type_resonance.baryon_number() ==
886  type_particle_a.baryon_number() + type_particle_b.baryon_number());
887 
888  // Calculate spin factor
889  const double spinfactor =
890  static_cast<double>(type_resonance.spin() + 1) /
891  ((type_particle_a.spin() + 1) * (type_particle_b.spin() + 1));
892  const int sym_factor =
893  (type_particle_a.pdgcode() == type_particle_b.pdgcode()) ? 2 : 1;
897  return spinfactor * sym_factor * 2. * M_PI * M_PI / cm_momentum_sqr *
898  type_resonance.spectral_function(sqrt_s_) * partial_width * hbarc *
899  hbarc / fm2_mb;
900 }
901 
902 CollisionBranchList CrossSections::two_to_two(
903  const ReactionsBitSet& included_2to2, const double KN_offset) const {
904  CollisionBranchList process_list;
905  const ParticleData& data_a = incoming_particles_[0];
906  const ParticleData& data_b = incoming_particles_[1];
907  const ParticleType& type_a = data_a.type();
908  const ParticleType& type_b = data_b.type();
909  const auto& pdg_a = data_a.pdgcode();
910  const auto& pdg_b = data_b.pdgcode();
911 
912  if (data_a.is_baryon() && data_b.is_baryon()) {
913  if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
914  pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
915  // Nucleon Nucleon Scattering
916  process_list = nn_xx(included_2to2);
917  } else {
918  // Baryon Baryon Scattering
919  process_list = bb_xx_except_nn(included_2to2);
920  }
921  } else if ((type_a.is_baryon() && type_b.is_meson()) ||
922  (type_a.is_meson() && type_b.is_baryon())) {
923  // Baryon Meson Scattering
924  if ((pdg_a.is_nucleon() && pdg_b.is_kaon()) ||
925  (pdg_b.is_nucleon() && pdg_a.is_kaon())) {
926  // Nucleon Kaon Scattering
927  process_list = nk_xx(included_2to2, KN_offset);
928  } else if ((pdg_a.is_hyperon() && pdg_b.is_pion()) ||
929  (pdg_b.is_hyperon() && pdg_a.is_pion())) {
930  // Hyperon Pion Scattering
931  process_list = ypi_xx(included_2to2);
932  } else if ((pdg_a.is_Delta() && pdg_b.is_kaon()) ||
933  (pdg_b.is_Delta() && pdg_a.is_kaon())) {
934  // Delta Kaon Scattering
935  process_list = deltak_xx(included_2to2);
936  }
937  } else if (type_a.is_nucleus() || type_b.is_nucleus()) {
938  if ((type_a.is_nucleon() && type_b.is_nucleus()) ||
939  (type_b.is_nucleon() && type_a.is_nucleus())) {
940  // Nucleon Deuteron and Nucleon d' Scattering
941  process_list = dn_xx(included_2to2);
942  } else if (((type_a.is_deuteron() || type_a.is_dprime()) &&
943  pdg_b.is_pion()) ||
944  ((type_b.is_deuteron() || type_b.is_dprime()) &&
945  pdg_a.is_pion())) {
946  // Pion Deuteron and Pion d' Scattering
947  process_list = dpi_xx(included_2to2);
948  }
949  }
950  return process_list;
951 }
952 
953 CollisionBranchList CrossSections::two_to_three() const {
954  CollisionBranchList process_list;
955  const ParticleType& type_a = incoming_particles_[0].type();
956  const ParticleType& type_b = incoming_particles_[1].type();
957 
958  if ((type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
959  (type_b.is_deuteron() && type_a.pdgcode().is_pion())) {
960  const ParticleType& type_pi = type_a.pdgcode().is_pion() ? type_a : type_b;
961  const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
962 
963  if (type_nucleus.baryon_number() > 0) {
964  // πd → πpn
965  const auto& type_p = ParticleType::find(pdg::p);
966  const auto& type_n = ParticleType::find(pdg::n);
967 
968  process_list.push_back(std::make_unique<CollisionBranch>(
969  type_pi, type_p, type_n, two_to_three_xs(type_a, type_b, sqrt_s_),
970  ProcessType::TwoToThree));
971  } else {
972  // πd̅ → πp̅n̅
973  const auto& type_anti_p = ParticleType::find(-pdg::p);
974  const auto& type_anti_n = ParticleType::find(-pdg::n);
975 
976  process_list.push_back(std::make_unique<CollisionBranch>(
977  type_pi, type_anti_p, type_anti_n,
978  two_to_three_xs(type_a, type_b, sqrt_s_), ProcessType::TwoToThree));
979  }
980  }
981 
982  if ((type_a.is_nucleon() && type_b.is_deuteron()) ||
983  (type_b.is_nucleon() && type_a.is_deuteron())) {
984  const ParticleType& type_N = type_a.is_nucleon() ? type_a : type_b;
985  const ParticleType& type_nucleus = type_a.is_deuteron() ? type_a : type_b;
986 
987  if (type_nucleus.baryon_number() > 0) {
988  // Nd → Nnp, N̅d → N̅np
989  const auto& type_p = ParticleType::find(pdg::p);
990  const auto& type_n = ParticleType::find(pdg::n);
991 
992  process_list.push_back(std::make_unique<CollisionBranch>(
993  type_N, type_p, type_n, two_to_three_xs(type_a, type_b, sqrt_s_),
994  ProcessType::TwoToThree));
995  } else {
996  // Nd̅ → Np̅n̅, N̅d̅ → N̅p̅n̅
997  const auto& type_anti_p = ParticleType::find(-pdg::p);
998  const auto& type_anti_n = ParticleType::find(-pdg::n);
999 
1000  process_list.push_back(std::make_unique<CollisionBranch>(
1001  type_N, type_anti_p, type_anti_n,
1002  two_to_three_xs(type_a, type_b, sqrt_s_), ProcessType::TwoToThree));
1003  }
1004  }
1005  return process_list;
1006 }
1007 
1008 CollisionBranchList CrossSections::two_to_four() const {
1009  CollisionBranchList process_list;
1010  ParticleTypePtr type_nucleus = &(incoming_particles_[0].type());
1011  ParticleTypePtr type_catalyzer = &(incoming_particles_[1].type());
1012  if (!type_nucleus->is_nucleus()) {
1013  type_nucleus = &(incoming_particles_[1].type());
1014  type_catalyzer = &(incoming_particles_[0].type());
1015  }
1016 
1017  if (type_nucleus->is_nucleus() &&
1018  std::abs(type_nucleus->baryon_number()) == 3 &&
1019  (type_catalyzer->is_pion() || type_catalyzer->is_nucleon())) {
1020  const ParticleTypePtr type_p = ParticleType::try_find(pdg::p);
1021  const ParticleTypePtr type_n = ParticleType::try_find(pdg::n);
1022  const ParticleTypePtr type_anti_p = ParticleType::try_find(-pdg::p);
1023  const ParticleTypePtr type_anti_n = ParticleType::try_find(-pdg::n);
1024  const ParticleTypePtr type_la = ParticleType::try_find(pdg::Lambda);
1025  const ParticleTypePtr type_anti_la = ParticleType::try_find(-pdg::Lambda);
1026 
1027  // Find nucleus components
1028  ParticleTypePtrList components;
1029  components.reserve(3);
1030  const PdgCode nucleus_pdg = type_nucleus->pdgcode();
1031  for (int i = 0; i < nucleus_pdg.nucleus_p(); i++) {
1032  components.push_back(type_p);
1033  }
1034  for (int i = 0; i < nucleus_pdg.nucleus_n(); i++) {
1035  components.push_back(type_n);
1036  }
1037  for (int i = 0; i < nucleus_pdg.nucleus_ap(); i++) {
1038  components.push_back(type_anti_p);
1039  }
1040  for (int i = 0; i < nucleus_pdg.nucleus_an(); i++) {
1041  components.push_back(type_anti_n);
1042  }
1043  for (int i = 0; i < nucleus_pdg.nucleus_La(); i++) {
1044  components.push_back(type_la);
1045  }
1046  for (int i = 0; i < nucleus_pdg.nucleus_aLa(); i++) {
1047  components.push_back(type_anti_la);
1048  }
1049  if (sqrt_s_ > type_catalyzer->mass() + components[0]->mass() +
1050  components[1]->mass() + components[2]->mass()) {
1051  process_list.push_back(std::make_unique<CollisionBranch>(
1052  *type_catalyzer, *(components[0]), *(components[1]), *(components[2]),
1053  two_to_four_xs(*type_nucleus, *type_catalyzer, sqrt_s_),
1054  ProcessType::TwoToFour));
1055  }
1056  }
1057  return process_list;
1058 }
1059 
1060 double CrossSections::d_pi_inelastic_xs(double pion_kinetic_energy) {
1061  const double x = pion_kinetic_energy;
1062  return x * (4.3 + 10.0 * x) / ((x - 0.16) * (x - 0.16) + 0.007);
1063 }
1064 
1065 double CrossSections::d_N_inelastic_xs(double N_kinetic_energy) {
1066  const double x = N_kinetic_energy;
1067  return x * (1.0 + 50 * x) / (x * x + 0.01) +
1068  4 * x / ((x - 0.008) * (x - 0.008) + 0.0004);
1069 }
1070 
1071 double CrossSections::d_aN_inelastic_xs(double aN_kinetic_energy) {
1072  return 55.0 / (aN_kinetic_energy + 0.17);
1073 }
1074 
1075 double CrossSections::two_to_three_xs(const ParticleType& type_a,
1076  const ParticleType& type_b,
1077  double sqrts) {
1078  double xs = 0.0;
1079  ParticleTypePtr type_nucleus = &type_a, type_catalyzer = &type_b;
1080  if (!type_nucleus->is_nucleus()) {
1081  type_nucleus = &type_b;
1082  type_catalyzer = &type_a;
1083  }
1084 
1085  bool nonzero_xs = type_nucleus->is_nucleus() &&
1086  (type_catalyzer->is_pion() || type_catalyzer->is_nucleon());
1087  if (!nonzero_xs) {
1088  return 0.0;
1089  }
1090 
1091  const double md = type_nucleus->mass(), mcat = type_catalyzer->mass();
1092  const double Tkin = (sqrts * sqrts - (md + mcat) * (md + mcat)) / (2.0 * md);
1093 
1094  // Should normally never happen, but may be a useful safeguard
1095  if (Tkin <= 0.0) {
1096  return 0.0;
1097  }
1098 
1099  if (type_catalyzer->is_pion()) {
1100  xs = d_pi_inelastic_xs(Tkin);
1101  } else if (type_catalyzer->is_nucleon()) {
1102  if (type_nucleus->pdgcode().antiparticle_sign() ==
1103  type_catalyzer->pdgcode().antiparticle_sign()) {
1104  // Nd and N̅d̅
1105  xs = d_N_inelastic_xs(Tkin);
1106  } else {
1107  // N̅d and Nd̅
1108  xs = d_aN_inelastic_xs(Tkin);
1109  }
1110  }
1111  return xs;
1112 }
1113 
1114 double CrossSections::two_to_four_xs(const ParticleType& type_a,
1115  const ParticleType& type_b, double sqrts) {
1116  double xs = 0.0;
1117  ParticleTypePtr type_nucleus = &type_a, type_catalyzer = &type_b;
1118  if (!type_nucleus->is_nucleus()) {
1119  type_nucleus = &type_b;
1120  type_catalyzer = &type_a;
1121  }
1122  bool nonzero_xs = type_nucleus->is_nucleus() &&
1123  (type_catalyzer->is_pion() || type_catalyzer->is_nucleon());
1124  if (!nonzero_xs) {
1125  return 0.0;
1126  }
1127 
1128  const double mA = type_nucleus->mass(), mcat = type_catalyzer->mass();
1129  const double Tkin = (sqrts * sqrts - (mA + mcat) * (mA + mcat)) / (2.0 * mA);
1130  const int A = type_nucleus->pdgcode().nucleus_A();
1131  // Should normally never happen, but may be a useful safeguard
1132  if (A != 3 || Tkin <= 0.0) {
1133  return 0.0;
1134  }
1135 
1136  if (type_catalyzer->is_pion()) {
1137  xs = A / 2. * d_pi_inelastic_xs(Tkin);
1138  } else if (type_catalyzer->is_nucleon()) {
1139  if (type_nucleus->pdgcode().antiparticle_sign() ==
1140  type_catalyzer->pdgcode().antiparticle_sign()) {
1141  // N + A, anti-N + anti-A
1142  xs = A / 2. * d_N_inelastic_xs(Tkin);
1143  } else {
1144  // N̅ + A and N + anti-A
1145  xs = A / 2. * d_aN_inelastic_xs(Tkin);
1146  }
1147  }
1148  return xs;
1149 }
1150 
1151 CollisionBranchList CrossSections::bb_xx_except_nn(
1152  const ReactionsBitSet& included_2to2) const {
1153  CollisionBranchList process_list;
1154  const ParticleType& type_a = incoming_particles_[0].type();
1155  const ParticleType& type_b = incoming_particles_[1].type();
1156 
1157  bool same_sign = type_a.antiparticle_sign() == type_b.antiparticle_sign();
1158  bool any_nucleus = type_a.is_nucleus() || type_b.is_nucleus();
1159  if (!same_sign && !any_nucleus) {
1160  return process_list;
1161  }
1162  bool anti_particles = type_a.antiparticle_sign() == -1;
1163  if (type_a.is_nucleon() || type_b.is_nucleon()) {
1164  // N R → N N, N̅ R → N̅ N̅
1165  if (included_2to2[IncludedReactions::NN_to_NR] == 1) {
1166  process_list = bar_bar_to_nuc_nuc(anti_particles);
1167  }
1168  } else if (type_a.is_Delta() || type_b.is_Delta()) {
1169  // Δ R → N N, Δ̅ R → N̅ N̅
1170  if (included_2to2[IncludedReactions::NN_to_DR] == 1) {
1171  process_list = bar_bar_to_nuc_nuc(anti_particles);
1172  }
1173  }
1174 
1175  return process_list;
1176 }
1177 
1178 CollisionBranchList CrossSections::nn_xx(
1179  const ReactionsBitSet& included_2to2) const {
1180  CollisionBranchList process_list, channel_list;
1181 
1182  const double sqrts = sqrt_s_;
1183 
1184  /* Find whether colliding particles are nucleons or anti-nucleons;
1185  * adjust lists of produced particles. */
1186  bool both_antinucleons =
1187  (incoming_particles_[0].type().antiparticle_sign() == -1) &&
1188  (incoming_particles_[1].type().antiparticle_sign() == -1);
1189  const ParticleTypePtrList& nuc_or_anti_nuc =
1190  both_antinucleons ? ParticleType::list_anti_nucleons()
1191  : ParticleType::list_nucleons();
1192  const ParticleTypePtrList& delta_or_anti_delta =
1193  both_antinucleons ? ParticleType::list_anti_Deltas()
1194  : ParticleType::list_Deltas();
1195  // Find N N → N R channels.
1196  if (included_2to2[IncludedReactions::NN_to_NR] == 1) {
1197  channel_list = find_nn_xsection_from_type(
1198  ParticleType::list_baryon_resonances(), nuc_or_anti_nuc,
1199  [&sqrts](const ParticleType& type_res_1, const ParticleType&) {
1200  return type_res_1.iso_multiplet()->get_integral_NR(sqrts);
1201  });
1202  process_list.reserve(process_list.size() + channel_list.size());
1203  std::move(channel_list.begin(), channel_list.end(),
1204  std::inserter(process_list, process_list.end()));
1205  channel_list.clear();
1206  }
1207 
1208  // Find N N → Δ R channels.
1209  if (included_2to2[IncludedReactions::NN_to_DR] == 1) {
1210  channel_list = find_nn_xsection_from_type(
1211  ParticleType::list_baryon_resonances(), delta_or_anti_delta,
1212  [&sqrts](const ParticleType& type_res_1,
1213  const ParticleType& type_res_2) {
1214  return type_res_1.iso_multiplet()->get_integral_RR(
1215  type_res_2.iso_multiplet(), sqrts);
1216  });
1217  process_list.reserve(process_list.size() + channel_list.size());
1218  std::move(channel_list.begin(), channel_list.end(),
1219  std::inserter(process_list, process_list.end()));
1220  channel_list.clear();
1221  }
1222 
1223  // Find N N → dπ and N̅ N̅→ d̅π channels.
1224  ParticleTypePtr deuteron = ParticleType::try_find(pdg::deuteron);
1225  ParticleTypePtr antideutron = ParticleType::try_find(pdg::antideuteron);
1226  ParticleTypePtr pim = ParticleType::try_find(pdg::pi_m);
1227  ParticleTypePtr pi0 = ParticleType::try_find(pdg::pi_z);
1228  ParticleTypePtr pip = ParticleType::try_find(pdg::pi_p);
1229  // Make sure all the necessary particle types are found
1230  if (deuteron && antideutron && pim && pi0 && pip &&
1231  included_2to2[IncludedReactions::PiDeuteron_to_NN] == 1) {
1232  const ParticleTypePtrList deutron_list = {deuteron};
1233  const ParticleTypePtrList antideutron_list = {antideutron};
1234  const ParticleTypePtrList pion_list = {pim, pi0, pip};
1235  channel_list = find_nn_xsection_from_type(
1236  (both_antinucleons ? antideutron_list : deutron_list), pion_list,
1237  [&sqrts](const ParticleType& type_res_1,
1238  const ParticleType& type_res_2) {
1239  return pCM(sqrts, type_res_1.mass(), type_res_2.mass());
1240  });
1241  process_list.reserve(process_list.size() + channel_list.size());
1242  std::move(channel_list.begin(), channel_list.end(),
1243  std::inserter(process_list, process_list.end()));
1244  channel_list.clear();
1245  }
1246 
1247  return process_list;
1248 }
1249 
1250 CollisionBranchList CrossSections::nk_xx(const ReactionsBitSet& included_2to2,
1251  const double KN_offset) const {
1252  const ParticleType& a = incoming_particles_[0].type();
1253  const ParticleType& b = incoming_particles_[1].type();
1254  const ParticleType& type_nucleon = a.pdgcode().is_nucleon() ? a : b;
1255  const ParticleType& type_kaon = a.pdgcode().is_nucleon() ? b : a;
1256 
1257  const auto pdg_nucleon = type_nucleon.pdgcode().code();
1258  const auto pdg_kaon = type_kaon.pdgcode().code();
1259 
1260  const double s = sqrt_s_ * sqrt_s_;
1261 
1262  // Some variable declarations for frequently used quantities
1263  const auto sigma_kplusp = kplusp_inelastic_background(s);
1264  const auto sigma_kplusn = kplusn_inelastic_background(s);
1265 
1266  /* At high energy, the parametrization we use diverges from experimental
1267  * data. This cutoff represents the point where the AQM cross section
1268  * becomes smaller than this parametrization, so we cut it here, and fully
1269  * switch to AQM beyond this point. */
1270  const double KN_to_KDelta_cutoff = KN_offset +
1271  incoming_particles_[0].pole_mass() +
1272  incoming_particles_[1].pole_mass();
1273 
1274  bool incl_KN_to_KN = included_2to2[IncludedReactions::KN_to_KN] == 1;
1275  bool incl_KN_to_KDelta =
1276  included_2to2[IncludedReactions::KN_to_KDelta] == 1 &&
1277  sqrt_s_ < KN_to_KDelta_cutoff;
1278  bool incl_Strangeness_exchange =
1279  included_2to2[IncludedReactions::Strangeness_exchange] == 1;
1280 
1281  CollisionBranchList process_list;
1282  switch (pdg_kaon) {
1283  case pdg::K_m: {
1284  /* All inelastic K- N channels here are strangeness exchange, plus one
1285  * charge exchange. */
1286  switch (pdg_nucleon) {
1287  case pdg::p: {
1288  if (incl_Strangeness_exchange) {
1289  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1290  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
1291  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
1292  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
1293  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
1294  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
1295  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
1296  add_channel(
1297  process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
1298  sqrt_s_, type_pi_m, type_Sigma_p);
1299  add_channel(
1300  process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
1301  sqrt_s_, type_pi_p, type_Sigma_m);
1302  add_channel(
1303  process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
1304  sqrt_s_, type_pi_z, type_Sigma_z);
1305  add_channel(
1306  process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
1307  sqrt_s_, type_pi_z, type_Lambda);
1308  }
1309  if (incl_KN_to_KN) {
1310  const auto& type_n = ParticleType::find(pdg::n);
1311  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1312  add_channel(
1313  process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
1314  type_Kbar_z, type_n);
1315  }
1316  break;
1317  }
1318  case pdg::n: {
1319  if (incl_Strangeness_exchange) {
1320  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1321  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
1322  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
1323  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
1324  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
1325  add_channel(
1326  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1327  sqrt_s_, type_pi_m, type_Sigma_z);
1328  add_channel(
1329  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1330  sqrt_s_, type_pi_z, type_Sigma_m);
1331  add_channel(
1332  process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
1333  sqrt_s_, type_pi_m, type_Lambda);
1334  }
1335  break;
1336  }
1337  case -pdg::p: {
1338  if (incl_KN_to_KDelta) {
1339  const auto& type_K_m = ParticleType::find(pdg::K_m);
1340  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1341  const auto& type_Delta_pp_bar = ParticleType::find(-pdg::Delta_pp);
1342  const auto& type_Delta_p_bar = ParticleType::find(-pdg::Delta_p);
1343  add_channel(
1344  process_list,
1345  [&] {
1346  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
1347  type_nucleon, type_kaon,
1348  type_Kbar_z, type_Delta_pp_bar);
1349  },
1350  sqrt_s_, type_Kbar_z, type_Delta_pp_bar);
1351  add_channel(
1352  process_list,
1353  [&] {
1354  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
1355  type_nucleon, type_kaon, type_K_m,
1356  type_Delta_p_bar);
1357  },
1358  sqrt_s_, type_K_m, type_Delta_p_bar);
1359  }
1360  break;
1361  }
1362  case -pdg::n: {
1363  if (incl_KN_to_KDelta) {
1364  const auto& type_K_m = ParticleType::find(pdg::K_m);
1365  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1366  const auto& type_Delta_p_bar = ParticleType::find(-pdg::Delta_p);
1367  const auto& type_Delta_z_bar = ParticleType::find(-pdg::Delta_z);
1368  add_channel(
1369  process_list,
1370  [&] {
1371  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
1372  type_nucleon, type_kaon,
1373  type_Kbar_z, type_Delta_p_bar);
1374  },
1375  sqrt_s_, type_Kbar_z, type_Delta_p_bar);
1376  add_channel(
1377  process_list,
1378  [&] {
1379  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
1380  type_nucleon, type_kaon, type_K_m,
1381  type_Delta_z_bar);
1382  },
1383  sqrt_s_, type_K_m, type_Delta_z_bar);
1384  }
1385  if (incl_KN_to_KN) {
1386  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1387  const auto& type_p_bar = ParticleType::find(-pdg::p);
1388  add_channel(
1389  process_list, [&] { return kplusn_k0p(s); }, sqrt_s_,
1390  type_Kbar_z, type_p_bar);
1391  }
1392  break;
1393  }
1394  }
1395  break;
1396  }
1397  case pdg::K_p: {
1398  /* All inelastic channels are K+ N -> K Delta -> K pi N, with identical
1399  * cross section, weighted by the isospin factor. */
1400  switch (pdg_nucleon) {
1401  case pdg::p: {
1402  if (incl_KN_to_KDelta) {
1403  const auto& type_K_p = ParticleType::find(pdg::K_p);
1404  const auto& type_K_z = ParticleType::find(pdg::K_z);
1405  const auto& type_Delta_pp = ParticleType::find(pdg::Delta_pp);
1406  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
1407  add_channel(
1408  process_list,
1409  [&] {
1410  return sigma_kplusp *
1411  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
1412  type_K_z, type_Delta_pp);
1413  },
1414  sqrt_s_, type_K_z, type_Delta_pp);
1415  add_channel(
1416  process_list,
1417  [&] {
1418  return sigma_kplusp *
1419  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
1420  type_K_p, type_Delta_p);
1421  },
1422  sqrt_s_, type_K_p, type_Delta_p);
1423  }
1424  break;
1425  }
1426  case pdg::n: {
1427  if (incl_KN_to_KDelta) {
1428  const auto& type_K_p = ParticleType::find(pdg::K_p);
1429  const auto& type_K_z = ParticleType::find(pdg::K_z);
1430  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
1431  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
1432  add_channel(
1433  process_list,
1434  [&] {
1435  return sigma_kplusn *
1436  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
1437  type_K_z, type_Delta_p);
1438  },
1439  sqrt_s_, type_K_z, type_Delta_p);
1440  add_channel(
1441  process_list,
1442  [&] {
1443  return sigma_kplusn *
1444  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
1445  type_K_p, type_Delta_z);
1446  },
1447  sqrt_s_, type_K_p, type_Delta_z);
1448  }
1449  if (incl_KN_to_KN) {
1450  const auto& type_K_z = ParticleType::find(pdg::K_z);
1451  const auto& type_p = ParticleType::find(pdg::p);
1452  add_channel(
1453  process_list, [&] { return kplusn_k0p(s); }, sqrt_s_, type_K_z,
1454  type_p);
1455  }
1456  break;
1457  }
1458  case -pdg::p: {
1459  if (incl_Strangeness_exchange) {
1460  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1461  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
1462  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
1463  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
1464  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
1465  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
1466  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
1467  add_channel(
1468  process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
1469  sqrt_s_, type_pi_p, type_Sigma_p_bar);
1470  add_channel(
1471  process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
1472  sqrt_s_, type_pi_m, type_Sigma_m_bar);
1473  add_channel(
1474  process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
1475  sqrt_s_, type_pi_z, type_Sigma_z_bar);
1476  add_channel(
1477  process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
1478  sqrt_s_, type_pi_z, type_Lambda_bar);
1479  }
1480  if (incl_KN_to_KN) {
1481  const auto& type_n_bar = ParticleType::find(-pdg::n);
1482  const auto& type_K_z = ParticleType::find(pdg::K_z);
1483  add_channel(
1484  process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
1485  type_K_z, type_n_bar);
1486  }
1487  break;
1488  }
1489  case -pdg::n: {
1490  if (incl_Strangeness_exchange) {
1491  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1492  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
1493  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
1494  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
1495  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
1496  add_channel(
1497  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1498  sqrt_s_, type_pi_p, type_Sigma_z_bar);
1499  add_channel(
1500  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1501  sqrt_s_, type_pi_z, type_Sigma_m_bar);
1502  add_channel(
1503  process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
1504  sqrt_s_, type_pi_p, type_Lambda_bar);
1505  }
1506  break;
1507  }
1508  }
1509  break;
1510  }
1511  case pdg::K_z: {
1512  // K+ and K0 have the same mass and spin, so their cross sections are
1513  // assumed to only differ in isospin factors. For the initial state, we
1514  // assume that K0 p is equivalent to K+ n and K0 n equivalent to K+ p,
1515  // like for the elastic background.
1516 
1517  switch (pdg_nucleon) {
1518  case pdg::p: {
1519  if (incl_KN_to_KDelta) {
1520  const auto& type_K_p = ParticleType::find(pdg::K_p);
1521  const auto& type_K_z = ParticleType::find(pdg::K_z);
1522  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
1523  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
1524  add_channel(
1525  process_list,
1526  [&] {
1527  return sigma_kplusn *
1528  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
1529  type_K_z, type_Delta_p);
1530  },
1531  sqrt_s_, type_K_z, type_Delta_p);
1532  add_channel(
1533  process_list,
1534  [&] {
1535  return sigma_kplusn *
1536  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
1537  type_K_p, type_Delta_z);
1538  },
1539  sqrt_s_, type_K_p, type_Delta_z);
1540  }
1541  if (incl_KN_to_KN) {
1542  const auto& type_K_p = ParticleType::find(pdg::K_p);
1543  const auto& type_n = ParticleType::find(pdg::n);
1544  add_channel(
1545  process_list,
1546  [&] {
1547  // The isospin factor is 1, see the parametrizations
1548  // tests.
1549  return kplusn_k0p(s);
1550  },
1551  sqrt_s_, type_K_p, type_n);
1552  }
1553  break;
1554  }
1555  case pdg::n: {
1556  if (incl_KN_to_KDelta) {
1557  const auto& type_K_p = ParticleType::find(pdg::K_p);
1558  const auto& type_K_z = ParticleType::find(pdg::K_z);
1559  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
1560  const auto& type_Delta_m = ParticleType::find(pdg::Delta_m);
1561  add_channel(
1562  process_list,
1563  [&] {
1564  return sigma_kplusp *
1565  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
1566  type_K_z, type_Delta_z);
1567  },
1568  sqrt_s_, type_K_z, type_Delta_z);
1569  add_channel(
1570  process_list,
1571  [&] {
1572  return sigma_kplusp *
1573  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
1574  type_K_p, type_Delta_m);
1575  },
1576  sqrt_s_, type_K_p, type_Delta_m);
1577  }
1578  break;
1579  }
1580  case -pdg::p: {
1581  if (incl_Strangeness_exchange) {
1582  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1583  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
1584  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
1585  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
1586  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
1587  add_channel(
1588  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1589  sqrt_s_, type_pi_m, type_Sigma_z_bar);
1590  add_channel(
1591  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1592  sqrt_s_, type_pi_z, type_Sigma_p_bar);
1593  add_channel(
1594  process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
1595  sqrt_s_, type_pi_m, type_Lambda_bar);
1596  }
1597  break;
1598  }
1599  case -pdg::n: {
1600  if (incl_Strangeness_exchange) {
1601  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1602  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
1603  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
1604  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
1605  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
1606  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
1607  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
1608  add_channel(
1609  process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
1610  sqrt_s_, type_pi_m, type_Sigma_m_bar);
1611  add_channel(
1612  process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
1613  sqrt_s_, type_pi_p, type_Sigma_p_bar);
1614  add_channel(
1615  process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
1616  sqrt_s_, type_pi_z, type_Sigma_z_bar);
1617  add_channel(
1618  process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
1619  sqrt_s_, type_pi_z, type_Lambda_bar);
1620  }
1621  if (incl_KN_to_KN) {
1622  const auto& type_K_p = ParticleType::find(pdg::K_p);
1623  const auto& type_p_bar = ParticleType::find(-pdg::p);
1624  add_channel(
1625  process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
1626  type_K_p, type_p_bar);
1627  }
1628  break;
1629  }
1630  }
1631  break;
1632  }
1633  case pdg::Kbar_z:
1634  switch (pdg_nucleon) {
1635  case pdg::p: {
1636  if (incl_Strangeness_exchange) {
1637  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1638  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
1639  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
1640  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
1641  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
1642  add_channel(
1643  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1644  sqrt_s_, type_pi_z, type_Sigma_p);
1645  add_channel(
1646  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1647  sqrt_s_, type_pi_p, type_Sigma_z);
1648  add_channel(
1649  process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
1650  sqrt_s_, type_pi_p, type_Lambda);
1651  }
1652  break;
1653  }
1654  case pdg::n: {
1655  if (incl_Strangeness_exchange) {
1656  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1657  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
1658  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
1659  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
1660  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
1661  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
1662  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
1663  add_channel(
1664  process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
1665  sqrt_s_, type_pi_p, type_Sigma_m);
1666  add_channel(
1667  process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
1668  sqrt_s_, type_pi_m, type_Sigma_p);
1669  add_channel(
1670  process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
1671  sqrt_s_, type_pi_z, type_Sigma_z);
1672  add_channel(
1673  process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
1674  sqrt_s_, type_pi_z, type_Lambda);
1675  }
1676  if (incl_KN_to_KN) {
1677  const auto& type_p = ParticleType::find(pdg::p);
1678  const auto& type_K_m = ParticleType::find(pdg::K_m);
1679  add_channel(
1680  process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
1681  type_K_m, type_p);
1682  }
1683  break;
1684  }
1685  case -pdg::p: {
1686  if (incl_KN_to_KDelta) {
1687  const auto& type_K_m = ParticleType::find(pdg::K_m);
1688  const auto& type_Kbar_z = type_kaon;
1689  const auto& type_Delta_bar_m = ParticleType::find(-pdg::Delta_p);
1690  const auto& type_Delta_bar_z = ParticleType::find(-pdg::Delta_z);
1691  add_channel(
1692  process_list,
1693  [&] {
1694  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
1695  type_nucleon, type_kaon,
1696  type_Kbar_z, type_Delta_bar_m);
1697  },
1698  sqrt_s_, type_Kbar_z, type_Delta_bar_m);
1699  add_channel(
1700  process_list,
1701  [&] {
1702  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
1703  type_nucleon, type_kaon, type_K_m,
1704  type_Delta_bar_z);
1705  },
1706  sqrt_s_, type_K_m, type_Delta_bar_z);
1707  }
1708  if (incl_KN_to_KN) {
1709  const auto& type_K_m = ParticleType::find(pdg::K_m);
1710  const auto& type_n_bar = ParticleType::find(-pdg::n);
1711  add_channel(
1712  process_list,
1713  [&] {
1714  // The isospin factor is 1, see the parametrizations
1715  // tests.
1716  return kplusn_k0p(s);
1717  },
1718  sqrt_s_, type_K_m, type_n_bar);
1719  }
1720  break;
1721  }
1722  case -pdg::n: {
1723  if (incl_KN_to_KDelta) {
1724  const auto& type_K_m = ParticleType::find(pdg::K_m);
1725  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1726  const auto& type_Delta_z_bar = ParticleType::find(-pdg::Delta_z);
1727  const auto& type_Delta_m_bar = ParticleType::find(-pdg::Delta_m);
1728  add_channel(
1729  process_list,
1730  [&] {
1731  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
1732  type_nucleon, type_kaon,
1733  type_Kbar_z, type_Delta_z_bar);
1734  },
1735  sqrt_s_, type_Kbar_z, type_Delta_z_bar);
1736  add_channel(
1737  process_list,
1738  [&] {
1739  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
1740  type_nucleon, type_kaon, type_K_m,
1741  type_Delta_m_bar);
1742  },
1743  sqrt_s_, type_K_m, type_Delta_m_bar);
1744  }
1745  break;
1746  }
1747  }
1748  break;
1749  }
1750 
1751  return process_list;
1752 }
1753 
1754 CollisionBranchList CrossSections::deltak_xx(
1755  const ReactionsBitSet& included_2to2) const {
1756  CollisionBranchList process_list;
1757  if (included_2to2[IncludedReactions::KN_to_KDelta] == 0) {
1758  return process_list;
1759  }
1760  const ParticleType& a = incoming_particles_[0].type();
1761  const ParticleType& b = incoming_particles_[1].type();
1762  const ParticleType& type_delta = a.pdgcode().is_Delta() ? a : b;
1763  const ParticleType& type_kaon = a.pdgcode().is_Delta() ? b : a;
1764 
1765  const auto pdg_delta = type_delta.pdgcode().code();
1766  const auto pdg_kaon = type_kaon.pdgcode().code();
1767 
1768  const double s = sqrt_s_ * sqrt_s_;
1769  const double pcm = cm_momentum();
1770  /* The cross sections are determined from the backward reactions via
1771  * detailed balance. The same isospin factors as for the backward reaction
1772  * are used. */
1773  switch (pack(pdg_delta, pdg_kaon)) {
1774  case pack(pdg::Delta_pp, pdg::K_z):
1775  case pack(pdg::Delta_p, pdg::K_p): {
1776  const auto& type_p = ParticleType::find(pdg::p);
1777  const auto& type_K_p = ParticleType::find(pdg::K_p);
1778  add_channel(
1779  process_list,
1780  [&] {
1781  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
1782  type_kaon, type_p, type_K_p) *
1783  kaon_nucleon_ratios.get_ratio(type_p, type_K_p, type_kaon,
1784  type_delta) *
1785  kplusp_inelastic_background(s);
1786  },
1787  sqrt_s_, type_p, type_K_p);
1788  break;
1789  }
1790  case pack(-pdg::Delta_pp, pdg::Kbar_z):
1791  case pack(-pdg::Delta_p, pdg::K_m): {
1792  const auto& type_p_bar = ParticleType::find(-pdg::p);
1793  const auto& type_K_m = ParticleType::find(pdg::K_m);
1794  add_channel(
1795  process_list,
1796  [&] {
1797  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
1798  type_kaon, type_p_bar, type_K_m) *
1799  kaon_nucleon_ratios.get_ratio(type_p_bar, type_K_m,
1800  type_kaon, type_delta) *
1801  kplusp_inelastic_background(s);
1802  },
1803  sqrt_s_, type_p_bar, type_K_m);
1804  break;
1805  }
1806  case pack(pdg::Delta_p, pdg::K_z):
1807  case pack(pdg::Delta_z, pdg::K_p): {
1808  const auto& type_n = ParticleType::find(pdg::n);
1809  const auto& type_p = ParticleType::find(pdg::p);
1810  const auto& type_K_p = ParticleType::find(pdg::K_p);
1811  const auto& type_K_z = ParticleType::find(pdg::K_z);
1812  add_channel(
1813  process_list,
1814  [&] {
1815  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
1816  type_kaon, type_n, type_K_p) *
1817  kaon_nucleon_ratios.get_ratio(type_n, type_K_p, type_kaon,
1818  type_delta) *
1819  kplusn_inelastic_background(s);
1820  },
1821  sqrt_s_, type_n, type_K_p);
1822 
1823  add_channel(
1824  process_list,
1825  [&] {
1826  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
1827  type_kaon, type_p, type_K_z) *
1828  kaon_nucleon_ratios.get_ratio(type_p, type_K_z, type_kaon,
1829  type_delta) *
1830  kplusn_inelastic_background(s);
1831  },
1832  sqrt_s_, type_p, type_K_z);
1833  break;
1834  }
1835  case pack(-pdg::Delta_p, pdg::Kbar_z):
1836  case pack(-pdg::Delta_z, pdg::K_m): {
1837  const auto& type_n_bar = ParticleType::find(-pdg::n);
1838  const auto& type_p_bar = ParticleType::find(-pdg::p);
1839  const auto& type_K_m = ParticleType::find(pdg::K_m);
1840  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1841  add_channel(
1842  process_list,
1843  [&] {
1844  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
1845  type_kaon, type_n_bar, type_K_m) *
1846  kaon_nucleon_ratios.get_ratio(type_n_bar, type_K_m,
1847  type_kaon, type_delta) *
1848  kplusn_inelastic_background(s);
1849  },
1850  sqrt_s_, type_n_bar, type_K_m);
1851 
1852  add_channel(
1853  process_list,
1854  [&] {
1855  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
1856  type_kaon, type_p_bar,
1857  type_Kbar_z) *
1858  kaon_nucleon_ratios.get_ratio(type_p_bar, type_Kbar_z,
1859  type_kaon, type_delta) *
1860  kplusn_inelastic_background(s);
1861  },
1862  sqrt_s_, type_p_bar, type_Kbar_z);
1863  break;
1864  }
1865  case pack(pdg::Delta_z, pdg::K_z):
1866  case pack(pdg::Delta_m, pdg::K_p): {
1867  const auto& type_n = ParticleType::find(pdg::n);
1868  const auto& type_K_z = ParticleType::find(pdg::K_z);
1869  add_channel(
1870  process_list,
1871  [&] {
1872  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
1873  type_kaon, type_n, type_K_z) *
1874  kaon_nucleon_ratios.get_ratio(type_n, type_K_z, type_kaon,
1875  type_delta) *
1876  kplusp_inelastic_background(s);
1877  },
1878  sqrt_s_, type_n, type_K_z);
1879  break;
1880  }
1881  case pack(-pdg::Delta_z, pdg::Kbar_z):
1882  case pack(-pdg::Delta_m, pdg::K_m): {
1883  const auto& type_n_bar = ParticleType::find(-pdg::n);
1884  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1885  add_channel(
1886  process_list,
1887  [&] {
1888  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
1889  type_kaon, type_n_bar,
1890  type_Kbar_z) *
1891  kaon_nucleon_ratios.get_ratio(type_n_bar, type_Kbar_z,
1892  type_kaon, type_delta) *
1893  kplusp_inelastic_background(s);
1894  },
1895  sqrt_s_, type_n_bar, type_Kbar_z);
1896  break;
1897  }
1898  default:
1899  break;
1900  }
1901 
1902  return process_list;
1903 }
1904 
1905 CollisionBranchList CrossSections::ypi_xx(
1906  const ReactionsBitSet& included_2to2) const {
1907  CollisionBranchList process_list;
1908  if (included_2to2[IncludedReactions::Strangeness_exchange] == 0) {
1909  return process_list;
1910  }
1911  const ParticleType& a = incoming_particles_[0].type();
1912  const ParticleType& b = incoming_particles_[1].type();
1913  const ParticleType& type_hyperon = a.pdgcode().is_hyperon() ? a : b;
1914  const ParticleType& type_pion = a.pdgcode().is_hyperon() ? b : a;
1915 
1916  const auto pdg_hyperon = type_hyperon.pdgcode().code();
1917  const auto pdg_pion = type_pion.pdgcode().code();
1918 
1919  const double s = sqrt_s_ * sqrt_s_;
1920 
1921  switch (pack(pdg_hyperon, pdg_pion)) {
1922  case pack(pdg::Sigma_z, pdg::pi_m): {
1923  const auto& type_n = ParticleType::find(pdg::n);
1924  const auto& type_K_m = ParticleType::find(pdg::K_m);
1925  add_channel(
1926  process_list,
1927  [&] {
1928  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
1929  type_n, type_K_m) *
1930  kminusn_piminussigma0(sqrt_s_);
1931  },
1932  sqrt_s_, type_n, type_K_m);
1933  break;
1934  }
1935  case pack(pdg::Sigma_z, pdg::pi_p): {
1936  const auto& type_p = ParticleType::find(pdg::p);
1937  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1938  add_channel(
1939  process_list,
1940  [&] {
1941  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
1942  type_p, type_Kbar_z) *
1943  kminusn_piminussigma0(sqrt_s_);
1944  },
1945  sqrt_s_, type_p, type_Kbar_z);
1946  break;
1947  }
1948  case pack(-pdg::Sigma_z, pdg::pi_p): {
1949  const auto& type_n_bar = ParticleType::find(-pdg::n);
1950  const auto& type_K_p = ParticleType::find(pdg::K_p);
1951  add_channel(
1952  process_list,
1953  [&] {
1954  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
1955  type_n_bar, type_K_p) *
1956  kminusn_piminussigma0(sqrt_s_);
1957  },
1958  sqrt_s_, type_n_bar, type_K_p);
1959  break;
1960  }
1961  case pack(-pdg::Sigma_z, pdg::pi_m): {
1962  const auto& type_p_bar = ParticleType::find(-pdg::p);
1963  const auto& type_K_z = ParticleType::find(pdg::K_z);
1964  add_channel(
1965  process_list,
1966  [&] {
1967  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
1968  type_p_bar, type_K_z) *
1969  kminusn_piminussigma0(sqrt_s_);
1970  },
1971  sqrt_s_, type_p_bar, type_K_z);
1972  break;
1973  }
1974  case pack(pdg::Sigma_m, pdg::pi_z): {
1975  const auto& type_n = ParticleType::find(pdg::n);
1976  const auto& type_K_m = ParticleType::find(pdg::K_m);
1977  add_channel(
1978  process_list,
1979  [&] {
1980  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
1981  type_n, type_K_m) *
1982  kminusn_piminussigma0(sqrt_s_);
1983  },
1984  sqrt_s_, type_n, type_K_m);
1985  break;
1986  }
1987  case pack(pdg::Sigma_p, pdg::pi_z): {
1988  const auto& type_p = ParticleType::find(pdg::p);
1989  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1990  add_channel(
1991  process_list,
1992  [&] {
1993  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
1994  type_p, type_Kbar_z) *
1995  kminusn_piminussigma0(sqrt_s_);
1996  },
1997  sqrt_s_, type_p, type_Kbar_z);
1998  break;
1999  }
2000  case pack(-pdg::Sigma_m, pdg::pi_z): {
2001  const auto& type_n_bar = ParticleType::find(-pdg::n);
2002  const auto& type_K_p = ParticleType::find(pdg::K_p);
2003  add_channel(
2004  process_list,
2005  [&] {
2006  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2007  type_n_bar, type_K_p) *
2008  kminusn_piminussigma0(sqrt_s_);
2009  },
2010  sqrt_s_, type_n_bar, type_K_p);
2011  break;
2012  }
2013  case pack(-pdg::Sigma_p, pdg::pi_z): {
2014  const auto& type_p_bar = ParticleType::find(-pdg::p);
2015  const auto& type_K_z = ParticleType::find(pdg::K_z);
2016  add_channel(
2017  process_list,
2018  [&] {
2019  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2020  type_p_bar, type_K_z) *
2021  kminusn_piminussigma0(sqrt_s_);
2022  },
2023  sqrt_s_, type_p_bar, type_K_z);
2024  break;
2025  }
2026  case pack(pdg::Lambda, pdg::pi_m): {
2027  const auto& type_n = ParticleType::find(pdg::n);
2028  const auto& type_K_m = ParticleType::find(pdg::K_m);
2029  add_channel(
2030  process_list,
2031  [&] {
2032  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2033  type_n, type_K_m) *
2034  kminusn_piminuslambda(sqrt_s_);
2035  },
2036  sqrt_s_, type_n, type_K_m);
2037  break;
2038  }
2039  case pack(pdg::Lambda, pdg::pi_p): {
2040  const auto& type_p = ParticleType::find(pdg::p);
2041  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2042  add_channel(
2043  process_list,
2044  [&] {
2045  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2046  type_p, type_Kbar_z) *
2047  kminusn_piminuslambda(sqrt_s_);
2048  },
2049  sqrt_s_, type_p, type_Kbar_z);
2050  break;
2051  }
2052  case pack(-pdg::Lambda, pdg::pi_p): {
2053  const auto& type_n_bar = ParticleType::find(-pdg::n);
2054  const auto& type_K_p = ParticleType::find(pdg::K_p);
2055  add_channel(
2056  process_list,
2057  [&] {
2058  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2059  type_n_bar, type_K_p) *
2060  kminusn_piminuslambda(sqrt_s_);
2061  },
2062  sqrt_s_, type_n_bar, type_K_p);
2063  break;
2064  }
2065  case pack(-pdg::Lambda, pdg::pi_m): {
2066  const auto& type_p_bar = ParticleType::find(-pdg::p);
2067  const auto& type_K_z = ParticleType::find(pdg::K_z);
2068  add_channel(
2069  process_list,
2070  [&] {
2071  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2072  type_p_bar, type_K_z) *
2073  kminusn_piminuslambda(sqrt_s_);
2074  },
2075  sqrt_s_, type_p_bar, type_K_z);
2076  break;
2077  }
2078  case pack(pdg::Sigma_z, pdg::pi_z): {
2079  const auto& type_p = ParticleType::find(pdg::p);
2080  const auto& type_n = ParticleType::find(pdg::n);
2081  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2082  const auto& type_K_m = ParticleType::find(pdg::K_m);
2083  add_channel(
2084  process_list,
2085  [&] {
2086  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2087  type_p, type_K_m) *
2088  kminusp_pi0sigma0(sqrt_s_);
2089  },
2090  sqrt_s_, type_p, type_K_m);
2091  add_channel(
2092  process_list,
2093  [&] {
2094  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2095  type_n, type_Kbar_z) *
2096  kminusp_pi0sigma0(sqrt_s_);
2097  },
2098  sqrt_s_, type_n, type_Kbar_z);
2099  break;
2100  }
2101  case pack(-pdg::Sigma_z, pdg::pi_z): {
2102  const auto& type_p_bar = ParticleType::find(-pdg::p);
2103  const auto& type_n_bar = ParticleType::find(-pdg::n);
2104  const auto& type_K_z = ParticleType::find(pdg::K_z);
2105  const auto& type_K_p = ParticleType::find(pdg::K_p);
2106  add_channel(
2107  process_list,
2108  [&] {
2109  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2110  type_p_bar, type_K_p) *
2111  kminusp_pi0sigma0(sqrt_s_);
2112  },
2113  sqrt_s_, type_p_bar, type_K_p);
2114  add_channel(
2115  process_list,
2116  [&] {
2117  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2118  type_n_bar, type_K_z) *
2119  kminusp_pi0sigma0(sqrt_s_);
2120  },
2121  sqrt_s_, type_n_bar, type_K_z);
2122  break;
2123  }
2124  case pack(pdg::Sigma_m, pdg::pi_p): {
2125  const auto& type_p = ParticleType::find(pdg::p);
2126  const auto& type_n = ParticleType::find(pdg::n);
2127  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2128  const auto& type_K_m = ParticleType::find(pdg::K_m);
2129  add_channel(
2130  process_list,
2131  [&] {
2132  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2133  type_p, type_K_m) *
2134  kminusp_piplussigmaminus(sqrt_s_);
2135  },
2136  sqrt_s_, type_p, type_K_m);
2137  add_channel(
2138  process_list,
2139  [&] {
2140  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2141  type_n, type_Kbar_z) *
2142  kminusp_piminussigmaplus(sqrt_s_);
2143  },
2144  sqrt_s_, type_n, type_Kbar_z);
2145  break;
2146  }
2147  case pack(-pdg::Sigma_m, pdg::pi_m): {
2148  const auto& type_p_bar = ParticleType::find(-pdg::p);
2149  const auto& type_n_bar = ParticleType::find(-pdg::n);
2150  const auto& type_K_z = ParticleType::find(pdg::K_z);
2151  const auto& type_K_p = ParticleType::find(pdg::K_p);
2152  add_channel(
2153  process_list,
2154  [&] {
2155  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2156  type_p_bar, type_K_p) *
2157  kminusp_piplussigmaminus(sqrt_s_);
2158  },
2159  sqrt_s_, type_p_bar, type_K_p);
2160  add_channel(
2161  process_list,
2162  [&] {
2163  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2164  type_n_bar, type_K_z) *
2165  kminusp_piminussigmaplus(sqrt_s_);
2166  },
2167  sqrt_s_, type_n_bar, type_K_z);
2168  break;
2169  }
2170  case pack(pdg::Lambda, pdg::pi_z): {
2171  const auto& type_p = ParticleType::find(pdg::p);
2172  const auto& type_n = ParticleType::find(pdg::n);
2173  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2174  const auto& type_K_m = ParticleType::find(pdg::K_m);
2175  add_channel(
2176  process_list,
2177  [&] {
2178  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2179  type_p, type_K_m) *
2180  kminusp_pi0lambda(sqrt_s_);
2181  },
2182  sqrt_s_, type_p, type_K_m);
2183  add_channel(
2184  process_list,
2185  [&] {
2186  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2187  type_n, type_Kbar_z) *
2188  kminusp_pi0lambda(sqrt_s_);
2189  },
2190  sqrt_s_, type_n, type_Kbar_z);
2191  break;
2192  }
2193  case pack(-pdg::Lambda, pdg::pi_z): {
2194  const auto& type_p_bar = ParticleType::find(-pdg::p);
2195  const auto& type_n_bar = ParticleType::find(-pdg::n);
2196  const auto& type_K_z = ParticleType::find(pdg::K_z);
2197  const auto& type_K_p = ParticleType::find(pdg::K_p);
2198  add_channel(
2199  process_list,
2200  [&] {
2201  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2202  type_p_bar, type_K_p) *
2203  kminusp_pi0lambda(sqrt_s_);
2204  },
2205  sqrt_s_, type_p_bar, type_K_p);
2206  add_channel(
2207  process_list,
2208  [&] {
2209  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2210  type_n_bar, type_K_z) *
2211  kminusp_pi0lambda(sqrt_s_);
2212  },
2213  sqrt_s_, type_n_bar, type_K_z);
2214  break;
2215  }
2216  case pack(pdg::Sigma_p, pdg::pi_m): {
2217  const auto& type_p = ParticleType::find(pdg::p);
2218  const auto& type_n = ParticleType::find(pdg::n);
2219  const auto& type_K_m = ParticleType::find(pdg::K_m);
2220  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2221  add_channel(
2222  process_list,
2223  [&] {
2224  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2225  type_p, type_K_m) *
2226  kminusp_piminussigmaplus(sqrt_s_);
2227  },
2228  sqrt_s_, type_p, type_K_m);
2229  add_channel(
2230  process_list,
2231  [&] {
2232  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2233  type_n, type_Kbar_z) *
2234  kminusp_piplussigmaminus(sqrt_s_);
2235  },
2236  sqrt_s_, type_n, type_Kbar_z);
2237  break;
2238  }
2239  case pack(-pdg::Sigma_p, pdg::pi_p): {
2240  const auto& type_p_bar = ParticleType::find(-pdg::p);
2241  const auto& type_n_bar = ParticleType::find(-pdg::n);
2242  const auto& type_K_p = ParticleType::find(pdg::K_p);
2243  const auto& type_K_z = ParticleType::find(pdg::K_z);
2244  add_channel(
2245  process_list,
2246  [&] {
2247  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2248  type_p_bar, type_K_p) *
2249  kminusp_piminussigmaplus(sqrt_s_);
2250  },
2251  sqrt_s_, type_p_bar, type_K_p);
2252  add_channel(
2253  process_list,
2254  [&] {
2255  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2256  type_n_bar, type_K_z) *
2257  kminusp_piplussigmaminus(sqrt_s_);
2258  },
2259  sqrt_s_, type_n_bar, type_K_z);
2260  break;
2261  }
2262  default:
2263  break;
2264  }
2265 
2266  return process_list;
2267 }
2268 
2269 double CrossSections::xs_dpi_dprimepi(const double sqrts, const double cm_mom,
2270  ParticleTypePtr produced_nucleus,
2271  const ParticleType& type_pi) {
2272  const double s = sqrts * sqrts;
2273  // same matrix element for πd and πd̅
2274  const double tmp = sqrts - pion_mass - deuteron_mass;
2275  // Matrix element is fit to match the inelastic pi+ d -> pi+ n p
2276  // cross-section from the Fig. 5 of [\iref{Arndt:1994bs}].
2277  const double matrix_element =
2278  295.5 + 2.862 / (0.00283735 + pow_int(sqrts - 2.181, 2)) +
2279  0.0672 / pow_int(tmp, 2) - 6.61753 / tmp;
2280 
2281  const double spin_factor =
2282  (produced_nucleus->spin() + 1) * (type_pi.spin() + 1);
2283  /* Isospin factor is always the same, so it is included into the
2284  * matrix element.
2285  * Symmetry factor is always 1 here.
2286  * The (hbarc)^2/16 pi factor is absorbed into matrix element. */
2287  double xsection = matrix_element * spin_factor / (s * cm_mom);
2288  if (produced_nucleus->is_stable()) {
2289  xsection *= pCM_from_s(s, type_pi.mass(), produced_nucleus->mass());
2290  } else {
2291  const double resonance_integral =
2292  produced_nucleus->iso_multiplet()->get_integral_piR(sqrts);
2293  xsection *= resonance_integral;
2294  logg[LScatterAction].debug("Resonance integral ", resonance_integral,
2295  ", matrix element: ", matrix_element,
2296  ", cm_momentum: ", cm_mom);
2297  }
2298  return xsection;
2299 }
2300 
2301 CollisionBranchList CrossSections::dpi_xx(
2302  const ReactionsBitSet& included_2to2) const {
2303  CollisionBranchList process_list;
2304  const double sqrts = sqrt_s_;
2305  const ParticleType& type_a = incoming_particles_[0].type();
2306  const ParticleType& type_b = incoming_particles_[1].type();
2307 
2308  // pi d -> N N
2309  bool is_pid = (type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
2310  (type_b.is_deuteron() && type_a.pdgcode().is_pion());
2311  if (is_pid && included_2to2[IncludedReactions::PiDeuteron_to_NN] == 1) {
2312  const int baryon_number = type_a.baryon_number() + type_b.baryon_number();
2313  ParticleTypePtrList nuc = (baryon_number > 0)
2314  ? ParticleType::list_nucleons()
2315  : ParticleType::list_anti_nucleons();
2316  const double s = sqrt_s_ * sqrt_s_;
2317  for (ParticleTypePtr nuc_a : nuc) {
2318  for (ParticleTypePtr nuc_b : nuc) {
2319  if (type_a.charge() + type_b.charge() !=
2320  nuc_a->charge() + nuc_b->charge()) {
2321  continue;
2322  }
2323  // loop over total isospin
2324  for (const int twoI : I_tot_range(*nuc_a, *nuc_b)) {
2325  const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
2326  type_a, type_b, *nuc_a, *nuc_b, twoI);
2327  // If Clebsch-Gordan coefficient = 0, don't bother with the rest.
2328  if (std::abs(isospin_factor) < really_small) {
2329  continue;
2330  }
2331 
2332  // Calculate matrix element for inverse process.
2333  const double matrix_element =
2334  nn_to_resonance_matrix_element(sqrts, type_a, type_b, twoI);
2335  if (matrix_element <= 0.) {
2336  continue;
2337  }
2338 
2339  const double spin_factor = (nuc_a->spin() + 1) * (nuc_b->spin() + 1);
2340  const int sym_fac_in =
2341  (type_a.iso_multiplet() == type_b.iso_multiplet()) ? 2 : 1;
2342  const int sym_fac_out =
2343  (nuc_a->iso_multiplet() == nuc_b->iso_multiplet()) ? 2 : 1;
2344  double p_cm_final = pCM_from_s(s, nuc_a->mass(), nuc_b->mass());
2345  const double xsection = isospin_factor * spin_factor * sym_fac_in /
2346  sym_fac_out * p_cm_final * matrix_element /
2347  (s * cm_momentum());
2348 
2349  if (xsection > really_small) {
2350  process_list.push_back(std::make_unique<CollisionBranch>(
2351  *nuc_a, *nuc_b, xsection, ProcessType::TwoToTwo));
2352  logg[LScatterAction].debug(type_a.name(), type_b.name(), "->",
2353  nuc_a->name(), nuc_b->name(),
2354  " at sqrts [GeV] = ", sqrts,
2355  " with cs[mb] = ", xsection);
2356  }
2357  }
2358  }
2359  }
2360  }
2361 
2362  // pi d -> pi d' (effectively pi d -> pi p n) AND reverse, pi d' -> pi d
2363  bool is_pid_or_pidprime = ((type_a.is_deuteron() || type_a.is_dprime()) &&
2364  type_b.pdgcode().is_pion()) ||
2365  ((type_b.is_deuteron() || type_b.is_dprime()) &&
2366  type_a.pdgcode().is_pion());
2367  if (is_pid_or_pidprime &&
2368  included_2to2[IncludedReactions::PiDeuteron_to_pidprime] == 1) {
2369  const ParticleType& type_pi = type_a.pdgcode().is_pion() ? type_a : type_b;
2370  const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
2371  ParticleTypePtrList nuclei = ParticleType::list_light_nuclei();
2372  for (ParticleTypePtr produced_nucleus : nuclei) {
2373  // Elastic collisions are treated in a different function
2374  if (produced_nucleus == &type_nucleus ||
2375  produced_nucleus->charge() != type_nucleus.charge() ||
2376  produced_nucleus->baryon_number() != type_nucleus.baryon_number()) {
2377  continue;
2378  }
2379  const double xsection =
2380  xs_dpi_dprimepi(sqrts, cm_momentum(), produced_nucleus, type_pi);
2381  process_list.push_back(std::make_unique<CollisionBranch>(
2382  type_pi, *produced_nucleus, xsection, ProcessType::TwoToTwo));
2383  logg[LScatterAction].debug(type_pi.name(), type_nucleus.name(), "→ ",
2384  type_pi.name(), produced_nucleus->name(),
2385  " at ", sqrts, " GeV, xs[mb] = ", xsection);
2386  }
2387  }
2388  return process_list;
2389 }
2390 
2391 double CrossSections::xs_dn_dprimen(const double sqrts, const double cm_mom,
2392  ParticleTypePtr produced_nucleus,
2393  const ParticleType& type_nucleus,
2394  const ParticleType& type_N) {
2395  const double s = sqrts * sqrts;
2396  double matrix_element = 0.0;
2397  double tmp = sqrts - nucleon_mass - deuteron_mass;
2398  assert(tmp >= 0.0);
2399  if (std::signbit(type_N.baryon_number()) ==
2400  std::signbit(type_nucleus.baryon_number())) {
2404  matrix_element = 79.0474 / std::pow(tmp, 0.7897) + 654.596 * tmp;
2405  } else {
2409  matrix_element = 342.572 / std::pow(tmp, 0.6);
2410  }
2411  const double spin_factor =
2412  (produced_nucleus->spin() + 1) * (type_N.spin() + 1);
2413  /* Isospin factor is always the same, so it is included into matrix element
2414  * Symmetry factor is always 1 here
2415  * Absorb (hbarc)^2/16 pi factor into matrix element */
2416  double xsection = matrix_element * spin_factor / (s * cm_mom);
2417  if (produced_nucleus->is_stable()) {
2418  assert(!type_nucleus.is_stable());
2419  xsection *= pCM_from_s(s, type_N.mass(), produced_nucleus->mass());
2420  } else {
2421  assert(type_nucleus.is_stable());
2422  const double resonance_integral =
2423  produced_nucleus->iso_multiplet()->get_integral_NR(sqrts);
2424  xsection *= resonance_integral;
2425  }
2426  return xsection;
2427 }
2428 
2429 CollisionBranchList CrossSections::dn_xx(
2430  const ReactionsBitSet& included_2to2) const {
2431  const ParticleType& type_a = incoming_particles_[0].type();
2432  const ParticleType& type_b = incoming_particles_[1].type();
2433  const ParticleType& type_N = type_a.is_nucleon() ? type_a : type_b;
2434  const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
2435  CollisionBranchList process_list;
2436  if (included_2to2[IncludedReactions::NDeuteron_to_Ndprime] == 0) {
2437  return process_list;
2438  }
2439  ParticleTypePtrList nuclei = ParticleType::list_light_nuclei();
2440 
2441  for (ParticleTypePtr produced_nucleus : nuclei) {
2442  // No elastic collisions for now, respect conservation laws
2443  if (produced_nucleus == &type_nucleus ||
2444  produced_nucleus->charge() != type_nucleus.charge() ||
2445  produced_nucleus->baryon_number() != type_nucleus.baryon_number()) {
2446  continue;
2447  }
2448  const double xsection = xs_dn_dprimen(
2449  sqrt_s_, cm_momentum(), produced_nucleus, type_nucleus, type_N);
2450  process_list.push_back(std::make_unique<CollisionBranch>(
2451  type_N, *produced_nucleus, xsection, ProcessType::TwoToTwo));
2452  logg[LScatterAction].debug(type_N.name(), type_nucleus.name(), "→ ",
2453  type_N.name(), produced_nucleus->name(), " at ",
2454  sqrt_s_, " GeV, xs[mb] = ", xsection);
2455  }
2456  return process_list;
2457 }
2458 
2459 CollisionBranchList CrossSections::string_excitation(
2460  double total_string_xs, StringProcess* string_process, bool use_AQM) const {
2461  if (!string_process) {
2462  throw std::runtime_error("string_process should be initialized.");
2463  }
2464 
2465  CollisionBranchList channel_list;
2466  if (total_string_xs <= 0.) {
2467  return channel_list;
2468  }
2469 
2470  /* Get mapped PDG id for evaluation of the parametrized cross sections
2471  * for diffractive processes.
2472  * This must be rescaled according to additive quark model
2473  * in the case of exotic hadrons.
2474  * Also calculate the multiplicative factor for AQM
2475  * based on the quark contents. */
2476  std::array<int, 2> pdgid;
2477  double AQM_factor = 1.;
2478  for (int i = 0; i < 2; i++) {
2479  PdgCode pdg = incoming_particles_[i].type().pdgcode();
2480  pdgid[i] = StringProcess::pdg_map_for_pythia(pdg);
2481  AQM_factor *= (1. - 0.4 * pdg.frac_strange());
2482  }
2483 
2484  /* Determine if the initial state is a baryon-antibaryon pair,
2485  * which can annihilate. */
2486  bool can_annihilate = false;
2487  if (is_BBbar_pair_) {
2488  int n_q_types = 5; // u, d, s, c, b
2489  for (int iq = 1; iq <= n_q_types; iq++) {
2490  std::array<int, 2> nquark;
2491  for (int i = 0; i < 2; i++) {
2492  nquark[i] =
2493  incoming_particles_[i].type().pdgcode().net_quark_number(iq);
2494  }
2495  if (nquark[0] != 0 && nquark[1] != 0) {
2496  can_annihilate = true;
2497  break;
2498  }
2499  }
2500  }
2501 
2502  /* Total parametrized cross-section (I) and pythia-produced total
2503  * cross-section (II) do not necessarily coincide. If I > II then
2504  * non-diffractive cross-section is reinforced to get I == II.
2505  * If I < II then partial cross-sections are drained one-by-one
2506  * to reduce II until I == II:
2507  * first non-diffractive, then double-diffractive, then
2508  * single-diffractive AB->AX and AB->XB in equal proportion.
2509  * The way it is done here is not unique. I (ryu) think that at high energy
2510  * collision this is not an issue, but at sqrt_s < 10 GeV it may
2511  * matter. */
2512  std::array<double, 3> xs =
2513  string_process->cross_sections_diffractive(pdgid[0], pdgid[1], sqrt_s_);
2514  if (use_AQM) {
2515  for (int ip = 0; ip < 3; ip++) {
2516  xs[ip] *= AQM_factor;
2517  }
2518  }
2519  double single_diffr_AX = xs[0], single_diffr_XB = xs[1], double_diffr = xs[2];
2520  double single_diffr = single_diffr_AX + single_diffr_XB;
2521  double diffractive = single_diffr + double_diffr;
2522 
2523  /* The case for baryon/anti-baryon annihilation is treated separately,
2524  * as in this case we use only one way to break up the particles, namely
2525  * into 2 mesonic strings of equal mass after annihilating one quark-
2526  * anti-quark pair. See StringProcess::next_BBbarAnn() */
2527  double sig_annihilation = 0.0;
2528  if (can_annihilate) {
2529  /* In the case of baryon-antibaryon pair,
2530  * the parametrized cross section for annihilation will be added.
2531  * See xs_ppbar_annihilation(). */
2532  double xs_param = xs_ppbar_annihilation(sqrt_s_ * sqrt_s_);
2533  if (use_AQM) {
2534  xs_param *= AQM_factor;
2535  }
2536  sig_annihilation = std::min(total_string_xs, xs_param);
2537  }
2538 
2539  const double nondiffractive_all =
2540  std::max(0., total_string_xs - sig_annihilation - diffractive);
2541  diffractive = total_string_xs - sig_annihilation - nondiffractive_all;
2542  double_diffr = std::max(0., diffractive - single_diffr);
2543  const double a = (diffractive - double_diffr) / single_diffr;
2544  single_diffr_AX *= a;
2545  single_diffr_XB *= a;
2546  assert(std::abs(single_diffr_AX + single_diffr_XB + double_diffr +
2547  sig_annihilation + nondiffractive_all - total_string_xs) <
2548  1.e-6);
2549 
2550  double nondiffractive_soft = 0.;
2551  double nondiffractive_hard = 0.;
2552  if (nondiffractive_all > 0.) {
2553  /* Hard string process is added by hard cross section
2554  * in conjunction with multipartion interaction picture
2555  * \iref{Sjostrand:1987su}. */
2556  const double hard_xsec = AQM_factor * string_hard_cross_section();
2557  nondiffractive_soft =
2558  nondiffractive_all * std::exp(-hard_xsec / nondiffractive_all);
2559  nondiffractive_hard = nondiffractive_all - nondiffractive_soft;
2560  }
2561  logg[LCrossSections].debug("String cross sections [mb] are");
2562  logg[LCrossSections].debug("Single-diffractive AB->AX: ", single_diffr_AX);
2563  logg[LCrossSections].debug("Single-diffractive AB->XB: ", single_diffr_XB);
2564  logg[LCrossSections].debug("Double-diffractive AB->XX: ", double_diffr);
2565  logg[LCrossSections].debug("Soft non-diffractive: ", nondiffractive_soft);
2566  logg[LCrossSections].debug("Hard non-diffractive: ", nondiffractive_hard);
2567  logg[LCrossSections].debug("B-Bbar annihilation: ", sig_annihilation);
2568 
2569  /* cross section of soft string excitation including annihilation */
2570  const double sig_string_soft = total_string_xs - nondiffractive_hard;
2571 
2572  /* fill the list of process channels */
2573  if (sig_string_soft > 0.) {
2574  channel_list.push_back(std::make_unique<CollisionBranch>(
2575  single_diffr_AX, ProcessType::StringSoftSingleDiffractiveAX));
2576  channel_list.push_back(std::make_unique<CollisionBranch>(
2577  single_diffr_XB, ProcessType::StringSoftSingleDiffractiveXB));
2578  channel_list.push_back(std::make_unique<CollisionBranch>(
2579  double_diffr, ProcessType::StringSoftDoubleDiffractive));
2580  channel_list.push_back(std::make_unique<CollisionBranch>(
2581  nondiffractive_soft, ProcessType::StringSoftNonDiffractive));
2582  if (can_annihilate) {
2583  channel_list.push_back(std::make_unique<CollisionBranch>(
2584  sig_annihilation, ProcessType::StringSoftAnnihilation));
2585  }
2586  }
2587  if (nondiffractive_hard > 0.) {
2588  channel_list.push_back(std::make_unique<CollisionBranch>(
2589  nondiffractive_hard, ProcessType::StringHard));
2590  }
2591  return channel_list;
2592 }
2593 
2594 double CrossSections::high_energy(
2595  const StringTransitionParameters& transition_high_energy) const {
2596  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
2597  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
2598 
2599  const double s = sqrt_s_ * sqrt_s_;
2600  double xs = 0.;
2601 
2602  // Currently all BB collisions use the nucleon-nucleon parametrizations.
2603  if (pdg_a.is_baryon() && pdg_b.is_baryon()) {
2604  if (pdg_a == pdg_b) {
2605  xs = pp_high_energy(s); // pp, nn
2606  } else if (pdg_a.antiparticle_sign() * pdg_b.antiparticle_sign() == 1) {
2607  xs = np_high_energy(s); // np, nbarpbar
2608  } else if (pdg_a.antiparticle_sign() * pdg_b.antiparticle_sign() == -1) {
2609  /* In the case of baryon-antibaryon interactions,
2610  * the low-energy cross section must be involved
2611  * due to annihilation processes (via strings). */
2612  double xs_l = ppbar_total(s);
2613  double xs_h = 0.;
2614  if (pdg_a.is_antiparticle_of(pdg_b)) {
2615  xs_h = ppbar_high_energy(s); // ppbar, nnbar
2616  } else {
2617  xs_h = npbar_high_energy(s); // npbar, nbarp
2618  }
2619  /* Transition between low and high energy is set to be consistent with
2620  * that defined in string_probability(). */
2621  auto [region_lower, region_upper] = transition_high_energy.sqrts_range_NN;
2622  double prob_high = probability_transit_high(region_lower, region_upper);
2623  xs = xs_l * (1. - prob_high) + xs_h * prob_high;
2624  }
2625  }
2626 
2627  // Pion nucleon interaction / baryon-meson
2628  if ((pdg_a == pdg::pi_p && pdg_b == pdg::p) ||
2629  (pdg_b == pdg::pi_p && pdg_a == pdg::p) ||
2630  (pdg_a == pdg::pi_m && pdg_b == pdg::n) ||
2631  (pdg_b == pdg::pi_m && pdg_a == pdg::n)) {
2632  xs = piplusp_high_energy(s); // pi+ p, pi- n
2633  } else if ((pdg_a == pdg::pi_m && pdg_b == pdg::p) ||
2634  (pdg_b == pdg::pi_m && pdg_a == pdg::p) ||
2635  (pdg_a == pdg::pi_p && pdg_b == pdg::n) ||
2636  (pdg_b == pdg::pi_p && pdg_a == pdg::n)) {
2637  xs = piminusp_high_energy(s); // pi- p, pi+ n
2638  } else if ((pdg_a.is_meson() && pdg_b.is_baryon()) ||
2639  (pdg_b.is_meson() && pdg_a.is_baryon())) {
2640  xs = piminusp_high_energy(s); // default for baryon-meson
2641  }
2642 
2643  /* Meson-meson interaction goes through AQM from pi+p,
2644  * see user guide "Use_AQM"*/
2645  if (pdg_a.is_meson() && pdg_b.is_meson()) {
2646  /* 2/3 factor since difference of 1 meson between meson-meson
2647  * and baryon-meson */
2648  xs = 2. / 3. * piplusp_high_energy(s);
2649  }
2650 
2651  // AQM scaling for cross-sections
2652  xs *= (1. - 0.4 * pdg_a.frac_strange()) * (1. - 0.4 * pdg_b.frac_strange());
2653 
2654  return xs;
2655 }
2656 
2657 double CrossSections::string_hard_cross_section() const {
2658  double cross_sec = 0.;
2659  /* Hard strings can only be excited if the lower cutoff by
2660  * Pythia is fulfilled */
2661  if (sqrt_s_ <= minimum_sqrts_pythia_can_handle) {
2662  return cross_sec;
2663  }
2664  const ParticleData& data_a = incoming_particles_[0];
2665  const ParticleData& data_b = incoming_particles_[1];
2666 
2667  if (data_a.is_baryon() && data_b.is_baryon()) {
2668  // Nucleon-nucleon cross section is used for all baryon-baryon cases.
2669  cross_sec = NN_string_hard(sqrt_s_ * sqrt_s_);
2670  } else if (data_a.is_baryon() || data_b.is_baryon()) {
2671  // Nucleon-pion cross section is used for all baryon-meson cases.
2672  cross_sec = Npi_string_hard(sqrt_s_ * sqrt_s_);
2673  } else {
2674  // Pion-pion cross section is used for all meson-meson cases.
2675  cross_sec = pipi_string_hard(sqrt_s_ * sqrt_s_);
2676  }
2677 
2678  return cross_sec;
2679 }
2680 
2681 CollisionBranchPtr CrossSections::NNbar_to_5pi(const double scale_xs) const {
2682  const double s = sqrt_s_ * sqrt_s_;
2683  /* Use difference between total and elastic in order to conserve detailed
2684  * balance for all inelastoc NNbar processes. */
2685  const double nnbar_xsec = std::max(0., ppbar_total(s) - ppbar_elastic(s));
2686  logg[LCrossSections].debug("NNbar cross section for 2-to-5 is: ", nnbar_xsec);
2687 
2688  /* Make collision channel NNbar -> 5π (with same final state as resonance
2689  * approach). */
2690  const auto& type_piz = ParticleType::find(pdg::pi_z);
2691  const auto& type_pip = ParticleType::find(pdg::pi_p);
2692  const auto& type_pim = ParticleType::find(pdg::pi_m);
2693  return std::make_unique<CollisionBranch>(
2694  type_pip, type_pim, type_pip, type_pim, type_piz, nnbar_xsec * scale_xs,
2695  ProcessType::TwoToFive);
2696 }
2697 
2698 CollisionBranchPtr CrossSections::NNbar_annihilation(
2699  const double current_xs, const double scale_xs) const {
2700  /* Calculate NNbar cross section:
2701  * Parametrized total minus all other present channels.*/
2702  const double s = sqrt_s_ * sqrt_s_;
2703  double nnbar_xsec = std::max(0., ppbar_total(s) * scale_xs - current_xs);
2704  logg[LCrossSections].debug("NNbar cross section is: ", nnbar_xsec);
2705  // Make collision channel NNbar -> ρh₁(1170); eventually decays into 5π
2706  return std::make_unique<CollisionBranch>(ParticleType::find(pdg::h1),
2707  ParticleType::find(pdg::rho_z),
2708  nnbar_xsec, ProcessType::TwoToTwo);
2709 }
2710 
2711 CollisionBranchList CrossSections::NNbar_creation() const {
2712  CollisionBranchList channel_list;
2713  const ParticleType& type_a = incoming_particles_[0].type();
2714  const ParticleType& type_b = incoming_particles_[1].type();
2715  if ((type_a.pdgcode() == pdg::rho_z && type_b.pdgcode() == pdg::h1) ||
2716  (type_a.pdgcode() == pdg::h1 && type_b.pdgcode() == pdg::rho_z)) {
2717  /* Calculate NNbar reverse cross section:
2718  * from reverse reaction (see NNbar_annihilation_cross_section).*/
2719  const double s = sqrt_s_ * sqrt_s_;
2720  const double pcm = cm_momentum();
2721 
2722  const auto& type_N = ParticleType::find(pdg::p);
2723  const auto& type_Nbar = ParticleType::find(-pdg::p);
2724 
2725  // Check available energy
2726  if (sqrt_s_ - 2 * type_N.mass() < 0) {
2727  return channel_list;
2728  }
2729 
2730  double xsection = detailed_balance_factor_RR(sqrt_s_, pcm, type_a, type_b,
2731  type_N, type_Nbar) *
2732  std::max(0., ppbar_total(s) - ppbar_elastic(s));
2733  logg[LCrossSections].debug("NNbar reverse cross section is: ", xsection);
2734  channel_list.push_back(std::make_unique<CollisionBranch>(
2735  type_N, type_Nbar, xsection, ProcessType::TwoToTwo));
2736  channel_list.push_back(std::make_unique<CollisionBranch>(
2737  ParticleType::find(pdg::n), ParticleType::find(-pdg::n), xsection,
2738  ProcessType::TwoToTwo));
2739  }
2740  return channel_list;
2741 }
2742 
2743 CollisionBranchList CrossSections::bar_bar_to_nuc_nuc(
2744  const bool is_anti_particles) const {
2745  const ParticleType& type_a = incoming_particles_[0].type();
2746  const ParticleType& type_b = incoming_particles_[1].type();
2747  CollisionBranchList process_list;
2748 
2749  const double s = sqrt_s_ * sqrt_s_;
2750  // CM momentum in final state
2751  double p_cm_final = std::sqrt(s - 4. * nucleon_mass * nucleon_mass) / 2.;
2752 
2753  ParticleTypePtrList nuc_or_anti_nuc;
2754  if (is_anti_particles) {
2755  nuc_or_anti_nuc = ParticleType::list_anti_nucleons();
2756  } else {
2757  nuc_or_anti_nuc = ParticleType::list_nucleons();
2758  }
2759 
2760  // Loop over all nucleon or anti-nucleon charge states.
2761  for (ParticleTypePtr nuc_a : nuc_or_anti_nuc) {
2762  for (ParticleTypePtr nuc_b : nuc_or_anti_nuc) {
2763  /* Check for charge conservation. */
2764  if (type_a.charge() + type_b.charge() !=
2765  nuc_a->charge() + nuc_b->charge()) {
2766  continue;
2767  }
2768  // loop over total isospin
2769  for (const int twoI : I_tot_range(*nuc_a, *nuc_b)) {
2770  const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
2771  type_a, type_b, *nuc_a, *nuc_b, twoI);
2772  // If Clebsch-Gordan coefficient is zero, don't bother with the rest
2773  if (std::abs(isospin_factor) < really_small) {
2774  continue;
2775  }
2776 
2777  // Calculate matrix element for inverse process.
2778  const double matrix_element =
2779  nn_to_resonance_matrix_element(sqrt_s_, type_a, type_b, twoI);
2780  if (matrix_element <= 0.) {
2781  continue;
2782  }
2783 
2788  const double spin_factor = (nuc_a->spin() + 1) * (nuc_b->spin() + 1);
2789  const int sym_fac_in =
2790  (type_a.iso_multiplet() == type_b.iso_multiplet()) ? 2 : 1;
2791  const int sym_fac_out =
2792  (nuc_a->iso_multiplet() == nuc_b->iso_multiplet()) ? 2 : 1;
2793  const double xsection = isospin_factor * spin_factor * sym_fac_in /
2794  sym_fac_out * p_cm_final * matrix_element /
2795  (s * cm_momentum());
2796 
2797  if (xsection > really_small) {
2798  process_list.push_back(std::make_unique<CollisionBranch>(
2799  *nuc_a, *nuc_b, xsection, ProcessType::TwoToTwo));
2800  logg[LCrossSections].debug(
2801  "2->2 absorption with original particles: ", type_a, type_b);
2802  }
2803  }
2804  }
2805  }
2806  return process_list;
2807 }
2808 
2809 double CrossSections::nn_to_resonance_matrix_element(double sqrts,
2810  const ParticleType& type_a,
2811  const ParticleType& type_b,
2812  const int twoI) {
2813  const double m_a = type_a.mass();
2814  const double m_b = type_b.mass();
2815  const double msqr = 2. * (m_a * m_a + m_b * m_b);
2816  /* If the c.m. energy is larger than the sum of the pole masses of the
2817  * outgoing particles plus three times of the sum of the widths plus 3 GeV,
2818  * the collision will be neglected.
2819  *
2820  * This can be problematic for some final-state cross sections, but at
2821  * energies that high strings are used anyway.
2822  */
2823  const double w_a = type_a.width_at_pole();
2824  const double w_b = type_b.width_at_pole();
2825  const double uplmt = m_a + m_b + 3.0 * (w_a + w_b) + 3.0;
2826  if (sqrts > uplmt) {
2827  return 0.;
2828  }
2830  if (((type_a.is_Delta() && type_b.is_nucleon()) ||
2831  (type_b.is_Delta() && type_a.is_nucleon())) &&
2832  (type_a.antiparticle_sign() == type_b.antiparticle_sign())) {
2833  return 68. / std::pow(sqrts - 1.104, 1.951);
2836  } else if (((type_a.is_Nstar() && type_b.is_nucleon()) ||
2837  (type_b.is_Nstar() && type_a.is_nucleon())) &&
2838  type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
2839  // NN → NN*
2840  if (twoI == 2) {
2841  return 4.5 / msqr;
2842  } else if (twoI == 0) {
2843  const double parametrization = 14. / msqr;
2849  if (type_a.is_Nstar1535() || type_b.is_Nstar1535()) {
2850  return 6.5 * parametrization;
2851  } else {
2852  return parametrization;
2853  }
2854  }
2855  } else if (((type_a.is_Deltastar() && type_b.is_nucleon()) ||
2856  (type_b.is_Deltastar() && type_a.is_nucleon())) &&
2857  type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
2858  // NN → NΔ*
2859  return 15. / msqr;
2860  } else if ((type_a.is_Delta() && type_b.is_Delta()) &&
2861  (type_a.antiparticle_sign() == type_b.antiparticle_sign())) {
2862  // NN → ΔΔ
2863  if (twoI == 2) {
2864  return 45. / msqr;
2865  } else if (twoI == 0) {
2866  return 120. / msqr;
2867  }
2868  } else if (((type_a.is_Nstar() && type_b.is_Delta()) ||
2869  (type_b.is_Nstar() && type_a.is_Delta())) &&
2870  type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
2871  // NN → ΔN*
2872  return 7. / msqr;
2873  } else if (((type_a.is_Deltastar() && type_b.is_Delta()) ||
2874  (type_b.is_Deltastar() && type_a.is_Delta())) &&
2875  type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
2876  // NN → ΔΔ*
2877  if (twoI == 2) {
2878  return 15. / msqr;
2879  } else if (twoI == 0) {
2880  return 25. / msqr;
2881  }
2882  } else if ((type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
2883  (type_b.is_deuteron() && type_a.pdgcode().is_pion())) {
2884  /* This parametrization is the result of fitting d+pi->NN cross-section.
2885  * Already Breit-Wigner-like part provides a good fit, exponential fixes
2886  * behaviour around the treshold. The d+pi experimental cross-section
2887  * was taken from Fig. 2 of [\iref{Tanabe:1987vg}]. */
2888  return 0.055 / (pow_int(sqrts - 2.145, 2) + pow_int(0.065, 2)) *
2889  (1.0 - std::exp(-(sqrts - 2.0) * 20.0));
2890  }
2891 
2892  // all cases not listed: zero!
2893  return 0.;
2894 }
2895 
2896 template <class IntegrationMethod>
2897 CollisionBranchList CrossSections::find_nn_xsection_from_type(
2898  const ParticleTypePtrList& list_res_1,
2899  const ParticleTypePtrList& list_res_2,
2900  const IntegrationMethod integrator) const {
2901  const ParticleType& type_particle_a = incoming_particles_[0].type();
2902  const ParticleType& type_particle_b = incoming_particles_[1].type();
2903 
2904  CollisionBranchList channel_list;
2905  const double s = sqrt_s_ * sqrt_s_;
2906 
2907  // Loop over specified first resonance list
2908  for (ParticleTypePtr type_res_1 : list_res_1) {
2909  // Loop over specified second resonance list
2910  for (ParticleTypePtr type_res_2 : list_res_2) {
2911  // Check for charge conservation.
2912  if (type_res_1->charge() + type_res_2->charge() !=
2913  type_particle_a.charge() + type_particle_b.charge()) {
2914  continue;
2915  }
2916 
2917  // loop over total isospin
2918  for (const int twoI : I_tot_range(type_particle_a, type_particle_b)) {
2919  const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
2920  type_particle_a, type_particle_b, *type_res_1, *type_res_2, twoI);
2921  // If Clebsch-Gordan coefficient is zero, don't bother with the rest.
2922  if (std::abs(isospin_factor) < really_small) {
2923  continue;
2924  }
2925 
2926  // Integration limits.
2927  const double lower_limit = type_res_1->min_mass_kinematic();
2928  const double upper_limit = sqrt_s_ - type_res_2->mass();
2929  /* Check the available energy (requiring it to be a little above the
2930  * threshold, because the integration will not work if it's too close).
2931  */
2932  if (upper_limit - lower_limit < 1E-3) {
2933  continue;
2934  }
2935 
2936  // Calculate matrix element.
2937  const double matrix_element = nn_to_resonance_matrix_element(
2938  sqrt_s_, *type_res_1, *type_res_2, twoI);
2939  if (matrix_element <= 0.) {
2940  continue;
2941  }
2942 
2943  /* Calculate resonance production cross section
2944  * using the Breit-Wigner distribution as probability amplitude.
2945  * Integrate over the allowed resonance mass range. */
2946  const double resonance_integral = integrator(*type_res_1, *type_res_2);
2947 
2951  const double spin_factor =
2952  (type_res_1->spin() + 1) * (type_res_2->spin() + 1);
2953  const double xsection = isospin_factor * spin_factor * matrix_element *
2954  resonance_integral / (s * cm_momentum());
2955 
2956  if (xsection > really_small) {
2957  channel_list.push_back(std::make_unique<CollisionBranch>(
2958  *type_res_1, *type_res_2, xsection, ProcessType::TwoToTwo));
2959  logg[LCrossSections].debug(
2960  "Found 2->2 creation process for resonance ", type_res_1, ", ",
2961  type_res_2);
2962  logg[LCrossSections].debug("2->2 with original particles: ",
2963  type_particle_a, type_particle_b);
2964  }
2965  }
2966  }
2967  }
2968  return channel_list;
2969 }
2970 
2971 double CrossSections::string_probability(
2972  const ScatterActionsFinderParameters& finder_parameters) const {
2973  /* string fragmentation is enabled when strings_switch is on and the process
2974  * is included in pythia. */
2975  if (!finder_parameters.strings_switch) {
2976  return 0.;
2977  }
2978 
2979  const ParticleType& t1 = incoming_particles_[0].type();
2980  const ParticleType& t2 = incoming_particles_[1].type();
2981  const bool treat_BBbar_with_strings =
2982  (finder_parameters.nnbar_treatment == NNbarTreatment::Strings);
2983  const bool is_NN_scattering =
2984  t1.is_nucleon() && t2.is_nucleon() &&
2985  t1.antiparticle_sign() == t2.antiparticle_sign();
2986  const bool is_BBbar_scattering =
2987  (treat_BBbar_with_strings && is_BBbar_pair_ &&
2988  finder_parameters.use_AQM) ||
2989  (t1.is_nucleon() && t2.is_nucleon() &&
2990  t1.antiparticle_sign() != t2.antiparticle_sign());
2991  const bool is_Npi_scattering = (t1.pdgcode().is_pion() && t2.is_nucleon()) ||
2992  (t1.is_nucleon() && t2.pdgcode().is_pion());
2993  /* True for baryon-baryon, anti-baryon-anti-baryon, baryon-meson,
2994  * anti-baryon-meson and meson-meson*/
2995  const bool is_AQM_scattering =
2996  finder_parameters.use_AQM &&
2997  ((t1.is_baryon() && t2.is_baryon() &&
2998  t1.antiparticle_sign() == t2.antiparticle_sign()) ||
2999  ((t1.is_baryon() && t2.is_meson()) ||
3000  (t2.is_baryon() && t1.is_meson())) ||
3001  (t1.is_meson() && t2.is_meson()));
3002  const double mass_sum =
3003  incoming_particles_[0].pole_mass() + incoming_particles_[1].pole_mass();
3004 
3005  if (!is_NN_scattering && !is_BBbar_scattering && !is_Npi_scattering &&
3006  !is_AQM_scattering) {
3007  return 0.;
3008  } else if (is_NNbar_pair_ && !treat_BBbar_with_strings) {
3009  return 0.;
3010  } else if (is_BBbar_scattering) {
3011  // BBbar only goes through strings, so there are no "window" considerations
3012  return 1.;
3013  } else {
3014  /* true for K+ p and K0 p (+ antiparticles), which have special treatment
3015  * to fit data */
3016  const PdgCode pdg1 = t1.pdgcode(), pdg2 = t2.pdgcode();
3017  const bool is_KplusP =
3018  ((pdg1 == pdg::K_p || pdg1 == pdg::K_z) && (pdg2 == pdg::p)) ||
3019  ((pdg2 == pdg::K_p || pdg2 == pdg::K_z) && (pdg1 == pdg::p)) ||
3020  ((pdg1 == -pdg::K_p || pdg1 == -pdg::K_z) && (pdg2 == -pdg::p)) ||
3021  ((pdg2 == -pdg::K_p || pdg2 == -pdg::K_z) && (pdg1 == -pdg::p));
3022  // where to start the AQM strings above mass sum
3023  double aqm_offset =
3024  finder_parameters.transition_high_energy.sqrts_add_lower;
3025  if (is_KplusP) {
3026  /* for this specific case we have data. This corresponds to the point
3027  * where the AQM parametrization is smaller than the current 2to2
3028  * parametrization, which starts growing and diverges from exp. data */
3029  aqm_offset = finder_parameters.transition_high_energy.KN_offset;
3030  } else if (pdg1.is_pion() && pdg2.is_pion()) {
3031  aqm_offset = finder_parameters.transition_high_energy.pipi_offset;
3032  }
3033  /* if we do not use the probability transition algorithm, this is always a
3034  * string contribution if the energy is large enough */
3035  if (!finder_parameters.strings_with_probability) {
3036  return static_cast<double>(sqrt_s_ > mass_sum + aqm_offset);
3037  }
3038  /* No strings at low energy, only strings at high energy and
3039  * a transition region in the middle. Determine transition region: */
3040  double region_lower, region_upper;
3041  if (is_Npi_scattering) {
3042  std::tie(region_lower, region_upper) =
3043  finder_parameters.transition_high_energy.sqrts_range_Npi;
3044  } else if (is_NN_scattering) {
3045  std::tie(region_lower, region_upper) =
3046  finder_parameters.transition_high_energy.sqrts_range_NN;
3047  } else { // AQM - Additive Quark Model
3048  /* Transition region around 0.9 larger than the sum of pole masses;
3049  * highly arbitrary, feel free to improve */
3050  region_lower = mass_sum + aqm_offset;
3051  region_upper = mass_sum + aqm_offset +
3052  finder_parameters.transition_high_energy.sqrts_range_width;
3053  }
3054 
3055  if (sqrt_s_ > region_upper) {
3056  return 1.;
3057  } else if (sqrt_s_ < region_lower) {
3058  return 0.;
3059  } else {
3060  // Rescale transition region to [-1, 1]
3061  return probability_transit_high(region_lower, region_upper);
3062  }
3063  }
3064 }
3065 
3066 double CrossSections::probability_transit_high(
3067  const double region_lower, const double region_upper) const {
3068  if (sqrt_s_ < region_lower) {
3069  return 0.;
3070  }
3071 
3072  if (sqrt_s_ > region_upper) {
3073  return 1.;
3074  }
3075 
3076  double x = (sqrt_s_ - 0.5 * (region_lower + region_upper)) /
3077  (region_upper - region_lower);
3078  assert(x >= -0.5 && x <= 0.5);
3079  double prob = 0.5 * (std::sin(M_PI * x) + 1.0);
3080  assert(prob >= 0. && prob <= 1.);
3081 
3082  return prob;
3083 }
3084 
3085 } // namespace smash
smash::CrossSections::sqrt_s_
const double sqrt_s_
Total energy in the center-of-mass frame.
Definition: crosssections.h:585
smash::CrossSections::d_N_inelastic_xs
static double d_N_inelastic_xs(double N_kinetic_energy)
Parametrization of deuteron-nucleon inelastic cross section.
Definition: crosssections.cc:1065
smash::CrossSections::xs_dpi_dprimepi
static double xs_dpi_dprimepi(double sqrts, double cm_mom, ParticleTypePtr produced_nucleus, const ParticleType &type_pi)
Parametrized cross section for πd→ πd' (mockup for πd→ πnp), πd̅→ πd̅' and reverse,...
Definition: crosssections.cc:2269
smash::CrossSections::two_to_two
CollisionBranchList two_to_two(const ReactionsBitSet &included_2to2, double KN_offset) const
Find all inelastic 2->2 processes for the given scattering.
Definition: crosssections.cc:902
smash::CrossSections::NNbar_to_5pi
CollisionBranchPtr NNbar_to_5pi(double scale_xs) const
Create collision branch for NNbar annihilation going directly into 5 pions.
Definition: crosssections.cc:2681
smash::CrossSections::probability_transit_high
double probability_transit_high(double region_lower, double region_upper) const
Definition: crosssections.cc:3066
smash::CrossSections::NNbar_creation
CollisionBranchList NNbar_creation() const
Determine the cross section for NNbar creation, which is given by detailed balance from the reverse r...
Definition: crosssections.cc:2711
smash::CrossSections::dpi_xx
CollisionBranchList dpi_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes involving Pion and (anti-) Deuteron (dpi), specifically dπ→ NN,...
Definition: crosssections.cc:2301
smash::CrossSections::bb_xx_except_nn
CollisionBranchList bb_xx_except_nn(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes for Baryon-Baryon (BB) Scattering except the more specific Nucleon-...
Definition: crosssections.cc:1151
smash::CrossSections::npi_yk
CollisionBranchList npi_yk() const
Find all processes for Nucleon-Pion to Hyperon-Kaon Scattering.
Definition: crosssections.cc:516
smash::CrossSections::find_nn_xsection_from_type
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().
Definition: crosssections.cc:2897
smash::CrossSections::cm_momentum
double cm_momentum() const
Determine the momenta of the incoming particles in the center-of-mass system.
Definition: crosssections.h:575
smash::CrossSections::bar_bar_to_nuc_nuc
CollisionBranchList bar_bar_to_nuc_nuc(const bool is_anti_particles) const
Calculate cross sections for resonance absorption (i.e.
Definition: crosssections.cc:2743
smash::CrossSections::string_probability
double string_probability(const ScatterActionsFinderParameters &finder_parameters) const
Definition: crosssections.cc:2971
smash::CrossSections::deltak_xx
CollisionBranchList deltak_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes for Delta-Kaon (DeltaK) Scattering.
Definition: crosssections.cc:1754
smash::CrossSections::string_excitation
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 para...
Definition: crosssections.cc:2459
smash::CrossSections::CrossSections
CrossSections(const ParticleList &incoming_particles, double sqrt_s, const std::pair< FourVector, FourVector > potentials)
Construct CrossSections instance.
Definition: crosssections.cc:100
smash::CrossSections::xs_dn_dprimen
static double xs_dn_dprimen(double sqrts, double cm_mom, ParticleTypePtr produced_nucleus, const ParticleType &type_nucleus, const ParticleType &type_N)
Parametrized cross section for Nd → Nd', N̅d → N̅d', N̅d̅→ N̅d̅', Nd̅→ Nd̅' and reverse (e....
Definition: crosssections.cc:2391
smash::CrossSections::two_to_four
CollisionBranchList two_to_four() const
Find all 2->4 processes for the given scattering.
Definition: crosssections.cc:1008
smash::CrossSections::elastic
CollisionBranchPtr elastic(const ScatterActionsFinderParameters &finder_parameters) const
Determine the elastic cross section for this collision.
Definition: crosssections.cc:306
smash::CrossSections::formation
double formation(const ParticleType &type_resonance, double cm_momentum_sqr) const
Return the 2-to-1 resonance production cross section for a given resonance.
Definition: crosssections.cc:871
smash::CrossSections::add_channel
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.
Definition: crosssections.h:610
smash::CrossSections::string_hard_cross_section
double string_hard_cross_section() const
Determine the (parametrized) hard non-diffractive string cross section for this collision.
Definition: crosssections.cc:2657
smash::CrossSections::two_to_three_xs
static double two_to_three_xs(const ParticleType &type_in1, const ParticleType &type_in2, double sqrts)
Determine 2->3 cross section for the scattering of the given particle types.
Definition: crosssections.cc:1075
smash::CrossSections::nk_el
double nk_el() const
Determine the elastic cross section for a nucleon-kaon (NK) collision.
Definition: crosssections.cc:746
smash::CrossSections::nk_xx
CollisionBranchList nk_xx(const ReactionsBitSet &included_2to2, double KN_offset) const
Find all inelastic 2->2 background processes for Nucleon-Kaon (NK) Scattering.
Definition: crosssections.cc:1250
smash::CrossSections::two_to_one
CollisionBranchList two_to_one() const
Find all resonances that can be produced in a 2->1 collision of the two input particles and the produ...
Definition: crosssections.cc:841
smash::CrossSections::incoming_particles_
const ParticleList incoming_particles_
List with data of scattering particles.
Definition: crosssections.h:582
smash::CrossSections::npi_el
double npi_el() const
Determine the elastic cross section for a nucleon-pion (Npi) collision.
Definition: crosssections.cc:433
smash::CrossSections::nn_el
double nn_el() const
Determine the (parametrized) elastic cross section for a nucleon-nucleon (NN) collision.
Definition: crosssections.cc:402
smash::CrossSections::is_BBbar_pair_
const bool is_BBbar_pair_
Whether incoming particles are a pair of a baryon and an antibaryon (could be different baryon types)
Definition: crosssections.h:597
smash::CrossSections::dn_xx
CollisionBranchList dn_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes involving Nucleon and (anti-) Deuteron (dN), specifically Nd → Nd',...
Definition: crosssections.cc:2429
smash::CrossSections::elastic_parametrization
double elastic_parametrization(bool use_AQM, double pipi_offset) const
Choose the appropriate parametrizations for given incoming particles and return the (parametrized) el...
Definition: crosssections.cc:341
smash::CrossSections::two_to_four_xs
static double two_to_four_xs(const ParticleType &type_in1, const ParticleType &type_in2, double sqrts)
Determine 2->4 cross section for the scattering of the given particle types.
Definition: crosssections.cc:1114
smash::CrossSections::ypi_xx
CollisionBranchList ypi_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes for Hyperon-Pion (Ypi) Scattering.
Definition: crosssections.cc:1905
smash::CrossSections::parametrized_total
double parametrized_total(const ScatterActionsFinderParameters &finder_parameters) const
Select the parametrization for the total cross section, given the types of incoming particles.
Definition: crosssections.cc:205
smash::CrossSections::is_NNbar_pair_
const bool is_NNbar_pair_
Whether incoming particles are a nulecon-antinucleon pair (same isospin)
Definition: crosssections.h:600
smash::CrossSections::nn_to_resonance_matrix_element
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 ...
Definition: crosssections.cc:2809
smash::CrossSections::rare_two_to_two
CollisionBranchList rare_two_to_two() const
Find all 2->2 processes which are suppressed at high energies when strings are turned on with probabi...
Definition: crosssections.cc:328
smash::CrossSections::sum_xs_of
static double sum_xs_of(const CollisionBranchList &list)
Helper function: Sum all cross sections of the given process list.
Definition: crosssections.h:76
smash::CrossSections::two_to_three
CollisionBranchList two_to_three() const
Find all 2->3 processes for the given scattering.
Definition: crosssections.cc:953
smash::CrossSections::d_aN_inelastic_xs
static double d_aN_inelastic_xs(double aN_kinetic_energy)
Parametrization of deuteron-antinucleon inelastic cross section.
Definition: crosssections.cc:1071
smash::CrossSections::NNbar_annihilation
CollisionBranchPtr NNbar_annihilation(double current_xs, double scale_xs) const
Determine the cross section for NNbar annihilation, which is given by the difference between the para...
Definition: crosssections.cc:2698
smash::CrossSections::generate_collision_list
CollisionBranchList generate_collision_list(const ScatterActionsFinderParameters &finder_parameters, StringProcess *string_process) const
Generate a list of all possible collisions between the incoming particles with the given c....
Definition: crosssections.cc:115
smash::CrossSections::d_pi_inelastic_xs
static double d_pi_inelastic_xs(double pion_kinetic_energy)
Parametrization of deuteron-pion inelastic cross section.
Definition: crosssections.cc:1060
smash::CrossSections::high_energy
double high_energy(const StringTransitionParameters &transition_high_energy) const
Determine the parametrized total cross section at high energies for the given collision,...
Definition: crosssections.cc:2594
smash::CrossSections::nn_xx
CollisionBranchList nn_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes for Nucelon-Nucelon Scattering.
Definition: crosssections.cc:1178
smash::I_tot_range
Range of total isospin for reaction of particle a with particle b.
Definition: clebschgordan.h:68
smash::IsoParticleType::get_integral_RR
double get_integral_RR(IsoParticleType *type_res_2, double sqrts)
Look up the tabulated resonance integral for the XX -> RR cross section.
Definition: isoparticletype.cc:340
smash::IsoParticleType::get_integral_NR
double get_integral_NR(double sqrts)
Look up the tabulated resonance integral for the XX -> NR cross section.
Definition: isoparticletype.cc:308
smash::IsoParticleType::get_integral_RK
double get_integral_RK(double sqrts)
Look up the tabulated resonance integral for the XX -> RK cross section.
Definition: isoparticletype.cc:324
smash::IsoParticleType::get_integral_piR
double get_integral_piR(double sqrts)
Look up the tabulated resonance integral for the XX -> piR cross section.
Definition: isoparticletype.cc:316
smash::KaonNucleonRatios::get_ratio
double get_ratio(const ParticleType &a, const ParticleType &b, const ParticleType &c, const ParticleType &d) const
Return the isospin ratio of the given K N -> K Delta cross section.
Definition: parametrizations.cc:743
smash::ParticleData
ParticleData contains the dynamic information of a certain particle.
Definition: particledata.h:58
smash::ParticleData::pdgcode
PdgCode pdgcode() const
Get the pdgcode of the particle.
Definition: particledata.h:87
smash::ParticleData::type
const ParticleType & type() const
Get the type of the particle.
Definition: particledata.h:128
smash::ParticleData::is_baryon
bool is_baryon() const
Definition: particledata.h:94
smash::ParticleTypePtr
A pointer-like interface to global references to ParticleType objects.
Definition: particletype.h:676
smash::ParticleType
Particle type contains the static properties of a particle species.
Definition: particletype.h:98
smash::ParticleType::is_pion
bool is_pion() const
Definition: particletype.h:219
smash::ParticleType::is_baryon
bool is_baryon() const
Definition: particletype.h:204
smash::ParticleType::spectral_function
double spectral_function(double m) const
Full spectral function of the resonance (relativistic Breit-Wigner distribution with mass-dependent ...
Definition: particletype.cc:569
smash::ParticleType::is_Nstar
bool is_Nstar() const
Definition: particletype.h:231
smash::ParticleType::try_find
static const ParticleTypePtr try_find(PdgCode pdgcode)
Returns the ParticleTypePtr for the given pdgcode.
Definition: particletype.cc:89
smash::ParticleType::find
static const ParticleType & find(PdgCode pdgcode)
Returns the ParticleType object for the given pdgcode.
Definition: particletype.cc:99
smash::ParticleType::is_nucleus
bool is_nucleus() const
Definition: particletype.h:249
smash::ParticleType::pdgcode
PdgCode pdgcode() const
Definition: particletype.h:157
smash::ParticleType::name
const std::string & name() const
Definition: particletype.h:142
smash::ParticleType::charge
int32_t charge() const
The charge of the particle.
Definition: particletype.h:189
smash::ParticleType::antiparticle_sign
int antiparticle_sign() const
Definition: particletype.h:166
smash::ParticleType::list_nucleons
static ParticleTypePtrList & list_nucleons()
Definition: particletype.cc:69
smash::ParticleType::list_anti_nucleons
static ParticleTypePtrList & list_anti_nucleons()
Definition: particletype.cc:71
smash::ParticleType::is_Nstar1535
bool is_Nstar1535() const
Definition: particletype.h:237
smash::ParticleType::is_stable
bool is_stable() const
Definition: particletype.h:246
smash::ParticleType::is_dprime
bool is_dprime() const
Definition: particletype.h:258
smash::ParticleType::is_Delta
bool is_Delta() const
Definition: particletype.h:225
smash::ParticleType::width_at_pole
double width_at_pole() const
Definition: particletype.h:151
smash::ParticleType::list_anti_Deltas
static ParticleTypePtrList & list_anti_Deltas()
Definition: particletype.cc:77
smash::ParticleType::is_nucleon
bool is_nucleon() const
Definition: particletype.h:216
smash::ParticleType::mass
double mass() const
Definition: particletype.h:145
smash::ParticleType::list_baryon_resonances
static ParticleTypePtrList & list_baryon_resonances()
Definition: particletype.cc:81
smash::ParticleType::list_Deltas
static ParticleTypePtrList & list_Deltas()
Definition: particletype.cc:75
smash::ParticleType::is_deuteron
bool is_deuteron() const
Definition: particletype.h:252
smash::ParticleType::is_meson
bool is_meson() const
Definition: particletype.h:207
smash::ParticleType::get_partial_in_width
double get_partial_in_width(const double m, const ParticleData &p_a, const ParticleData &p_b) const
Get the mass-dependent partial in-width of a resonance with mass m, decaying into two given daughter ...
Definition: particletype.cc:550
smash::ParticleType::spin
unsigned int spin() const
Definition: particletype.h:192
smash::ParticleType::baryon_number
int baryon_number() const
Definition: particletype.h:210
smash::ParticleType::is_Deltastar
bool is_Deltastar() const
Definition: particletype.h:240
smash::ParticleType::iso_multiplet
IsoParticleType * iso_multiplet() const
Definition: particletype.h:186
smash::ParticleType::list_light_nuclei
static ParticleTypePtrList & list_light_nuclei()
Definition: particletype.cc:85
smash::PdgCode
PdgCode stores a Particle Data Group Particle Numbering Scheme particle type number.
Definition: pdgcode.h:111
smash::PdgCode::code
std::int32_t code() const
Definition: pdgcode.h:322
smash::PdgCode::antiparticle_sign
int antiparticle_sign() const
Definition: pdgcode.h:651
smash::PdgCode::nucleus_n
int nucleus_n() const
Number of neutrons in nucleus.
Definition: pdgcode.h:814
smash::PdgCode::is_meson
bool is_meson() const
Definition: pdgcode.h:394
smash::PdgCode::nucleus_ap
int nucleus_ap() const
Number of antiprotons in nucleus.
Definition: pdgcode.h:824
smash::PdgCode::nucleus_an
int nucleus_an() const
Number of antineutrons in nucleus.
Definition: pdgcode.h:828
smash::PdgCode::spin
unsigned int spin() const
Definition: pdgcode.h:608
smash::PdgCode::is_pion
bool is_pion() const
Definition: pdgcode.h:457
smash::PdgCode::is_kaon
bool is_kaon() const
Definition: pdgcode.h:451
smash::PdgCode::is_hyperon
bool is_hyperon() const
Definition: pdgcode.h:428
smash::PdgCode::is_nucleus
bool is_nucleus() const
Definition: pdgcode.h:364
smash::PdgCode::is_deuteron
bool is_deuteron() const
Definition: pdgcode.h:475
smash::PdgCode::is_antiparticle_of
bool is_antiparticle_of(const PdgCode rhs) const
Definition: pdgcode.h:746
smash::PdgCode::nucleus_p
int nucleus_p() const
Number of protons in nucleus.
Definition: pdgcode.h:810
smash::PdgCode::is_baryon
bool is_baryon() const
Definition: pdgcode.h:391
smash::PdgCode::isospin3
int isospin3() const
Definition: pdgcode.h:506
smash::PdgCode::is_nucleon
bool is_nucleon() const
Definition: pdgcode.h:397
smash::PdgCode::nucleus_La
int nucleus_La() const
Number of Lambdas in nucleus.
Definition: pdgcode.h:820
smash::PdgCode::nucleus_A
int nucleus_A() const
Nucleus mass number.
Definition: pdgcode.h:838
smash::PdgCode::nucleus_aLa
int nucleus_aLa() const
Number of anti-Lambdas in nucleus.
Definition: pdgcode.h:834
smash::PdgCode::is_Delta
bool is_Delta() const
Definition: pdgcode.h:421
smash::PdgCode::frac_strange
double frac_strange() const
Definition: pdgcode.h:519
smash::StringProcess
String excitation processes used in SMASH.
Definition: stringprocess.h:45
smash::StringProcess::pdg_map_for_pythia
static int pdg_map_for_pythia(PdgCode &pdg)
Take pdg code and map onto particle specie which can be handled by PYTHIA.
Definition: stringprocess.cc:3259
smash::StringProcess::cross_sections_diffractive
std::array< double, 3 > cross_sections_diffractive(int pdg_a, int pdg_b, double sqrt_s)
Interface to pythia_sigmatot_ to compute cross-sections of A+B-> different final states Schuler:1993w...
Definition: stringprocess.h:319
constants.h
Collection of useful constants that are known at compile time.
ReactionsBitSet
std::bitset< 10 > ReactionsBitSet
Container for the 2 to 2 reactions in the code.
Definition: forwarddeclarations.h:269
NNbarTreatment::TwoToFive
@ TwoToFive
Directly create 5 pions, use with multi-particle reactions.
NNbarTreatment::Resonances
@ Resonances
Use intermediate Resonances.
NNbarTreatment::Strings
@ Strings
Use string fragmentation.
A3_Nuclei_4to2
@ A3_Nuclei_4to2
Definition: forwarddeclarations.h:276
Deuteron_3to2
@ Deuteron_3to2
Definition: forwarddeclarations.h:274
KN_to_KDelta
@ KN_to_KDelta
Definition: forwarddeclarations.h:260
KN_to_KN
@ KN_to_KN
Definition: forwarddeclarations.h:259
NN_to_NR
@ NN_to_NR
Definition: forwarddeclarations.h:257
PiDeuteron_to_pidprime
@ PiDeuteron_to_pidprime
Definition: forwarddeclarations.h:264
NDeuteron_to_Ndprime
@ NDeuteron_to_Ndprime
Definition: forwarddeclarations.h:265
Strangeness_exchange
@ Strangeness_exchange
Definition: forwarddeclarations.h:261
PiDeuteron_to_NN
@ PiDeuteron_to_NN
Definition: forwarddeclarations.h:263
NN_to_DR
@ NN_to_DR
Definition: forwarddeclarations.h:258
smash::logg
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:39
smash::pdg::pi_p
constexpr int pi_p
π⁺.
Definition: pdgcode_constants.h:67
smash::pdg::Delta_p
constexpr int Delta_p
Δ⁺.
Definition: pdgcode_constants.h:45
smash::pdg::Delta_pp
constexpr int Delta_pp
Δ⁺⁺.
Definition: pdgcode_constants.h:43
smash::pdg::antideuteron
constexpr int64_t antideuteron
Anti-deuteron in decimal digits.
Definition: pdgcode_constants.h:104
smash::pdg::Sigma_m
constexpr int Sigma_m
Σ⁻.
Definition: pdgcode_constants.h:58
smash::pdg::K_p
constexpr int K_p
K⁺.
Definition: pdgcode_constants.h:74
smash::pdg::K_z
constexpr int K_z
K⁰.
Definition: pdgcode_constants.h:76
smash::pdg::p
constexpr int p
Proton.
Definition: pdgcode_constants.h:28
smash::pdg::h1
constexpr int h1
h₁(1170).
Definition: pdgcode_constants.h:95
smash::pdg::K_m
constexpr int K_m
K̄⁻.
Definition: pdgcode_constants.h:80
smash::pdg::pi_z
constexpr int pi_z
π⁰.
Definition: pdgcode_constants.h:69
smash::pdg::n
constexpr int n
Neutron.
Definition: pdgcode_constants.h:30
smash::pdg::deuteron
constexpr int64_t deuteron
Deuteron.
Definition: pdgcode_constants.h:102
smash::pdg::Delta_m
constexpr int Delta_m
Δ⁻.
Definition: pdgcode_constants.h:49
smash::pdg::Delta_z
constexpr int Delta_z
Δ⁰.
Definition: pdgcode_constants.h:47
smash::pdg::rho_z
constexpr int rho_z
ρ⁰.
Definition: pdgcode_constants.h:90
smash::pdg::Lambda
constexpr int Lambda
Λ.
Definition: pdgcode_constants.h:52
smash::pdg::pi_m
constexpr int pi_m
π⁻.
Definition: pdgcode_constants.h:71
smash::pdg::Sigma_p
constexpr int Sigma_p
Σ⁺.
Definition: pdgcode_constants.h:54
smash::pdg::Kbar_z
constexpr int Kbar_z
K̄⁰.
Definition: pdgcode_constants.h:78
smash::pdg::Sigma_z
constexpr int Sigma_z
Σ⁰.
Definition: pdgcode_constants.h:56
smash
Definition: action.h:24
smash::kplusn_k0p
double kplusn_k0p(double mandelstam_s)
K+ n charge exchange cross section parametrization.
Definition: parametrizations.cc:552
smash::kplusp_total
double kplusp_total(double mandelstam_s)
K+ p total cross section parametrization.
Definition: parametrizations.cc:486
smash::pizeropizero_total
double pizeropizero_total(double sqrts)
pi0 pi0 total cross section parametrized from PDG2018, smoothed using the LOWESS algorithm.
Definition: parametrizations.cc:128
smash::kminusp_pi0lambda
double kminusp_pi0lambda(double sqrts)
K- p <-> pi0 Lambda cross section parametrization Fit to Landolt-Börnstein instead of UrQMD values.
Definition: parametrizations.cc:796
smash::pCM
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
smash::pipluspiminus_total
double pipluspiminus_total(double sqrts)
pi+ pi- total cross section parametrized from PDG2018, smoothed using the LOWESS algorithm.
Definition: parametrizations.cc:113
smash::list_possible_resonances
ParticleTypePtrList list_possible_resonances(const ParticleTypePtr type_a, const ParticleTypePtr type_b)
Lists the possible resonances that decay into two particles.
Definition: particletype.cc:775
smash::piminusp_sigma0k0_res
double piminusp_sigma0k0_res(double mandelstam_s)
pi- p -> Sigma0 K0 cross section parametrization, resonance contribution.
Definition: parametrizations.cc:351
smash::pCM_sqr
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:91
smash::ppbar_total
double ppbar_total(double mandelstam_s)
ppbar total cross section parametrization Source: Bass:1998ca
Definition: parametrizations.cc:458
smash::np_total
double np_total(double mandelstam_s)
np total cross section parametrization Sources: low-p: Cugnon:1996kh highest-p: Buss:2011mx
Definition: parametrizations.cc:425
smash::piminusp_elastic
double piminusp_elastic(double mandelstam_s)
pi-p elastic cross section parametrization Source: GiBUU:parametrizationBarMes_HighEnergy....
Definition: parametrizations.cc:273
smash::npbar_high_energy
double npbar_high_energy(double mandelstam_s)
npbar total cross section at high energies
Definition: parametrizations.cc:66
smash::detailed_balance_factor_RR
static double detailed_balance_factor_RR(double sqrts, double pcm, const ParticleType &a, const ParticleType &b, const ParticleType &c, const ParticleType &d)
Helper function: Calculate the detailed balance factor R such that.
Definition: crosssections.cc:70
smash::kminusn_piminussigma0
double kminusn_piminussigma0(double sqrts)
K- n <-> pi- Sigma0 cross section parametrization Follow from the parametrization with the same stran...
Definition: parametrizations.cc:800
smash::kbar0p_elastic_background
double kbar0p_elastic_background(double mandelstam_s)
Kbar0 p elastic background cross section parametrization Source: Buss:2011mx , B.3....
Definition: parametrizations.cc:629
smash::kaon_nucleon_ratios
KaonNucleonRatios kaon_nucleon_ratios
Definition: parametrizations.cc:769
smash::ppbar_elastic
double ppbar_elastic(double mandelstam_s)
ppbar elastic cross section parametrization Source: Bass:1998ca
Definition: parametrizations.cc:441
smash::kminusp_elastic_background
double kminusp_elastic_background(double mandelstam_s)
K- p elastic background cross section parametrization Source: Buss:2011mx , B.3.9.
Definition: parametrizations.cc:573
smash::np_high_energy
double np_high_energy(double mandelstam_s)
np total cross section at high energies
Definition: parametrizations.cc:62
smash::pp_elastic_high_energy
double pp_elastic_high_energy(double mandelstam_s, double m1, double m2)
pp elastic cross section parametrization, with only the high energy part generalized to all energy re...
Definition: parametrizations.cc:382
smash::Npi_string_hard
double Npi_string_hard(double mandelstam_s)
nucleon-pion hard scattering cross section (with partonic scattering)
Definition: parametrizations.cc:105
smash::kminusn_piminuslambda
double kminusn_piminuslambda(double sqrts)
K- n <-> pi- Lambda cross section parametrization Follow from the parametrization with the same stran...
Definition: parametrizations.cc:805
smash::isospin_clebsch_gordan_sqr_2to2
double isospin_clebsch_gordan_sqr_2to2(const ParticleType &p_a, const ParticleType &p_b, const ParticleType &p_c, const ParticleType &p_d, const int I=-1)
Calculate the squared isospin Clebsch-Gordan coefficient for a 2-to-2 reaction A + B -> C + D.
Definition: clebschgordan.cc:71
smash::ProcessType::TwoToOne
@ TwoToOne
See here for a short description.
smash::ProcessType::StringSoftDoubleDiffractive
@ StringSoftDoubleDiffractive
See here for a short description.
smash::ProcessType::TwoToFive
@ TwoToFive
See here for a short description.
smash::ProcessType::StringSoftSingleDiffractiveXB
@ StringSoftSingleDiffractiveXB
See here for a short description.
smash::ProcessType::TwoToTwo
@ TwoToTwo
See here for a short description.
smash::ProcessType::Elastic
@ Elastic
See here for a short description.
smash::ProcessType::TwoToFour
@ TwoToFour
See here for a short description.
smash::ProcessType::StringSoftAnnihilation
@ StringSoftAnnihilation
See here for a short description.
smash::ProcessType::StringSoftNonDiffractive
@ StringSoftNonDiffractive
See here for a short description.
smash::ProcessType::StringSoftSingleDiffractiveAX
@ StringSoftSingleDiffractiveAX
See here for a short description.
smash::ProcessType::StringHard
@ StringHard
See here for a short description.
smash::ProcessType::TwoToThree
@ TwoToThree
See here for a short description.
smash::minimum_sqrts_pythia_can_handle
constexpr double minimum_sqrts_pythia_can_handle
Energy in GeV, below which hard reactions via pythia are impossible.
Definition: constants.h:111
smash::ppbar_high_energy
double ppbar_high_energy(double mandelstam_s)
ppbar total cross section at high energies
Definition: parametrizations.cc:58
smash::pp_high_energy
double pp_high_energy(double mandelstam_s)
pp total cross section at high energies
Definition: parametrizations.cc:54
smash::pCM_sqr_from_s
T pCM_sqr_from_s(const T s, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:52
smash::pipi_string_hard
double pipi_string_hard(double mandelstam_s)
pion-pion hard scattering cross section (with partonic scattering)
Definition: parametrizations.cc:109
smash::piplusp_high_energy
double piplusp_high_energy(double mandelstam_s)
pi+p total cross section at high energies
Definition: parametrizations.cc:70
smash::kminusp_piminussigmaplus
double kminusp_piminussigmaplus(double sqrts)
K- p <-> pi- Sigma+ cross section parametrization Taken from UrQMD (Graef:2014mra ).
Definition: parametrizations.cc:784
smash::piplusp_elastic_high_energy
double piplusp_elastic_high_energy(double mandelstam_s, double m1, double m2)
pi+p elactic cross section parametrization.
Definition: parametrizations.cc:175
smash::piplusp_sigmapluskplus_pdg
double piplusp_sigmapluskplus_pdg(double mandelstam_s)
pi+ p to Sigma+ K+ cross section parametrization, PDG data.
Definition: parametrizations.cc:226
smash::LCrossSections
static constexpr int LCrossSections
Definition: crosssections.cc:19
smash::piminusp_total
double piminusp_total(double sqrts)
pi- p total cross section parametrized from PDG2018, smoothed using the LOWESS algorithm.
Definition: parametrizations.cc:241
smash::deuteron_mass
constexpr double deuteron_mass
Deuteron mass in GeV.
Definition: constants.h:92
smash::deuteron_nucleon_elastic
double deuteron_nucleon_elastic(double mandelstam_s)
Deuteron nucleon elastic cross-section [mb] parametrized by Oh:2009gx .
Definition: parametrizations.cc:480
smash::nucleon_mass
constexpr double nucleon_mass
Nucleon mass in GeV.
Definition: constants.h:58
smash::pow_int
constexpr T pow_int(const T base, unsigned const exponent)
Efficient template for calculating integer powers using squaring.
Definition: pow.h:23
smash::piminusp_sigmaminuskplus_pdg
double piminusp_sigmaminuskplus_pdg(double mandelstam_s)
pi- p -> Sigma- K+ cross section parametrization, PDG data.
Definition: parametrizations.cc:334
smash::piminusp_lambdak0_pdg
double piminusp_lambdak0_pdg(double mandelstam_s)
pi- p -> Lambda K0 cross section parametrization, PDG data.
Definition: parametrizations.cc:317
smash::append_list
static void append_list(CollisionBranchList &main_list, CollisionBranchList in_list, double weight=1.)
Helper function: Append a list of processes to another (main) list of processes.
Definition: crosssections.cc:91
smash::detailed_balance_factor_stable
static double detailed_balance_factor_stable(double s, const ParticleType &a, const ParticleType &b, const ParticleType &c, const ParticleType &d)
Helper function: Calculate the detailed balance factor R such that.
Definition: crosssections.cc:28
smash::k0p_elastic_background
double k0p_elastic_background(double mandelstam_s)
K0 p elastic background cross section parametrization Source: Buss:2011mx , B.3.9.
Definition: parametrizations.cc:619
smash::pack
constexpr uint64_t pack(int32_t x, int32_t y)
Pack two int32_t into an uint64_t.
Definition: pdgcode_constants.h:130
smash::deuteron_pion_elastic
double deuteron_pion_elastic(double mandelstam_s)
Deuteron pion elastic cross-section [mb] parametrized to fit pi-d elastic scattering data (the data c...
Definition: parametrizations.cc:475
smash::NN_string_hard
double NN_string_hard(double mandelstam_s)
nucleon-nucleon hard scattering cross section (with partonic scattering)
Definition: parametrizations.cc:101
smash::xs_ppbar_annihilation
double xs_ppbar_annihilation(double mandelstam_s)
parametrized cross-section for proton-antiproton annihilation used in the UrQMD model
Definition: parametrizations.cc:78
smash::kplusp_inelastic_background
double kplusp_inelastic_background(double mandelstam_s)
K+ p inelastic background cross section parametrization Source: Buss:2011mx , B.3....
Definition: parametrizations.cc:639
smash::pion_mass
constexpr double pion_mass
Pion mass in GeV.
Definition: constants.h:65
smash::hbarc
constexpr double hbarc
GeV <-> fm conversion factor.
Definition: constants.h:25
smash::kminusp_pi0sigma0
double kminusp_pi0sigma0(double sqrts)
K- p <-> pi0 Sigma0 cross section parametrization Fit to Landolt-Börnstein instead of UrQMD values.
Definition: parametrizations.cc:792
smash::kplusn_elastic_background
double kplusn_elastic_background(double mandelstam_s)
K+ n elastic background cross section parametrization sigma(K+n->K+n) = sigma(K+n->K0p) = 0....
Definition: parametrizations.cc:548
smash::pp_total
double pp_total(double mandelstam_s)
pp total cross section parametrization Sources: low-p: Cugnon:1996kh highest-p: Buss:2011mx
Definition: parametrizations.cc:390
smash::pCM_from_s
T pCM_from_s(const T s, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:66
smash::kplusn_total
double kplusn_total(double mandelstam_s)
K+ n total cross section parametrization.
Definition: parametrizations.cc:498
smash::really_small
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:37
smash::np_elastic
double np_elastic(double mandelstam_s)
np elastic cross section parametrization Source: Weil:2013mya , eq.
Definition: parametrizations.cc:406
smash::detailed_balance_factor_RK
static double detailed_balance_factor_RK(double sqrts, double pcm, const ParticleType &a, const ParticleType &b, const ParticleType &c, const ParticleType &d)
Helper function: Calculate the detailed balance factor R such that.
Definition: crosssections.cc:47
smash::k0n_elastic_background
double k0n_elastic_background(double mandelstam_s)
K0 n elastic background cross section parametrization Source: Buss:2011mx , B.3.9.
Definition: parametrizations.cc:624
smash::LScatterAction
static constexpr int LScatterAction
Definition: bremsstrahlungaction.cc:17
smash::kminusp_piplussigmaminus
double kminusp_piplussigmaminus(double sqrts)
K- p <-> pi+ Sigma- cross section parametrization Taken from UrQMD (Graef:2014mra ).
Definition: parametrizations.cc:788
smash::kminusp_total
double kminusp_total(double mandelstam_s)
K- p total cross section parametrization.
Definition: parametrizations.cc:510
smash::kbar0n_elastic_background
double kbar0n_elastic_background(double mandelstam_s)
Kbar0 n elastic background cross section parametrization Source: Buss:2011mx , B.3....
Definition: parametrizations.cc:634
smash::kminusp_kbar0n
double kminusp_kbar0n(double mandelstam_s)
K- p <-> Kbar0 n cross section parametrization.
Definition: parametrizations.cc:771
smash::piplusp_total
double piplusp_total(double sqrts)
pi+ p total cross section parametrized from PDG2018, smoothed using the LOWESS algorithm.
Definition: parametrizations.cc:143
smash::piminusp_high_energy
double piminusp_high_energy(double mandelstam_s)
pi-p total cross section at high energies
Definition: parametrizations.cc:74
smash::kminusn_total
double kminusn_total(double mandelstam_s)
K- n total cross section parametrization.
Definition: parametrizations.cc:523
smash::kminusn_elastic_background
double kminusn_elastic_background(double mandelstam_s)
K- n elastic background cross section parametrization Source: Buss:2011mx , B.3.9.
Definition: parametrizations.cc:617
smash::pp_elastic
double pp_elastic(double mandelstam_s)
pp elastic cross section parametrization Source: Weil:2013mya , eq.
Definition: parametrizations.cc:363
smash::kplusn_inelastic_background
double kplusn_inelastic_background(double mandelstam_s)
K+ n inelastic background cross section parametrization Source: Buss:2011mx , B.3....
Definition: parametrizations.cc:652
smash::kplusp_elastic_background
double kplusp_elastic_background(double mandelstam_s)
K+ p elastic background cross section parametrization.
Definition: parametrizations.cc:535
smash::piplusp_elastic_AQM
double piplusp_elastic_AQM(double mandelstam_s, double m1, double m2)
An overload of piplusp_elastic_high_energy in which the very low part is replaced by a flat 5 mb cros...
Definition: parametrizations.cc:182
smash::piplusp_elastic
double piplusp_elastic(double mandelstam_s)
pi+p elastic cross section parametrization, PDG data.
Definition: parametrizations.cc:192
smash::fm2_mb
constexpr double fm2_mb
mb <-> fm^2 conversion factor.
Definition: constants.h:28
smash::ScatterActionsFinderParameters
Helper structure for ScatterActionsFinder.
Definition: scatteractionsfinderparameters.h:65
smash::ScatterActionsFinderParameters::strings_with_probability
const bool strings_with_probability
This indicates whether the string fragmentation is swiched on with a probability smoothly increasing ...
Definition: scatteractionsfinderparameters.h:110
smash::ScatterActionsFinderParameters::included_2to2
const ReactionsBitSet included_2to2
List of included 2<->2 reactions.
Definition: scatteractionsfinderparameters.h:89
smash::ScatterActionsFinderParameters::use_AQM
const bool use_AQM
Switch to control whether to use AQM or not.
Definition: scatteractionsfinderparameters.h:102
smash::ScatterActionsFinderParameters::scale_xs
const double scale_xs
Factor by which all (partial) cross sections are scaled.
Definition: scatteractionsfinderparameters.h:74
smash::ScatterActionsFinderParameters::elastic_parameter
const double elastic_parameter
Elastic cross section parameter (in mb).
Definition: scatteractionsfinderparameters.h:67
smash::ScatterActionsFinderParameters::strings_switch
const bool strings_switch
Indicates whether string fragmentation is switched on.
Definition: scatteractionsfinderparameters.h:100
smash::ScatterActionsFinderParameters::transition_high_energy
const StringTransitionParameters transition_high_energy
Constants related to transition between low collision energies - mediated via resonances - and high c...
Definition: scatteractionsfinderparameters.h:121
smash::ScatterActionsFinderParameters::additional_el_xs
const double additional_el_xs
Additional constant contribution (in mb) to the elastic cross sections.
Definition: scatteractionsfinderparameters.h:81
smash::ScatterActionsFinderParameters::nnbar_treatment
const NNbarTreatment nnbar_treatment
Switch for NNbar reactions.
Definition: scatteractionsfinderparameters.h:87
smash::ScatterActionsFinderParameters::low_snn_cut
const double low_snn_cut
Elastic collsions between two nucleons with sqrt_s below low_snn_cut_ are excluded.
Definition: scatteractionsfinderparameters.h:72
smash::ScatterActionsFinderParameters::two_to_one
const bool two_to_one
Enables resonance production.
Definition: scatteractionsfinderparameters.h:95
smash::StringTransitionParameters
Constants related to transition between low and high collision energies.
Definition: scatteractionsfinderparameters.h:19
smash::StringTransitionParameters::sqrts_range_Npi
const std::pair< double, double > sqrts_range_Npi
Transition range in N collisions.
Definition: scatteractionsfinderparameters.h:21
smash::StringTransitionParameters::pipi_offset
const double pipi_offset
Constant offset as to where to turn on the strings and elastic processes for reactions (this is an e...
Definition: scatteractionsfinderparameters.h:49
smash::StringTransitionParameters::sqrts_add_lower
const double sqrts_add_lower
Constant for the lower end of transition region in the case of AQM this is added to the sum of masses...
Definition: scatteractionsfinderparameters.h:35
smash::StringTransitionParameters::KN_offset
const double KN_offset
Constant offset as to where to shift from 2to2 to string processes (in GeV) in the case of KN reactio...
Definition: scatteractionsfinderparameters.h:55
smash::StringTransitionParameters::sqrts_range_NN
const std::pair< double, double > sqrts_range_NN
Transition range in NN collisions.
Definition: scatteractionsfinderparameters.h:29
smash::StringTransitionParameters::sqrts_range_width
const double sqrts_range_width
Constant for the range of transition region, in the case of AQM this is added to the sum of masses + ...
Definition: scatteractionsfinderparameters.h:41