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