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