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