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