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