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