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