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