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