Version: SMASH-3.4
crosssections.cc
Go to the documentation of this file.
1 /*
2  *
3  * Copyright (c) 2018-2026
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"
15 #include "smash/logging.h"
16 #include "smash/parametrizations.h"
17 #include "smash/pow.h"
18 
19 namespace smash {
20 static constexpr int LCrossSections = LogArea::CrossSections::id;
21 static constexpr int LScatterAction = LogArea::ScatterAction::id;
22 
23 /**
24  * Helper function:
25  * Calculate the detailed balance factor R such that
26  * \f[ R = \sigma(AB \to CD) / \sigma(CD \to AB) \f]
27  * where \f$ A, B, C, D \f$ are stable.
28  */
29 static double detailed_balance_factor_stable(double s, const ParticleType& a,
30  const ParticleType& b,
31  const ParticleType& c,
32  const ParticleType& d) {
33  double spin_factor = (c.spin() + 1) * (d.spin() + 1);
34  spin_factor /= (a.spin() + 1) * (b.spin() + 1);
35  double symmetry_factor = (1 + (a == b));
36  symmetry_factor /= (1 + (c == d));
37  const double momentum_factor = pCM_sqr_from_s(s, c.mass(), d.mass()) /
38  pCM_sqr_from_s(s, a.mass(), b.mass());
39  return spin_factor * symmetry_factor * momentum_factor;
40 }
41 
42 /**
43  * Helper function:
44  * Calculate the detailed balance factor R such that
45  * \f[ R = \sigma(AB \to CD) / \sigma(CD \to AB) \f]
46  * where \f$A\f$ is unstable, \f$B\f$ is a kaon and \f$C, D\f$ are stable.
47  */
48 static double detailed_balance_factor_RK(double sqrts, double pcm,
49  const ParticleType& a,
50  const ParticleType& b,
51  const ParticleType& c,
52  const ParticleType& d) {
53  assert(!a.is_stable());
54  assert(b.pdgcode().is_kaon());
55  double spin_factor = (c.spin() + 1) * (d.spin() + 1);
56  spin_factor /= (a.spin() + 1) * (b.spin() + 1);
57  double symmetry_factor = (1 + (a == b));
58  symmetry_factor /= (1 + (c == d));
59  const double momentum_factor =
60  pCM_sqr(sqrts, c.mass(), d.mass()) /
61  (pcm * a.iso_multiplet()->get_integral_RK(sqrts));
62  return spin_factor * symmetry_factor * momentum_factor;
63 }
64 
65 /**
66  * Helper function:
67  * Calculate the detailed balance factor R such that
68  * \f[ R = \sigma(AB \to CD) / \sigma(CD \to AB) \f]
69  * where \f$A\f$ and \f$B\f$ are unstable, and \f$C\f$ and \f$D\f$ are stable.
70  */
71 static double detailed_balance_factor_RR(double sqrts, double pcm,
72  const ParticleType& a,
73  const ParticleType& b,
74  const ParticleType& c,
75  const ParticleType& d) {
76  assert(!a.is_stable());
77  assert(!b.is_stable());
78  double spin_factor = (c.spin() + 1) * (d.spin() + 1);
79  spin_factor /= (a.spin() + 1) * (b.spin() + 1);
80  double symmetry_factor = (1 + (a == b));
81  symmetry_factor /= (1 + (c == d));
82  const double momentum_factor =
83  pCM_sqr(sqrts, c.mass(), d.mass()) /
84  (pcm * a.iso_multiplet()->get_integral_RR(b.iso_multiplet(), sqrts));
85  return spin_factor * symmetry_factor * momentum_factor;
86 }
87 
88 /**
89  * Helper function:
90  * Append a list of processes to another (main) list of processes.
91  */
92 static void append_list(CollisionBranchList& main_list,
93  CollisionBranchList in_list, double weight = 1.) {
94  main_list.reserve(main_list.size() + in_list.size());
95  for (auto& proc : in_list) {
96  proc->set_weight(proc->weight() * weight);
97  main_list.emplace_back(std::move(proc));
98  }
99 }
100 
101 /**
102  * Helper function:
103  * Throw if cross section between two particles is not implemented.
104  *
105  * \param[in] data_a incoming particle a
106  * \param[in] data_b incoming particle b
107  * \param[in] func_name name of the function that encountered the throw
108  */
109 [[noreturn]] static void throw_xsec_is_not_implemented(
110  const ParticleData& data_a, const ParticleData& data_b,
111  const std::string func_name) {
112  std::stringstream ss{};
113  const ParticleType& a = data_a.type();
114  const ParticleType& b = data_b.type();
115  ss << "Cross section for scattering of " << a.name() << " and " << b.name()
116  << " is not implemented in function 'CrossSections::" << func_name << "'.";
117  throw std::runtime_error(ss.str());
118 }
119 
120 /**
121  * Helper function:
122  * Throw if cross section is negative.
123  *
124  * \param[in] sqrts center of mass energy of incoming particles
125  * \param[in] xsec cross section
126  * \param[in] data_a incoming particle a
127  * \param[in] data_b incoming particle b
128  * \param[in] func_name name of the function that encountered the throw
129  */
130 [[noreturn]] static void throw_xsec_is_negative(const double sqrts,
131  const double xsec,
132  const ParticleData& data_a,
133  const ParticleData& data_b,
134  std::string func_name) {
135  const ParticleType& a = data_a.type();
136  const ParticleType& b = data_b.type();
137  const PdgCode& pdg_a = a.pdgcode();
138  const PdgCode& pdg_b = b.pdgcode();
139  std::stringstream ss{};
140  ss << "Negative cross section encountered in function 'CrossSections::"
141  << func_name << "':\na=" << a.name() << " b=" << b.name()
142  << " j_a=" << pdg_a.spin() << " j_b=" << pdg_b.spin() << " sigma=" << xsec
143  << " s=" << sqrts * sqrts << " sqrt(s)=" << sqrts;
144  throw std::runtime_error(ss.str());
145 }
146 
147 /**
148  * Helper function:
149  * Shift the energy of a collision for AQM rescaled cross sections.
150  *
151  * \param[in] mandelstam_s the rest frame total energy squared
152  * \param[in] m1 effective mass of incoming first particle
153  * \param[in] m2 effective mass of incoming second particle
154  * \param[in] m1_ref mass of the first AQM reference
155  * \param[in] m2_ref mass of the second AQM reference
156  * \return the shifted center of mass energy squared
157  */
158 static double effective_AQM_s(const double mandelstam_s, const double m1,
159  const double m2, const double m1_ref,
160  const double m2_ref) {
161  const double eff_sqrt_s = std::sqrt(mandelstam_s) - m1 - m2 + m1_ref + m2_ref;
162  return eff_sqrt_s * eff_sqrt_s;
163 }
164 
165 /**
166  * Helper function:
167  * Approximate cross section using AQM based on function `piminusp_high_energy`.
168  *
169  * \param[in] sqrts center of mass energy of incoming particles
170  * \param[in] pdg_a PDG code of incoming particle a
171  * \param[in] pdg_b PDG code of incoming particle b
172  * \param[in] AQM_scaling_factor_a AQM scaling factor of incoming particle a
173  * \param[in] AQM_scaling_factor_b AQM scaling factor of incoming particle b
174  * \return the approximated cross section
175  */
177  const double sqrts, const PdgCode& pdg_a, const PdgCode& pdg_b,
178  const double AQM_scaling_factor_a, const double AQM_scaling_factor_b) {
179  int n_mesons = 0;
180  if (pdg_a.is_meson()) {
181  n_mesons += 1;
182  }
183  if (pdg_b.is_meson()) {
184  n_mesons += 1;
185  }
186  return std::pow(2. / 3., n_mesons - 1) * piminusp_high_energy(sqrts * sqrts) *
187  AQM_scaling_factor_a * AQM_scaling_factor_b;
188 }
189 
190 /**
191  * Helper function:
192  * Print a warning message if `Charm_Rescattering_Method` is not set to `none`
193  * and AQM is disabled but should be used.
194  *
195  * \param[in] sqrts center of mass energy of incoming particles
196  * \param[in] type_a type of incoming particle a
197  * \param[in] type_b type of incoming particle b
198  * \param[in] charm_rescattering type of charm rescattering
199  */
201  const double sqrts, const ParticleType& type_a, const ParticleType& type_b,
202  const CharmRescattering charm_rescattering) {
203  std::ostringstream warn_msg{
204  "AQM is disabled and 'Charm_Rescattering_Method' is set to ",
205  std::ios::ate};
206  if (charm_rescattering == CharmRescattering::T_Matrix) {
207  warn_msg << "'T-matrix' with sqrt(s) = " << sqrts
208  << " GeV out of bounds of the underlying data";
209  } else if (charm_rescattering == CharmRescattering::Resonances) {
210  warn_msg << "'resonances'";
211  }
212  warn_msg << ".\nElastic interactions of " << type_a.name() << " and "
213  << type_b.name()
214  << " are disabled under these circumstances.\nPlease enable AQM "
215  "if these interactions should occur.";
216  logg[LCrossSections].warn(warn_msg.str());
217 }
218 
219 CrossSections::CrossSections(const ParticleList& incoming_particles,
220  const double sqrt_s,
221  const std::pair<FourVector, FourVector> potentials)
222  : incoming_particles_(incoming_particles),
223  sqrt_s_(sqrt_s),
224  potentials_(potentials),
225  is_BBbar_pair_(incoming_particles_[0].type().is_baryon() &&
226  incoming_particles_[1].type().is_baryon() &&
227  incoming_particles_[0].type().antiparticle_sign() ==
228  -incoming_particles_[1].type().antiparticle_sign()),
229  is_NNbar_pair_(
230  incoming_particles_[0].type().is_nucleon() &&
231  incoming_particles_[1].pdgcode() ==
232  incoming_particles_[0].type().get_antiparticle()->pdgcode()) {}
233 
235  const ScatterActionsFinderParameters& finder_parameters,
236  StringProcess* string_process) const {
237  CollisionBranchList process_list;
238  const ParticleType& t1 = incoming_particles_[0].type();
239  const ParticleType& t2 = incoming_particles_[1].type();
240 
241  double p_pythia = 0.;
242  if (finder_parameters.strings_with_probability) {
243  p_pythia = string_probability(finder_parameters);
244  }
245 
246  /* Elastic collisions between two nucleons with sqrt_s below
247  * low_snn_cut can not happen. */
248  const bool reject_by_nucleon_elastic_cutoff =
249  t1.is_nucleon() && t2.is_nucleon() &&
250  t1.antiparticle_sign() == t2.antiparticle_sign() &&
251  sqrt_s_ < finder_parameters.low_snn_cut;
252  bool incl_elastic =
253  finder_parameters.included_2to2[IncludedReactions::Elastic];
254  if (incl_elastic && !reject_by_nucleon_elastic_cutoff) {
255  process_list.emplace_back(elastic(finder_parameters));
256  }
257  if (incoming_particles_[0].is_core() != incoming_particles_[1].is_core()) {
258  return process_list;
259  }
260  if (p_pythia > 0.) {
261  /* String-excitation cross section =
262  * Parametrized total cross - the contributions
263  * from all other present channels. */
264  const double sig_current = sum_xs_of(process_list);
265  const double sig_string = std::max(
266  0., finder_parameters.scale_xs * high_energy(finder_parameters) -
267  sig_current);
268  append_list(
269  process_list,
270  string_excitation(sig_string, string_process, finder_parameters),
271  p_pythia);
272  append_list(process_list, rare_two_to_two(),
273  p_pythia * finder_parameters.scale_xs);
274  }
275  if (p_pythia < 1.) {
276  if (finder_parameters.two_to_one) {
277  // resonance formation (2->1)
278  append_list(process_list,
279  two_to_one(finder_parameters.charm_rescattering),
280  (1. - p_pythia) * finder_parameters.scale_xs);
281  }
282  if (finder_parameters.included_2to2.any()) {
283  // 2->2 (inelastic)
284  append_list(process_list,
285  two_to_two(finder_parameters.included_2to2,
286  finder_parameters.transition_high_energy.KN_offset,
287  finder_parameters.charm_rescattering),
288  (1. - p_pythia) * finder_parameters.scale_xs);
289  }
290  if (finder_parameters
292  1) {
293  // 2->3 (deuterons only 2-to-3 reaction at the moment)
294  append_list(process_list, two_to_three(),
295  (1. - p_pythia) * finder_parameters.scale_xs);
296  }
297  if (finder_parameters
299  1) {
300  // 2->4
301  append_list(process_list, two_to_four(),
302  (1. - p_pythia) * finder_parameters.scale_xs);
303  }
304  }
305  if (finder_parameters.nnbar_treatment == NNbarTreatment::TwoToFive &&
306  is_NNbar_pair_) {
307  // NNbar directly to 5 pions (2-to-5)
308  process_list.emplace_back(NNbar_to_5pi(finder_parameters.scale_xs));
309  }
310 
311  /* NNbar annihilation thru NNbar → ρh₁(1170); combined with the decays
312  * ρ → ππ and h₁(1170) → πρ, this gives a final state of 5 pions.
313  * Only use in cases when detailed balance MUST happen, i.e. in a box! */
314  if (finder_parameters.nnbar_treatment == NNbarTreatment::Resonances) {
315  if (is_NNbar_pair_) {
316  /* Has to be called after the other processes are already determined,
317  * so that the sum of the cross sections includes all other processes. */
318  process_list.emplace_back(NNbar_annihilation(sum_xs_of(process_list),
319  finder_parameters.scale_xs));
320  } else {
321  append_list(process_list, NNbar_creation(), finder_parameters.scale_xs);
322  }
323  }
324  return process_list;
325 }
326 
328  const ScatterActionsFinderParameters& finder_parameters) const {
329  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
330  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
331  double total_xs = 0.;
332  if (pdg_a.is_baryon() && pdg_b.is_baryon() &&
333  sqrt_s_ > finder_parameters.low_snn_cut) {
334  if (pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
335  // NN
336  total_xs = (pdg_a == pdg_b) ? pp_total(sqrt_s_ * sqrt_s_)
337  : np_total(sqrt_s_ * sqrt_s_);
338  } else {
339  // NNbar
340  total_xs = ppbar_total(sqrt_s_ * sqrt_s_);
341  }
342  total_xs *= finder_parameters.AQM_scaling_factor(pdg_a) *
343  finder_parameters.AQM_scaling_factor(pdg_b);
344  } else if ((pdg_a.is_baryon() && pdg_b.is_meson()) ||
345  (pdg_a.is_meson() && pdg_b.is_baryon())) {
346  const PdgCode& meson = pdg_a.is_meson() ? pdg_a : pdg_b;
347  const PdgCode& baryon = pdg_a.is_meson() ? pdg_b : pdg_a;
348  if (meson.is_kaon() && baryon.is_nucleon()) {
349  if ((meson.code() == pdg::K_p && baryon.code() == pdg::p) ||
350  (meson.code() == pdg::K_z && baryon.code() == pdg::n) ||
351  (meson.code() == pdg::K_m && baryon.code() == -pdg::p) ||
352  (meson.code() == pdg::Kbar_z && baryon.code() == -pdg::n)) {
353  // K⁺p, K⁰n, and anti-processes
354  total_xs = kplusp_total(sqrt_s_ * sqrt_s_);
355  } else if ((meson.code() == pdg::K_p && baryon.code() == -pdg::p) ||
356  (meson.code() == pdg::K_z && baryon.code() == -pdg::n) ||
357  (meson.code() == pdg::K_m && baryon.code() == pdg::p) ||
358  (meson.code() == pdg::Kbar_z && baryon.code() == pdg::n)) {
359  // K⁻p, K̅⁰n, and anti-processes
360  total_xs = kminusp_total(sqrt_s_ * sqrt_s_);
361  } else if ((meson.code() == pdg::K_p && baryon.code() == pdg::n) ||
362  (meson.code() == pdg::K_z && baryon.code() == pdg::p) ||
363  (meson.code() == pdg::K_m && baryon.code() == -pdg::n) ||
364  (meson.code() == pdg::Kbar_z && baryon.code() == -pdg::p)) {
365  // K⁺n, K⁰p, and anti-processes
366  total_xs = kplusn_total(sqrt_s_ * sqrt_s_);
367  } else if ((meson.code() == pdg::K_p && baryon.code() == -pdg::n) ||
368  (meson.code() == pdg::K_z && baryon.code() == -pdg::p) ||
369  (meson.code() == pdg::K_m && baryon.code() == pdg::n) ||
370  (meson.code() == pdg::Kbar_z && baryon.code() == pdg::p)) {
371  // K⁻n, K̅⁰p and anti-processes
372  total_xs = kminusn_total(sqrt_s_ * sqrt_s_);
373  }
374  } else if (meson.is_pion() && baryon.is_nucleon()) {
375  // π⁺(p,nbar), π⁻(n,pbar)
376  if ((meson.code() == pdg::pi_p &&
377  (baryon.code() == pdg::p || baryon.code() == -pdg::n)) ||
378  (meson.code() == pdg::pi_m &&
379  (baryon.code() == pdg::n || baryon.code() == -pdg::p))) {
380  total_xs = piplusp_total(sqrt_s_);
381  } else if (meson.code() == pdg::pi_z) {
382  // π⁰N
383  total_xs = 0.5 * (piplusp_total(sqrt_s_) + piminusp_total(sqrt_s_));
384  } else {
385  // π⁻(p,nbar), π⁺(n,pbar)
386  total_xs = piminusp_total(sqrt_s_);
387  }
388  } else if (meson.is_Dmeson()) {
389  const CharmRescattering& charm_rescattering =
390  finder_parameters.charm_rescattering;
391  std::optional<double> elastic_xs = std::nullopt;
392  double inelastic_xs = 0.;
393  if (baryon.is_nucleon()) {
394  elastic_xs = DN_elastic();
395  inelastic_xs = DN_inelastic();
396  } else if (baryon.is_Delta()) {
397  elastic_xs = DDelta_elastic();
398  inelastic_xs = DDelta_inelastic();
399  }
400  if ((charm_rescattering == CharmRescattering::T_Matrix) &&
401  elastic_xs.has_value()) {
402  total_xs = elastic_xs.value() + inelastic_xs;
403  } else {
404  /* use AQM if charm_rescattering == CharmRescattering::Resonances or if
405  * tmp_elastic_xs has no value, which happens either when sqrts is above
406  * the upper bound of the energy range of the underlying cross section
407  * data or there is no underlying data for the two colliding particles
408  */
410  sqrt_s_, pdg_a, pdg_b, finder_parameters.AQM_scaling_factor(pdg_a),
411  finder_parameters.AQM_scaling_factor(pdg_b));
412  }
413  } else {
414  // M*+B* goes to AQM high energy π⁻p
416  sqrt_s_, pdg_a, pdg_b, finder_parameters.AQM_scaling_factor(pdg_a),
417  finder_parameters.AQM_scaling_factor(pdg_b));
418  }
419  } else if (pdg_a.is_meson() && pdg_b.is_meson()) {
420  if (pdg_a.is_pion() && pdg_b.is_pion()) {
421  switch (pdg_a.isospin3() * pdg_b.isospin3() / 4) {
422  // π⁺π⁻
423  case -1:
424  total_xs = pipluspiminus_total(sqrt_s_);
425  break;
426  case 0:
427  // π⁰π⁰
428  if (pdg_a.isospin3() + pdg_b.isospin3() == 0) {
429  total_xs = pizeropizero_total(sqrt_s_);
430  } else {
431  // π⁺π⁰: similar to π⁺π⁻
432  total_xs = pipluspiminus_total(sqrt_s_);
433  }
434  break;
435  // π⁺π⁺ goes to π⁻p AQM
436  case 1:
437  total_xs = (2. / 3.) * piminusp_high_energy(sqrt_s_ * sqrt_s_);
438  break;
439  default:
440  throw std::runtime_error("wrong isospin in ππ scattering");
441  }
442  } else if ((pdg_a.is_Dmeson() || pdg_b.is_Dmeson()) ||
443  (pdg_a.is_Dstar2007() || pdg_b.is_Dstar2007())) {
444  const CharmRescattering& charm_rescattering =
445  finder_parameters.charm_rescattering;
446  std::optional<double> elastic_xs = std::nullopt;
447  double inelastic_xs = 0.;
448  if (pdg_a.is_pion() || pdg_b.is_pion()) {
449  elastic_xs = Dpi_and_Dstarpi_elastic();
450  inelastic_xs = Dpi_and_Dstarpi_inelastic();
451  } else if (pdg_a.is_eta() || pdg_b.is_eta()) {
452  elastic_xs = Deta_and_Dstareta_elastic();
453  // no inelastic scattering for Deta or D*eta
454  } else if (pdg_a.is_kaon() || pdg_b.is_kaon()) {
455  elastic_xs = DK_and_DstarK_elastic();
456  inelastic_xs = DK_and_DstarK_inelastic();
457  }
458  if ((charm_rescattering == CharmRescattering::T_Matrix) &&
459  elastic_xs.has_value()) {
460  total_xs = elastic_xs.value() + inelastic_xs;
461  } else {
462  /* use AQM if charm_rescattering == CharmRescattering::Resonances or if
463  * tmp_elastic_xs has no value, which happens either when sqrts is above
464  * the upper bound of the energy range of the underlying cross section
465  * data or there is no underlying data for the two colliding particles
466  */
468  sqrt_s_, pdg_a, pdg_b, finder_parameters.AQM_scaling_factor(pdg_a),
469  finder_parameters.AQM_scaling_factor(pdg_b));
470  }
471  } else {
472  // M*+M* goes to AQM high energy π⁻p
474  sqrt_s_, pdg_a, pdg_b, finder_parameters.AQM_scaling_factor(pdg_a),
475  finder_parameters.AQM_scaling_factor(pdg_b));
476  }
477  }
478  return (total_xs + finder_parameters.additional_el_xs) *
479  finder_parameters.scale_xs;
480 }
481 
482 CollisionBranchPtr CrossSections::elastic(
483  const ScatterActionsFinderParameters& finder_parameters) const {
484  double elastic_xs = 0.;
485 
486  if (finder_parameters.elastic_parameter >= 0.) {
487  // use constant elastic cross section from config file
488  elastic_xs = finder_parameters.elastic_parameter;
489  } else {
490  // use parametrization
491  elastic_xs = elastic_parametrization(finder_parameters);
492  }
493  /* when using a factor to scale the cross section and an additional
494  * contribution to the elastic cross section, the contribution is added first
495  * and then everything is scaled */
496  return std::make_unique<CollisionBranch>(
497  incoming_particles_[0].type(), incoming_particles_[1].type(),
498  (elastic_xs + finder_parameters.additional_el_xs) *
499  finder_parameters.scale_xs,
501 }
502 
503 CollisionBranchList CrossSections::rare_two_to_two() const {
504  CollisionBranchList process_list;
505  const ParticleData& data_a = incoming_particles_[0];
506  const ParticleData& data_b = incoming_particles_[1];
507  const auto& pdg_a = data_a.pdgcode();
508  const auto& pdg_b = data_b.pdgcode();
509  if ((pdg_a.is_nucleon() && pdg_b.is_pion()) ||
510  (pdg_b.is_nucleon() && pdg_a.is_pion())) {
511  process_list = npi_yk();
512  }
513  return process_list;
514 }
515 
517  const ScatterActionsFinderParameters& finder_parameters) const {
518  const bool use_AQM = finder_parameters.use_AQM;
519  const double pipi_offset =
520  finder_parameters.transition_high_energy.pipi_offset;
521  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
522  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
523  double elastic_xs = 0.0;
524  if ((pdg_a.is_nucleon() && pdg_b.is_pion()) ||
525  (pdg_b.is_nucleon() && pdg_a.is_pion())) {
526  // Elastic Nucleon Pion Scattering
527  elastic_xs = npi_el();
528  } else if ((pdg_a.is_nucleon() && pdg_b.is_kaon()) ||
529  (pdg_b.is_nucleon() && pdg_a.is_kaon())) {
530  // Elastic Nucleon Kaon Scattering
531  elastic_xs = nk_el();
532  } else if (pdg_a.is_Dmeson() || pdg_b.is_Dmeson()) {
533  const CharmRescattering& charm_rescattering =
534  finder_parameters.charm_rescattering;
535  std::optional<double> tmp_elastic_xs = std::nullopt;
536  if (pdg_a.is_nucleon() || pdg_b.is_nucleon()) {
537  tmp_elastic_xs = DN_elastic();
538  } else if (pdg_a.is_Delta() || pdg_b.is_Delta()) {
539  tmp_elastic_xs = DDelta_elastic();
540  }
541  if ((charm_rescattering == CharmRescattering::T_Matrix) &&
542  tmp_elastic_xs.has_value()) {
543  elastic_xs = tmp_elastic_xs.value();
544  } else if (use_AQM) {
545  /* use AQM if charm_rescattering == CharmRescattering::Resonances or if
546  * tmp_elastic_xs has no value, which happens either when sqrts is above
547  * the upper bound of the energy range of the underlying cross section
548  * data or there is no underlying data for the two colliding particles */
549  const double m1 = incoming_particles_[0].effective_mass();
550  const double m2 = incoming_particles_[1].effective_mass();
551  const double s = sqrt_s_ * sqrt_s_;
552  elastic_xs = 2. / 3. * piplusp_elastic_AQM(s, m1, m2) *
553  finder_parameters.AQM_scaling_factor(pdg_a) *
554  finder_parameters.AQM_scaling_factor(pdg_b);
555  } else {
556  const ParticleType& a = incoming_particles_[0].type();
557  const ParticleType& b = incoming_particles_[1].type();
559  charm_rescattering);
560  return 0.;
561  }
562  } else if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
563  pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
564  // Elastic Nucleon Nucleon Scattering
565  elastic_xs = nn_el();
566  } else if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
567  pdg_a.antiparticle_sign() == -pdg_b.antiparticle_sign()) {
568  // Elastic Nucleon anti-Nucleon Scattering
569  elastic_xs = ppbar_elastic(sqrt_s_ * sqrt_s_);
570  } else if (pdg_a.is_nucleus() || pdg_b.is_nucleus()) {
571  const PdgCode& pdg_nucleus = pdg_a.is_nucleus() ? pdg_a : pdg_b;
572  const PdgCode& pdg_other = pdg_a.is_nucleus() ? pdg_b : pdg_a;
573  const bool is_deuteron = pdg_nucleus.is_deuteron(); // d or anti-d
574  if (is_deuteron && pdg_other.is_pion()) {
575  // Elastic (Anti-)deuteron Pion Scattering
576  elastic_xs = deuteron_pion_elastic(sqrt_s_ * sqrt_s_);
577  } else if (is_deuteron && pdg_other.is_nucleon()) {
578  // Elastic (Anti-)deuteron (Anti-)Nucleon Scattering
579  elastic_xs = deuteron_nucleon_elastic(sqrt_s_ * sqrt_s_);
580  }
581  } else if ((pdg_a.is_Dmeson() || pdg_b.is_Dmeson()) ||
582  (pdg_a.is_Dstar2007() || pdg_b.is_Dstar2007())) {
583  const CharmRescattering& charm_rescattering =
584  finder_parameters.charm_rescattering;
585  std::optional<double> tmp_elastic_xs = std::nullopt;
586  if (pdg_a.is_pion() || pdg_b.is_pion()) {
587  tmp_elastic_xs = Dpi_and_Dstarpi_elastic();
588  } else if (pdg_a.is_eta() || pdg_b.is_eta()) {
589  tmp_elastic_xs = Deta_and_Dstareta_elastic();
590  } else if (pdg_a.is_kaon() || pdg_b.is_kaon()) {
591  tmp_elastic_xs = DK_and_DstarK_elastic();
592  }
593  if ((charm_rescattering == CharmRescattering::T_Matrix) &&
594  tmp_elastic_xs.has_value()) {
595  elastic_xs = tmp_elastic_xs.value();
596  } else if (use_AQM) {
597  /* use AQM if charm_rescattering == CharmRescattering::Resonances or if
598  * tmp_elastic_xs has no value, which happens either when sqrts is above
599  * the upper bound of the energy range of the underlying cross section
600  * data or there is no underlying data for the two colliding particles */
601  const double m1 = incoming_particles_[0].effective_mass();
602  const double m2 = incoming_particles_[1].effective_mass();
603  const double s = sqrt_s_ * sqrt_s_;
604  elastic_xs = 2. / 3. * piplusp_elastic_AQM(s, m1, m2) *
605  finder_parameters.AQM_scaling_factor(pdg_a) *
606  finder_parameters.AQM_scaling_factor(pdg_b);
607  } else {
608  const ParticleType& a = incoming_particles_[0].type();
609  const ParticleType& b = incoming_particles_[1].type();
611  charm_rescattering);
612  return 0.;
613  }
614  } else if (use_AQM) {
615  const double m1 = incoming_particles_[0].effective_mass();
616  const double m2 = incoming_particles_[1].effective_mass();
617  const double s = sqrt_s_ * sqrt_s_;
618  if (pdg_a.is_baryon() && pdg_b.is_baryon()) {
619  elastic_xs = nn_el(); // valid also for annihilation
620  } else if ((pdg_a.is_meson() && pdg_b.is_baryon()) ||
621  (pdg_b.is_meson() && pdg_a.is_baryon())) {
622  elastic_xs = piplusp_elastic_AQM(s, m1, m2);
623  } else if (pdg_a.is_meson() && pdg_b.is_meson()) {
624  /* Special case: the pi+pi- elastic cross-section goes through resonances
625  * at low sqrt_s, so we turn it off for this region so as not to destroy
626  * the agreement with experimental data; this does not
627  * apply to other pi pi cross-sections, which do not have any data */
628  if (((pdg_a == pdg::pi_p && pdg_b == pdg::pi_m) ||
629  (pdg_a == pdg::pi_m && pdg_b == pdg::pi_p)) &&
630  (m1 + m2 + pipi_offset) > sqrt_s_) {
631  elastic_xs = 0.0;
632  } else {
633  // meson-meson goes through scaling from π+p parametrization
634  elastic_xs = 2. / 3. * piplusp_elastic_AQM(s, m1, m2);
635  }
636  }
637  elastic_xs *= finder_parameters.AQM_scaling_factor(pdg_a) *
638  finder_parameters.AQM_scaling_factor(pdg_b);
639  }
640  if (elastic_xs < 0.) {
642  incoming_particles_[1], __func__);
643  } else {
644  return elastic_xs;
645  }
646 }
647 
648 double CrossSections::nn_el() const {
649  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
650  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
651 
652  // Use parametrized cross sections.
653  double sig_el = -1.;
654  const double s = sqrt_s_ * sqrt_s_;
655  const bool is_NN_pair = pdg_a.is_nucleon() && pdg_b.is_nucleon();
656  if (is_NN_pair) {
657  if (is_BBbar_pair_) {
658  // npbar and ppbar
659  sig_el = ppbar_elastic(s);
660  } else {
661  sig_el = (pdg_a == pdg_b) ? pp_elastic(s) : np_elastic(s);
662  }
663  } else {
664  // AQM - Additive Quark Model
665  const double m1 = incoming_particles_[0].effective_mass();
666  const double m2 = incoming_particles_[1].effective_mass();
667  if (is_BBbar_pair_) {
668  sig_el =
670  } else {
671  sig_el = pp_elastic_high_energy(s, m1, m2);
672  }
673  }
674 
675  if (sig_el > 0.) {
676  return sig_el;
677  } else {
679  incoming_particles_[1], __func__);
680  }
681 }
682 
683 double CrossSections::npi_el() const {
684  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
685  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
686 
687  const PdgCode& nucleon = pdg_a.is_nucleon() ? pdg_a : pdg_b;
688  const PdgCode& pion = pdg_a.is_nucleon() ? pdg_b : pdg_a;
689  assert(pion != nucleon);
690 
691  const double s = sqrt_s_ * sqrt_s_;
692 
693  double sig_el = 0.;
694  switch (nucleon.code()) {
695  case pdg::p:
696  switch (pion.code()) {
697  case pdg::pi_p:
698  sig_el = piplusp_elastic(s);
699  break;
700  case pdg::pi_m:
701  sig_el = piminusp_elastic(s);
702  break;
703  case pdg::pi_z:
704  sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
705  break;
706  }
707  break;
708  case pdg::n:
709  switch (pion.code()) {
710  case pdg::pi_p:
711  sig_el = piminusp_elastic(s);
712  break;
713  case pdg::pi_m:
714  sig_el = piplusp_elastic(s);
715  break;
716  case pdg::pi_z:
717  sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
718  break;
719  }
720  break;
721  case -pdg::p:
722  switch (pion.code()) {
723  case pdg::pi_p:
724  sig_el = piminusp_elastic(s);
725  break;
726  case pdg::pi_m:
727  sig_el = piplusp_elastic(s);
728  break;
729  case pdg::pi_z:
730  sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
731  break;
732  }
733  break;
734  case -pdg::n:
735  switch (pion.code()) {
736  case pdg::pi_p:
737  sig_el = piplusp_elastic(s);
738  break;
739  case pdg::pi_m:
740  sig_el = piminusp_elastic(s);
741  break;
742  case pdg::pi_z:
743  sig_el = 0.5 * (piplusp_elastic(s) + piminusp_elastic(s));
744  break;
745  }
746  break;
747  default:
749  incoming_particles_[1], __func__);
750  }
751 
752  if (sig_el > 0) {
753  return sig_el;
754  } else {
756  incoming_particles_[1], __func__);
757  }
758 }
759 
760 CollisionBranchList CrossSections::npi_yk() const {
761  const ParticleType& a = incoming_particles_[0].type();
762  const ParticleType& b = incoming_particles_[1].type();
763  const ParticleType& type_nucleon = a.pdgcode().is_nucleon() ? a : b;
764  const ParticleType& type_pion = a.pdgcode().is_nucleon() ? b : a;
765 
766  const auto pdg_nucleon = type_nucleon.pdgcode().code();
767  const auto pdg_pion = type_pion.pdgcode().code();
768 
769  const double s = sqrt_s_ * sqrt_s_;
770 
771  /* The cross sections are paramectrized for four isospin channels. The
772  * cross sections of the rest isospin channels are obtained using
773  * Clebsch-Gordan coefficients */
774 
775  CollisionBranchList process_list;
776  switch (pdg_nucleon) {
777  case pdg::p: {
778  switch (pdg_pion) {
779  case pdg::pi_p: {
780  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
781  const auto& type_K_p = ParticleType::find(pdg::K_p);
782  add_channel(
783  process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
784  sqrt_s_, type_K_p, type_Sigma_p);
785  break;
786  }
787  case pdg::pi_m: {
788  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
789  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
790  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
791  const auto& type_K_p = ParticleType::find(pdg::K_p);
792  const auto& type_K_z = ParticleType::find(pdg::K_z);
793  add_channel(
794  process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
795  sqrt_s_, type_K_p, type_Sigma_m);
796  add_channel(
797  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
798  type_K_z, type_Sigma_z);
799  add_channel(
800  process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
801  type_K_z, type_Lambda);
802  break;
803  }
804  case pdg::pi_z: {
805  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
806  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
807  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
808  const auto& type_K_p = ParticleType::find(pdg::K_p);
809  const auto& type_K_z = ParticleType::find(pdg::K_z);
810  add_channel(
811  process_list,
812  [&] {
813  return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
816  },
817  sqrt_s_, type_K_p, type_Sigma_z);
818  add_channel(
819  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
820  type_K_z, type_Sigma_p);
821  add_channel(
822  process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
823  sqrt_s_, type_K_p, type_Lambda);
824  break;
825  }
826  }
827  break;
828  }
829  case pdg::n: {
830  switch (pdg_pion) {
831  case pdg::pi_p: {
832  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
833  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
834  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
835  const auto& type_K_p = ParticleType::find(pdg::K_p);
836  const auto& type_K_z = ParticleType::find(pdg::K_z);
837  add_channel(
838  process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
839  sqrt_s_, type_K_z, type_Sigma_p);
840  add_channel(
841  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
842  type_K_p, type_Sigma_z);
843  add_channel(
844  process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
845  type_K_p, type_Lambda);
846  break;
847  }
848  case pdg::pi_m: {
849  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
850  const auto& type_K_z = ParticleType::find(pdg::K_z);
851  add_channel(
852  process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
853  sqrt_s_, type_K_z, type_Sigma_m);
854  break;
855  }
856  case pdg::pi_z: {
857  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
858  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
859  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
860  const auto& type_K_p = ParticleType::find(pdg::K_p);
861  const auto& type_K_z = ParticleType::find(pdg::K_z);
862  add_channel(
863  process_list,
864  [&] {
865  return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
868  },
869  sqrt_s_, type_K_z, type_Sigma_z);
870  add_channel(
871  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
872  type_K_p, type_Sigma_m);
873  add_channel(
874  process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
875  sqrt_s_, type_K_z, type_Lambda);
876  break;
877  }
878  }
879  break;
880  }
881  case -pdg::p: {
882  switch (pdg_pion) {
883  case pdg::pi_p: {
884  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
885  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
886  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
887  const auto& type_K_m = ParticleType::find(-pdg::K_p);
888  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
889  add_channel(
890  process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
891  sqrt_s_, type_K_m, type_Sigma_m_bar);
892  add_channel(
893  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
894  type_Kbar_z, type_Sigma_z_bar);
895  add_channel(
896  process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
897  type_Kbar_z, type_Lambda_bar);
898  break;
899  }
900  case pdg::pi_m: {
901  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
902  const auto& type_K_m = ParticleType::find(-pdg::K_p);
903  add_channel(
904  process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
905  sqrt_s_, type_K_m, type_Sigma_p_bar);
906  break;
907  }
908  case pdg::pi_z: {
909  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
910  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
911  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
912  const auto& type_K_m = ParticleType::find(-pdg::K_p);
913  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
914  add_channel(
915  process_list,
916  [&] {
917  return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
920  },
921  sqrt_s_, type_K_m, type_Sigma_z_bar);
922  add_channel(
923  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
924  type_Kbar_z, type_Sigma_p_bar);
925  add_channel(
926  process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
927  sqrt_s_, type_K_m, type_Lambda_bar);
928  break;
929  }
930  }
931  break;
932  }
933  case -pdg::n: {
934  switch (pdg_pion) {
935  case pdg::pi_p: {
936  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
937  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
938  add_channel(
939  process_list, [&] { return piplusp_sigmapluskplus_pdg(s); },
940  sqrt_s_, type_Kbar_z, type_Sigma_m_bar);
941  break;
942  }
943  case pdg::pi_m: {
944  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
945  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
946  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
947  const auto& type_K_m = ParticleType::find(-pdg::K_p);
948  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
949  add_channel(
950  process_list, [&] { return piminusp_sigmaminuskplus_pdg(s); },
951  sqrt_s_, type_Kbar_z, type_Sigma_p_bar);
952  add_channel(
953  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
954  type_K_m, type_Sigma_z_bar);
955  add_channel(
956  process_list, [&] { return piminusp_lambdak0_pdg(s); }, sqrt_s_,
957  type_K_m, type_Lambda_bar);
958  break;
959  }
960  case pdg::pi_z: {
961  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
962  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
963  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
964  const auto& type_K_m = ParticleType::find(-pdg::K_p);
965  const auto& type_Kbar_z = ParticleType::find(-pdg::K_z);
966  add_channel(
967  process_list,
968  [&] {
969  return 0.5 * (piplusp_sigmapluskplus_pdg(s) -
972  },
973  sqrt_s_, type_Kbar_z, type_Sigma_z_bar);
974  add_channel(
975  process_list, [&] { return piminusp_sigma0k0_res(s); }, sqrt_s_,
976  type_K_m, type_Sigma_m_bar);
977  add_channel(
978  process_list, [&] { return 0.5 * piminusp_lambdak0_pdg(s); },
979  sqrt_s_, type_Kbar_z, type_Lambda_bar);
980  break;
981  }
982  }
983  break;
984  }
985  }
986 
987  return process_list;
988 }
989 
990 double CrossSections::nk_el() const {
991  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
992  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
993 
994  const PdgCode& nucleon = pdg_a.is_nucleon() ? pdg_a : pdg_b;
995  const PdgCode& kaon = pdg_a.is_nucleon() ? pdg_b : pdg_a;
996  assert(kaon != nucleon);
997 
998  const double s = sqrt_s_ * sqrt_s_;
999 
1000  double sig_el = 0.;
1001  switch (nucleon.code()) {
1002  case pdg::p:
1003  switch (kaon.code()) {
1004  case pdg::K_p:
1005  sig_el = kplusp_elastic_background(s);
1006  break;
1007  case pdg::K_m:
1008  sig_el = kminusp_elastic_background(s);
1009  break;
1010  case pdg::K_z:
1011  sig_el = k0p_elastic_background(s);
1012  break;
1013  case pdg::Kbar_z:
1014  sig_el = kbar0p_elastic_background(s);
1015  break;
1016  }
1017  break;
1018  case pdg::n:
1019  switch (kaon.code()) {
1020  case pdg::K_p:
1021  sig_el = kplusn_elastic_background(s);
1022  break;
1023  case pdg::K_m:
1024  sig_el = kminusn_elastic_background(s);
1025  break;
1026  case pdg::K_z:
1027  sig_el = k0n_elastic_background(s);
1028  break;
1029  case pdg::Kbar_z:
1030  sig_el = kbar0n_elastic_background(s);
1031  break;
1032  }
1033  break;
1034  case -pdg::p:
1035  switch (kaon.code()) {
1036  case pdg::K_p:
1037  sig_el = kminusp_elastic_background(s);
1038  break;
1039  case pdg::K_m:
1040  sig_el = kplusp_elastic_background(s);
1041  break;
1042  case pdg::K_z:
1043  sig_el = kbar0p_elastic_background(s);
1044  break;
1045  case pdg::Kbar_z:
1046  sig_el = k0p_elastic_background(s);
1047  break;
1048  }
1049  break;
1050  case -pdg::n:
1051  switch (kaon.code()) {
1052  case pdg::K_p:
1053  sig_el = kminusn_elastic_background(s);
1054  break;
1055  case pdg::K_m:
1056  sig_el = kplusn_elastic_background(s);
1057  break;
1058  case pdg::K_z:
1059  sig_el = kbar0n_elastic_background(s);
1060  break;
1061  case pdg::Kbar_z:
1062  sig_el = k0n_elastic_background(s);
1063  break;
1064  }
1065  break;
1066  default:
1068  incoming_particles_[1], __func__);
1069  }
1070 
1071  if (sig_el > 0) {
1072  return sig_el;
1073  } else {
1075  incoming_particles_[1], __func__);
1076  }
1077 }
1078 
1079 std::optional<double> CrossSections::Dpi_and_Dstarpi_elastic() const {
1080  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
1081  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
1082  const auto pdg_D =
1083  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
1084  const auto pdg_pion =
1085  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
1086 
1087  std::optional<double> sig_el = std::nullopt;
1088  switch (pack(pdg_D, pdg_pion)) {
1089  // Checks for D mesons scatterings
1090  case pack(pdg::D_z, pdg::pi_p):
1091  case pack(pdg::Dbar_z, pdg::pi_m): { // Same xsec for charge conjugation.
1092  sig_el = Dzeropiplus_elastic(sqrt_s_);
1093  break;
1094  }
1095  case pack(pdg::D_z, pdg::pi_m):
1096  case pack(pdg::Dbar_z, pdg::pi_p): { // Same xsec for charge conjugation.
1097  sig_el = Dzeropiminus_elastic(sqrt_s_);
1098  break;
1099  }
1100  case pack(pdg::D_z, pdg::pi_z):
1101  case pack(pdg::Dbar_z, pdg::pi_z): { // Same xsec for charge conjugation.
1102  sig_el = Dzeropizero_elastic(sqrt_s_);
1103  break;
1104  }
1105  case pack(pdg::D_p, pdg::pi_p):
1106  case pack(pdg::D_m, pdg::pi_m): { // Same xsec for charge conjugation.
1107  sig_el = Dpluspiplus_elastic(sqrt_s_);
1108  break;
1109  }
1110  case pack(pdg::D_p, pdg::pi_m):
1111  case pack(pdg::D_m, pdg::pi_p): { // Same xsec for charge conjugation.
1112  sig_el = Dpluspiminus_elastic(sqrt_s_);
1113  break;
1114  }
1115  case pack(pdg::D_p, pdg::pi_z):
1116  case pack(pdg::D_m, pdg::pi_z): { // Same xsec for charge conjugation.
1117  sig_el = Dpluspizero_elastic(sqrt_s_);
1118  break;
1119  }
1120  // Checks for D* mesons scatterings
1121  case pack(pdg::Dstar_z, pdg::pi_p):
1122  case pack(pdg::Dstarbar_z, pdg::pi_m): { // Same xs for charge conjugation.
1123  sig_el = Dstarzeropiplus_elastic(sqrt_s_);
1124  break;
1125  }
1126  case pack(pdg::Dstar_z, pdg::pi_m):
1127  case pack(pdg::Dstarbar_z, pdg::pi_p): { // Same xs for charge conjugation.
1129  break;
1130  }
1131  case pack(pdg::Dstar_z, pdg::pi_z):
1132  case pack(pdg::Dstarbar_z, pdg::pi_z): { // Same xs for charge conjugation.
1133  sig_el = Dstarzeropizero_elastic(sqrt_s_);
1134  break;
1135  }
1136  case pack(pdg::Dstar_p, pdg::pi_p):
1137  case pack(pdg::Dstar_m, pdg::pi_m): { // Same xsec for charge conjugation.
1138  sig_el = Dstarpluspiplus_elastic(sqrt_s_);
1139  break;
1140  }
1141  case pack(pdg::Dstar_p, pdg::pi_m):
1142  case pack(pdg::Dstar_m, pdg::pi_p): { // Same xsec for charge conjugation.
1144  break;
1145  }
1146  case pack(pdg::Dstar_p, pdg::pi_z):
1147  case pack(pdg::Dstar_m, pdg::pi_z): { // Same xsec for charge conjugation.
1148  sig_el = Dstarpluspizero_elastic(sqrt_s_);
1149  break;
1150  }
1151  default:
1153  incoming_particles_[1], __func__);
1154  }
1155 
1156  if (sig_el.has_value() && sig_el.value() < 0.) {
1158  incoming_particles_[1], __func__);
1159  } else {
1160  return sig_el;
1161  }
1162 }
1163 
1164 std::optional<double> CrossSections::Deta_and_Dstareta_elastic() const {
1165  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
1166  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
1167  const auto pdg_D =
1168  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
1169  const auto pdg_eta =
1170  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
1171 
1172  std::optional<double> sig_el = std::nullopt;
1173  switch (pack(pdg_D, pdg_eta)) {
1174  // Checks for D mesons scatterings
1175  case pack(pdg::D_p, pdg::eta):
1176  case pack(pdg::D_m, pdg::eta): { // Same xsec for charge conjugation.
1177  sig_el = Dpluseta_elastic(sqrt_s_);
1178  break;
1179  }
1180  case pack(pdg::D_z, pdg::eta):
1181  case pack(pdg::Dbar_z, pdg::eta): { // Same xsec for charge conjugation.
1182  sig_el = Dzeroeta_elastic(sqrt_s_);
1183  break;
1184  }
1185  // Checks for D* mesons scatterings
1186  case pack(pdg::Dstar_p, pdg::eta):
1187  case pack(pdg::Dstar_m, pdg::eta): { // Same xsec for charge conjugation.
1188  sig_el = Dstarpluseta_elastic(sqrt_s_);
1189  break;
1190  }
1191  case pack(pdg::Dstar_z, pdg::eta):
1192  case pack(pdg::Dstarbar_z, pdg::eta): { // Same xs for charge conjugation.
1193  sig_el = Dstarzeroeta_elastic(sqrt_s_);
1194  break;
1195  }
1196  default:
1198  incoming_particles_[1], __func__);
1199  }
1200 
1201  if (sig_el.has_value() && sig_el.value() < 0.) {
1203  incoming_particles_[1], __func__);
1204  } else {
1205  return sig_el;
1206  }
1207 }
1208 
1209 std::optional<double> CrossSections::DK_and_DstarK_elastic() const {
1210  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
1211  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
1212  const auto pdg_D =
1213  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
1214  const auto pdg_kaon =
1215  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
1216 
1217  std::optional<double> sig_el = std::nullopt;
1218  switch (pack(pdg_D, pdg_kaon)) {
1219  // Checks for D mesons scatterings
1220  case pack(pdg::D_p, pdg::K_p):
1221  case pack(pdg::D_m, pdg::K_m): { // Same xsec for charge conjugation.
1222  sig_el = DplusKplus_elastic(sqrt_s_);
1223  break;
1224  }
1225  case pack(pdg::D_p, pdg::K_z):
1226  case pack(pdg::D_m, pdg::Kbar_z): { // Same xsec for charge conjugation.
1227  sig_el = DplusKzero_elastic(sqrt_s_);
1228  break;
1229  }
1230  case pack(pdg::D_z, pdg::K_p):
1231  case pack(pdg::Dbar_z, pdg::K_m): { // Same xsec for charge conjugation.
1232  sig_el = DzeroKplus_elastic(sqrt_s_);
1233  break;
1234  }
1235  case pack(pdg::D_z, pdg::K_z):
1236  case pack(pdg::Dbar_z, pdg::Kbar_z): { // Same xsec for charge conjugation.
1237  sig_el = DzeroKzero_elastic(sqrt_s_);
1238  break;
1239  }
1240  case pack(pdg::D_p, pdg::Kbar_z):
1241  case pack(pdg::D_m, pdg::K_z): { // Same xsec for charge conjugation.
1242  sig_el = DplusKbarzero_elastic(sqrt_s_);
1243  break;
1244  }
1245  case pack(pdg::D_p, pdg::K_m):
1246  case pack(pdg::D_m, pdg::K_p): { // Same xsec for charge conjugation.
1247  sig_el = DplusKminus_elastic(sqrt_s_);
1248  break;
1249  }
1250  case pack(pdg::D_z, pdg::Kbar_z):
1251  case pack(pdg::Dbar_z, pdg::K_z): { // Same xsec for charge conjugation.
1252  sig_el = DzeroKbarzero_elastic(sqrt_s_);
1253  break;
1254  }
1255  case pack(pdg::D_z, pdg::K_m):
1256  case pack(pdg::Dbar_z, pdg::K_p): { // Same xsec for charge conjugation.
1257  sig_el = DzeroKminus_elastic(sqrt_s_);
1258  break;
1259  }
1260  // Checks for D* mesons scatterings
1261  case pack(pdg::Dstar_p, pdg::K_p):
1262  case pack(pdg::Dstar_m, pdg::K_m): { // Same xsec for charge conjugation.
1263  sig_el = DstarplusKplus_elastic(sqrt_s_);
1264  break;
1265  }
1266  case pack(pdg::Dstar_p, pdg::K_z):
1267  case pack(pdg::Dstar_m, pdg::Kbar_z): { // Same xs for charge conjugation.
1268  sig_el = DstarplusKzero_elastic(sqrt_s_);
1269  break;
1270  }
1271  case pack(pdg::Dstar_z, pdg::K_p):
1272  case pack(pdg::Dstarbar_z, pdg::K_m): { // Same xs for charge conjugation.
1273  sig_el = DstarzeroKplus_elastic(sqrt_s_);
1274  break;
1275  }
1276  case pack(pdg::Dstar_z, pdg::K_z):
1277  case pack(pdg::Dstarbar_z, pdg::Kbar_z): { // Same xs for charge conjugat.
1278  sig_el = DstarzeroKzero_elastic(sqrt_s_);
1279  break;
1280  }
1281  case pack(pdg::Dstar_p, pdg::Kbar_z):
1282  case pack(pdg::Dstar_m, pdg::K_z): { // Same xsec for charge conjugation.
1284  break;
1285  }
1286  case pack(pdg::Dstar_p, pdg::K_m):
1287  case pack(pdg::Dstar_m, pdg::K_p): { // Same xsec for charge conjugation.
1288  sig_el = DstarplusKminus_elastic(sqrt_s_);
1289  break;
1290  }
1291  case pack(pdg::Dstar_z, pdg::Kbar_z):
1292  case pack(pdg::Dstarbar_z, pdg::K_z): { // Same xs for charge conjugation.
1294  break;
1295  }
1296  case pack(pdg::Dstar_z, pdg::K_m):
1297  case pack(pdg::Dstarbar_z, pdg::K_p): { // Same xs for charge conjugation.
1298  sig_el = DstarzeroKminus_elastic(sqrt_s_);
1299  break;
1300  }
1301  default:
1303  incoming_particles_[1], __func__);
1304  }
1305 
1306  if (sig_el.has_value() && sig_el.value() < 0.) {
1308  incoming_particles_[1], __func__);
1309  } else {
1310  return sig_el;
1311  }
1312 }
1313 
1314 std::optional<double> CrossSections::DN_elastic() const {
1315  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
1316  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
1317  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
1318  const auto pdg_nucleon = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
1319 
1320  std::optional<double> sig_el = std::nullopt;
1321  switch (pack(pdg_D, pdg_nucleon)) {
1322  case pack(pdg::D_p, pdg::n):
1323  case pack(pdg::D_m, -pdg::n): { // Same xsec for charge conjugation.
1324  sig_el = Dplusn_elastic(sqrt_s_);
1325  break;
1326  }
1327  case pack(pdg::D_p, pdg::p):
1328  case pack(pdg::D_m, -pdg::p): { // Same xsec for charge conjugation.
1329  sig_el = Dplusp_elastic(sqrt_s_);
1330  break;
1331  }
1332  case pack(pdg::D_z, pdg::n):
1333  case pack(pdg::Dbar_z, -pdg::n): { // Same xsec for charge conjugation.
1334  sig_el = Dzeron_elastic(sqrt_s_);
1335  break;
1336  }
1337  case pack(pdg::D_z, pdg::p):
1338  case pack(pdg::Dbar_z, -pdg::p): { // Same xsec for charge conjugation.
1339  sig_el = Dzerop_elastic(sqrt_s_);
1340  break;
1341  }
1342  case pack(pdg::D_m, pdg::n):
1343  case pack(pdg::D_p, -pdg::n): { // Same xsec for charge conjugation.
1344  sig_el = Dminusn_elastic(sqrt_s_);
1345  break;
1346  }
1347  case pack(pdg::D_m, pdg::p):
1348  case pack(pdg::D_p, -pdg::p): { // Same xsec for charge conjugation.
1349  sig_el = Dminusp_elastic(sqrt_s_);
1350  break;
1351  }
1352  case pack(pdg::Dbar_z, pdg::n):
1353  case pack(pdg::D_z, -pdg::n): { // Same xsec for charge conjugation.
1354  sig_el = Dbarzeron_elastic(sqrt_s_);
1355  break;
1356  }
1357  case pack(pdg::Dbar_z, pdg::p):
1358  case pack(pdg::D_z, -pdg::p): { // Same xsec for charge conjugation.
1359  sig_el = Dbarzerop_elastic(sqrt_s_);
1360  break;
1361  }
1362  default:
1364  incoming_particles_[1], __func__);
1365  }
1366 
1367  if (sig_el.has_value() && sig_el.value() < 0.) {
1369  incoming_particles_[1], __func__);
1370  } else {
1371  return sig_el;
1372  }
1373 }
1374 
1375 std::optional<double> CrossSections::DDelta_elastic() const {
1376  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
1377  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
1378  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
1379  const auto pdg_Delta = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
1380 
1381  std::optional<double> sig_el = std::nullopt;
1382  switch (pack(pdg_D, pdg_Delta)) {
1383  case pack(pdg::D_p, pdg::Delta_p):
1384  case pack(pdg::D_m, -pdg::Delta_p): { // Same xsec for charge conjugation.
1385  sig_el = DplusDeltaplus_elastic(sqrt_s_);
1386  break;
1387  }
1388  case pack(pdg::D_p, pdg::Delta_pp):
1389  case pack(pdg::D_m, -pdg::Delta_pp): { // Same xsec for charge conjugation.
1391  break;
1392  }
1393  case pack(pdg::D_p, pdg::Delta_m):
1394  case pack(pdg::D_m, -pdg::Delta_m): { // Same xsec for charge conjugation.
1395  sig_el = DplusDeltaminus_elastic(sqrt_s_);
1396  break;
1397  }
1398  case pack(pdg::D_p, pdg::Delta_z):
1399  case pack(pdg::D_m, -pdg::Delta_z): { // Same xsec for charge conjugation.
1400  sig_el = DplusDeltazero_elastic(sqrt_s_);
1401  break;
1402  }
1403  case pack(pdg::D_z, pdg::Delta_p):
1404  case pack(pdg::Dbar_z, -pdg::Delta_p): { // Same xs for charge conjugation.
1405  sig_el = DzeroDeltaplus_elastic(sqrt_s_);
1406  break;
1407  }
1408  case pack(pdg::D_z, pdg::Delta_pp):
1409  case pack(pdg::Dbar_z, -pdg::Delta_pp): { // Same xs for charge conjugate.
1411  break;
1412  }
1413  case pack(pdg::D_z, pdg::Delta_m):
1414  case pack(pdg::Dbar_z, -pdg::Delta_m): { // Same xs for charge conjugation.
1415  sig_el = DzeroDeltaminus_elastic(sqrt_s_);
1416  break;
1417  }
1418  case pack(pdg::D_z, pdg::Delta_z):
1419  case pack(pdg::Dbar_z, -pdg::Delta_z): { // Same xs for charge conjugation.
1420  sig_el = DzeroDeltazero_elastic(sqrt_s_);
1421  break;
1422  }
1423  case pack(pdg::D_m, pdg::Delta_p):
1424  case pack(pdg::D_p, -pdg::Delta_p): { // Same xsec for charge conjugation.
1425  sig_el = DminusDeltaplus_elastic(sqrt_s_);
1426  break;
1427  }
1428  case pack(pdg::D_m, pdg::Delta_pp):
1429  case pack(pdg::D_p, -pdg::Delta_pp): { // Same xsec for charge conjugation.
1431  break;
1432  }
1433  case pack(pdg::D_m, pdg::Delta_m):
1434  case pack(pdg::D_p, -pdg::Delta_m): { // Same xsec for charge conjugation.
1436  break;
1437  }
1438  case pack(pdg::D_m, pdg::Delta_z):
1439  case pack(pdg::D_p, -pdg::Delta_z): { // Same xsec for charge conjugation.
1440  sig_el = DminusDeltazero_elastic(sqrt_s_);
1441  break;
1442  }
1443  case pack(pdg::Dbar_z, pdg::Delta_p):
1444  case pack(pdg::D_z, -pdg::Delta_p): { // Same xsec for charge conjugation.
1446  break;
1447  }
1449  case pack(pdg::D_z, -pdg::Delta_pp): { // Same xsec for charge conjugate.
1451  break;
1452  }
1453  case pack(pdg::Dbar_z, pdg::Delta_m):
1454  case pack(pdg::D_z, -pdg::Delta_m): { // Same xsec for charge conjugation.
1456  break;
1457  }
1458  case pack(pdg::Dbar_z, pdg::Delta_z):
1459  case pack(pdg::D_z, -pdg::Delta_z): { // Same xsec for charge conjugation.
1461  break;
1462  }
1463  default:
1465  incoming_particles_[1], __func__);
1466  }
1467 
1468  if (sig_el.has_value() && sig_el.value() < 0.) {
1470  incoming_particles_[1], __func__);
1471  } else {
1472  return sig_el;
1473  }
1474 }
1475 
1476 CollisionBranchList CrossSections::two_to_one(
1477  const CharmRescattering charm_rescattering) const {
1478  CollisionBranchList resonance_process_list;
1479  const ParticleType& type_particle_a = incoming_particles_[0].type();
1480  const ParticleType& type_particle_b = incoming_particles_[1].type();
1481 
1482  if (charm_rescattering == CharmRescattering::T_Matrix) {
1483  const PdgCode& pdg_a = type_particle_a.pdgcode();
1484  const PdgCode& pdg_b = type_particle_b.pdgcode();
1485  const bool Dmeson_present = pdg_a.is_Dmeson() || pdg_b.is_Dmeson();
1486  const bool Dstar_present = pdg_a.is_Dstar2007() || pdg_b.is_Dstar2007();
1487  const bool light_meson_present = pdg_a.is_pion() || pdg_b.is_pion() ||
1488  pdg_a.is_eta() || pdg_b.is_eta() ||
1489  pdg_a.is_kaon() || pdg_b.is_kaon();
1490  const bool nucleon_or_Delta_present = pdg_a.is_nucleon() ||
1491  pdg_b.is_nucleon() ||
1492  pdg_a.is_Delta() || pdg_b.is_Delta();
1493  if ((Dmeson_present && (light_meson_present || nucleon_or_Delta_present)) ||
1494  (Dstar_present && light_meson_present)) {
1495  return resonance_process_list;
1496  }
1497  }
1498 
1499  const double m1 = incoming_particles_[0].effective_mass();
1500  const double m2 = incoming_particles_[1].effective_mass();
1501  const double p_cm_sqr = pCM_sqr(sqrt_s_, m1, m2);
1502 
1503  ParticleTypePtrList possible_resonances =
1504  list_possible_resonances(&type_particle_a, &type_particle_b);
1505 
1506  // Find all the possible resonances
1507  for (const ParticleTypePtr type_resonance : possible_resonances) {
1508  double resonance_xsection = formation(*type_resonance, p_cm_sqr);
1509 
1510  // If cross section is non-negligible, add resonance to the list
1511  if (resonance_xsection > really_small) {
1512  resonance_process_list.push_back(std::make_unique<CollisionBranch>(
1513  *type_resonance, resonance_xsection, ProcessType::TwoToOne));
1514  logg[LCrossSections].debug("Found resonance: ", *type_resonance);
1515  logg[LCrossSections].debug(type_particle_a.name(), type_particle_b.name(),
1516  "->", type_resonance->name(),
1517  " at sqrt(s)[GeV] = ", sqrt_s_,
1518  " with xs[mb] = ", resonance_xsection);
1519  }
1520  }
1521  return resonance_process_list;
1522 }
1523 
1524 double CrossSections::formation(const ParticleType& type_resonance,
1525  double cm_momentum_sqr) const {
1526  const ParticleType& type_particle_a = incoming_particles_[0].type();
1527  const ParticleType& type_particle_b = incoming_particles_[1].type();
1528 
1529  // Calculate partial in-width.
1530  const double partial_width = type_resonance.get_partial_in_width(
1532  if (partial_width <= 0.) {
1533  return 0.;
1534  }
1535 
1536  assert(type_resonance.charge() ==
1537  type_particle_a.charge() + type_particle_b.charge());
1538  assert(type_resonance.baryon_number() ==
1539  type_particle_a.baryon_number() + type_particle_b.baryon_number());
1540 
1541  const double spinfactor =
1542  static_cast<double>(type_resonance.spin() + 1) /
1543  ((type_particle_a.spin() + 1) * (type_particle_b.spin() + 1));
1544  const int sym_factor =
1545  (type_particle_a.pdgcode() == type_particle_b.pdgcode()) ? 2 : 1;
1546  return spinfactor * sym_factor * 2. * M_PI * M_PI / cm_momentum_sqr *
1547  type_resonance.full_spectral_function(sqrt_s_) * partial_width *
1548  hbarc * hbarc / fm2_mb;
1549 }
1550 
1551 CollisionBranchList CrossSections::two_to_two(
1552  const ReactionsBitSet& included_2to2, const double KN_offset,
1553  const CharmRescattering charm_rescattering) const {
1554  CollisionBranchList process_list;
1555  const ParticleData& data_a = incoming_particles_[0];
1556  const ParticleData& data_b = incoming_particles_[1];
1557  const ParticleType& type_a = data_a.type();
1558  const ParticleType& type_b = data_b.type();
1559  const auto& pdg_a = data_a.pdgcode();
1560  const auto& pdg_b = data_b.pdgcode();
1561 
1562  if (data_a.is_baryon() && data_b.is_baryon()) {
1563  if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
1564  pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
1565  // Nucleon Nucleon Scattering
1566  process_list = nn_xx(included_2to2);
1567  } else {
1568  // Baryon Baryon Scattering
1569  process_list = bb_xx_except_nn(included_2to2);
1570  }
1571  } else if ((type_a.is_baryon() && type_b.is_meson()) ||
1572  (type_a.is_meson() && type_b.is_baryon())) {
1573  // Baryon Meson Scattering
1574  if ((pdg_a.is_nucleon() && pdg_b.is_kaon()) ||
1575  (pdg_b.is_nucleon() && pdg_a.is_kaon())) {
1576  // Nucleon Kaon Scattering
1577  process_list = nk_xx(included_2to2, KN_offset);
1578  } else if ((pdg_a.is_hyperon() && pdg_b.is_pion()) ||
1579  (pdg_b.is_hyperon() && pdg_a.is_pion())) {
1580  // Hyperon Pion Scattering
1581  process_list = ypi_xx(included_2to2);
1582  } else if ((pdg_a.is_Delta() && pdg_b.is_kaon()) ||
1583  (pdg_b.is_Delta() && pdg_a.is_kaon())) {
1584  // Delta Kaon Scattering
1585  process_list = deltak_xx(included_2to2);
1586  } else if ((pdg_a.is_nucleon() && pdg_b.is_Dmeson()) ||
1587  (pdg_b.is_nucleon() && pdg_a.is_Dmeson())) {
1588  // Nucleon D meson Scattering
1589  process_list = DN_xx(included_2to2, charm_rescattering);
1590  } else if ((pdg_a.is_Delta() && pdg_b.is_Dmeson()) ||
1591  (pdg_b.is_Delta() && pdg_a.is_Dmeson())) {
1592  // Delta D meson Scattering
1593  process_list = DDelta_xx(included_2to2, charm_rescattering);
1594  }
1595  } else if (type_a.is_meson() && type_b.is_meson()) {
1596  if ((pdg_a.is_Dmeson() || pdg_b.is_Dmeson()) ||
1597  (pdg_a.is_Dstar2007() || pdg_b.is_Dstar2007())) {
1598  if (pdg_a.is_pion() || pdg_b.is_pion()) {
1599  // D or D* - Pion Scattering
1600  process_list = Dpi_and_Dstarpi_xx(included_2to2, charm_rescattering);
1601  } else if (pdg_a.is_eta() || pdg_b.is_eta()) {
1602  // D or D* - Eta Scattering (inelastic) not existent in T-matrix method
1603  return process_list;
1604  } else if (pdg_a.is_kaon() || pdg_b.is_kaon()) {
1605  // D or D* - Kaon Scattering
1606  process_list = DK_and_DstarK_xx(included_2to2, charm_rescattering);
1607  }
1608  }
1609  } else if (type_a.is_nucleus() || type_b.is_nucleus()) {
1610  if ((type_a.is_nucleon() && type_b.is_nucleus()) ||
1611  (type_b.is_nucleon() && type_a.is_nucleus())) {
1612  // Nucleon Deuteron and Nucleon d' Scattering
1613  process_list = dn_xx(included_2to2);
1614  } else if (((type_a.is_deuteron() || type_a.is_dprime()) &&
1615  pdg_b.is_pion()) ||
1616  ((type_b.is_deuteron() || type_b.is_dprime()) &&
1617  pdg_a.is_pion())) {
1618  // Pion Deuteron and Pion d' Scattering
1619  process_list = dpi_xx(included_2to2);
1620  }
1621  }
1622  return process_list;
1623 }
1624 
1625 CollisionBranchList CrossSections::two_to_three() const {
1626  CollisionBranchList process_list;
1627  const ParticleType& type_a = incoming_particles_[0].type();
1628  const ParticleType& type_b = incoming_particles_[1].type();
1629 
1630  if ((type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
1631  (type_b.is_deuteron() && type_a.pdgcode().is_pion())) {
1632  const ParticleType& type_pi = type_a.pdgcode().is_pion() ? type_a : type_b;
1633  const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
1634 
1635  if (type_nucleus.baryon_number() > 0) {
1636  // πd → πpn
1637  const auto& type_p = ParticleType::find(pdg::p);
1638  const auto& type_n = ParticleType::find(pdg::n);
1639 
1640  process_list.push_back(std::make_unique<CollisionBranch>(
1641  type_pi, type_p, type_n, two_to_three_xs(type_a, type_b, sqrt_s_),
1643  } else {
1644  // πd̅ → πp̅n̅
1645  const auto& type_anti_p = ParticleType::find(-pdg::p);
1646  const auto& type_anti_n = ParticleType::find(-pdg::n);
1647 
1648  process_list.push_back(std::make_unique<CollisionBranch>(
1649  type_pi, type_anti_p, type_anti_n,
1650  two_to_three_xs(type_a, type_b, sqrt_s_), ProcessType::TwoToThree));
1651  }
1652  }
1653 
1654  if ((type_a.is_nucleon() && type_b.is_deuteron()) ||
1655  (type_b.is_nucleon() && type_a.is_deuteron())) {
1656  const ParticleType& type_N = type_a.is_nucleon() ? type_a : type_b;
1657  const ParticleType& type_nucleus = type_a.is_deuteron() ? type_a : type_b;
1658 
1659  if (type_nucleus.baryon_number() > 0) {
1660  // Nd → Nnp, N̅d → N̅np
1661  const auto& type_p = ParticleType::find(pdg::p);
1662  const auto& type_n = ParticleType::find(pdg::n);
1663 
1664  process_list.push_back(std::make_unique<CollisionBranch>(
1665  type_N, type_p, type_n, two_to_three_xs(type_a, type_b, sqrt_s_),
1667  } else {
1668  // Nd̅ → Np̅n̅, N̅d̅ → N̅p̅n̅
1669  const auto& type_anti_p = ParticleType::find(-pdg::p);
1670  const auto& type_anti_n = ParticleType::find(-pdg::n);
1671 
1672  process_list.push_back(std::make_unique<CollisionBranch>(
1673  type_N, type_anti_p, type_anti_n,
1674  two_to_three_xs(type_a, type_b, sqrt_s_), ProcessType::TwoToThree));
1675  }
1676  }
1677  return process_list;
1678 }
1679 
1680 CollisionBranchList CrossSections::two_to_four() const {
1681  CollisionBranchList process_list;
1682  ParticleTypePtr type_nucleus = &(incoming_particles_[0].type());
1683  ParticleTypePtr type_catalyzer = &(incoming_particles_[1].type());
1684  if (!type_nucleus->is_nucleus()) {
1685  type_nucleus = &(incoming_particles_[1].type());
1686  type_catalyzer = &(incoming_particles_[0].type());
1687  }
1688 
1689  if (type_nucleus->is_nucleus() &&
1690  std::abs(type_nucleus->baryon_number()) == 3 &&
1691  (type_catalyzer->is_pion() || type_catalyzer->is_nucleon())) {
1694  const ParticleTypePtr type_anti_p = ParticleType::try_find(-pdg::p);
1695  const ParticleTypePtr type_anti_n = ParticleType::try_find(-pdg::n);
1697  const ParticleTypePtr type_anti_la = ParticleType::try_find(-pdg::Lambda);
1698 
1699  // Find nucleus components
1700  ParticleTypePtrList components;
1701  components.reserve(3);
1702  const PdgCode nucleus_pdg = type_nucleus->pdgcode();
1703  for (int i = 0; i < nucleus_pdg.nucleus_p(); i++) {
1704  components.push_back(type_p);
1705  }
1706  for (int i = 0; i < nucleus_pdg.nucleus_n(); i++) {
1707  components.push_back(type_n);
1708  }
1709  for (int i = 0; i < nucleus_pdg.nucleus_ap(); i++) {
1710  components.push_back(type_anti_p);
1711  }
1712  for (int i = 0; i < nucleus_pdg.nucleus_an(); i++) {
1713  components.push_back(type_anti_n);
1714  }
1715  for (int i = 0; i < nucleus_pdg.nucleus_La(); i++) {
1716  components.push_back(type_la);
1717  }
1718  for (int i = 0; i < nucleus_pdg.nucleus_aLa(); i++) {
1719  components.push_back(type_anti_la);
1720  }
1721  if (sqrt_s_ > type_catalyzer->mass() + components[0]->mass() +
1722  components[1]->mass() + components[2]->mass()) {
1723  process_list.push_back(std::make_unique<CollisionBranch>(
1724  *type_catalyzer, *(components[0]), *(components[1]), *(components[2]),
1725  two_to_four_xs(*type_nucleus, *type_catalyzer, sqrt_s_),
1727  }
1728  }
1729  return process_list;
1730 }
1731 
1733  const ParticleType& type_b,
1734  double sqrts) {
1735  double xs = 0.0;
1736  ParticleTypePtr type_nucleus = &type_a, type_catalyzer = &type_b;
1737  if (!type_nucleus->is_nucleus()) {
1738  type_nucleus = &type_b;
1739  type_catalyzer = &type_a;
1740  }
1741 
1742  bool nonzero_xs = type_nucleus->is_nucleus() &&
1743  (type_catalyzer->is_pion() || type_catalyzer->is_nucleon());
1744  if (!nonzero_xs) {
1745  return 0.0;
1746  }
1747 
1748  const double md = type_nucleus->mass(), mcat = type_catalyzer->mass();
1749  const double Tkin = (sqrts * sqrts - (md + mcat) * (md + mcat)) / (2.0 * md);
1750 
1751  // Should normally never happen, but may be a useful safeguard
1752  if (Tkin <= 0.0) {
1753  return 0.0;
1754  }
1755 
1756  if (type_catalyzer->is_pion()) {
1757  xs = deuteron_pion_inelastic(Tkin);
1758  } else if (type_catalyzer->is_nucleon()) {
1759  if (type_nucleus->pdgcode().antiparticle_sign() ==
1760  type_catalyzer->pdgcode().antiparticle_sign()) {
1761  // Nd and N̅d̅
1762  xs = deuteron_nucleon_inelastic(Tkin);
1763  } else {
1764  // N̅d and Nd̅
1765  xs = deuteron_antinucleon_inelastic(Tkin);
1766  }
1767  }
1768  return xs;
1769 }
1770 
1772  const ParticleType& type_b, double sqrts) {
1773  double xs = 0.0;
1774  ParticleTypePtr type_nucleus = &type_a, type_catalyzer = &type_b;
1775  if (!type_nucleus->is_nucleus()) {
1776  type_nucleus = &type_b;
1777  type_catalyzer = &type_a;
1778  }
1779  bool nonzero_xs = type_nucleus->is_nucleus() &&
1780  (type_catalyzer->is_pion() || type_catalyzer->is_nucleon());
1781  if (!nonzero_xs) {
1782  return 0.0;
1783  }
1784 
1785  const double mA = type_nucleus->mass(), mcat = type_catalyzer->mass();
1786  const double Tkin = (sqrts * sqrts - (mA + mcat) * (mA + mcat)) / (2.0 * mA);
1787  const int A = type_nucleus->pdgcode().nucleus_A();
1788  // Should normally never happen, but may be a useful safeguard
1789  if (A != 3 || Tkin <= 0.0) {
1790  return 0.0;
1791  }
1792 
1793  if (type_catalyzer->is_pion()) {
1794  xs = A / 2. * deuteron_pion_inelastic(Tkin);
1795  } else if (type_catalyzer->is_nucleon()) {
1796  if (type_nucleus->pdgcode().antiparticle_sign() ==
1797  type_catalyzer->pdgcode().antiparticle_sign()) {
1798  // N + A, anti-N + anti-A
1799  xs = A / 2. * deuteron_nucleon_inelastic(Tkin);
1800  } else {
1801  // N̅ + A and N + anti-A
1802  xs = A / 2. * deuteron_antinucleon_inelastic(Tkin);
1803  }
1804  }
1805  return xs;
1806 }
1807 
1809  const ReactionsBitSet& included_2to2) const {
1810  CollisionBranchList process_list;
1811  const ParticleType& type_a = incoming_particles_[0].type();
1812  const ParticleType& type_b = incoming_particles_[1].type();
1813 
1814  bool same_sign = type_a.antiparticle_sign() == type_b.antiparticle_sign();
1815  bool any_nucleus = type_a.is_nucleus() || type_b.is_nucleus();
1816  if (!same_sign && !any_nucleus) {
1817  return process_list;
1818  }
1819  bool anti_particles = type_a.antiparticle_sign() == -1;
1820  if (type_a.is_nucleon() || type_b.is_nucleon()) {
1821  // N R → N N, N̅ R → N̅ N̅
1822  if (included_2to2[IncludedReactions::NN_to_NR] == 1) {
1823  process_list = bar_bar_to_nuc_nuc(anti_particles);
1824  }
1825  } else if (type_a.is_Delta() || type_b.is_Delta()) {
1826  // Δ R → N N, Δ̅ R → N̅ N̅
1827  if (included_2to2[IncludedReactions::NN_to_DR] == 1) {
1828  process_list = bar_bar_to_nuc_nuc(anti_particles);
1829  }
1830  }
1831 
1832  return process_list;
1833 }
1834 
1835 CollisionBranchList CrossSections::nn_xx(
1836  const ReactionsBitSet& included_2to2) const {
1837  CollisionBranchList process_list, channel_list;
1838 
1839  const double sqrts = sqrt_s_;
1840 
1841  /* Find whether colliding particles are nucleons or anti-nucleons;
1842  * adjust lists of produced particles. */
1843  bool both_antinucleons =
1844  (incoming_particles_[0].type().antiparticle_sign() == -1) &&
1845  (incoming_particles_[1].type().antiparticle_sign() == -1);
1846  const ParticleTypePtrList& nuc_or_anti_nuc =
1847  both_antinucleons ? ParticleType::list_anti_nucleons()
1849  const ParticleTypePtrList& delta_or_anti_delta =
1850  both_antinucleons ? ParticleType::list_anti_Deltas()
1852  // Find N N → N R channels.
1853  if (included_2to2[IncludedReactions::NN_to_NR] == 1) {
1854  channel_list = find_nn_xsection_from_type(
1855  ParticleType::list_baryon_resonances(), nuc_or_anti_nuc,
1856  [&sqrts](const ParticleType& type_res_1, const ParticleType&) {
1857  return type_res_1.iso_multiplet()->get_integral_NR(sqrts);
1858  });
1859  process_list.reserve(process_list.size() + channel_list.size());
1860  std::move(channel_list.begin(), channel_list.end(),
1861  std::inserter(process_list, process_list.end()));
1862  channel_list.clear();
1863  }
1864 
1865  // Find N N → Δ R channels.
1866  if (included_2to2[IncludedReactions::NN_to_DR] == 1) {
1867  channel_list = find_nn_xsection_from_type(
1868  ParticleType::list_baryon_resonances(), delta_or_anti_delta,
1869  [&sqrts](const ParticleType& type_res_1,
1870  const ParticleType& type_res_2) {
1871  return type_res_1.iso_multiplet()->get_integral_RR(
1872  type_res_2.iso_multiplet(), sqrts);
1873  });
1874  process_list.reserve(process_list.size() + channel_list.size());
1875  std::move(channel_list.begin(), channel_list.end(),
1876  std::inserter(process_list, process_list.end()));
1877  channel_list.clear();
1878  }
1879 
1880  // Find N N → dπ and N̅ N̅→ d̅π channels.
1886  // Make sure all the necessary particle types are found
1887  if (deuteron && antideutron && pim && pi0 && pip &&
1888  included_2to2[IncludedReactions::PiDeuteron_to_NN] == 1) {
1889  const ParticleTypePtrList deutron_list = {deuteron};
1890  const ParticleTypePtrList antideutron_list = {antideutron};
1891  const ParticleTypePtrList pion_list = {pim, pi0, pip};
1892  channel_list = find_nn_xsection_from_type(
1893  (both_antinucleons ? antideutron_list : deutron_list), pion_list,
1894  [&sqrts](const ParticleType& type_res_1,
1895  const ParticleType& type_res_2) {
1896  return pCM(sqrts, type_res_1.mass(), type_res_2.mass());
1897  });
1898  process_list.reserve(process_list.size() + channel_list.size());
1899  std::move(channel_list.begin(), channel_list.end(),
1900  std::inserter(process_list, process_list.end()));
1901  channel_list.clear();
1902  }
1903 
1904  return process_list;
1905 }
1906 
1907 CollisionBranchList CrossSections::nk_xx(const ReactionsBitSet& included_2to2,
1908  const double KN_offset) const {
1909  const ParticleType& a = incoming_particles_[0].type();
1910  const ParticleType& b = incoming_particles_[1].type();
1911  const ParticleType& type_nucleon = a.pdgcode().is_nucleon() ? a : b;
1912  const ParticleType& type_kaon = a.pdgcode().is_nucleon() ? b : a;
1913 
1914  const auto pdg_nucleon = type_nucleon.pdgcode().code();
1915  const auto pdg_kaon = type_kaon.pdgcode().code();
1916 
1917  const double s = sqrt_s_ * sqrt_s_;
1918 
1919  // Some variable declarations for frequently used quantities
1920  const auto sigma_kplusp = kplusp_inelastic_background(s);
1921  const auto sigma_kplusn = kplusn_inelastic_background(s);
1922 
1923  /* At high energy, the parametrization we use diverges from experimental
1924  * data. This cutoff represents the point where the AQM cross section
1925  * becomes smaller than this parametrization, so we cut it here, and fully
1926  * switch to AQM beyond this point. */
1927  const double KN_to_KDelta_cutoff = KN_offset +
1928  incoming_particles_[0].pole_mass() +
1929  incoming_particles_[1].pole_mass();
1930 
1931  bool incl_KN_to_KN = included_2to2[IncludedReactions::KN_to_KN] == 1;
1932  bool incl_KN_to_KDelta =
1933  included_2to2[IncludedReactions::KN_to_KDelta] == 1 &&
1934  sqrt_s_ < KN_to_KDelta_cutoff;
1935  bool incl_Strangeness_exchange =
1936  included_2to2[IncludedReactions::Strangeness_exchange] == 1;
1937 
1938  CollisionBranchList process_list;
1939  switch (pdg_kaon) {
1940  case pdg::K_m: {
1941  /* All inelastic K- N channels here are strangeness exchange, plus one
1942  * charge exchange. */
1943  switch (pdg_nucleon) {
1944  case pdg::p: {
1945  if (incl_Strangeness_exchange) {
1946  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1947  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
1948  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
1949  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
1950  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
1951  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
1952  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
1953  add_channel(
1954  process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
1955  sqrt_s_, type_pi_m, type_Sigma_p);
1956  add_channel(
1957  process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
1958  sqrt_s_, type_pi_p, type_Sigma_m);
1959  add_channel(
1960  process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
1961  sqrt_s_, type_pi_z, type_Sigma_z);
1962  add_channel(
1963  process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
1964  sqrt_s_, type_pi_z, type_Lambda);
1965  }
1966  if (incl_KN_to_KN) {
1967  const auto& type_n = ParticleType::find(pdg::n);
1968  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1969  add_channel(
1970  process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
1971  type_Kbar_z, type_n);
1972  }
1973  break;
1974  }
1975  case pdg::n: {
1976  if (incl_Strangeness_exchange) {
1977  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
1978  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
1979  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
1980  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
1981  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
1982  add_channel(
1983  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1984  sqrt_s_, type_pi_m, type_Sigma_z);
1985  add_channel(
1986  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
1987  sqrt_s_, type_pi_z, type_Sigma_m);
1988  add_channel(
1989  process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
1990  sqrt_s_, type_pi_m, type_Lambda);
1991  }
1992  break;
1993  }
1994  case -pdg::p: {
1995  if (incl_KN_to_KDelta) {
1996  const auto& type_K_m = ParticleType::find(pdg::K_m);
1997  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
1998  const auto& type_Delta_pp_bar = ParticleType::find(-pdg::Delta_pp);
1999  const auto& type_Delta_p_bar = ParticleType::find(-pdg::Delta_p);
2000  add_channel(
2001  process_list,
2002  [&] {
2003  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
2004  type_nucleon, type_kaon,
2005  type_Kbar_z, type_Delta_pp_bar);
2006  },
2007  sqrt_s_, type_Kbar_z, type_Delta_pp_bar);
2008  add_channel(
2009  process_list,
2010  [&] {
2011  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
2012  type_nucleon, type_kaon, type_K_m,
2013  type_Delta_p_bar);
2014  },
2015  sqrt_s_, type_K_m, type_Delta_p_bar);
2016  }
2017  break;
2018  }
2019  case -pdg::n: {
2020  if (incl_KN_to_KDelta) {
2021  const auto& type_K_m = ParticleType::find(pdg::K_m);
2022  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2023  const auto& type_Delta_p_bar = ParticleType::find(-pdg::Delta_p);
2024  const auto& type_Delta_z_bar = ParticleType::find(-pdg::Delta_z);
2025  add_channel(
2026  process_list,
2027  [&] {
2028  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
2029  type_nucleon, type_kaon,
2030  type_Kbar_z, type_Delta_p_bar);
2031  },
2032  sqrt_s_, type_Kbar_z, type_Delta_p_bar);
2033  add_channel(
2034  process_list,
2035  [&] {
2036  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
2037  type_nucleon, type_kaon, type_K_m,
2038  type_Delta_z_bar);
2039  },
2040  sqrt_s_, type_K_m, type_Delta_z_bar);
2041  }
2042  if (incl_KN_to_KN) {
2043  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2044  const auto& type_p_bar = ParticleType::find(-pdg::p);
2045  add_channel(
2046  process_list, [&] { return kplusn_k0p(s); }, sqrt_s_,
2047  type_Kbar_z, type_p_bar);
2048  }
2049  break;
2050  }
2051  }
2052  break;
2053  }
2054  case pdg::K_p: {
2055  /* All inelastic channels are K+ N -> K Delta -> K pi N, with identical
2056  * cross section, weighted by the isospin factor. */
2057  switch (pdg_nucleon) {
2058  case pdg::p: {
2059  if (incl_KN_to_KDelta) {
2060  const auto& type_K_p = ParticleType::find(pdg::K_p);
2061  const auto& type_K_z = ParticleType::find(pdg::K_z);
2062  const auto& type_Delta_pp = ParticleType::find(pdg::Delta_pp);
2063  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
2064  add_channel(
2065  process_list,
2066  [&] {
2067  return sigma_kplusp *
2068  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
2069  type_K_z, type_Delta_pp);
2070  },
2071  sqrt_s_, type_K_z, type_Delta_pp);
2072  add_channel(
2073  process_list,
2074  [&] {
2075  return sigma_kplusp *
2076  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
2077  type_K_p, type_Delta_p);
2078  },
2079  sqrt_s_, type_K_p, type_Delta_p);
2080  }
2081  break;
2082  }
2083  case pdg::n: {
2084  if (incl_KN_to_KDelta) {
2085  const auto& type_K_p = ParticleType::find(pdg::K_p);
2086  const auto& type_K_z = ParticleType::find(pdg::K_z);
2087  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
2088  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
2089  add_channel(
2090  process_list,
2091  [&] {
2092  return sigma_kplusn *
2093  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
2094  type_K_z, type_Delta_p);
2095  },
2096  sqrt_s_, type_K_z, type_Delta_p);
2097  add_channel(
2098  process_list,
2099  [&] {
2100  return sigma_kplusn *
2101  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
2102  type_K_p, type_Delta_z);
2103  },
2104  sqrt_s_, type_K_p, type_Delta_z);
2105  }
2106  if (incl_KN_to_KN) {
2107  const auto& type_K_z = ParticleType::find(pdg::K_z);
2108  const auto& type_p = ParticleType::find(pdg::p);
2109  add_channel(
2110  process_list, [&] { return kplusn_k0p(s); }, sqrt_s_, type_K_z,
2111  type_p);
2112  }
2113  break;
2114  }
2115  case -pdg::p: {
2116  if (incl_Strangeness_exchange) {
2117  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
2118  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
2119  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
2120  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
2121  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
2122  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
2123  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
2124  add_channel(
2125  process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
2126  sqrt_s_, type_pi_p, type_Sigma_p_bar);
2127  add_channel(
2128  process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
2129  sqrt_s_, type_pi_m, type_Sigma_m_bar);
2130  add_channel(
2131  process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
2132  sqrt_s_, type_pi_z, type_Sigma_z_bar);
2133  add_channel(
2134  process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
2135  sqrt_s_, type_pi_z, type_Lambda_bar);
2136  }
2137  if (incl_KN_to_KN) {
2138  const auto& type_n_bar = ParticleType::find(-pdg::n);
2139  const auto& type_K_z = ParticleType::find(pdg::K_z);
2140  add_channel(
2141  process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
2142  type_K_z, type_n_bar);
2143  }
2144  break;
2145  }
2146  case -pdg::n: {
2147  if (incl_Strangeness_exchange) {
2148  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
2149  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
2150  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
2151  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
2152  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
2153  add_channel(
2154  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
2155  sqrt_s_, type_pi_p, type_Sigma_z_bar);
2156  add_channel(
2157  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
2158  sqrt_s_, type_pi_z, type_Sigma_m_bar);
2159  add_channel(
2160  process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
2161  sqrt_s_, type_pi_p, type_Lambda_bar);
2162  }
2163  break;
2164  }
2165  }
2166  break;
2167  }
2168  case pdg::K_z: {
2169  /* K+ and K0 have the same mass and spin, so their cross sections are
2170  * assumed to only differ in isospin factors. For the initial state, we
2171  * assume that K0 p is equivalent to K+ n and K0 n equivalent to K+ p,
2172  * like for the elastic background. */
2173  switch (pdg_nucleon) {
2174  case pdg::p: {
2175  if (incl_KN_to_KDelta) {
2176  const auto& type_K_p = ParticleType::find(pdg::K_p);
2177  const auto& type_K_z = ParticleType::find(pdg::K_z);
2178  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
2179  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
2180  add_channel(
2181  process_list,
2182  [&] {
2183  return sigma_kplusn *
2184  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
2185  type_K_z, type_Delta_p);
2186  },
2187  sqrt_s_, type_K_z, type_Delta_p);
2188  add_channel(
2189  process_list,
2190  [&] {
2191  return sigma_kplusn *
2192  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
2193  type_K_p, type_Delta_z);
2194  },
2195  sqrt_s_, type_K_p, type_Delta_z);
2196  }
2197  if (incl_KN_to_KN) {
2198  const auto& type_K_p = ParticleType::find(pdg::K_p);
2199  const auto& type_n = ParticleType::find(pdg::n);
2200  add_channel(
2201  process_list,
2202  [&] {
2203  // The isospin factor is 1, see the parametrizations tests.
2204  return kplusn_k0p(s);
2205  },
2206  sqrt_s_, type_K_p, type_n);
2207  }
2208  break;
2209  }
2210  case pdg::n: {
2211  if (incl_KN_to_KDelta) {
2212  const auto& type_K_p = ParticleType::find(pdg::K_p);
2213  const auto& type_K_z = ParticleType::find(pdg::K_z);
2214  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
2215  const auto& type_Delta_m = ParticleType::find(pdg::Delta_m);
2216  add_channel(
2217  process_list,
2218  [&] {
2219  return sigma_kplusp *
2220  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
2221  type_K_z, type_Delta_z);
2222  },
2223  sqrt_s_, type_K_z, type_Delta_z);
2224  add_channel(
2225  process_list,
2226  [&] {
2227  return sigma_kplusp *
2228  kaon_nucleon_ratios.get_ratio(type_nucleon, type_kaon,
2229  type_K_p, type_Delta_m);
2230  },
2231  sqrt_s_, type_K_p, type_Delta_m);
2232  }
2233  break;
2234  }
2235  case -pdg::p: {
2236  if (incl_Strangeness_exchange) {
2237  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
2238  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
2239  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
2240  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
2241  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
2242  add_channel(
2243  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
2244  sqrt_s_, type_pi_m, type_Sigma_z_bar);
2245  add_channel(
2246  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
2247  sqrt_s_, type_pi_z, type_Sigma_p_bar);
2248  add_channel(
2249  process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
2250  sqrt_s_, type_pi_m, type_Lambda_bar);
2251  }
2252  break;
2253  }
2254  case -pdg::n: {
2255  if (incl_Strangeness_exchange) {
2256  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
2257  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
2258  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
2259  const auto& type_Sigma_p_bar = ParticleType::find(-pdg::Sigma_p);
2260  const auto& type_Sigma_m_bar = ParticleType::find(-pdg::Sigma_m);
2261  const auto& type_Sigma_z_bar = ParticleType::find(-pdg::Sigma_z);
2262  const auto& type_Lambda_bar = ParticleType::find(-pdg::Lambda);
2263  add_channel(
2264  process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
2265  sqrt_s_, type_pi_m, type_Sigma_m_bar);
2266  add_channel(
2267  process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
2268  sqrt_s_, type_pi_p, type_Sigma_p_bar);
2269  add_channel(
2270  process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
2271  sqrt_s_, type_pi_z, type_Sigma_z_bar);
2272  add_channel(
2273  process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
2274  sqrt_s_, type_pi_z, type_Lambda_bar);
2275  }
2276  if (incl_KN_to_KN) {
2277  const auto& type_K_p = ParticleType::find(pdg::K_p);
2278  const auto& type_p_bar = ParticleType::find(-pdg::p);
2279  add_channel(
2280  process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
2281  type_K_p, type_p_bar);
2282  }
2283  break;
2284  }
2285  }
2286  break;
2287  }
2288  case pdg::Kbar_z:
2289  switch (pdg_nucleon) {
2290  case pdg::p: {
2291  if (incl_Strangeness_exchange) {
2292  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
2293  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
2294  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
2295  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
2296  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
2297  add_channel(
2298  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
2299  sqrt_s_, type_pi_z, type_Sigma_p);
2300  add_channel(
2301  process_list, [&] { return kminusn_piminussigma0(sqrt_s_); },
2302  sqrt_s_, type_pi_p, type_Sigma_z);
2303  add_channel(
2304  process_list, [&] { return kminusn_piminuslambda(sqrt_s_); },
2305  sqrt_s_, type_pi_p, type_Lambda);
2306  }
2307  break;
2308  }
2309  case pdg::n: {
2310  if (incl_Strangeness_exchange) {
2311  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
2312  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
2313  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
2314  const auto& type_Sigma_p = ParticleType::find(pdg::Sigma_p);
2315  const auto& type_Sigma_m = ParticleType::find(pdg::Sigma_m);
2316  const auto& type_Sigma_z = ParticleType::find(pdg::Sigma_z);
2317  const auto& type_Lambda = ParticleType::find(pdg::Lambda);
2318  add_channel(
2319  process_list, [&] { return kminusp_piminussigmaplus(sqrt_s_); },
2320  sqrt_s_, type_pi_p, type_Sigma_m);
2321  add_channel(
2322  process_list, [&] { return kminusp_piplussigmaminus(sqrt_s_); },
2323  sqrt_s_, type_pi_m, type_Sigma_p);
2324  add_channel(
2325  process_list, [&] { return kminusp_pi0sigma0(sqrt_s_); },
2326  sqrt_s_, type_pi_z, type_Sigma_z);
2327  add_channel(
2328  process_list, [&] { return kminusp_pi0lambda(sqrt_s_); },
2329  sqrt_s_, type_pi_z, type_Lambda);
2330  }
2331  if (incl_KN_to_KN) {
2332  const auto& type_p = ParticleType::find(pdg::p);
2333  const auto& type_K_m = ParticleType::find(pdg::K_m);
2334  add_channel(
2335  process_list, [&] { return kminusp_kbar0n(s); }, sqrt_s_,
2336  type_K_m, type_p);
2337  }
2338  break;
2339  }
2340  case -pdg::p: {
2341  if (incl_KN_to_KDelta) {
2342  const auto& type_K_m = ParticleType::find(pdg::K_m);
2343  const auto& type_Kbar_z = type_kaon;
2344  const auto& type_Delta_bar_m = ParticleType::find(-pdg::Delta_p);
2345  const auto& type_Delta_bar_z = ParticleType::find(-pdg::Delta_z);
2346  add_channel(
2347  process_list,
2348  [&] {
2349  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
2350  type_nucleon, type_kaon,
2351  type_Kbar_z, type_Delta_bar_m);
2352  },
2353  sqrt_s_, type_Kbar_z, type_Delta_bar_m);
2354  add_channel(
2355  process_list,
2356  [&] {
2357  return sigma_kplusn * kaon_nucleon_ratios.get_ratio(
2358  type_nucleon, type_kaon, type_K_m,
2359  type_Delta_bar_z);
2360  },
2361  sqrt_s_, type_K_m, type_Delta_bar_z);
2362  }
2363  if (incl_KN_to_KN) {
2364  const auto& type_K_m = ParticleType::find(pdg::K_m);
2365  const auto& type_n_bar = ParticleType::find(-pdg::n);
2366  add_channel(
2367  process_list,
2368  [&] {
2369  // The isospin factor is 1, see the parametrizations tests.
2370  return kplusn_k0p(s);
2371  },
2372  sqrt_s_, type_K_m, type_n_bar);
2373  }
2374  break;
2375  }
2376  case -pdg::n: {
2377  if (incl_KN_to_KDelta) {
2378  const auto& type_K_m = ParticleType::find(pdg::K_m);
2379  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2380  const auto& type_Delta_z_bar = ParticleType::find(-pdg::Delta_z);
2381  const auto& type_Delta_m_bar = ParticleType::find(-pdg::Delta_m);
2382  add_channel(
2383  process_list,
2384  [&] {
2385  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
2386  type_nucleon, type_kaon,
2387  type_Kbar_z, type_Delta_z_bar);
2388  },
2389  sqrt_s_, type_Kbar_z, type_Delta_z_bar);
2390  add_channel(
2391  process_list,
2392  [&] {
2393  return sigma_kplusp * kaon_nucleon_ratios.get_ratio(
2394  type_nucleon, type_kaon, type_K_m,
2395  type_Delta_m_bar);
2396  },
2397  sqrt_s_, type_K_m, type_Delta_m_bar);
2398  }
2399  break;
2400  }
2401  }
2402  break;
2403  }
2404 
2405  return process_list;
2406 }
2407 
2408 CollisionBranchList CrossSections::deltak_xx(
2409  const ReactionsBitSet& included_2to2) const {
2410  CollisionBranchList process_list;
2411  if (included_2to2[IncludedReactions::KN_to_KDelta] == 0) {
2412  return process_list;
2413  }
2414  const ParticleType& a = incoming_particles_[0].type();
2415  const ParticleType& b = incoming_particles_[1].type();
2416  const ParticleType& type_delta = a.pdgcode().is_Delta() ? a : b;
2417  const ParticleType& type_kaon = a.pdgcode().is_Delta() ? b : a;
2418 
2419  const auto pdg_delta = type_delta.pdgcode().code();
2420  const auto pdg_kaon = type_kaon.pdgcode().code();
2421 
2422  const double s = sqrt_s_ * sqrt_s_;
2423  const double pcm = cm_momentum();
2424  /* The cross sections are determined from the backward reactions via detailed
2425  * balance. The same isospin factors as for the backward reaction are used. */
2426  switch (pack(pdg_delta, pdg_kaon)) {
2427  case pack(pdg::Delta_pp, pdg::K_z):
2428  case pack(pdg::Delta_p, pdg::K_p): {
2429  const auto& type_p = ParticleType::find(pdg::p);
2430  const auto& type_K_p = ParticleType::find(pdg::K_p);
2431  add_channel(
2432  process_list,
2433  [&] {
2434  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
2435  type_kaon, type_p, type_K_p) *
2436  kaon_nucleon_ratios.get_ratio(type_p, type_K_p, type_kaon,
2437  type_delta) *
2439  },
2440  sqrt_s_, type_p, type_K_p);
2441  break;
2442  }
2443  case pack(-pdg::Delta_pp, pdg::Kbar_z):
2444  case pack(-pdg::Delta_p, pdg::K_m): {
2445  const auto& type_p_bar = ParticleType::find(-pdg::p);
2446  const auto& type_K_m = ParticleType::find(pdg::K_m);
2447  add_channel(
2448  process_list,
2449  [&] {
2450  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
2451  type_kaon, type_p_bar, type_K_m) *
2452  kaon_nucleon_ratios.get_ratio(type_p_bar, type_K_m,
2453  type_kaon, type_delta) *
2455  },
2456  sqrt_s_, type_p_bar, type_K_m);
2457  break;
2458  }
2459  case pack(pdg::Delta_p, pdg::K_z):
2460  case pack(pdg::Delta_z, pdg::K_p): {
2461  const auto& type_n = ParticleType::find(pdg::n);
2462  const auto& type_p = ParticleType::find(pdg::p);
2463  const auto& type_K_p = ParticleType::find(pdg::K_p);
2464  const auto& type_K_z = ParticleType::find(pdg::K_z);
2465  add_channel(
2466  process_list,
2467  [&] {
2468  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
2469  type_kaon, type_n, type_K_p) *
2470  kaon_nucleon_ratios.get_ratio(type_n, type_K_p, type_kaon,
2471  type_delta) *
2473  },
2474  sqrt_s_, type_n, type_K_p);
2475 
2476  add_channel(
2477  process_list,
2478  [&] {
2479  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
2480  type_kaon, type_p, type_K_z) *
2481  kaon_nucleon_ratios.get_ratio(type_p, type_K_z, type_kaon,
2482  type_delta) *
2484  },
2485  sqrt_s_, type_p, type_K_z);
2486  break;
2487  }
2488  case pack(-pdg::Delta_p, pdg::Kbar_z):
2489  case pack(-pdg::Delta_z, pdg::K_m): {
2490  const auto& type_n_bar = ParticleType::find(-pdg::n);
2491  const auto& type_p_bar = ParticleType::find(-pdg::p);
2492  const auto& type_K_m = ParticleType::find(pdg::K_m);
2493  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2494  add_channel(
2495  process_list,
2496  [&] {
2497  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
2498  type_kaon, type_n_bar, type_K_m) *
2499  kaon_nucleon_ratios.get_ratio(type_n_bar, type_K_m,
2500  type_kaon, type_delta) *
2502  },
2503  sqrt_s_, type_n_bar, type_K_m);
2504 
2505  add_channel(
2506  process_list,
2507  [&] {
2508  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
2509  type_kaon, type_p_bar,
2510  type_Kbar_z) *
2511  kaon_nucleon_ratios.get_ratio(type_p_bar, type_Kbar_z,
2512  type_kaon, type_delta) *
2514  },
2515  sqrt_s_, type_p_bar, type_Kbar_z);
2516  break;
2517  }
2518  case pack(pdg::Delta_z, pdg::K_z):
2519  case pack(pdg::Delta_m, pdg::K_p): {
2520  const auto& type_n = ParticleType::find(pdg::n);
2521  const auto& type_K_z = ParticleType::find(pdg::K_z);
2522  add_channel(
2523  process_list,
2524  [&] {
2525  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
2526  type_kaon, type_n, type_K_z) *
2527  kaon_nucleon_ratios.get_ratio(type_n, type_K_z, type_kaon,
2528  type_delta) *
2530  },
2531  sqrt_s_, type_n, type_K_z);
2532  break;
2533  }
2534  case pack(-pdg::Delta_z, pdg::Kbar_z):
2535  case pack(-pdg::Delta_m, pdg::K_m): {
2536  const auto& type_n_bar = ParticleType::find(-pdg::n);
2537  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2538  add_channel(
2539  process_list,
2540  [&] {
2541  return detailed_balance_factor_RK(sqrt_s_, pcm, type_delta,
2542  type_kaon, type_n_bar,
2543  type_Kbar_z) *
2544  kaon_nucleon_ratios.get_ratio(type_n_bar, type_Kbar_z,
2545  type_kaon, type_delta) *
2547  },
2548  sqrt_s_, type_n_bar, type_Kbar_z);
2549  break;
2550  }
2551  default:
2552  break;
2553  }
2554 
2555  return process_list;
2556 }
2557 
2558 CollisionBranchList CrossSections::ypi_xx(
2559  const ReactionsBitSet& included_2to2) const {
2560  CollisionBranchList process_list;
2561  if (included_2to2[IncludedReactions::Strangeness_exchange] == 0) {
2562  return process_list;
2563  }
2564  const ParticleType& a = incoming_particles_[0].type();
2565  const ParticleType& b = incoming_particles_[1].type();
2566  const ParticleType& type_hyperon = a.pdgcode().is_hyperon() ? a : b;
2567  const ParticleType& type_pion = a.pdgcode().is_hyperon() ? b : a;
2568 
2569  const auto pdg_hyperon = type_hyperon.pdgcode().code();
2570  const auto pdg_pion = type_pion.pdgcode().code();
2571 
2572  const double s = sqrt_s_ * sqrt_s_;
2573 
2574  switch (pack(pdg_hyperon, pdg_pion)) {
2575  case pack(pdg::Sigma_z, pdg::pi_m): {
2576  const auto& type_n = ParticleType::find(pdg::n);
2577  const auto& type_K_m = ParticleType::find(pdg::K_m);
2578  add_channel(
2579  process_list,
2580  [&] {
2581  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2582  type_n, type_K_m) *
2584  },
2585  sqrt_s_, type_n, type_K_m);
2586  break;
2587  }
2588  case pack(pdg::Sigma_z, pdg::pi_p): {
2589  const auto& type_p = ParticleType::find(pdg::p);
2590  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2591  add_channel(
2592  process_list,
2593  [&] {
2594  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2595  type_p, type_Kbar_z) *
2597  },
2598  sqrt_s_, type_p, type_Kbar_z);
2599  break;
2600  }
2601  case pack(-pdg::Sigma_z, pdg::pi_p): {
2602  const auto& type_n_bar = ParticleType::find(-pdg::n);
2603  const auto& type_K_p = ParticleType::find(pdg::K_p);
2604  add_channel(
2605  process_list,
2606  [&] {
2607  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2608  type_n_bar, type_K_p) *
2610  },
2611  sqrt_s_, type_n_bar, type_K_p);
2612  break;
2613  }
2614  case pack(-pdg::Sigma_z, pdg::pi_m): {
2615  const auto& type_p_bar = ParticleType::find(-pdg::p);
2616  const auto& type_K_z = ParticleType::find(pdg::K_z);
2617  add_channel(
2618  process_list,
2619  [&] {
2620  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2621  type_p_bar, type_K_z) *
2623  },
2624  sqrt_s_, type_p_bar, type_K_z);
2625  break;
2626  }
2627  case pack(pdg::Sigma_m, pdg::pi_z): {
2628  const auto& type_n = ParticleType::find(pdg::n);
2629  const auto& type_K_m = ParticleType::find(pdg::K_m);
2630  add_channel(
2631  process_list,
2632  [&] {
2633  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2634  type_n, type_K_m) *
2636  },
2637  sqrt_s_, type_n, type_K_m);
2638  break;
2639  }
2640  case pack(pdg::Sigma_p, pdg::pi_z): {
2641  const auto& type_p = ParticleType::find(pdg::p);
2642  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2643  add_channel(
2644  process_list,
2645  [&] {
2646  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2647  type_p, type_Kbar_z) *
2649  },
2650  sqrt_s_, type_p, type_Kbar_z);
2651  break;
2652  }
2653  case pack(-pdg::Sigma_m, pdg::pi_z): {
2654  const auto& type_n_bar = ParticleType::find(-pdg::n);
2655  const auto& type_K_p = ParticleType::find(pdg::K_p);
2656  add_channel(
2657  process_list,
2658  [&] {
2659  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2660  type_n_bar, type_K_p) *
2662  },
2663  sqrt_s_, type_n_bar, type_K_p);
2664  break;
2665  }
2666  case pack(-pdg::Sigma_p, pdg::pi_z): {
2667  const auto& type_p_bar = ParticleType::find(-pdg::p);
2668  const auto& type_K_z = ParticleType::find(pdg::K_z);
2669  add_channel(
2670  process_list,
2671  [&] {
2672  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2673  type_p_bar, type_K_z) *
2675  },
2676  sqrt_s_, type_p_bar, type_K_z);
2677  break;
2678  }
2679  case pack(pdg::Lambda, pdg::pi_m): {
2680  const auto& type_n = ParticleType::find(pdg::n);
2681  const auto& type_K_m = ParticleType::find(pdg::K_m);
2682  add_channel(
2683  process_list,
2684  [&] {
2685  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2686  type_n, type_K_m) *
2688  },
2689  sqrt_s_, type_n, type_K_m);
2690  break;
2691  }
2692  case pack(pdg::Lambda, pdg::pi_p): {
2693  const auto& type_p = ParticleType::find(pdg::p);
2694  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2695  add_channel(
2696  process_list,
2697  [&] {
2698  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2699  type_p, type_Kbar_z) *
2701  },
2702  sqrt_s_, type_p, type_Kbar_z);
2703  break;
2704  }
2705  case pack(-pdg::Lambda, pdg::pi_p): {
2706  const auto& type_n_bar = ParticleType::find(-pdg::n);
2707  const auto& type_K_p = ParticleType::find(pdg::K_p);
2708  add_channel(
2709  process_list,
2710  [&] {
2711  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2712  type_n_bar, type_K_p) *
2714  },
2715  sqrt_s_, type_n_bar, type_K_p);
2716  break;
2717  }
2718  case pack(-pdg::Lambda, pdg::pi_m): {
2719  const auto& type_p_bar = ParticleType::find(-pdg::p);
2720  const auto& type_K_z = ParticleType::find(pdg::K_z);
2721  add_channel(
2722  process_list,
2723  [&] {
2724  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2725  type_p_bar, type_K_z) *
2727  },
2728  sqrt_s_, type_p_bar, type_K_z);
2729  break;
2730  }
2731  case pack(pdg::Sigma_z, pdg::pi_z): {
2732  const auto& type_p = ParticleType::find(pdg::p);
2733  const auto& type_n = ParticleType::find(pdg::n);
2734  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2735  const auto& type_K_m = ParticleType::find(pdg::K_m);
2736  add_channel(
2737  process_list,
2738  [&] {
2739  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2740  type_p, type_K_m) *
2742  },
2743  sqrt_s_, type_p, type_K_m);
2744  add_channel(
2745  process_list,
2746  [&] {
2747  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2748  type_n, type_Kbar_z) *
2750  },
2751  sqrt_s_, type_n, type_Kbar_z);
2752  break;
2753  }
2754  case pack(-pdg::Sigma_z, pdg::pi_z): {
2755  const auto& type_p_bar = ParticleType::find(-pdg::p);
2756  const auto& type_n_bar = ParticleType::find(-pdg::n);
2757  const auto& type_K_z = ParticleType::find(pdg::K_z);
2758  const auto& type_K_p = ParticleType::find(pdg::K_p);
2759  add_channel(
2760  process_list,
2761  [&] {
2762  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2763  type_p_bar, type_K_p) *
2765  },
2766  sqrt_s_, type_p_bar, type_K_p);
2767  add_channel(
2768  process_list,
2769  [&] {
2770  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2771  type_n_bar, type_K_z) *
2773  },
2774  sqrt_s_, type_n_bar, type_K_z);
2775  break;
2776  }
2777  case pack(pdg::Sigma_m, pdg::pi_p): {
2778  const auto& type_p = ParticleType::find(pdg::p);
2779  const auto& type_n = ParticleType::find(pdg::n);
2780  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2781  const auto& type_K_m = ParticleType::find(pdg::K_m);
2782  add_channel(
2783  process_list,
2784  [&] {
2785  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2786  type_p, type_K_m) *
2788  },
2789  sqrt_s_, type_p, type_K_m);
2790  add_channel(
2791  process_list,
2792  [&] {
2793  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2794  type_n, type_Kbar_z) *
2796  },
2797  sqrt_s_, type_n, type_Kbar_z);
2798  break;
2799  }
2800  case pack(-pdg::Sigma_m, pdg::pi_m): {
2801  const auto& type_p_bar = ParticleType::find(-pdg::p);
2802  const auto& type_n_bar = ParticleType::find(-pdg::n);
2803  const auto& type_K_z = ParticleType::find(pdg::K_z);
2804  const auto& type_K_p = ParticleType::find(pdg::K_p);
2805  add_channel(
2806  process_list,
2807  [&] {
2808  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2809  type_p_bar, type_K_p) *
2811  },
2812  sqrt_s_, type_p_bar, type_K_p);
2813  add_channel(
2814  process_list,
2815  [&] {
2816  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2817  type_n_bar, type_K_z) *
2819  },
2820  sqrt_s_, type_n_bar, type_K_z);
2821  break;
2822  }
2823  case pack(pdg::Lambda, pdg::pi_z): {
2824  const auto& type_p = ParticleType::find(pdg::p);
2825  const auto& type_n = ParticleType::find(pdg::n);
2826  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2827  const auto& type_K_m = ParticleType::find(pdg::K_m);
2828  add_channel(
2829  process_list,
2830  [&] {
2831  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2832  type_p, type_K_m) *
2834  },
2835  sqrt_s_, type_p, type_K_m);
2836  add_channel(
2837  process_list,
2838  [&] {
2839  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2840  type_n, type_Kbar_z) *
2842  },
2843  sqrt_s_, type_n, type_Kbar_z);
2844  break;
2845  }
2846  case pack(-pdg::Lambda, pdg::pi_z): {
2847  const auto& type_p_bar = ParticleType::find(-pdg::p);
2848  const auto& type_n_bar = ParticleType::find(-pdg::n);
2849  const auto& type_K_z = ParticleType::find(pdg::K_z);
2850  const auto& type_K_p = ParticleType::find(pdg::K_p);
2851  add_channel(
2852  process_list,
2853  [&] {
2854  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2855  type_p_bar, type_K_p) *
2857  },
2858  sqrt_s_, type_p_bar, type_K_p);
2859  add_channel(
2860  process_list,
2861  [&] {
2862  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2863  type_n_bar, type_K_z) *
2865  },
2866  sqrt_s_, type_n_bar, type_K_z);
2867  break;
2868  }
2869  case pack(pdg::Sigma_p, pdg::pi_m): {
2870  const auto& type_p = ParticleType::find(pdg::p);
2871  const auto& type_n = ParticleType::find(pdg::n);
2872  const auto& type_K_m = ParticleType::find(pdg::K_m);
2873  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
2874  add_channel(
2875  process_list,
2876  [&] {
2877  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2878  type_p, type_K_m) *
2880  },
2881  sqrt_s_, type_p, type_K_m);
2882  add_channel(
2883  process_list,
2884  [&] {
2885  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2886  type_n, type_Kbar_z) *
2888  },
2889  sqrt_s_, type_n, type_Kbar_z);
2890  break;
2891  }
2892  case pack(-pdg::Sigma_p, pdg::pi_p): {
2893  const auto& type_p_bar = ParticleType::find(-pdg::p);
2894  const auto& type_n_bar = ParticleType::find(-pdg::n);
2895  const auto& type_K_p = ParticleType::find(pdg::K_p);
2896  const auto& type_K_z = ParticleType::find(pdg::K_z);
2897  add_channel(
2898  process_list,
2899  [&] {
2900  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2901  type_p_bar, type_K_p) *
2903  },
2904  sqrt_s_, type_p_bar, type_K_p);
2905  add_channel(
2906  process_list,
2907  [&] {
2908  return detailed_balance_factor_stable(s, type_hyperon, type_pion,
2909  type_n_bar, type_K_z) *
2911  },
2912  sqrt_s_, type_n_bar, type_K_z);
2913  break;
2914  }
2915  default:
2916  break;
2917  }
2918 
2919  return process_list;
2920 }
2921 
2922 double CrossSections::xs_dpi_dprimepi(const double sqrts, const double cm_mom,
2923  ParticleTypePtr produced_nucleus,
2924  const ParticleType& type_pi) {
2925  const double s = sqrts * sqrts;
2926  // same matrix element for πd and πd̅
2927  const double tmp = sqrts - pion_mass - deuteron_mass;
2928  /**
2929  * Matrix element is fit to match the inelastic pi+ d -> pi+ n p cross-section
2930  * from the Fig. 5 of \iref{Arndt:1994bs}.
2931  */
2932  const double matrix_element =
2933  295.5 + 2.862 / (0.00283735 + pow_int(sqrts - 2.181, 2)) +
2934  0.0672 / pow_int(tmp, 2) - 6.61753 / tmp;
2935 
2936  const double spin_factor =
2937  (produced_nucleus->spin() + 1) * (type_pi.spin() + 1);
2938  /* Isospin factor is always the same, so it is included into the matrix
2939  * element. Symmetry factor is always 1 here. The (hbarc)^2/16 pi factor is
2940  * absorbed into matrix element. */
2941  double xsection = matrix_element * spin_factor / (s * cm_mom);
2942  if (produced_nucleus->is_stable()) {
2943  xsection *= pCM_from_s(s, type_pi.mass(), produced_nucleus->mass());
2944  } else {
2945  const double resonance_integral =
2946  produced_nucleus->iso_multiplet()->get_integral_piR(sqrts);
2947  xsection *= resonance_integral;
2948  logg[LScatterAction].debug("Resonance integral ", resonance_integral,
2949  ", matrix element: ", matrix_element,
2950  ", cm_momentum: ", cm_mom);
2951  }
2952  return xsection;
2953 }
2954 
2955 CollisionBranchList CrossSections::dpi_xx(
2956  const ReactionsBitSet& included_2to2) const {
2957  CollisionBranchList process_list;
2958  const double sqrts = sqrt_s_;
2959  const ParticleType& type_a = incoming_particles_[0].type();
2960  const ParticleType& type_b = incoming_particles_[1].type();
2961 
2962  // pi d -> N N
2963  bool is_pid = (type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
2964  (type_b.is_deuteron() && type_a.pdgcode().is_pion());
2965  if (is_pid && included_2to2[IncludedReactions::PiDeuteron_to_NN] == 1) {
2966  const int baryon_number = type_a.baryon_number() + type_b.baryon_number();
2967  ParticleTypePtrList nuc = (baryon_number > 0)
2970  const double s = sqrt_s_ * sqrt_s_;
2971  for (ParticleTypePtr nuc_a : nuc) {
2972  for (ParticleTypePtr nuc_b : nuc) {
2973  if (type_a.charge() + type_b.charge() !=
2974  nuc_a->charge() + nuc_b->charge()) {
2975  continue;
2976  }
2977  // loop over total isospin
2978  for (const int twoI : I_tot_range(*nuc_a, *nuc_b)) {
2979  const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
2980  type_a, type_b, *nuc_a, *nuc_b, twoI);
2981  // If Clebsch-Gordan coefficient = 0, don't bother with the rest.
2982  if (std::abs(isospin_factor) < really_small) {
2983  continue;
2984  }
2985 
2986  // Calculate matrix element for inverse process.
2987  const double matrix_element =
2988  nn_to_resonance_matrix_element(sqrts, type_a, type_b, twoI);
2989  if (matrix_element <= 0.) {
2990  continue;
2991  }
2992 
2993  const double spin_factor = (nuc_a->spin() + 1) * (nuc_b->spin() + 1);
2994  const int sym_fac_in =
2995  (type_a.iso_multiplet() == type_b.iso_multiplet()) ? 2 : 1;
2996  const int sym_fac_out =
2997  (nuc_a->iso_multiplet() == nuc_b->iso_multiplet()) ? 2 : 1;
2998  double p_cm_final = pCM_from_s(s, nuc_a->mass(), nuc_b->mass());
2999  const double xsection = isospin_factor * spin_factor * sym_fac_in /
3000  sym_fac_out * p_cm_final * matrix_element /
3001  (s * cm_momentum());
3002 
3003  if (xsection > really_small) {
3004  process_list.push_back(std::make_unique<CollisionBranch>(
3005  *nuc_a, *nuc_b, xsection, ProcessType::TwoToTwo));
3006  logg[LScatterAction].debug(type_a.name(), type_b.name(), "->",
3007  nuc_a->name(), nuc_b->name(),
3008  " at sqrts [GeV] = ", sqrts,
3009  " with cs[mb] = ", xsection);
3010  }
3011  }
3012  }
3013  }
3014  }
3015 
3016  // pi d -> pi d' (effectively pi d -> pi p n) AND reverse, pi d' -> pi d
3017  bool is_pid_or_pidprime = ((type_a.is_deuteron() || type_a.is_dprime()) &&
3018  type_b.pdgcode().is_pion()) ||
3019  ((type_b.is_deuteron() || type_b.is_dprime()) &&
3020  type_a.pdgcode().is_pion());
3021  if (is_pid_or_pidprime &&
3022  included_2to2[IncludedReactions::PiDeuteron_to_pidprime] == 1) {
3023  const ParticleType& type_pi = type_a.pdgcode().is_pion() ? type_a : type_b;
3024  const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
3025  ParticleTypePtrList nuclei = ParticleType::list_light_nuclei();
3026  for (ParticleTypePtr produced_nucleus : nuclei) {
3027  // Elastic collisions are treated in a different function
3028  if (produced_nucleus == &type_nucleus ||
3029  produced_nucleus->charge() != type_nucleus.charge() ||
3030  produced_nucleus->baryon_number() != type_nucleus.baryon_number()) {
3031  continue;
3032  }
3033  const double xsection =
3034  xs_dpi_dprimepi(sqrts, cm_momentum(), produced_nucleus, type_pi);
3035  process_list.push_back(std::make_unique<CollisionBranch>(
3036  type_pi, *produced_nucleus, xsection, ProcessType::TwoToTwo));
3037  logg[LScatterAction].debug(type_pi.name(), type_nucleus.name(), "→ ",
3038  type_pi.name(), produced_nucleus->name(),
3039  " at ", sqrts, " GeV, xs[mb] = ", xsection);
3040  }
3041  }
3042  return process_list;
3043 }
3044 
3045 double CrossSections::xs_dn_dprimen(const double sqrts, const double cm_mom,
3046  ParticleTypePtr produced_nucleus,
3047  const ParticleType& type_nucleus,
3048  const ParticleType& type_N) {
3049  const double s = sqrts * sqrts;
3050  double matrix_element = 0.0;
3051  double tmp = sqrts - nucleon_mass - deuteron_mass;
3052  assert(tmp >= 0.0);
3053  if (std::signbit(type_N.baryon_number()) ==
3054  std::signbit(type_nucleus.baryon_number())) {
3055  /**
3056  * Nd → Nd', N̅d̅→ N̅d̅' and reverse:
3057  * Fit to match experimental cross-section Nd -> Nnp from \cite Carlson1973.
3058  */
3059  matrix_element = 79.0474 / std::pow(tmp, 0.7897) + 654.596 * tmp;
3060  } else {
3061  /**
3062  * N̅d → N̅d', Nd̅→ Nd̅' and reverse:
3063  * Fit to roughly match experimental cross-section N̅d -> N̅ np from
3064  * \iref{Bizzarri:1973sp}.
3065  */
3066  matrix_element = 342.572 / std::pow(tmp, 0.6);
3067  }
3068  const double spin_factor =
3069  (produced_nucleus->spin() + 1) * (type_N.spin() + 1);
3070  /* Isospin factor is always the same, so it is included into matrix element
3071  * Symmetry factor is always 1 here. Absorb (hbarc)^2/16 pi factor into matrix
3072  * element. */
3073  double xsection = matrix_element * spin_factor / (s * cm_mom);
3074  if (produced_nucleus->is_stable()) {
3075  assert(!type_nucleus.is_stable());
3076  xsection *= pCM_from_s(s, type_N.mass(), produced_nucleus->mass());
3077  } else {
3078  assert(type_nucleus.is_stable());
3079  const double resonance_integral =
3080  produced_nucleus->iso_multiplet()->get_integral_NR(sqrts);
3081  xsection *= resonance_integral;
3082  }
3083  return xsection;
3084 }
3085 
3086 CollisionBranchList CrossSections::dn_xx(
3087  const ReactionsBitSet& included_2to2) const {
3088  const ParticleType& type_a = incoming_particles_[0].type();
3089  const ParticleType& type_b = incoming_particles_[1].type();
3090  const ParticleType& type_N = type_a.is_nucleon() ? type_a : type_b;
3091  const ParticleType& type_nucleus = type_a.is_nucleus() ? type_a : type_b;
3092  CollisionBranchList process_list;
3093  if (included_2to2[IncludedReactions::NDeuteron_to_Ndprime] == 0) {
3094  return process_list;
3095  }
3096  ParticleTypePtrList nuclei = ParticleType::list_light_nuclei();
3097 
3098  for (ParticleTypePtr produced_nucleus : nuclei) {
3099  // No elastic collisions for now, respect conservation laws
3100  if (produced_nucleus == &type_nucleus ||
3101  produced_nucleus->charge() != type_nucleus.charge() ||
3102  produced_nucleus->baryon_number() != type_nucleus.baryon_number()) {
3103  continue;
3104  }
3105  const double xsection = xs_dn_dprimen(
3106  sqrt_s_, cm_momentum(), produced_nucleus, type_nucleus, type_N);
3107  process_list.push_back(std::make_unique<CollisionBranch>(
3108  type_N, *produced_nucleus, xsection, ProcessType::TwoToTwo));
3109  logg[LScatterAction].debug(type_N.name(), type_nucleus.name(), "→ ",
3110  type_N.name(), produced_nucleus->name(), " at ",
3111  sqrt_s_, " GeV, xs[mb] = ", xsection);
3112  }
3113  return process_list;
3114 }
3115 
3117  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
3118  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
3119  const auto pdg_D =
3120  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
3121  const auto pdg_pion =
3122  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
3123 
3124  double sig_inel = -1.;
3125  switch (pack(pdg_D, pdg_pion)) {
3126  // Checks for D mesons scatterings
3127  case pack(pdg::D_z, pdg::pi_p):
3128  case pack(pdg::Dbar_z, pdg::pi_m): { // Same xsec for charge conjugation.
3129  sig_inel = Dzeropiplus_Dpluspizero(sqrt_s_);
3130  break;
3131  }
3132  case pack(pdg::D_z, pdg::pi_z):
3133  case pack(pdg::Dbar_z, pdg::pi_z): { // Same xsec for charge conjugation.
3134  sig_inel = Dzeropizero_Dpluspiminus(sqrt_s_);
3135  break;
3136  }
3137  case pack(pdg::D_p, pdg::pi_m):
3138  case pack(pdg::D_m, pdg::pi_p): { // Same xsec for charge conjugation.
3139  sig_inel = Dpluspiminus_Dzeropizero(sqrt_s_);
3140  break;
3141  }
3142  case pack(pdg::D_p, pdg::pi_z):
3143  case pack(pdg::D_m, pdg::pi_z): { // Same xsec for charge conjugation.
3144  sig_inel = Dpluspizero_Dzeropiplus(sqrt_s_);
3145  break;
3146  }
3147  // Checks for D* mesons scatterings
3148  case pack(pdg::Dstar_z, pdg::pi_p):
3149  case pack(pdg::Dstarbar_z, pdg::pi_m): { // Same xs for charge conjugation.
3151  break;
3152  }
3153  case pack(pdg::Dstar_z, pdg::pi_z):
3154  case pack(pdg::Dstarbar_z, pdg::pi_z): { // Same xs for charge conjugation.
3156  break;
3157  }
3158  case pack(pdg::Dstar_p, pdg::pi_m):
3159  case pack(pdg::Dstar_m, pdg::pi_p): { // Same xsec for charge conjugation.
3161  break;
3162  }
3163  case pack(pdg::Dstar_p, pdg::pi_z):
3164  case pack(pdg::Dstar_m, pdg::pi_z): { // Same xsec for charge conjugation.
3166  break;
3167  }
3168  case pack(pdg::D_z, pdg::pi_m):
3169  case pack(pdg::Dbar_z, pdg::pi_p):
3170  case pack(pdg::D_p, pdg::pi_p):
3171  case pack(pdg::D_m, pdg::pi_m):
3172  case pack(pdg::Dstar_z, pdg::pi_m):
3174  case pack(pdg::Dstar_p, pdg::pi_p):
3175  case pack(pdg::Dstar_m, pdg::pi_m): {
3176  // These combinations can only scatter elastically.
3177  return 0.;
3178  }
3179  default:
3181  incoming_particles_[1], __func__);
3182  }
3183 
3184  if (sig_inel < 0.) {
3186  incoming_particles_[1], __func__);
3187  } else {
3188  return sig_inel;
3189  }
3190 }
3191 
3193  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
3194  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
3195  const auto pdg_D =
3196  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
3197  const auto pdg_kaon =
3198  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
3199 
3200  double sig_inel = -1.;
3201  switch (pack(pdg_D, pdg_kaon)) {
3202  // Checks for D mesons scatterings
3203  case pack(pdg::D_p, pdg::K_z):
3204  case pack(pdg::D_m, pdg::Kbar_z): { // Same xsec for charge conjugation.
3205  sig_inel = DplusKzero_DzeroKplus(sqrt_s_);
3206  break;
3207  }
3208  case pack(pdg::D_z, pdg::K_p):
3209  case pack(pdg::Dbar_z, pdg::K_m): { // Same xsec for charge conjugation.
3210  sig_inel = DzeroKplus_DplusKzero(sqrt_s_);
3211  break;
3212  }
3213  case pack(pdg::D_p, pdg::K_m):
3214  case pack(pdg::D_m, pdg::K_p): { // Same xsec for charge conjugation.
3215  sig_inel = DplusKminus_DzeroKbarzero(sqrt_s_);
3216  break;
3217  }
3218  case pack(pdg::D_z, pdg::Kbar_z):
3219  case pack(pdg::Dbar_z, pdg::K_z): { // Same xsec for charge conjugation.
3220  sig_inel = DzeroKbarzero_DplusKminus(sqrt_s_);
3221  break;
3222  }
3223  // Checks for D* mesons scatterings
3224  case pack(pdg::Dstar_p, pdg::K_z):
3225  case pack(pdg::Dstar_m, pdg::Kbar_z): { // Same xs for charge conjugation.
3227  break;
3228  }
3229  case pack(pdg::Dstar_z, pdg::K_p):
3230  case pack(pdg::Dstarbar_z, pdg::K_m): { // Same xs for charge conjugation.
3232  break;
3233  }
3234  case pack(pdg::Dstar_p, pdg::K_m):
3235  case pack(pdg::Dstar_m, pdg::K_p): { // Same xsec for charge conjugation.
3237  break;
3238  }
3239  case pack(pdg::Dstar_z, pdg::Kbar_z):
3240  case pack(pdg::Dstarbar_z, pdg::K_z): { // Same xs for charge conjugation.
3242  break;
3243  }
3244  case pack(pdg::D_p, pdg::K_p):
3245  case pack(pdg::D_p, pdg::Kbar_z):
3246  case pack(pdg::D_z, pdg::K_z):
3247  case pack(pdg::D_z, pdg::K_m):
3248  case pack(pdg::D_m, pdg::K_z):
3249  case pack(pdg::D_m, pdg::K_m):
3250  case pack(pdg::Dbar_z, pdg::K_p):
3251  case pack(pdg::Dbar_z, pdg::Kbar_z):
3252  case pack(pdg::Dstar_p, pdg::K_p):
3253  case pack(pdg::Dstar_p, pdg::Kbar_z):
3254  case pack(pdg::Dstar_z, pdg::K_z):
3255  case pack(pdg::Dstar_z, pdg::K_m):
3256  case pack(pdg::Dstar_m, pdg::K_z):
3257  case pack(pdg::Dstar_m, pdg::K_m):
3258  case pack(pdg::Dstarbar_z, pdg::K_p):
3259  case pack(pdg::Dstarbar_z, pdg::Kbar_z): {
3260  // These combinations can only scatter elastically.
3261  return 0.;
3262  break;
3263  }
3264  default:
3266  incoming_particles_[1], __func__);
3267  }
3268 
3269  if (sig_inel < 0.) {
3271  incoming_particles_[1], __func__);
3272  } else {
3273  return sig_inel;
3274  }
3275 }
3276 
3278  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
3279  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
3280  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
3281  const auto pdg_nucleon = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
3282 
3283  double sig_inel = -1.;
3284  switch (pack(pdg_D, pdg_nucleon)) {
3285  case pack(pdg::D_p, pdg::n):
3286  case pack(pdg::D_m, -pdg::n): { // Same xsec for charge conjugation.
3287  sig_inel = Dplusn_Dzerop(sqrt_s_);
3288  break;
3289  }
3290  case pack(pdg::D_z, pdg::p):
3291  case pack(pdg::Dbar_z, -pdg::p): { // Same xsec for charge conjugation.
3292  sig_inel = Dzerop_Dplusn(sqrt_s_);
3293  break;
3294  }
3295  case pack(pdg::D_m, pdg::p):
3296  case pack(pdg::D_p, -pdg::p): { // Same xsec for charge conjugation.
3297  sig_inel = Dminusp_Dbarzeron(sqrt_s_);
3298  break;
3299  }
3300  case pack(pdg::Dbar_z, pdg::n):
3301  case pack(pdg::D_z, -pdg::n): { // Same xsec for charge conjugation.
3302  sig_inel = Dbarzeron_Dminusp(sqrt_s_);
3303  break;
3304  }
3305  case pack(pdg::D_p, -pdg::n):
3306  case pack(pdg::D_p, pdg::p):
3307  case pack(pdg::D_z, -pdg::p):
3308  case pack(pdg::D_z, pdg::n):
3309  case pack(pdg::D_m, pdg::n):
3310  case pack(pdg::D_m, -pdg::p):
3311  case pack(pdg::Dbar_z, pdg::p):
3312  case pack(pdg::Dbar_z, -pdg::n): {
3313  // These combinations can only scatter elastically.
3314  return 0.;
3315  break;
3316  }
3317  default:
3319  incoming_particles_[1], __func__);
3320  }
3321 
3322  if (sig_inel < 0.) {
3324  incoming_particles_[1], __func__);
3325  } else {
3326  return sig_inel;
3327  }
3328 }
3329 
3331  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
3332  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
3333  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
3334  const auto pdg_Delta = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
3335 
3336  double sig_inel = -1.;
3337  switch (pack(pdg_D, pdg_Delta)) {
3338  case pack(pdg::D_p, pdg::Delta_p):
3339  case pack(pdg::D_m, -pdg::Delta_p): { // Same xsec for charge conjugation.
3341  break;
3342  }
3343  case pack(pdg::D_p, pdg::Delta_m):
3344  case pack(pdg::D_m, -pdg::Delta_m): { // Same xsec for charge conjugation.
3346  break;
3347  }
3348  case pack(pdg::D_p, pdg::Delta_z):
3349  case pack(pdg::D_m, -pdg::Delta_z): { // Same xsec for charge conjugation.
3351  break;
3352  }
3353  case pack(pdg::D_z, pdg::Delta_p):
3354  case pack(pdg::Dbar_z, -pdg::Delta_p): { // Same xs for charge conjugation.
3356  break;
3357  }
3358  case pack(pdg::D_z, pdg::Delta_pp):
3359  case pack(pdg::Dbar_z, -pdg::Delta_pp): { // Same xs for charge conjugate.
3361  break;
3362  }
3363  case pack(pdg::D_z, pdg::Delta_z):
3364  case pack(pdg::Dbar_z, -pdg::Delta_z): { // Same xs for charge conjugation.
3366  break;
3367  }
3368  case pack(pdg::D_m, pdg::Delta_p):
3369  case pack(pdg::D_p, -pdg::Delta_p): { // Same xsec for charge conjugation.
3371  break;
3372  }
3373  case pack(pdg::D_m, pdg::Delta_pp):
3374  case pack(pdg::D_p, -pdg::Delta_pp): { // Same xsec for charge conjugation.
3376  break;
3377  }
3378  case pack(pdg::D_m, pdg::Delta_z):
3379  case pack(pdg::D_p, -pdg::Delta_z): { // Same xsec for charge conjugation.
3381  break;
3382  }
3383  case pack(pdg::Dbar_z, pdg::Delta_p):
3384  case pack(pdg::D_z, -pdg::Delta_p): { // Same xsec for charge conjugation.
3386  break;
3387  }
3388  case pack(pdg::Dbar_z, pdg::Delta_m):
3389  case pack(pdg::D_z, -pdg::Delta_m): { // Same xsec for charge conjugation.
3391  break;
3392  }
3393  case pack(pdg::Dbar_z, pdg::Delta_z):
3394  case pack(pdg::D_z, -pdg::Delta_z): { // Same xsec for charge conjugation.
3396  break;
3397  }
3398  case pack(pdg::D_p, pdg::Delta_pp):
3399  case pack(pdg::D_p, -pdg::Delta_m):
3400  case pack(pdg::D_z, -pdg::Delta_pp):
3401  case pack(pdg::D_z, pdg::Delta_m):
3402  case pack(pdg::D_m, -pdg::Delta_pp):
3403  case pack(pdg::D_m, pdg::Delta_m):
3405  case pack(pdg::Dbar_z, -pdg::Delta_m): {
3406  // These combinations can only scatter elastically.
3407  return 0.;
3408  break;
3409  }
3410  default:
3412  incoming_particles_[1], __func__);
3413  }
3414 
3415  if (sig_inel < 0.) {
3417  incoming_particles_[1], __func__);
3418  } else {
3419  return sig_inel;
3420  }
3421 }
3422 
3424  const ReactionsBitSet& included_2to2,
3425  const CharmRescattering charm_rescattering) const {
3426  CollisionBranchList process_list;
3427  if ((included_2to2[IncludedReactions::Charm_T_matrix] == 0) ||
3428  !(charm_rescattering == CharmRescattering::T_Matrix)) {
3429  return process_list;
3430  }
3431  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
3432  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
3433  const auto pdg_D =
3434  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
3435  const auto pdg_pion =
3436  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
3437 
3438  /* Adding the following channels, a check for detailed balance is not needed
3439  * because the parametrization is explicit. */
3440  switch (pack(pdg_D, pdg_pion)) {
3441  // Channels for Dpi scatterings
3442  case pack(pdg::D_z, pdg::pi_p): {
3443  const auto& type_D_p = ParticleType::find(pdg::D_p);
3444  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
3445  // D0pi+ -> D+pi0
3446  add_channel(
3447  process_list, [&] { return Dzeropiplus_Dpluspizero(sqrt_s_); },
3448  sqrt_s_, type_D_p, type_pi_z);
3449  break;
3450  }
3451  case pack(pdg::Dbar_z, pdg::pi_m): {
3452  const auto& type_D_m = ParticleType::find(pdg::D_m);
3453  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
3454  // D0barpi- -> D-pi0 (charge conjugation of D0pi+ -> D+pi0)
3455  add_channel(
3456  process_list, [&] { return Dzeropiplus_Dpluspizero(sqrt_s_); },
3457  sqrt_s_, type_D_m, type_pi_z);
3458  break;
3459  }
3460  case pack(pdg::D_z, pdg::pi_z): {
3461  const auto& type_D_p = ParticleType::find(pdg::D_p);
3462  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
3463  // D0pi0 -> D+pi-
3464  add_channel(
3465  process_list, [&] { return Dzeropizero_Dpluspiminus(sqrt_s_); },
3466  sqrt_s_, type_D_p, type_pi_m);
3467  break;
3468  }
3469  case pack(pdg::Dbar_z, pdg::pi_z): {
3470  const auto& type_D_m = ParticleType::find(pdg::D_m);
3471  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
3472  // Dbar0pi0 -> D-pi+ (charge conjugation of D0pi0 -> D+pi-)
3473  add_channel(
3474  process_list, [&] { return Dzeropizero_Dpluspiminus(sqrt_s_); },
3475  sqrt_s_, type_D_m, type_pi_p);
3476  break;
3477  }
3478  case pack(pdg::D_p, pdg::pi_m): {
3479  const auto& type_D_z = ParticleType::find(pdg::D_z);
3480  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
3481  // D+pi- -> D0pi0
3482  add_channel(
3483  process_list, [&] { return Dpluspiminus_Dzeropizero(sqrt_s_); },
3484  sqrt_s_, type_D_z, type_pi_z);
3485  break;
3486  }
3487  case pack(pdg::D_m, pdg::pi_p): {
3488  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3489  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
3490  // D-pi+ -> Dbar0pi0 (charge conjugation of D+pi- -> D0pi0)
3491  add_channel(
3492  process_list, [&] { return Dpluspiminus_Dzeropizero(sqrt_s_); },
3493  sqrt_s_, type_Dbar_z, type_pi_z);
3494  break;
3495  }
3496  case pack(pdg::D_p, pdg::pi_z): {
3497  const auto& type_D_z = ParticleType::find(pdg::D_z);
3498  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
3499  // D+pi0 -> D0pi+
3500  add_channel(
3501  process_list, [&] { return Dpluspizero_Dzeropiplus(sqrt_s_); },
3502  sqrt_s_, type_D_z, type_pi_p);
3503  break;
3504  }
3505  case pack(pdg::D_m, pdg::pi_z): {
3506  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3507  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
3508  // D-pi0 -> Dbar0pi- (charge conjugation of D+pi0 -> D0pi+)
3509  add_channel(
3510  process_list, [&] { return Dpluspizero_Dzeropiplus(sqrt_s_); },
3511  sqrt_s_, type_Dbar_z, type_pi_m);
3512  break;
3513  }
3514  // Channels for D*pi scatterings
3515  case pack(pdg::Dstar_z, pdg::pi_p): {
3516  const auto& type_Dstar_p = ParticleType::find(pdg::Dstar_p);
3517  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
3518  // D*(2007)0pi+ -> D*(2010)+pi0
3519  add_channel(
3520  process_list,
3522  type_Dstar_p, type_pi_z);
3523  break;
3524  }
3525  case pack(pdg::Dstarbar_z, pdg::pi_m): {
3526  const auto& type_Dstar_m = ParticleType::find(pdg::Dstar_m);
3527  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
3528  /* D*(2007)bar0pi- -> D*(2010)-pi0
3529  * (charge conjugation of D*(2007)0pi+ -> D*(2010)+pi0) */
3530  add_channel(
3531  process_list,
3533  type_Dstar_m, type_pi_z);
3534  break;
3535  }
3536  case pack(pdg::Dstar_z, pdg::pi_z): {
3537  const auto& type_Dstar_p = ParticleType::find(pdg::Dstar_p);
3538  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
3539  // D*(2007)0pi0 -> D*(2010)+pi-
3540  add_channel(
3541  process_list,
3543  type_Dstar_p, type_pi_m);
3544  break;
3545  }
3546  case pack(pdg::Dstarbar_z, pdg::pi_z): {
3547  const auto& type_Dstar_m = ParticleType::find(pdg::Dstar_m);
3548  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
3549  /* D*(2007)bar0pi0 -> D*(2010)-pi+
3550  * (charge conjugation of D*(2007)0pi0 -> D*(2010)+pi-) */
3551  add_channel(
3552  process_list,
3554  type_Dstar_m, type_pi_p);
3555  break;
3556  }
3557  case pack(pdg::Dstar_p, pdg::pi_m): {
3558  const auto& type_Dstar_z = ParticleType::find(pdg::Dstar_z);
3559  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
3560  // D*(2010)+pi- -> D*(2007)0pi0
3561  add_channel(
3562  process_list,
3564  type_Dstar_z, type_pi_z);
3565  break;
3566  }
3567  case pack(pdg::Dstar_m, pdg::pi_p): {
3568  const auto& type_Dstarbar_z = ParticleType::find(pdg::Dstarbar_z);
3569  const auto& type_pi_z = ParticleType::find(pdg::pi_z);
3570  /* D*(2010)-pi+ -> D*(2007)bar0pi0
3571  * (charge conjugation of D*(2010)+pi- -> D*(2007)0pi0) */
3572  add_channel(
3573  process_list,
3575  type_Dstarbar_z, type_pi_z);
3576  break;
3577  }
3578  case pack(pdg::Dstar_p, pdg::pi_z): {
3579  const auto& type_Dstar_z = ParticleType::find(pdg::Dstar_z);
3580  const auto& type_pi_p = ParticleType::find(pdg::pi_p);
3581  // D*(2010)+pi0 -> D*(2007)0pi+
3582  add_channel(
3583  process_list,
3585  type_Dstar_z, type_pi_p);
3586  break;
3587  }
3588  case pack(pdg::Dstar_m, pdg::pi_z): {
3589  const auto& type_Dstarbar_z = ParticleType::find(pdg::Dstarbar_z);
3590  const auto& type_pi_m = ParticleType::find(pdg::pi_m);
3591  /* D*(2010)-pi0 -> D*(2007)bar0pi-
3592  * (charge conjugation of D*(2010)+pi0 -> D*(2007)0pi+) */
3593  add_channel(
3594  process_list,
3596  type_Dstarbar_z, type_pi_m);
3597  break;
3598  }
3599  default:
3600  break;
3601  }
3602 
3603  return process_list;
3604 }
3605 
3607  const ReactionsBitSet& included_2to2,
3608  const CharmRescattering charm_rescattering) const {
3609  CollisionBranchList process_list;
3610  if ((included_2to2[IncludedReactions::Charm_T_matrix] == 0) ||
3611  !(charm_rescattering == CharmRescattering::T_Matrix)) {
3612  return process_list;
3613  }
3614  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
3615  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
3616  const auto pdg_D =
3617  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_a.code() : pdg_b.code();
3618  const auto pdg_kaon =
3619  (pdg_a.is_Dmeson() || pdg_a.is_Dstar2007()) ? pdg_b.code() : pdg_a.code();
3620 
3621  /* Adding the following channels, a check for detailed balance is not needed
3622  * because the parametrization is explicit. */
3623  switch (pack(pdg_D, pdg_kaon)) {
3624  // Channels for DK scatterings
3625  case pack(pdg::D_p, pdg::K_z): {
3626  const auto& type_D_z = ParticleType::find(pdg::D_z);
3627  const auto& type_K_p = ParticleType::find(pdg::K_p);
3628  // D+K0 -> D0K+
3629  add_channel(
3630  process_list, [&] { return DplusKzero_DzeroKplus(sqrt_s_); }, sqrt_s_,
3631  type_D_z, type_K_p);
3632  break;
3633  }
3634  case pack(pdg::D_m, pdg::Kbar_z): {
3635  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3636  const auto& type_K_m = ParticleType::find(pdg::K_m);
3637  // D-Kbar0 -> Dbar0K- (charge conjugation of D+K0 -> D0K+)
3638  add_channel(
3639  process_list, [&] { return DplusKzero_DzeroKplus(sqrt_s_); }, sqrt_s_,
3640  type_Dbar_z, type_K_m);
3641  break;
3642  }
3643  case pack(pdg::D_z, pdg::K_p): {
3644  const auto& type_D_p = ParticleType::find(pdg::D_p);
3645  const auto& type_K_z = ParticleType::find(pdg::K_z);
3646  // D0K+ -> D+K0
3647  add_channel(
3648  process_list, [&] { return DzeroKplus_DplusKzero(sqrt_s_); }, sqrt_s_,
3649  type_D_p, type_K_z);
3650  break;
3651  }
3652  case pack(pdg::Dbar_z, pdg::K_m): {
3653  const auto& type_D_m = ParticleType::find(pdg::D_m);
3654  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
3655  // Dbar0K- -> D-Kbar0 (charge conjugation of D0K+ -> D+K0)
3656  add_channel(
3657  process_list, [&] { return DzeroKplus_DplusKzero(sqrt_s_); }, sqrt_s_,
3658  type_D_m, type_Kbar_z);
3659  break;
3660  }
3661  case pack(pdg::D_p, pdg::K_m): {
3662  const auto& type_D_z = ParticleType::find(pdg::D_z);
3663  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
3664  // D+K- -> D0Kbar0
3665  add_channel(
3666  process_list, [&] { return DplusKminus_DzeroKbarzero(sqrt_s_); },
3667  sqrt_s_, type_D_z, type_Kbar_z);
3668  break;
3669  }
3670  case pack(pdg::D_m, pdg::K_p): {
3671  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3672  const auto& type_K_z = ParticleType::find(pdg::K_z);
3673  // D-K+ -> Dbar0K0 (charge conjugation of D+K- -> D0Kbar0)
3674  add_channel(
3675  process_list, [&] { return DplusKminus_DzeroKbarzero(sqrt_s_); },
3676  sqrt_s_, type_Dbar_z, type_K_z);
3677  break;
3678  }
3679  case pack(pdg::D_z, pdg::Kbar_z): {
3680  const auto& type_D_p = ParticleType::find(pdg::D_p);
3681  const auto& type_K_m = ParticleType::find(pdg::K_m);
3682  // D0Kbar0 -> D+K-
3683  add_channel(
3684  process_list, [&] { return DzeroKbarzero_DplusKminus(sqrt_s_); },
3685  sqrt_s_, type_D_p, type_K_m);
3686  break;
3687  }
3688  case pack(pdg::Dbar_z, pdg::K_z): {
3689  const auto& type_D_m = ParticleType::find(pdg::D_m);
3690  const auto& type_K_p = ParticleType::find(pdg::K_p);
3691  // Dbar0K0 -> D-K+ (charge conjugation of D0Kbar0 -> D+K-)
3692  add_channel(
3693  process_list, [&] { return DzeroKbarzero_DplusKminus(sqrt_s_); },
3694  sqrt_s_, type_D_m, type_K_p);
3695  break;
3696  }
3697  // Channels for D*K scatterings
3698  case pack(pdg::Dstar_p, pdg::K_z): {
3699  const auto& type_Dstar_z = ParticleType::find(pdg::Dstar_z);
3700  const auto& type_K_p = ParticleType::find(pdg::K_p);
3701  // D*(2010)+K0 -> D*(2007)0K+
3702  add_channel(
3703  process_list, [&] { return DstarplusKzero_DstarzeroKplus(sqrt_s_); },
3704  sqrt_s_, type_Dstar_z, type_K_p);
3705  break;
3706  }
3707  case pack(pdg::Dstar_m, pdg::Kbar_z): {
3708  const auto& type_Dstarbar_z = ParticleType::find(pdg::Dstarbar_z);
3709  const auto& type_K_m = ParticleType::find(pdg::K_m);
3710  /* D*(2010)-Kbar0 -> D*(2007)bar0K-
3711  * (charge conjugation of D*(2010)+K0 -> D*(2007)0K+) */
3712  add_channel(
3713  process_list, [&] { return DstarplusKzero_DstarzeroKplus(sqrt_s_); },
3714  sqrt_s_, type_Dstarbar_z, type_K_m);
3715  break;
3716  }
3717  case pack(pdg::Dstar_z, pdg::K_p): {
3718  const auto& type_Dstar_p = ParticleType::find(pdg::Dstar_p);
3719  const auto& type_K_z = ParticleType::find(pdg::K_z);
3720  // D*(2007)0K+ -> D*(2010)+K0
3721  add_channel(
3722  process_list, [&] { return DstarzeroKplus_DstarplusKzero(sqrt_s_); },
3723  sqrt_s_, type_Dstar_p, type_K_z);
3724  break;
3725  }
3726  case pack(pdg::Dstarbar_z, pdg::K_m): {
3727  const auto& type_Dstar_m = ParticleType::find(pdg::Dstar_m);
3728  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
3729  /* D*(2007)bar0K- -> D*(2010)-Kbar0
3730  * (charge conjugation of D*(2007)0K+ -> D*(2010)+K0) */
3731  add_channel(
3732  process_list, [&] { return DstarzeroKplus_DstarplusKzero(sqrt_s_); },
3733  sqrt_s_, type_Dstar_m, type_Kbar_z);
3734  break;
3735  }
3736  case pack(pdg::Dstar_p, pdg::K_m): {
3737  const auto& type_Dstar_z = ParticleType::find(pdg::Dstar_z);
3738  const auto& type_Kbar_z = ParticleType::find(pdg::Kbar_z);
3739  // D*(2010)+K- -> D*(2007)0Kbar0
3740  add_channel(
3741  process_list,
3743  type_Dstar_z, type_Kbar_z);
3744  break;
3745  }
3746  case pack(pdg::Dstar_m, pdg::K_p): {
3747  const auto& type_Dstarbar_z = ParticleType::find(pdg::Dstarbar_z);
3748  const auto& type_K_z = ParticleType::find(pdg::K_z);
3749  /* D*(2010)-K+ -> D*(2007)bar0K0
3750  * (charge conjugation of D*(2010)+K- -> D*(2007)0Kbar0) */
3751  add_channel(
3752  process_list,
3754  type_Dstarbar_z, type_K_z);
3755  break;
3756  }
3757  case pack(pdg::Dstar_z, pdg::Kbar_z): {
3758  const auto& type_Dstar_p = ParticleType::find(pdg::Dstar_p);
3759  const auto& type_K_m = ParticleType::find(pdg::K_m);
3760  // D*(2007)0Kbar0 -> D*(2010)+K-
3761  add_channel(
3762  process_list,
3764  type_Dstar_p, type_K_m);
3765  break;
3766  }
3767  case pack(pdg::Dstarbar_z, pdg::K_z): {
3768  const auto& type_Dstar_m = ParticleType::find(pdg::Dstar_m);
3769  const auto& type_K_p = ParticleType::find(pdg::K_p);
3770  /* D*(2007)bar0K0 -> D*(2010)-K+
3771  * (charge conjugation of D*(2007)0Kbar0 -> D*(2010)+K-) */
3772  add_channel(
3773  process_list,
3775  type_Dstar_m, type_K_p);
3776  break;
3777  }
3778  default:
3779  break;
3780  }
3781 
3782  return process_list;
3783 }
3784 
3785 CollisionBranchList CrossSections::DN_xx(
3786  const ReactionsBitSet& included_2to2,
3787  const CharmRescattering charm_rescattering) const {
3788  CollisionBranchList process_list;
3789  if ((included_2to2[IncludedReactions::Charm_T_matrix] == 0) ||
3790  !(charm_rescattering == CharmRescattering::T_Matrix)) {
3791  return process_list;
3792  }
3793  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
3794  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
3795  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
3796  const auto pdg_nucleon = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
3797 
3798  /* Adding the following channels, a check for detailed balance is not needed
3799  * because the parametrization is explicit. */
3800  switch (pack(pdg_D, pdg_nucleon)) {
3801  case pack(pdg::D_p, pdg::n): {
3802  const auto& type_D_z = ParticleType::find(pdg::D_z);
3803  const auto& type_p = ParticleType::find(pdg::p);
3804  // D+n -> D0p
3805  add_channel(
3806  process_list, [&] { return Dplusn_Dzerop(sqrt_s_); }, sqrt_s_,
3807  type_D_z, type_p);
3808  break;
3809  }
3810  case pack(pdg::D_m, -pdg::n): {
3811  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3812  const auto& type_pbar = ParticleType::find(-pdg::p);
3813  // D-nbar -> Dbar0pbar (charge conjugation of D+n -> D0p)
3814  add_channel(
3815  process_list, [&] { return Dplusn_Dzerop(sqrt_s_); }, sqrt_s_,
3816  type_Dbar_z, type_pbar);
3817  break;
3818  }
3819  case pack(pdg::D_z, pdg::p): {
3820  const auto& type_D_p = ParticleType::find(pdg::D_p);
3821  const auto& type_n = ParticleType::find(pdg::n);
3822  // D0p -> D+n
3823  add_channel(
3824  process_list, [&] { return Dzerop_Dplusn(sqrt_s_); }, sqrt_s_,
3825  type_D_p, type_n);
3826  break;
3827  }
3828  case pack(pdg::Dbar_z, -pdg::p): {
3829  const auto& type_D_m = ParticleType::find(pdg::D_m);
3830  const auto& type_nbar = ParticleType::find(-pdg::n);
3831  // Dbar0pbar -> D-nbar (charge conjugation of D0p -> D+n)
3832  add_channel(
3833  process_list, [&] { return Dzerop_Dplusn(sqrt_s_); }, sqrt_s_,
3834  type_D_m, type_nbar);
3835  break;
3836  }
3837  case pack(pdg::D_m, pdg::p): {
3838  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3839  const auto& type_n = ParticleType::find(pdg::n);
3840  // D-p -> Dbar0n
3841  add_channel(
3842  process_list, [&] { return Dminusp_Dbarzeron(sqrt_s_); }, sqrt_s_,
3843  type_Dbar_z, type_n);
3844  break;
3845  }
3846  case pack(pdg::D_p, -pdg::p): {
3847  const auto& type_D_z = ParticleType::find(pdg::D_z);
3848  const auto& type_nbar = ParticleType::find(-pdg::n);
3849  // D+pbar -> D0nbar (charge conjugation of D-p -> Dbar0n)
3850  add_channel(
3851  process_list, [&] { return Dminusp_Dbarzeron(sqrt_s_); }, sqrt_s_,
3852  type_D_z, type_nbar);
3853  break;
3854  }
3855  case pack(pdg::Dbar_z, pdg::n): {
3856  const auto& type_D_m = ParticleType::find(pdg::D_m);
3857  const auto& type_p = ParticleType::find(pdg::p);
3858  // Dbar0n -> D-p
3859  add_channel(
3860  process_list, [&] { return Dbarzeron_Dminusp(sqrt_s_); }, sqrt_s_,
3861  type_D_m, type_p);
3862  break;
3863  }
3864  case pack(pdg::D_z, -pdg::n): {
3865  const auto& type_D_p = ParticleType::find(pdg::D_p);
3866  const auto& type_pbar = ParticleType::find(-pdg::p);
3867  // D0nbar -> D+pbar (charge conjugation of Dbar0n -> D-p)
3868  add_channel(
3869  process_list, [&] { return Dbarzeron_Dminusp(sqrt_s_); }, sqrt_s_,
3870  type_D_p, type_pbar);
3871  break;
3872  }
3873  default:
3874  break;
3875  }
3876 
3877  return process_list;
3878 }
3879 
3880 CollisionBranchList CrossSections::DDelta_xx(
3881  const ReactionsBitSet& included_2to2,
3882  const CharmRescattering charm_rescattering) const {
3883  CollisionBranchList process_list;
3884  if ((included_2to2[IncludedReactions::Charm_T_matrix] == 0) ||
3885  !(charm_rescattering == CharmRescattering::T_Matrix)) {
3886  return process_list;
3887  }
3888  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
3889  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
3890  const auto pdg_D = pdg_a.is_Dmeson() ? pdg_a.code() : pdg_b.code();
3891  const auto pdg_Delta = pdg_a.is_Dmeson() ? pdg_b.code() : pdg_a.code();
3892 
3893  /* Adding the following channels, a check for detailed balance is not needed
3894  * because the parametrization is explicit. */
3895  switch (pack(pdg_D, pdg_Delta)) {
3896  case pack(pdg::D_p, pdg::Delta_p): {
3897  const auto& type_D_z = ParticleType::find(pdg::D_z);
3898  const auto& type_Delta_pp = ParticleType::find(pdg::Delta_pp);
3899  // D+Δ+ -> D0Δ++
3900  add_channel(
3901  process_list,
3903  type_D_z, type_Delta_pp);
3904  break;
3905  }
3906  case pack(pdg::D_m, -pdg::Delta_p): {
3907  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3908  const auto& type_Deltabar_mm = ParticleType::find(-pdg::Delta_pp);
3909  // D-Δbar- -> Dbar0Δbar-- (charge conjugation of D+Δ+ -> D0Δ++)
3910  add_channel(
3911  process_list,
3913  type_Dbar_z, type_Deltabar_mm);
3914  break;
3915  }
3916  case pack(pdg::D_p, pdg::Delta_m): {
3917  const auto& type_D_z = ParticleType::find(pdg::D_z);
3918  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
3919  // D+Δ- -> D0Δ0
3920  add_channel(
3921  process_list, [&] { return DplusDeltaminus_DzeroDeltazero(sqrt_s_); },
3922  sqrt_s_, type_D_z, type_Delta_z);
3923  break;
3924  }
3925  case pack(pdg::D_m, -pdg::Delta_m): {
3926  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3927  const auto& type_Deltabar_z = ParticleType::find(-pdg::Delta_z);
3928  // D-Δbar- -> Dbar0Δbar0 (charge conjugation of D+Δ- -> D0Δ0)
3929  add_channel(
3930  process_list, [&] { return DplusDeltaminus_DzeroDeltazero(sqrt_s_); },
3931  sqrt_s_, type_Dbar_z, type_Deltabar_z);
3932  break;
3933  }
3934  case pack(pdg::D_p, pdg::Delta_z): {
3935  const auto& type_D_z = ParticleType::find(pdg::D_z);
3936  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
3937  // D+Δ0 -> D0Δ+
3938  add_channel(
3939  process_list, [&] { return DplusDeltazero_DzeroDeltaplus(sqrt_s_); },
3940  sqrt_s_, type_D_z, type_Delta_p);
3941  break;
3942  }
3943  case pack(pdg::D_m, -pdg::Delta_z): {
3944  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
3945  const auto& type_Deltabar_m = ParticleType::find(-pdg::Delta_p);
3946  // D-Δbar0 -> Dbar0Δbar- (charge conjugation of D+Δ0 -> D0Δ+)
3947  add_channel(
3948  process_list, [&] { return DplusDeltazero_DzeroDeltaplus(sqrt_s_); },
3949  sqrt_s_, type_Dbar_z, type_Deltabar_m);
3950  break;
3951  }
3952  case pack(pdg::D_z, pdg::Delta_p): {
3953  const auto& type_D_p = ParticleType::find(pdg::D_p);
3954  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
3955  // D0Δ+ -> D+Δ0
3956  add_channel(
3957  process_list, [&] { return DzeroDeltaplus_DplusDeltazero(sqrt_s_); },
3958  sqrt_s_, type_D_p, type_Delta_z);
3959  break;
3960  }
3961  case pack(pdg::Dbar_z, -pdg::Delta_p): {
3962  const auto& type_D_m = ParticleType::find(pdg::D_m);
3963  const auto& type_Deltabar_z = ParticleType::find(-pdg::Delta_z);
3964  // Dbar0Δbar- -> D-Δbar0 (charge conjugation of D0Δ+ -> D+Δ0)
3965  add_channel(
3966  process_list, [&] { return DzeroDeltaplus_DplusDeltazero(sqrt_s_); },
3967  sqrt_s_, type_D_m, type_Deltabar_z);
3968  break;
3969  }
3970  case pack(pdg::D_z, pdg::Delta_pp): {
3971  const auto& type_D_p = ParticleType::find(pdg::D_p);
3972  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
3973  // D0Δ++ -> D+Δ+
3974  add_channel(
3975  process_list,
3977  type_D_p, type_Delta_p);
3978  break;
3979  }
3980  case pack(pdg::Dbar_z, -pdg::Delta_pp): {
3981  const auto& type_D_m = ParticleType::find(pdg::D_m);
3982  const auto& type_Deltabar_m = ParticleType::find(-pdg::Delta_p);
3983  // Dbar0Δbar-- -> D-Δbar- (charge conjugation of D0Δ++ -> D+Δ+)
3984  add_channel(
3985  process_list,
3987  type_D_m, type_Deltabar_m);
3988  break;
3989  }
3990  case pack(pdg::D_z, pdg::Delta_z): {
3991  const auto& type_D_p = ParticleType::find(pdg::D_p);
3992  const auto& type_Delta_m = ParticleType::find(pdg::Delta_m);
3993  // D0Δ0 -> D+Δ-
3994  add_channel(
3995  process_list, [&] { return DzeroDeltazero_DplusDeltaminus(sqrt_s_); },
3996  sqrt_s_, type_D_p, type_Delta_m);
3997  break;
3998  }
3999  case pack(pdg::Dbar_z, -pdg::Delta_z): {
4000  const auto& type_D_m = ParticleType::find(pdg::D_m);
4001  const auto& type_Deltabar_p = ParticleType::find(-pdg::Delta_m);
4002  // Dbar0Δbar0 -> D-Δbar+ (charge conjugation of D0Δ0 -> D+Δ-)
4003  add_channel(
4004  process_list, [&] { return DzeroDeltazero_DplusDeltaminus(sqrt_s_); },
4005  sqrt_s_, type_D_m, type_Deltabar_p);
4006  break;
4007  }
4008  case pack(pdg::D_m, pdg::Delta_p): {
4009  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
4010  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
4011  // D-Δ+ -> Dbar0Δ0
4012  add_channel(
4013  process_list,
4015  type_Dbar_z, type_Delta_z);
4016  break;
4017  }
4018  case pack(pdg::D_p, -pdg::Delta_p): {
4019  const auto& type_D_z = ParticleType::find(pdg::D_z);
4020  const auto& type_Deltabar_z = ParticleType::find(-pdg::Delta_z);
4021  // D+Δbar- -> D0Δbar0 (charge conjugation of D-Δ+ -> Dbar0Δ0)
4022  add_channel(
4023  process_list,
4025  type_D_z, type_Deltabar_z);
4026  break;
4027  }
4028  case pack(pdg::D_m, pdg::Delta_pp): {
4029  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
4030  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
4031  // D-Δ++ -> Dbar0Δ+
4032  add_channel(
4033  process_list,
4035  sqrt_s_, type_Dbar_z, type_Delta_p);
4036  break;
4037  }
4038  case pack(pdg::D_p, -pdg::Delta_pp): {
4039  const auto& type_D_z = ParticleType::find(pdg::D_z);
4040  const auto& type_Deltabar_m = ParticleType::find(-pdg::Delta_p);
4041  // D+Δbar-- -> D0Δbar- (charge conjugation of D-Δ++ -> Dbar0Δ+)
4042  add_channel(
4043  process_list,
4045  sqrt_s_, type_D_z, type_Deltabar_m);
4046  break;
4047  }
4048  case pack(pdg::D_m, pdg::Delta_z): {
4049  const auto& type_Dbar_z = ParticleType::find(pdg::Dbar_z);
4050  const auto& type_Delta_m = ParticleType::find(pdg::Delta_m);
4051  // D-Δ0 -> Dbar0Δ-
4052  add_channel(
4053  process_list,
4055  type_Dbar_z, type_Delta_m);
4056  break;
4057  }
4058  case pack(pdg::D_p, -pdg::Delta_z): {
4059  const auto& type_D_z = ParticleType::find(pdg::D_z);
4060  const auto& type_Deltabar_p = ParticleType::find(-pdg::Delta_m);
4061  // D+Δbar0 -> D0Δbar+ (charge conjugation of D-Δ0 -> Dbar0Δ-)
4062  add_channel(
4063  process_list,
4065  type_D_z, type_Deltabar_p);
4066  break;
4067  }
4068  case pack(pdg::Dbar_z, pdg::Delta_p): {
4069  const auto& type_D_m = ParticleType::find(pdg::D_m);
4070  const auto& type_Delta_pp = ParticleType::find(pdg::Delta_pp);
4071  // Dbar0Δ+ -> D-Δ++
4072  add_channel(
4073  process_list,
4075  sqrt_s_, type_D_m, type_Delta_pp);
4076  break;
4077  }
4078  case pack(pdg::D_z, -pdg::Delta_p): {
4079  const auto& type_D_p = ParticleType::find(pdg::D_p);
4080  const auto& type_Deltabar_mm = ParticleType::find(-pdg::Delta_pp);
4081  // D0Δbar- -> D+Δbar-- (charge conjugation of Dbar0Δ+ -> D-Δ++)
4082  add_channel(
4083  process_list,
4085  sqrt_s_, type_D_p, type_Deltabar_mm);
4086  break;
4087  }
4088  case pack(pdg::Dbar_z, pdg::Delta_m): {
4089  const auto& type_D_m = ParticleType::find(pdg::D_m);
4090  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
4091  // Dbar0Δ- -> D-Δ0
4092  add_channel(
4093  process_list,
4095  type_D_m, type_Delta_z);
4096  break;
4097  }
4098  case pack(pdg::D_z, -pdg::Delta_m): {
4099  const auto& type_D_p = ParticleType::find(pdg::D_p);
4100  const auto& type_Deltabar_z = ParticleType::find(-pdg::Delta_z);
4101  // D0Δbar+ -> D+Δbar0 (charge conjugation of Dbar0Δ- -> D-Δ0)
4102  add_channel(
4103  process_list,
4105  type_D_p, type_Deltabar_z);
4106  break;
4107  }
4108  case pack(pdg::Dbar_z, pdg::Delta_z): {
4109  const auto& type_D_m = ParticleType::find(pdg::D_m);
4110  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
4111  // Dbar0Δ0 -> D-Δ+
4112  add_channel(
4113  process_list,
4115  type_D_m, type_Delta_p);
4116  break;
4117  }
4118  case pack(pdg::D_z, -pdg::Delta_z): {
4119  const auto& type_D_p = ParticleType::find(pdg::D_p);
4120  const auto& type_Deltabar_m = ParticleType::find(-pdg::Delta_p);
4121  // D0Δbar0 -> D+Δbar- (charge conjugation of Dbar0Δ0 -> D-Δ+)
4122  add_channel(
4123  process_list,
4125  type_D_p, type_Deltabar_m);
4126  break;
4127  }
4128  default:
4129  break;
4130  }
4131 
4132  return process_list;
4133 }
4134 
4136  double total_string_xs, StringProcess* string_process,
4137  const ScatterActionsFinderParameters& finder_parameters) const {
4138  if (!string_process) {
4139  throw std::runtime_error("string_process should be initialized.");
4140  }
4141 
4142  CollisionBranchList channel_list;
4143  if (total_string_xs <= 0.) {
4144  return channel_list;
4145  }
4146 
4147  double mandelstam_s = sqrt_s_ * sqrt_s_;
4148  /* Get mapped PDG id for evaluation of the parametrized cross sections for
4149  * diffractive processes. This must be rescaled according to additive quark
4150  * model in the case of exotic hadrons. Also calculate the multiplicative
4151  * factor for AQM based on the quark contents. */
4152  std::array<int, 2> pdgid;
4153  double AQM_scaling = 1.;
4154  for (int i = 0; i < 2; i++) {
4155  PdgCode pdg = incoming_particles_[i].type().pdgcode();
4156  pdgid[i] = StringProcess::pdg_map_for_pythia(pdg);
4157  AQM_scaling *= finder_parameters.AQM_scaling_factor(pdg);
4158  }
4159 
4160  /* Determine if the initial state is a baryon-antibaryon pair,
4161  * which can annihilate. */
4162  bool can_annihilate = false;
4163  if (is_BBbar_pair_) {
4164  int n_q_types = 5; // u, d, s, c, b
4165  for (int iq = 1; iq <= n_q_types; iq++) {
4166  std::array<int, 2> nquark;
4167  for (int i = 0; i < 2; i++) {
4168  nquark[i] =
4169  incoming_particles_[i].type().pdgcode().net_quark_number(iq);
4170  }
4171  if (nquark[0] != 0 && nquark[1] != 0) {
4172  can_annihilate = true;
4173  break;
4174  }
4175  }
4176  }
4177 
4178  /* The case for baryon/anti-baryon annihilation is treated separately,
4179  * as in this case we use only one way to break up the particles, namely
4180  * into 2 mesonic strings of equal mass after annihilating one quark-
4181  * anti-quark pair. See StringProcess::next_BBbarAnn() */
4182  double sig_annihilation = 0.0;
4183  if (can_annihilate) {
4184  /* In the case of baryon-antibaryon pair, the parametrized cross section for
4185  * annihilation will be added. See xs_ppbar_annihilation(). */
4186  mandelstam_s = effective_AQM_s(
4187  mandelstam_s, incoming_particles_[0].effective_mass(),
4188  incoming_particles_[1].effective_mass(), nucleon_mass, nucleon_mass);
4189  double xs_param = xs_ppbar_annihilation(mandelstam_s);
4190  if (finder_parameters.use_AQM) {
4191  xs_param *= AQM_scaling;
4192  }
4193  sig_annihilation = std::min(total_string_xs, xs_param);
4194  }
4195 
4196  /* Total parametrized cross-section (I) and pythia-produced total
4197  * cross-section (II) do not necessarily coincide. If I > II then
4198  * non-diffractive cross-section is reinforced to get I == II.
4199  * If I < II then partial cross-sections are drained one-by-one
4200  * to reduce II until I == II:
4201  * first non-diffractive, then double-diffractive, then
4202  * single-diffractive AB->AX and AB->XB in equal proportion.
4203  * The way it is done here is not unique. I (ryu) think that at high energy
4204  * collision this is not an issue, but at sqrt_s < 10 GeV it may matter. */
4205  std::array<double, 3> xs_diffractive =
4206  string_process->cross_sections_diffractive(pdgid[0], pdgid[1],
4207  std::sqrt(mandelstam_s));
4208 
4209  if (finder_parameters.use_AQM) {
4210  for (double& x : xs_diffractive) {
4211  x *= AQM_scaling;
4212  }
4213  }
4214 
4215  double single_diffr_AX = xs_diffractive[0];
4216  double single_diffr_XB = xs_diffractive[1];
4217  double double_diffr = xs_diffractive[2];
4218 
4219  double single_diffr = single_diffr_AX + single_diffr_XB;
4220  double diffractive = single_diffr + double_diffr;
4221 
4222  const double nondiffractive =
4223  std::max(0., total_string_xs - sig_annihilation - diffractive);
4224 
4225  diffractive = total_string_xs - sig_annihilation - nondiffractive;
4226  double_diffr = std::max(0., diffractive - single_diffr);
4227 
4228  const double a =
4229  single_diffr > 0.0 ? (diffractive - double_diffr) / single_diffr : 0.0;
4230 
4231  single_diffr_AX *= a;
4232  single_diffr_XB *= a;
4233 
4234  assert(std::abs(single_diffr_AX + single_diffr_XB + double_diffr +
4235  sig_annihilation + nondiffractive - total_string_xs) < 1.e-6);
4236  enum class Proc { ND, SD_AX, SD_XB, DD, N };
4237  enum class Comp { Soft, Hard, N };
4238 
4239  constexpr std::size_t n_proc = static_cast<std::size_t>(Proc::N);
4240  constexpr std::size_t n_comp = static_cast<std::size_t>(Comp::N);
4241 
4242  std::array<std::array<double, n_comp>, n_proc> split_xs{};
4243 
4244  auto proc_idx = [](Proc p) { return static_cast<std::size_t>(p); };
4245 
4246  auto comp_idx = [](Comp c) { return static_cast<std::size_t>(c); };
4247 
4248  auto soft = [&](Proc p) -> double& {
4249  return split_xs[proc_idx(p)][comp_idx(Comp::Soft)];
4250  };
4251 
4252  auto hard = [&](Proc p) -> double& {
4253  return split_xs[proc_idx(p)][comp_idx(Comp::Hard)];
4254  };
4255 
4256  auto set_all_soft = [&](Proc p, double total) {
4257  soft(p) = total;
4258  hard(p) = 0.0;
4259  };
4260 
4261  auto split = [&](Proc p, double total, double weight_hard) {
4262  hard(p) = total * weight_hard;
4263  soft(p) = total - hard(p);
4264  };
4265  auto split_all_string_processes = [&](double weight_hard) {
4266  split(Proc::ND, nondiffractive, weight_hard);
4267  split(Proc::SD_AX, single_diffr_AX, weight_hard);
4268  split(Proc::SD_XB, single_diffr_XB, weight_hard);
4269  split(Proc::DD, double_diffr, weight_hard);
4270  };
4271 
4272  if (finder_parameters.hard_string_transition_mode ==
4274  const auto& [hard_transition_start, hard_transition_end] =
4275  finder_parameters.hard_string_transition_energy_range;
4276 
4277  const double weight_hard = transition_probability_at_sqrts(
4278  hard_transition_start, hard_transition_end);
4279  split_all_string_processes(weight_hard);
4280  } else if (nondiffractive > 0.0) {
4281  const double hard_xsec = AQM_scaling * string_hard_cross_section();
4282 
4283  /* Use the non-diffractive exponential transition probability for all
4284 
4285  * string-excitation channels, so that soft strings are suppressed at high
4286  * energies also for diffractive processes. */
4287  const double weight_soft = std::exp(-hard_xsec / nondiffractive);
4288  const double weight_hard = std::clamp(1.0 - weight_soft, 0.0, 1.0);
4289  split_all_string_processes(weight_hard);
4290  } else {
4291  set_all_soft(Proc::ND, nondiffractive);
4292  set_all_soft(Proc::SD_AX, single_diffr_AX);
4293  set_all_soft(Proc::SD_XB, single_diffr_XB);
4294  set_all_soft(Proc::DD, double_diffr);
4295  }
4296  logg[LCrossSections].debug("Soft string cross sections [mb] are");
4297  logg[LCrossSections].debug("Soft single-diffractive AB->AX: ",
4298  soft(Proc::SD_AX));
4299  logg[LCrossSections].debug("Soft single-diffractive AB->XB: ",
4300  soft(Proc::SD_XB));
4301  logg[LCrossSections].debug("Soft double-diffractive AB->XX: ",
4302  soft(Proc::DD));
4303  logg[LCrossSections].debug("Soft non-diffractive: ", soft(Proc::ND));
4304  logg[LCrossSections].debug("B-Bbar annihilation: ", sig_annihilation);
4305 
4306  logg[LCrossSections].debug("Hard string cross sections [mb] are");
4307  logg[LCrossSections].debug("Hard single-diffractive AB->AX: ",
4308  hard(Proc::SD_AX));
4309  logg[LCrossSections].debug("Hard single-diffractive AB->XB: ",
4310  hard(Proc::SD_XB));
4311  logg[LCrossSections].debug("Hard double-diffractive AB->XX: ",
4312  hard(Proc::DD));
4313  logg[LCrossSections].debug("Hard non-diffractive: ", hard(Proc::ND));
4314 
4315  // cross section of soft string excitation including annihilation
4316  const double sig_string_soft = soft(Proc::SD_AX) + soft(Proc::SD_XB) +
4317  soft(Proc::DD) + soft(Proc::ND) +
4318  sig_annihilation;
4319 
4320  // fill the list of process channels
4321  if (sig_string_soft > 0.) {
4322  channel_list.push_back(std::make_unique<CollisionBranch>(
4323  soft(Proc::SD_AX), ProcessType::StringSoftSingleDiffractiveAX));
4324  channel_list.push_back(std::make_unique<CollisionBranch>(
4325  soft(Proc::SD_XB), ProcessType::StringSoftSingleDiffractiveXB));
4326  channel_list.push_back(std::make_unique<CollisionBranch>(
4327  soft(Proc::DD), ProcessType::StringSoftDoubleDiffractive));
4328  channel_list.push_back(std::make_unique<CollisionBranch>(
4329  soft(Proc::ND), ProcessType::StringSoftNonDiffractive));
4330 
4331  if (can_annihilate) {
4332  channel_list.push_back(std::make_unique<CollisionBranch>(
4333  sig_annihilation, ProcessType::StringSoftAnnihilation));
4334  }
4335  }
4336 
4337  if (hard(Proc::SD_AX) > 0.) {
4338  channel_list.push_back(std::make_unique<CollisionBranch>(
4339  hard(Proc::SD_AX), ProcessType::StringHardSingleDiffractiveAX));
4340  }
4341  if (hard(Proc::SD_XB) > 0.) {
4342  channel_list.push_back(std::make_unique<CollisionBranch>(
4343  hard(Proc::SD_XB), ProcessType::StringHardSingleDiffractiveXB));
4344  }
4345  if (hard(Proc::DD) > 0.) {
4346  channel_list.push_back(std::make_unique<CollisionBranch>(
4347  hard(Proc::DD), ProcessType::StringHardDoubleDiffractive));
4348  }
4349  if (hard(Proc::ND) > 0.) {
4350  channel_list.push_back(std::make_unique<CollisionBranch>(
4351  hard(Proc::ND), ProcessType::StringHardNonDiffractive));
4352  }
4353 
4354  return channel_list;
4355 }
4357  const ScatterActionsFinderParameters& finder_parameters) const {
4358  const PdgCode& pdg_a = incoming_particles_[0].type().pdgcode();
4359  const PdgCode& pdg_b = incoming_particles_[1].type().pdgcode();
4360 
4361  const double s = sqrt_s_ * sqrt_s_;
4362  double xs = 0.;
4363 
4364  // Currently all BB collisions use the nucleon-nucleon parametrizations.
4365  if (pdg_a.is_baryon() && pdg_b.is_baryon()) {
4366  const double eff_s = effective_AQM_s(
4367  s, incoming_particles_[0].effective_mass(),
4368  incoming_particles_[1].effective_mass(), nucleon_mass, nucleon_mass);
4369  if (pdg_a == pdg_b) {
4370  xs = pp_high_energy(eff_s); // pp, nn
4371  } else if (pdg_a.antiparticle_sign() * pdg_b.antiparticle_sign() == 1) {
4372  xs = np_high_energy(eff_s); // np, nbarpbar
4373  } else if (pdg_a.antiparticle_sign() * pdg_b.antiparticle_sign() == -1) {
4374  /* In the case of baryon-antibaryon interactions,
4375  * the low-energy cross section must be involved
4376  * due to annihilation processes (via strings). */
4377  double xs_l = ppbar_total(eff_s);
4378  double xs_h = 0.;
4379  if (pdg_a.is_antiparticle_of(pdg_b)) {
4380  xs_h = ppbar_high_energy(eff_s); // ppbar, nnbar
4381  } else {
4382  xs_h = npbar_high_energy(eff_s); // npbar, nbarp
4383  }
4384  /* Transition between low and high energy is set to be consistent with
4385  * that defined in string_probability(). */
4386  auto [region_lower, region_upper] =
4387  finder_parameters.transition_high_energy.sqrts_range_NN;
4388  double prob_high =
4389  transition_probability_at_sqrts(region_lower, region_upper);
4390  xs = xs_l * (1. - prob_high) + xs_h * prob_high;
4391  }
4392  }
4393 
4394  // Pion nucleon interaction / baryon-meson
4395  if ((pdg_a == pdg::pi_p && pdg_b == pdg::p) ||
4396  (pdg_b == pdg::pi_p && pdg_a == pdg::p) ||
4397  (pdg_a == pdg::pi_m && pdg_b == pdg::n) ||
4398  (pdg_b == pdg::pi_m && pdg_a == pdg::n)) {
4399  xs = piplusp_high_energy(s); // pi+ p, pi- n
4400  } else if ((pdg_a == pdg::pi_m && pdg_b == pdg::p) ||
4401  (pdg_b == pdg::pi_m && pdg_a == pdg::p) ||
4402  (pdg_a == pdg::pi_p && pdg_b == pdg::n) ||
4403  (pdg_b == pdg::pi_p && pdg_a == pdg::n)) {
4404  xs = piminusp_high_energy(s); // pi- p, pi+ n
4405  } else if ((pdg_a.is_meson() && pdg_b.is_baryon()) ||
4406  (pdg_b.is_meson() && pdg_a.is_baryon())) {
4407  xs = piminusp_high_energy(s); // default for baryon-meson
4408  }
4409 
4410  /* Meson-meson interaction goes through AQM from pi+p,
4411  * see user guide "Use_AQM" */
4412  if (pdg_a.is_meson() && pdg_b.is_meson()) {
4413  /* 2/3 factor since difference of 1 meson between meson-meson
4414  * and baryon-meson */
4415  xs = 2. / 3. * piplusp_high_energy(s);
4416  }
4417 
4418  // AQM scaling for cross-sections
4419  xs *= finder_parameters.AQM_scaling_factor(pdg_a) *
4420  finder_parameters.AQM_scaling_factor(pdg_b);
4421 
4422  return xs;
4423 }
4424 
4426  double cross_sec = 0.;
4427  // Hard strings can only be excited if the lower cutoff by Pythia is fulfilled
4429  return cross_sec;
4430  }
4431  const ParticleData& data_a = incoming_particles_[0];
4432  const ParticleData& data_b = incoming_particles_[1];
4433 
4434  if (data_a.is_baryon() && data_b.is_baryon()) {
4435  // Nucleon-nucleon cross section is used for all baryon-baryon cases.
4436  const double eff_s =
4439  cross_sec = NN_string_hard(eff_s);
4440  } else if (data_a.is_baryon() || data_b.is_baryon()) {
4441  // Nucleon-pion cross section is used for all baryon-meson cases.
4442  cross_sec = Npi_string_hard(sqrt_s_ * sqrt_s_);
4443  } else {
4444  // Pion-pion cross section is used for all meson-meson cases.
4445  cross_sec = pipi_string_hard(sqrt_s_ * sqrt_s_);
4446  }
4447 
4448  return cross_sec;
4449 }
4450 
4451 CollisionBranchPtr CrossSections::NNbar_to_5pi(const double scale_xs) const {
4452  const double s = sqrt_s_ * sqrt_s_;
4453  /* Use difference between total and elastic in order to conserve detailed
4454  * balance for all inelastoc NNbar processes. */
4455  const double nnbar_xsec = std::max(0., ppbar_total(s) - ppbar_elastic(s));
4456  logg[LCrossSections].debug("NNbar cross section for 2-to-5 is: ", nnbar_xsec);
4457 
4458  /* Make collision channel NNbar -> 5π (with same final state as resonance
4459  * approach). */
4460  const auto& type_piz = ParticleType::find(pdg::pi_z);
4461  const auto& type_pip = ParticleType::find(pdg::pi_p);
4462  const auto& type_pim = ParticleType::find(pdg::pi_m);
4463  return std::make_unique<CollisionBranch>(
4464  type_pip, type_pim, type_pip, type_pim, type_piz, nnbar_xsec * scale_xs,
4466 }
4467 
4469  const double current_xs, const double scale_xs) const {
4470  /* Calculate NNbar cross section:
4471  * Parametrized total minus all other present channels.*/
4472  const double s = sqrt_s_ * sqrt_s_;
4473  double nnbar_xsec = std::max(0., ppbar_total(s) * scale_xs - current_xs);
4474  logg[LCrossSections].debug("NNbar cross section is: ", nnbar_xsec);
4475  // Make collision channel NNbar -> ρh₁(1170); eventually decays into 5π
4476  return std::make_unique<CollisionBranch>(ParticleType::find(pdg::h1),
4478  nnbar_xsec, ProcessType::TwoToTwo);
4479 }
4480 
4481 CollisionBranchList CrossSections::NNbar_creation() const {
4482  CollisionBranchList channel_list;
4483  const ParticleType& type_a = incoming_particles_[0].type();
4484  const ParticleType& type_b = incoming_particles_[1].type();
4485  if ((type_a.pdgcode() == pdg::rho_z && type_b.pdgcode() == pdg::h1) ||
4486  (type_a.pdgcode() == pdg::h1 && type_b.pdgcode() == pdg::rho_z)) {
4487  /* Calculate NNbar reverse cross section:
4488  * from reverse reaction (see NNbar_annihilation_cross_section).*/
4489  const double s = sqrt_s_ * sqrt_s_;
4490  const double pcm = cm_momentum();
4491 
4492  const auto& type_N = ParticleType::find(pdg::p);
4493  const auto& type_Nbar = ParticleType::find(-pdg::p);
4494 
4495  // Check available energy
4496  if (sqrt_s_ - 2 * type_N.mass() < 0) {
4497  return channel_list;
4498  }
4499 
4500  double xsection = detailed_balance_factor_RR(sqrt_s_, pcm, type_a, type_b,
4501  type_N, type_Nbar) *
4502  std::max(0., ppbar_total(s) - ppbar_elastic(s));
4503  logg[LCrossSections].debug("NNbar reverse cross section is: ", xsection);
4504  channel_list.push_back(std::make_unique<CollisionBranch>(
4505  type_N, type_Nbar, xsection, ProcessType::TwoToTwo));
4506  channel_list.push_back(std::make_unique<CollisionBranch>(
4509  }
4510  return channel_list;
4511 }
4512 
4514  const bool is_anti_particles) const {
4515  const ParticleType& type_a = incoming_particles_[0].type();
4516  const ParticleType& type_b = incoming_particles_[1].type();
4517  CollisionBranchList process_list;
4518 
4519  const double s = sqrt_s_ * sqrt_s_;
4520  // CM momentum in final state
4521  double p_cm_final = std::sqrt(s - 4. * nucleon_mass * nucleon_mass) / 2.;
4522 
4523  ParticleTypePtrList nuc_or_anti_nuc;
4524  if (is_anti_particles) {
4525  nuc_or_anti_nuc = ParticleType::list_anti_nucleons();
4526  } else {
4527  nuc_or_anti_nuc = ParticleType::list_nucleons();
4528  }
4529 
4530  // Loop over all nucleon or anti-nucleon charge states.
4531  for (ParticleTypePtr nuc_a : nuc_or_anti_nuc) {
4532  for (ParticleTypePtr nuc_b : nuc_or_anti_nuc) {
4533  /* Check for charge conservation. */
4534  if (type_a.charge() + type_b.charge() !=
4535  nuc_a->charge() + nuc_b->charge()) {
4536  continue;
4537  }
4538  // loop over total isospin
4539  for (const int twoI : I_tot_range(*nuc_a, *nuc_b)) {
4540  const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
4541  type_a, type_b, *nuc_a, *nuc_b, twoI);
4542  // If Clebsch-Gordan coefficient is zero, don't bother with the rest
4543  if (std::abs(isospin_factor) < really_small) {
4544  continue;
4545  }
4546 
4547  // Calculate matrix element for inverse process.
4548  const double matrix_element =
4549  nn_to_resonance_matrix_element(sqrt_s_, type_a, type_b, twoI);
4550  if (matrix_element <= 0.) {
4551  continue;
4552  }
4553 
4554  /* Cross section for 2->2 resonance absorption, obtained via detailed
4555  * balance from the inverse reaction. */
4556  const double spin_factor = (nuc_a->spin() + 1) * (nuc_b->spin() + 1);
4557  const int sym_fac_in =
4558  (type_a.iso_multiplet() == type_b.iso_multiplet()) ? 2 : 1;
4559  const int sym_fac_out =
4560  (nuc_a->iso_multiplet() == nuc_b->iso_multiplet()) ? 2 : 1;
4561  const double xsection = isospin_factor * spin_factor * sym_fac_in /
4562  sym_fac_out * p_cm_final * matrix_element /
4563  (s * cm_momentum());
4564 
4565  if (xsection > really_small) {
4566  process_list.push_back(std::make_unique<CollisionBranch>(
4567  *nuc_a, *nuc_b, xsection, ProcessType::TwoToTwo));
4568  logg[LCrossSections].debug(
4569  "2->2 absorption with original particles: ", type_a, type_b);
4570  }
4571  }
4572  }
4573  }
4574  return process_list;
4575 }
4576 
4578  const ParticleType& type_a,
4579  const ParticleType& type_b,
4580  const int twoI) {
4581  const double m_a = type_a.mass();
4582  const double m_b = type_b.mass();
4583  const double msqr = 2. * (m_a * m_a + m_b * m_b);
4584  /* If the c.m. energy is larger than the sum of the pole masses of the
4585  * outgoing particles plus three times the sum of the widths plus 3 GeV, the
4586  * collision will be neglected. This can be problematic for some final-state
4587  * cross sections, but at energies that high strings are used anyway. */
4588  const double w_a = type_a.width_at_pole();
4589  const double w_b = type_b.width_at_pole();
4590  const double uplmt = m_a + m_b + 3.0 * (w_a + w_b) + 3.0;
4591  if (sqrts > uplmt) {
4592  return 0.;
4593  }
4594  /// NN → NΔ: fit sqrt(s)-dependence to OBE model [\iref{Dmitriev:1986st}]
4595  if (((type_a.is_Delta() && type_b.is_nucleon()) ||
4596  (type_b.is_Delta() && type_a.is_nucleon())) &&
4597  (type_a.antiparticle_sign() == type_b.antiparticle_sign())) {
4598  return 68. / std::pow(sqrts - 1.104, 1.951);
4599  /**
4600  * All other processes use a constant matrix element, similar to
4601  * \iref{Bass:1998ca}, eq. (3.35).
4602  */
4603  } else if (((type_a.is_Nstar() && type_b.is_nucleon()) ||
4604  (type_b.is_Nstar() && type_a.is_nucleon())) &&
4605  type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
4606  // NN → NN*
4607  if (twoI == 2) {
4608  return 4.5 / msqr;
4609  } else if (twoI == 0) {
4610  const double parametrization = 14. / msqr;
4611  /**
4612  * pn → pnη cross section is known to be larger than the corresponding
4613  * pp → ppη cross section by a factor of 6.5 [\iref{Calen:1998vh}].
4614  * Since the eta is mainly produced by an intermediate N*(1535) we
4615  * introduce an explicit isospin asymmetry for the production of N*(1535)
4616  * produced in pn vs. pp similar to [\iref{Teis:1996kx}], eq. (29).
4617  */
4618  if (type_a.is_Nstar1535() || type_b.is_Nstar1535()) {
4619  return 6.5 * parametrization;
4620  } else {
4621  return parametrization;
4622  }
4623  }
4624  } else if (((type_a.is_Deltastar() && type_b.is_nucleon()) ||
4625  (type_b.is_Deltastar() && type_a.is_nucleon())) &&
4626  type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
4627  // NN → NΔ*
4628  return 15. / msqr;
4629  } else if ((type_a.is_Delta() && type_b.is_Delta()) &&
4630  (type_a.antiparticle_sign() == type_b.antiparticle_sign())) {
4631  // NN → ΔΔ
4632  if (twoI == 2) {
4633  return 45. / msqr;
4634  } else if (twoI == 0) {
4635  return 120. / msqr;
4636  }
4637  } else if (((type_a.is_Nstar() && type_b.is_Delta()) ||
4638  (type_b.is_Nstar() && type_a.is_Delta())) &&
4639  type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
4640  // NN → ΔN*
4641  return 7. / msqr;
4642  } else if (((type_a.is_Deltastar() && type_b.is_Delta()) ||
4643  (type_b.is_Deltastar() && type_a.is_Delta())) &&
4644  type_a.antiparticle_sign() == type_b.antiparticle_sign()) {
4645  // NN → ΔΔ*
4646  if (twoI == 2) {
4647  return 15. / msqr;
4648  } else if (twoI == 0) {
4649  return 25. / msqr;
4650  }
4651  } else if ((type_a.is_deuteron() && type_b.pdgcode().is_pion()) ||
4652  (type_b.is_deuteron() && type_a.pdgcode().is_pion())) {
4653  /* This parametrization is the result of fitting d+pi->NN cross-section.
4654  * Already Breit-Wigner-like part provides a good fit, exponential fixes
4655  * behaviour around the treshold. The d+pi experimental cross-section
4656  * was taken from Fig. 2 of [\iref{Tanabe:1987vg}]. */
4657  return 0.055 / (pow_int(sqrts - 2.145, 2) + pow_int(0.065, 2)) *
4658  (1.0 - std::exp(-(sqrts - 2.0) * 20.0));
4659  }
4660 
4661  // all cases not listed: zero!
4662  return 0.;
4663 }
4664 
4665 template <class IntegrationMethod>
4667  const ParticleTypePtrList& list_res_1,
4668  const ParticleTypePtrList& list_res_2,
4669  const IntegrationMethod integrator) const {
4670  const ParticleType& type_particle_a = incoming_particles_[0].type();
4671  const ParticleType& type_particle_b = incoming_particles_[1].type();
4672 
4673  CollisionBranchList channel_list;
4674  const double s = sqrt_s_ * sqrt_s_;
4675 
4676  // Loop over specified first resonance list
4677  for (ParticleTypePtr type_res_1 : list_res_1) {
4678  // Loop over specified second resonance list
4679  for (ParticleTypePtr type_res_2 : list_res_2) {
4680  // Check for charge conservation.
4681  if (type_res_1->charge() + type_res_2->charge() !=
4682  type_particle_a.charge() + type_particle_b.charge()) {
4683  continue;
4684  }
4685 
4686  // loop over total isospin
4687  for (const int twoI : I_tot_range(type_particle_a, type_particle_b)) {
4688  const double isospin_factor = isospin_clebsch_gordan_sqr_2to2(
4689  type_particle_a, type_particle_b, *type_res_1, *type_res_2, twoI);
4690  // If Clebsch-Gordan coefficient is zero, don't bother with the rest.
4691  if (std::abs(isospin_factor) < really_small) {
4692  continue;
4693  }
4694 
4695  // Integration limits.
4696  const double lower_limit = type_res_1->min_mass_kinematic();
4697  const double upper_limit = sqrt_s_ - type_res_2->mass();
4698  /* Check the available energy (requiring it to be a little above the
4699  * threshold, because the integration will not work if it's too close).
4700  */
4701  if (upper_limit - lower_limit < 1E-3) {
4702  continue;
4703  }
4704 
4705  // Calculate matrix element.
4706  const double matrix_element = nn_to_resonance_matrix_element(
4707  sqrt_s_, *type_res_1, *type_res_2, twoI);
4708  if (matrix_element <= 0.) {
4709  continue;
4710  }
4711 
4712  /* Calculate resonance production cross section
4713  * using the Breit-Wigner distribution as probability amplitude.
4714  * Integrate over the allowed resonance mass range. */
4715  const double resonance_integral = integrator(*type_res_1, *type_res_2);
4716 
4717  /**
4718  * Cross section for 2->2 process with 1/2 resonance(s) in final state.
4719  * Based on eq. (46) in \iref{Weil:2013mya} and eq. (3.29) in
4720  * \iref{Bass:1998ca}
4721  */
4722  const double spin_factor =
4723  (type_res_1->spin() + 1) * (type_res_2->spin() + 1);
4724  const double xsection = isospin_factor * spin_factor * matrix_element *
4725  resonance_integral / (s * cm_momentum());
4726 
4727  if (xsection > really_small) {
4728  channel_list.push_back(std::make_unique<CollisionBranch>(
4729  *type_res_1, *type_res_2, xsection, ProcessType::TwoToTwo));
4730  logg[LCrossSections].debug(
4731  "Found 2->2 creation process for resonance ", type_res_1, ", ",
4732  type_res_2);
4733  logg[LCrossSections].debug("2->2 with original particles: ",
4734  type_particle_a, type_particle_b);
4735  }
4736  }
4737  }
4738  }
4739  return channel_list;
4740 }
4741 
4743  const ScatterActionsFinderParameters& finder_parameters) const {
4744  /* string fragmentation is enabled when strings_switch is on and the process
4745  * is included in pythia. */
4746  if (!finder_parameters.strings_switch) {
4747  return 0.;
4748  }
4749 
4750  const ParticleType& t1 = incoming_particles_[0].type();
4751  const ParticleType& t2 = incoming_particles_[1].type();
4752  const bool treat_BBbar_with_strings =
4753  (finder_parameters.nnbar_treatment == NNbarTreatment::Strings);
4754  const bool is_NN_scattering =
4755  t1.is_nucleon() && t2.is_nucleon() &&
4756  t1.antiparticle_sign() == t2.antiparticle_sign();
4757  const bool is_BBbar_scattering =
4758  (treat_BBbar_with_strings && is_BBbar_pair_ &&
4759  finder_parameters.use_AQM) ||
4760  (t1.is_nucleon() && t2.is_nucleon() &&
4761  t1.antiparticle_sign() != t2.antiparticle_sign());
4762  const bool is_Npi_scattering = (t1.pdgcode().is_pion() && t2.is_nucleon()) ||
4763  (t1.is_nucleon() && t2.pdgcode().is_pion());
4764  /* True for baryon-baryon, anti-baryon-anti-baryon, baryon-meson,
4765  * anti-baryon-meson and meson-meson*/
4766  const bool is_AQM_scattering =
4767  finder_parameters.use_AQM &&
4768  ((t1.is_baryon() && t2.is_baryon() &&
4769  t1.antiparticle_sign() == t2.antiparticle_sign()) ||
4770  ((t1.is_baryon() && t2.is_meson()) ||
4771  (t2.is_baryon() && t1.is_meson())) ||
4772  (t1.is_meson() && t2.is_meson()));
4773  const double mass_sum =
4774  incoming_particles_[0].pole_mass() + incoming_particles_[1].pole_mass();
4775 
4776  if (!is_NN_scattering && !is_BBbar_scattering && !is_Npi_scattering &&
4777  !is_AQM_scattering) {
4778  return 0.;
4779  } else if (is_NNbar_pair_ && !treat_BBbar_with_strings) {
4780  return 0.;
4781  } else if (is_BBbar_scattering) {
4782  // BBbar only goes through strings, so there are no "window" considerations
4783  return 1.;
4784  } else {
4785  /* true for K+ p and K0 p (+ antiparticles), which have special treatment
4786  * to fit data */
4787  const PdgCode pdg1 = t1.pdgcode(), pdg2 = t2.pdgcode();
4788  const bool is_KplusP =
4789  ((pdg1 == pdg::K_p || pdg1 == pdg::K_z) && (pdg2 == pdg::p)) ||
4790  ((pdg2 == pdg::K_p || pdg2 == pdg::K_z) && (pdg1 == pdg::p)) ||
4791  ((pdg1 == -pdg::K_p || pdg1 == -pdg::K_z) && (pdg2 == -pdg::p)) ||
4792  ((pdg2 == -pdg::K_p || pdg2 == -pdg::K_z) && (pdg1 == -pdg::p));
4793  // where to start the AQM strings above mass sum
4794  double aqm_offset =
4795  finder_parameters.transition_high_energy.sqrts_add_lower;
4796  if (is_KplusP) {
4797  /* for this specific case we have data. This corresponds to the point
4798  * where the AQM parametrization is smaller than the current 2to2
4799  * parametrization, which starts growing and diverges from exp. data */
4800  aqm_offset = finder_parameters.transition_high_energy.KN_offset;
4801  } else if (pdg1.is_pion() && pdg2.is_pion()) {
4802  aqm_offset = finder_parameters.transition_high_energy.pipi_offset;
4803  }
4804  /* if we do not use the probability transition algorithm, this is always a
4805  * string contribution if the energy is large enough */
4806  if (!finder_parameters.strings_with_probability) {
4807  return static_cast<double>(sqrt_s_ > mass_sum + aqm_offset);
4808  }
4809  /* No strings at low energy, only strings at high energy and
4810  * a transition region in the middle. Determine transition region: */
4811  double region_lower, region_upper;
4812  if (is_Npi_scattering) {
4813  std::tie(region_lower, region_upper) =
4814  finder_parameters.transition_high_energy.sqrts_range_Npi;
4815  } else if (is_NN_scattering) {
4816  std::tie(region_lower, region_upper) =
4817  finder_parameters.transition_high_energy.sqrts_range_NN;
4818  } else { // AQM - Additive Quark Model
4819  /* Transition region around 0.9 larger than the sum of pole masses;
4820  * highly arbitrary, feel free to improve */
4821  region_lower = mass_sum + aqm_offset;
4822  region_upper = mass_sum + aqm_offset +
4823  finder_parameters.transition_high_energy.sqrts_range_width;
4824  }
4825 
4826  return transition_probability_at_sqrts(region_lower, region_upper);
4827  }
4828 }
4829 
4831  double region_lower, double region_upper) const {
4832  if (sqrt_s_ < region_lower) {
4833  return 0.;
4834  } else if (sqrt_s_ > region_upper) {
4835  return 1.;
4836  }
4837 
4838  /* Map sqrt_s_ from [region_lower, region_upper] to [-0.5, 0.5] that
4839  * sin(pi * x) goes from -1 to 1 leading to a probability within 0 and 1. */
4840  const double x = (sqrt_s_ - 0.5 * (region_lower + region_upper)) /
4841  (region_upper - region_lower);
4842  assert(x >= -0.5 && x <= 0.5);
4843  double prob = 0.5 * (std::sin(M_PI * x) + 1.0);
4844  assert(prob >= 0. && prob <= 1.);
4845 
4846  return prob;
4847 }
4848 } // namespace smash
const double sqrt_s_
Total energy in the center-of-mass frame.
CollisionBranchList two_to_two(const ReactionsBitSet &included_2to2, double KN_offset, CharmRescattering charm_rescattering) const
Find all inelastic 2->2 processes for the given scattering.
double transition_probability_at_sqrts(double region_lower, double region_upper) const
Computes a smooth transition probability as a function of sqrt(s).
CollisionBranchList DK_and_DstarK_xx(const ReactionsBitSet &included_2to2, CharmRescattering charm_rescattering) const
Find all inelastic 2->2 processes for D meson-kaon (DK) and D*-kaon (D*K) scattering.
static double xs_dpi_dprimepi(double sqrts, double cm_mom, ParticleTypePtr produced_nucleus, const ParticleType &type_pi)
Parametrized cross section for πd→ πd' (mockup for πd→ πnp), πd̅→ πd̅' and reverse,...
double DK_and_DstarK_inelastic() const
Determine the inelastic cross section for a D meson-kaon (DK) or a D*-kaon (D*K) collision.
CollisionBranchPtr NNbar_to_5pi(double scale_xs) const
Create collision branch for NNbar annihilation going directly into 5 pions.
std::optional< double > DDelta_elastic() const
Determine the elastic cross section for a D meson-Delta (DΔ) collision, If the center-of-mass energy ...
CollisionBranchList NNbar_creation() const
Determine the cross section for NNbar creation, which is given by detailed balance from the reverse r...
CollisionBranchList dpi_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes involving Pion and (anti-) Deuteron (dpi), specifically dπ→ NN,...
std::optional< double > Dpi_and_Dstarpi_elastic() const
Determine the elastic cross section for a D meson-pion (Dpi) or a D*-pion (D*pi) collision.
double high_energy(const ScatterActionsFinderParameters &finder_parameters) const
Determine the parametrized total cross section at high energies for the given collision,...
CollisionBranchList bb_xx_except_nn(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes for Baryon-Baryon (BB) Scattering except the more specific Nucleon-...
CollisionBranchList npi_yk() const
Find all processes for Nucleon-Pion to Hyperon-Kaon Scattering.
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().
double DDelta_inelastic() const
Determine the inelastic cross section for a D meson-Delta (DΔ) collision.
double cm_momentum() const
Determine the momenta of the incoming particles in the center-of-mass system.
double string_probability(const ScatterActionsFinderParameters &finder_parameters) const
CollisionBranchList deltak_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes for Delta-Kaon (DeltaK) Scattering.
static double nn_to_resonance_matrix_element(double sqrts, const ParticleType &type_a, const ParticleType &type_b, int twoI)
Scattering matrix amplitude squared (divided by 16π) for resonance production processes like NN → NR ...
double DN_inelastic() const
Determine the inelastic cross section for a D meson-nucleon (DN) collision.
std::optional< double > DK_and_DstarK_elastic() const
Determine the elastic cross section for a D meson-kaon (DK) or a D*-kaon (D*K) collision.
double Dpi_and_Dstarpi_inelastic() const
Determine the inelastic cross section for a D meson-pion (Dpi) or a D*-pion (D*pi) collision.
CrossSections(const ParticleList &incoming_particles, double sqrt_s, const std::pair< FourVector, FourVector > potentials)
Construct CrossSections instance.
CollisionBranchList Dpi_and_Dstarpi_xx(const ReactionsBitSet &included_2to2, CharmRescattering charm_rescattering) const
Find all inelastic 2->2 processes for D meson-pion (Dpi) and D*-pion (D*pi) scattering.
static double xs_dn_dprimen(double sqrts, double cm_mom, ParticleTypePtr produced_nucleus, const ParticleType &type_nucleus, const ParticleType &type_N)
Parametrized cross section for Nd → Nd', N̅d → N̅d', N̅d̅→ N̅d̅', Nd̅→ Nd̅' and reverse (e....
double elastic_parametrization(const ScatterActionsFinderParameters &finder_parameters) const
Choose the appropriate parametrizations for given incoming particles and return the (parametrized) el...
CollisionBranchList DN_xx(const ReactionsBitSet &included_2to2, CharmRescattering charm_rescattering) const
Find all inelastic 2->2 processes for D meson-nucleon (DN) scatterings.
CollisionBranchList two_to_four() const
Find all 2->4 processes for the given scattering.
CollisionBranchPtr elastic(const ScatterActionsFinderParameters &finder_parameters) const
Determine the elastic cross section for this collision.
double formation(const ParticleType &type_resonance, double cm_momentum_sqr) const
Calculates the 2-to-1 resonance production cross section for a given resonance using the Breit-Wigner...
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.
double string_hard_cross_section() const
Determine the (parametrized) hard non-diffractive string cross section for this collision.
static double two_to_three_xs(const ParticleType &type_in1, const ParticleType &type_in2, double sqrts)
Determine 2->3 cross section for the scattering of the given particle types.
double nk_el() const
Determine the elastic cross section for a nucleon-kaon (NK) collision.
CollisionBranchList nk_xx(const ReactionsBitSet &included_2to2, double KN_offset) const
Find all inelastic 2->2 background processes for Nucleon-Kaon (NK) Scattering.
const ParticleList incoming_particles_
List with data of scattering particles.
double npi_el() const
Determine the elastic cross section for a nucleon-pion (Npi) collision.
double nn_el() const
Determine the (parametrized) elastic cross section for a nucleon-nucleon (NN) collision.
const bool is_BBbar_pair_
Whether incoming particles are a pair of a baryon and an antibaryon (could be different baryon types)
CollisionBranchList dn_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes involving Nucleon and (anti-) Deuteron (dN), specifically Nd → Nd',...
static double two_to_four_xs(const ParticleType &type_in1, const ParticleType &type_in2, double sqrts)
Determine 2->4 cross section for the scattering of the given particle types.
CollisionBranchList ypi_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes for Hyperon-Pion (Ypi) Scattering.
CollisionBranchList DDelta_xx(const ReactionsBitSet &included_2to2, CharmRescattering charm_rescattering) const
Find all inelastic 2->2 processes for D meson-Delta (DΔ) scatterings.
CollisionBranchList string_excitation(double total_string_xs, StringProcess *string_process, const ScatterActionsFinderParameters &finder_parameters) const
Determine the cross section for string excitations, which is given by the difference between the para...
CollisionBranchList two_to_one(CharmRescattering charm_rescattering) const
Find all resonances that can be produced in a 2->1 collision of the two input particles and the produ...
double parametrized_total(const ScatterActionsFinderParameters &finder_parameters) const
Select the parametrization for the total cross section, given the types of incoming particles.
const bool is_NNbar_pair_
Whether incoming particles are a nulecon-antinucleon pair (same isospin)
CollisionBranchList rare_two_to_two() const
Find all 2->2 processes which are suppressed at high energies when strings are turned on with probabi...
static double sum_xs_of(const CollisionBranchList &list)
Helper function: Sum all cross sections of the given process list.
Definition: crosssections.h:77
CollisionBranchList two_to_three() const
Find all 2->3 processes for the given scattering.
CollisionBranchList bar_bar_to_nuc_nuc(bool is_anti_particles) const
Calculate cross sections for 2 → 2 resonance absorption (i.e.
CollisionBranchPtr NNbar_annihilation(double current_xs, double scale_xs) const
Determine the cross section for NNbar annihilation, which is given by the difference between the para...
CollisionBranchList generate_collision_list(const ScatterActionsFinderParameters &finder_parameters, StringProcess *string_process) const
Generate a list of all possible collisions between the incoming particles with the given c....
std::optional< double > DN_elastic() const
Determine the elastic cross section for a D meson-nucleon (DN) collision, If the center-of-mass energ...
std::optional< double > Deta_and_Dstareta_elastic() const
Determine the elastic cross section for a D meson-eta (Deta) or a D*-eta (D*eta) collision.
CollisionBranchList nn_xx(const ReactionsBitSet &included_2to2) const
Find all inelastic 2->2 processes for Nucelon-Nucelon Scattering.
Range of total isospin for reaction of particle a with particle b.
Definition: clebschgordan.h:68
double get_integral_RR(IsoParticleType *type_res_2, double sqrts)
Look up the tabulated resonance integral for the XX -> RR cross section.
double get_integral_NR(double sqrts)
Look up the tabulated resonance integral for the XX -> NR cross section.
double get_integral_RK(double sqrts)
Look up the tabulated resonance integral for the XX -> RK cross section.
double get_integral_piR(double sqrts)
Look up the tabulated resonance integral for the XX -> piR cross section.
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.
ParticleData contains the dynamic information of a certain particle.
Definition: particledata.h:59
PdgCode pdgcode() const
Get the pdgcode of the particle.
Definition: particledata.h:88
const ParticleType & type() const
Get the type of the particle.
Definition: particledata.h:132
bool is_baryon() const
Definition: particledata.h:95
double effective_mass() const
Get the particle's effective mass.
Definition: particledata.cc:25
A pointer-like interface to global references to ParticleType objects.
Definition: particletype.h:731
Particle type contains the static properties of a particle species.
Definition: particletype.h:100
bool is_pion() const
Definition: particletype.h:221
bool is_baryon() const
Definition: particletype.h:206
bool is_Nstar() const
Definition: particletype.h:233
static const ParticleTypePtr try_find(PdgCode pdgcode)
Returns the ParticleTypePtr for the given pdgcode.
Definition: particletype.cc:89
static const ParticleType & find(PdgCode pdgcode)
Returns the ParticleType object for the given pdgcode.
Definition: particletype.cc:99
bool is_nucleus() const
Definition: particletype.h:254
PdgCode pdgcode() const
Definition: particletype.h:159
double full_spectral_function(double m) const
Full spectral function of the resonance (relativistic Breit-Wigner distribution with mass-dependent ...
const std::string & name() const
Definition: particletype.h:144
int32_t charge() const
The charge of the particle.
Definition: particletype.h:191
int antiparticle_sign() const
Definition: particletype.h:168
static ParticleTypePtrList & list_nucleons()
Definition: particletype.cc:69
static ParticleTypePtrList & list_anti_nucleons()
Definition: particletype.cc:71
bool is_Nstar1535() const
Definition: particletype.h:239
bool is_stable() const
Definition: particletype.h:251
bool is_dprime() const
Definition: particletype.h:263
bool is_Delta() const
Definition: particletype.h:227
double width_at_pole() const
Definition: particletype.h:153
static ParticleTypePtrList & list_anti_Deltas()
Definition: particletype.cc:77
bool is_nucleon() const
Definition: particletype.h:218
double mass() const
Definition: particletype.h:147
static ParticleTypePtrList & list_baryon_resonances()
Definition: particletype.cc:81
static ParticleTypePtrList & list_Deltas()
Definition: particletype.cc:75
bool is_deuteron() const
Definition: particletype.h:257
bool is_meson() const
Definition: particletype.h:209
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 ...
unsigned int spin() const
Definition: particletype.h:194
int baryon_number() const
Definition: particletype.h:212
bool is_Deltastar() const
Definition: particletype.h:245
IsoParticleType * iso_multiplet() const
Definition: particletype.h:188
static ParticleTypePtrList & list_light_nuclei()
Definition: particletype.cc:85
PdgCode stores a Particle Data Group Particle Numbering Scheme particle type number.
Definition: pdgcode.h:108
std::int32_t code() const
Definition: pdgcode.h:319
bool is_Dmeson() const
Definition: pdgcode.h:489
int antiparticle_sign() const
Definition: pdgcode.h:734
int nucleus_n() const
Number of neutrons in nucleus.
Definition: pdgcode.h:919
bool is_meson() const
Definition: pdgcode.h:401
bool is_Dstar2007() const
Definition: pdgcode.h:498
int nucleus_ap() const
Number of antiprotons in nucleus.
Definition: pdgcode.h:929
int nucleus_an() const
Number of antineutrons in nucleus.
Definition: pdgcode.h:933
unsigned int spin() const
Definition: pdgcode.h:691
bool is_pion() const
Definition: pdgcode.h:471
bool is_kaon() const
Definition: pdgcode.h:465
bool is_hyperon() const
Definition: pdgcode.h:435
bool is_nucleus() const
Definition: pdgcode.h:361
bool is_deuteron() const
Definition: pdgcode.h:504
bool is_antiparticle_of(const PdgCode rhs) const
Definition: pdgcode.h:829
int nucleus_p() const
Number of protons in nucleus.
Definition: pdgcode.h:915
bool is_baryon() const
Definition: pdgcode.h:398
int isospin3() const
Definition: pdgcode.h:535
bool is_nucleon() const
Definition: pdgcode.h:404
int nucleus_La() const
Number of Lambdas in nucleus.
Definition: pdgcode.h:925
bool is_eta() const
Definition: pdgcode.h:477
int nucleus_A() const
Nucleus mass number.
Definition: pdgcode.h:943
int nucleus_aLa() const
Number of anti-Lambdas in nucleus.
Definition: pdgcode.h:939
bool is_Delta() const
Definition: pdgcode.h:428
Helper class for ScatterActionsFinder.
const bool strings_with_probability
This indicates whether the string fragmentation is swiched on with a probability smoothly increasing ...
const std::pair< double, double > hard_string_transition_energy_range
Invariant energy range for the soft-to-hard string transition (in measured in GeV).
const ReactionsBitSet included_2to2
List of included 2<->2 reactions.
const bool use_AQM
Switch to control whether to use AQM or not.
const double scale_xs
Factor by which all (partial) cross sections are scaled.
const double elastic_parameter
Elastic cross section parameter (in mb).
const bool strings_switch
Indicates whether string fragmentation is switched on.
const StringTransitionParameters transition_high_energy
Constants related to transition between low collision energies - mediated via resonances - and high c...
const double additional_el_xs
Additional constant contribution (in mb) to the elastic cross sections.
const HardStringTransitionMode hard_string_transition_mode
Mode used to control the transition from soft to hard string excitation.
const CharmRescattering charm_rescattering
Specifies kind of charm rescattering.
double AQM_scaling_factor(const PdgCode &pdg) const
AQM scaling factor for a hadron.
const NNbarTreatment nnbar_treatment
Switch for NNbar reactions.
const double low_snn_cut
Elastic collsions between two nucleons with sqrt_s below low_snn_cut_ are excluded.
const bool two_to_one
Enables resonance production.
String excitation processes used in SMASH.
Definition: stringprocess.h:46
static int pdg_map_for_pythia(PdgCode &pdg)
Take pdg code and map onto particle specie which can be handled by PYTHIA.
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...
Collection of useful constants that are known at compile time.
@ TwoToFive
Directly create 5 pions, use with multi-particle reactions.
@ Resonances
Use intermediate Resonances.
@ Strings
Use string fragmentation.
CharmRescattering
Possible charm scattering options.
@ T_Matrix
Charm interactions via T-matrix approach.
@ Resonances
Charm interactions via resonances.
@ Custom_Range
Smooth transition within a user-defined invariant energy range.
@ A3_Nuclei_4to2
@ Deuteron_3to2
@ KN_to_KDelta
@ KN_to_KN
@ NN_to_NR
@ PiDeuteron_to_pidprime
@ NDeuteron_to_Ndprime
@ Strangeness_exchange
@ Charm_T_matrix
@ PiDeuteron_to_NN
@ NN_to_DR
std::bitset< 11 > ReactionsBitSet
Container for the 2 to 2 reactions in the code.
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > & logg
An array that stores all pre-configured Logger objects.
Definition: logging.h:245
constexpr Section potentials
Section for the potentials information.
Definition: input_keys.h:228
constexpr int pi_p
π⁺.
constexpr int Delta_p
Δ⁺.
constexpr int Delta_pp
Δ⁺⁺.
constexpr int64_t antideuteron
Anti-deuteron in decimal digits.
constexpr int Sigma_m
Σ⁻.
constexpr int Dstar_p
D*(2010)⁺.
constexpr int D_z
D⁰.
constexpr int Dbar_z
D̄⁰.
constexpr int K_p
K⁺.
constexpr int K_z
K⁰.
constexpr int p
Proton.
constexpr int h1
h₁(1170).
constexpr int K_m
K̄⁻.
constexpr int eta
η.
constexpr int Dstarbar_z
D̄*(2007)⁰.
constexpr int pi_z
π⁰.
constexpr int n
Neutron.
constexpr int64_t deuteron
Deuteron.
constexpr int Delta_m
Δ⁻.
constexpr int Delta_z
Δ⁰.
constexpr int rho_z
ρ⁰.
constexpr int Lambda
Λ.
constexpr int pi_m
π⁻.
constexpr int Dstar_m
D*(2010)⁻.
constexpr int Sigma_p
Σ⁺.
constexpr int Dstar_z
D*(2007)⁰.
constexpr int Kbar_z
K̄⁰.
constexpr int Sigma_z
Σ⁰.
constexpr int D_m
D⁻.
constexpr int D_p
D⁺.
Definition: action.h:24
double kplusn_k0p(double mandelstam_s)
K+ n charge exchange cross section parametrization.
std::optional< double > Dzeron_elastic(double sqrts)
D⁰n elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > Dplusn_elastic(double sqrts)
D⁺n elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double kplusp_total(double mandelstam_s)
K+ p total cross section parametrization.
std::optional< double > DplusKzero_elastic(double sqrts)
D⁺K⁰ elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double deuteron_pion_inelastic(double pion_kinetic_energy)
Parametrization of deuteron-pion inelastic cross section.
double DzeroDeltazero_DplusDeltaminus(double sqrts)
D⁰Δ⁰ -> D⁺Δ⁻ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DstarzeroKzero_elastic(double sqrts)
D*(2007)⁰K⁰ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
double pizeropizero_total(double sqrts)
pi0 pi0 total cross section parametrized from PDG2018, smoothed using the LOWESS algorithm.
double kminusp_pi0lambda(double sqrts)
K- p <-> pi0 Lambda cross section parametrization Fit to Landolt-Börnstein instead of UrQMD values.
double Dpluspiminus_Dzeropizero(double sqrts)
D⁺π⁻ -> D⁰π⁰ cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
double Dstarzeropiplus_Dstarpluspizero(double sqrts)
D*(2007)⁰π⁺ -> D*(2010)⁺π⁰ cross section (closest reference Song:2015sfa , data provided by Juan Torr...
double DplusDeltaplus_DzeroDeltaplusplus(double sqrts)
D⁺Δ⁺ -> D⁰Δ⁺⁺ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double Dminusp_Dbarzeron(double sqrts)
D⁻p -> D̄⁰n cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
static void throw_xsec_is_not_implemented(const ParticleData &data_a, const ParticleData &data_b, const std::string func_name)
Helper function: Throw if cross section between two particles is not implemented.
double pipluspiminus_total(double sqrts)
pi+ pi- total cross section parametrized from PDG2018, smoothed using the LOWESS algorithm.
ParticleTypePtrList list_possible_resonances(const ParticleTypePtr type_a, const ParticleTypePtr type_b)
Lists the possible resonances that decay into two particles.
double piminusp_sigma0k0_res(double mandelstam_s)
pi- p -> Sigma0 K0 cross section parametrization, resonance contribution.
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:91
double ppbar_total(double mandelstam_s)
ppbar total cross section parametrization Source: Bass:1998ca
double np_total(double mandelstam_s)
np total cross section parametrization Sources: low-p: Cugnon:1996kh highest-p: Buss:2011mx
double Dstarzeropizero_Dstarpluspiminus(double sqrts)
D*(2007)⁰π⁰ -> D*(2010)⁺π⁻ cross section (closest reference Song:2015sfa , data provided by Juan Torr...
double DstarplusKzero_DstarzeroKplus(double sqrts)
D*(2010)⁺K⁰ -> D*(2007)⁰K⁺ cross section (closest reference Song:2015sfa , data provided by Juan Torr...
double DstarplusKminus_DstarzeroKbarzero(double sqrts)
D*(2010)⁺K⁻ -> D*(2007)⁰K̄⁰ cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double piminusp_elastic(double mandelstam_s)
pi-p elastic cross section parametrization Source: GiBUU:parametrizationBarMes_HighEnergy....
std::optional< double > Dbarzeron_elastic(double sqrts)
D̄⁰n elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DzeroDeltaplus_elastic(double sqrts)
D⁰Δ⁺ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double Dpluspizero_Dzeropiplus(double sqrts)
D⁺π⁰ -> D⁰π⁺ cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
std::vector< std::string > split(const std::string &s, char delim)
Split string by delimiter.
double npbar_high_energy(double mandelstam_s)
npbar total cross section at high energies
double Dstarpluspiminus_Dstarzeropizero(double sqrts)
D*(2010)⁺π⁻ -> D*(2007)⁰π⁰ cross section (closest reference Song:2015sfa , data provided by Juan Torr...
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.
double kminusn_piminussigma0(double sqrts)
K- n <-> pi- Sigma0 cross section parametrization Follow from the parametrization with the same stran...
double DzeroKbarzero_DplusKminus(double sqrts)
D⁰K̄⁰ -> D⁺K⁻ cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double kbar0p_elastic_background(double mandelstam_s)
Kbar0 p elastic background cross section parametrization Source: Buss:2011mx , B.3....
std::optional< double > DplusKplus_elastic(double sqrts)
D⁺K⁺ elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double DstarzeroKbarzero_DstarplusKminus(double sqrts)
D*(2007)⁰K̄⁰ -> D*(2010)⁺K⁻ cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
std::optional< double > DplusDeltazero_elastic(double sqrts)
D⁺Δ⁰ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
KaonNucleonRatios kaon_nucleon_ratios
double ppbar_elastic(double mandelstam_s)
ppbar elastic cross section parametrization Source: Bass:1998ca
std::optional< double > Dbarzerop_elastic(double sqrts)
D̄⁰p elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DzeroKplus_elastic(double sqrts)
D⁰K⁺ elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
std::optional< double > Dstarzeropiminus_elastic(double sqrts)
D*(2007)⁰π- elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
std::optional< double > Dstarzeroeta_elastic(double sqrts)
D*(2007)⁰η elastic cross section (data provided by Juan Torres-Rincon).
std::optional< double > Dstarpluspizero_elastic(double sqrts)
D*(2010)⁺π⁰ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
double kminusp_elastic_background(double mandelstam_s)
K- p elastic background cross section parametrization Source: Buss:2011mx , B.3.9.
double np_high_energy(double mandelstam_s)
np total cross section at high energies
std::optional< double > Dminusn_elastic(double sqrts)
D⁻n elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
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...
double DstarzeroKplus_DstarplusKzero(double sqrts)
D*(2007)⁰K⁺ -> D*(2010)⁺K⁰ cross section (closest reference Song:2015sfa , data provided by Juan Torr...
double Npi_string_hard(double mandelstam_s)
nucleon-pion hard scattering cross section (with partonic scattering)
double Dstarpluspizero_Dstarzeropiplus(double sqrts)
D*(2010)⁺π⁰ -> D*(2007)⁰π⁺ cross section (closest reference Song:2015sfa , data provided by Juan Torr...
std::optional< double > DzeroDeltaminus_elastic(double sqrts)
D⁰Δ⁻ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double kminusn_piminuslambda(double sqrts)
K- n <-> pi- Lambda cross section parametrization Follow from the parametrization with the same stran...
std::optional< double > Dpluspiplus_elastic(double sqrts)
D⁺π⁺ elastic cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
std::optional< double > DplusKbarzero_elastic(double sqrts)
D⁺K̄⁰ elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
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.
@ TwoToOne
See here for a short description.
@ StringHardSingleDiffractiveAX
See here for a short description.
@ StringSoftDoubleDiffractive
See here for a short description.
@ TwoToFive
See here for a short description.
@ StringSoftSingleDiffractiveXB
See here for a short description.
@ TwoToTwo
See here for a short description.
@ Elastic
See here for a short description.
@ TwoToFour
See here for a short description.
@ StringHardNonDiffractive
See here for a short description.
@ StringSoftAnnihilation
See here for a short description.
@ StringSoftNonDiffractive
See here for a short description.
@ StringSoftSingleDiffractiveAX
See here for a short description.
@ StringHardSingleDiffractiveXB
See here for a short description.
@ StringHardDoubleDiffractive
See here for a short description.
@ TwoToThree
See here for a short description.
constexpr double minimum_sqrts_pythia_can_handle
Energy in GeV, below which hard reactions via pythia are impossible.
Definition: constants.h:122
std::optional< double > Dplusp_elastic(double sqrts)
D⁺p elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double ppbar_high_energy(double mandelstam_s)
ppbar total cross section at high energies
std::optional< double > DstarzeroKplus_elastic(double sqrts)
D*(2007)⁰K⁺ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
double pp_high_energy(double mandelstam_s)
pp total cross section at high energies
std::optional< double > Dminusp_elastic(double sqrts)
D⁻p elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double Dbarzeron_Dminusp(double sqrts)
D̄⁰n -> D⁻p cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DzeroKzero_elastic(double sqrts)
D⁰K⁰ elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
T pCM_sqr_from_s(const T s, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:52
double pipi_string_hard(double mandelstam_s)
pion-pion hard scattering cross section (with partonic scattering)
std::optional< double > Dstarpluspiplus_elastic(double sqrts)
D*(2010)⁺π⁺ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
double piplusp_high_energy(double mandelstam_s)
pi+p total cross section at high energies
std::optional< double > Dpluseta_elastic(double sqrts)
D⁺η elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double kminusp_piminussigmaplus(double sqrts)
K- p <-> pi- Sigma+ cross section parametrization Taken from UrQMD (Graef:2014mra ).
std::optional< double > Dpluspiminus_elastic(double sqrts)
D⁺π⁻ elastic cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
double DplusKminus_DzeroKbarzero(double sqrts)
D⁺K⁻ -> D⁰K̄⁰ cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double DbarzeroDeltaminus_DminusDeltazero(double sqrts)
D̄⁰Δ⁻ -> D⁻Δ⁰ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double DminusDeltaplusplus_DbarzeroDeltaplus(double sqrts)
D⁻Δ⁺⁺ -> D̄⁰Δ⁺ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DminusDeltazero_elastic(double sqrts)
D⁻Δ⁰ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double piplusp_sigmapluskplus_pdg(double mandelstam_s)
pi+ p to Sigma+ K+ cross section parametrization, PDG data.
static void throw_xsec_is_negative(const double sqrts, const double xsec, const ParticleData &data_a, const ParticleData &data_b, std::string func_name)
Helper function: Throw if cross section is negative.
double DplusKzero_DzeroKplus(double sqrts)
D⁺K⁰ -> D⁰K⁺ cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
std::optional< double > DzeroDeltaplusplus_elastic(double sqrts)
D⁰Δ⁺⁺ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DminusDeltaplus_elastic(double sqrts)
D⁻Δ⁺ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double Dplusn_Dzerop(double sqrts)
D⁺n -> D⁰p cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DstarzeroKminus_elastic(double sqrts)
D*(2007)⁰K⁻ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
std::optional< double > Dzeropizero_elastic(double sqrts)
D⁰π⁰ elastic cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
static double AQM_based_on_piminusp_high_energy(const double sqrts, const PdgCode &pdg_a, const PdgCode &pdg_b, const double AQM_scaling_factor_a, const double AQM_scaling_factor_b)
Helper function: Approximate cross section using AQM based on function piminusp_high_energy.
static constexpr int LCrossSections
double piminusp_total(double sqrts)
pi- p total cross section parametrized from PDG2018, smoothed using the LOWESS algorithm.
static void warn_if_charm_rescattering_enabled_and_AQM_disabled(const double sqrts, const ParticleType &type_a, const ParticleType &type_b, const CharmRescattering charm_rescattering)
Helper function: Print a warning message if Charm_Rescattering_Method is not set to none and AQM is d...
constexpr double deuteron_mass
Deuteron mass in GeV.
Definition: constants.h:103
double DzeroDeltaplus_DplusDeltazero(double sqrts)
D⁰Δ⁺ -> D⁺Δ⁰ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double deuteron_nucleon_elastic(double mandelstam_s)
Deuteron nucleon elastic cross-section [mb] parametrized by Oh:2009gx .
constexpr double nucleon_mass
Nucleon mass in GeV.
Definition: constants.h:69
constexpr T pow_int(const T base, unsigned const exponent)
Efficient template for calculating integer powers using squaring.
Definition: pow.h:23
std::optional< double > DplusDeltaminus_elastic(double sqrts)
D⁺Δ⁻ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DzeroKminus_elastic(double sqrts)
D⁰K⁻ elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double piminusp_sigmaminuskplus_pdg(double mandelstam_s)
pi- p -> Sigma- K+ cross section parametrization, PDG data.
double piminusp_lambdak0_pdg(double mandelstam_s)
pi- p -> Lambda K0 cross section parametrization, PDG data.
std::optional< double > Dzeropiminus_elastic(double sqrts)
D⁰π⁻ elastic cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
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.
std::optional< double > Dstarzeropiplus_elastic(double sqrts)
D*(2007)⁰π⁺ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
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.
std::optional< double > DplusDeltaplusplus_elastic(double sqrts)
D⁺Δ⁺⁺ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DminusDeltaplusplus_elastic(double sqrts)
D⁻Δ⁺⁺ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double DplusDeltaminus_DzeroDeltazero(double sqrts)
D⁺Δ⁻ -> D⁰Δ⁰ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DbarzeroDeltaplusplus_elastic(double sqrts)
D̄⁰Δ⁺⁺ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double deuteron_antinucleon_inelastic(double aN_kinetic_energy)
Parametrization of deuteron-antinucleon inelastic cross section.
double k0p_elastic_background(double mandelstam_s)
K0 p elastic background cross section parametrization Source: Buss:2011mx , B.3.9.
std::optional< double > Dpluspizero_elastic(double sqrts)
D⁺π⁰ elastic cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
constexpr uint64_t pack(int32_t x, int32_t y)
Pack two int32_t into an uint64_t.
double deuteron_pion_elastic(double mandelstam_s)
Deuteron pion elastic cross-section [mb] parametrized to fit pi-d elastic scattering data (the data c...
std::optional< double > DzeroKbarzero_elastic(double sqrts)
D⁰K̄⁰ elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double DminusDeltazero_DbarzeroDeltaminus(double sqrts)
D⁻Δ⁰ -> D̄⁰Δ⁻ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double NN_string_hard(double mandelstam_s)
nucleon-nucleon hard scattering cross section (with partonic scattering)
double Dzerop_Dplusn(double sqrts)
D⁰p -> D⁺n cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double xs_ppbar_annihilation(double mandelstam_s)
parametrized cross-section for proton-antiproton annihilation used in the UrQMD model
double DbarzeroDeltaplus_DminusDeltaplusplus(double sqrts)
D̄⁰Δ⁺ -> D⁻Δ⁺⁺ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DplusDeltaplus_elastic(double sqrts)
D⁺Δ⁺ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double kplusp_inelastic_background(double mandelstam_s)
K+ p inelastic background cross section parametrization Source: Buss:2011mx , B.3....
std::optional< double > Dzerop_elastic(double sqrts)
D⁰p elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double DminusDeltaplus_DbarzeroDeltazero(double sqrts)
D⁻Δ⁺ -> D̄⁰Δ⁰ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
constexpr double pion_mass
Pion mass in GeV.
Definition: constants.h:76
constexpr double hbarc
GeV <-> fm conversion factor.
Definition: constants.h:29
double kminusp_pi0sigma0(double sqrts)
K- p <-> pi0 Sigma0 cross section parametrization Fit to Landolt-Börnstein instead of UrQMD values.
double kplusn_elastic_background(double mandelstam_s)
K+ n elastic background cross section parametrization sigma(K+n->K+n) = sigma(K+n->K0p) = 0....
double pp_total(double mandelstam_s)
pp total cross section parametrization Sources: low-p: Cugnon:1996kh highest-p: Buss:2011mx
double DplusDeltazero_DzeroDeltaplus(double sqrts)
D⁺Δ⁰ -> D⁰Δ⁺ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
T pCM_from_s(const T s, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:66
std::optional< double > Dstarpluseta_elastic(double sqrts)
D*(2010)⁺η elastic cross section (data provided by Juan Torres-Rincon).
std::optional< double > DplusKminus_elastic(double sqrts)
D⁺K⁻ elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double kplusn_total(double mandelstam_s)
K+ n total cross section parametrization.
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:41
std::optional< double > DstarplusKzero_elastic(double sqrts)
D*(2010)⁺K⁰ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
double np_elastic(double mandelstam_s)
np elastic cross section parametrization Source: Weil:2013mya , eq.
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.
std::optional< double > DbarzeroDeltaminus_elastic(double sqrts)
D̄⁰Δ⁻ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > DbarzeroDeltazero_elastic(double sqrts)
D̄⁰Δ⁰ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double k0n_elastic_background(double mandelstam_s)
K0 n elastic background cross section parametrization Source: Buss:2011mx , B.3.9.
static constexpr int LScatterAction
double kminusp_piplussigmaminus(double sqrts)
K- p <-> pi+ Sigma- cross section parametrization Taken from UrQMD (Graef:2014mra ).
static double effective_AQM_s(const double mandelstam_s, const double m1, const double m2, const double m1_ref, const double m2_ref)
Helper function: Shift the energy of a collision for AQM rescaled cross sections.
double kminusp_total(double mandelstam_s)
K- p total cross section parametrization.
double Dzeropiplus_Dpluspizero(double sqrts)
D⁰π⁺ -> D⁺π⁰ cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
double kbar0n_elastic_background(double mandelstam_s)
Kbar0 n elastic background cross section parametrization Source: Buss:2011mx , B.3....
double kminusp_kbar0n(double mandelstam_s)
K- p <-> Kbar0 n cross section parametrization.
std::optional< double > DstarplusKplus_elastic(double sqrts)
D*(2010)⁺K⁺ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
double DzeroKplus_DplusKzero(double sqrts)
D⁰K⁺ -> D⁺K⁰ cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
double deuteron_nucleon_inelastic(double N_kinetic_energy)
Parametrization of deuteron-nucleon inelastic cross section.
std::optional< double > DminusDeltaminus_elastic(double sqrts)
D⁻Δ⁻ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double piplusp_total(double sqrts)
pi+ p total cross section parametrized from PDG2018, smoothed using the LOWESS algorithm.
double piminusp_high_energy(double mandelstam_s)
pi-p total cross section at high energies
std::optional< double > DstarzeroKbarzero_elastic(double sqrts)
D*(2007)⁰K̄⁰ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rin...
double kminusn_total(double mandelstam_s)
K- n total cross section parametrization.
double kminusn_elastic_background(double mandelstam_s)
K- n elastic background cross section parametrization Source: Buss:2011mx , B.3.9.
std::optional< double > Dzeropiplus_elastic(double sqrts)
D⁰π⁺ elastic cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
double pp_elastic(double mandelstam_s)
pp elastic cross section parametrization Source: Weil:2013mya , eq.
double DzeroDeltaplusplus_DplusDeltaplus(double sqrts)
D⁰Δ⁺⁺ -> D⁺Δ⁺ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > Dstarpluspiminus_elastic(double sqrts)
D*(2010)⁺π⁻ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
std::optional< double > DstarplusKbarzero_elastic(double sqrts)
D*(2010)⁺K̄⁰ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rin...
double Dzeropizero_Dpluspiminus(double sqrts)
D⁰π⁰ -> D⁺π⁻ cross section (Abreu:2011ic , data provided by Juan Torres-Rincon).
double kplusn_inelastic_background(double mandelstam_s)
K+ n inelastic background cross section parametrization Source: Buss:2011mx , B.3....
std::optional< double > DstarplusKminus_elastic(double sqrts)
D*(2010)⁺K⁻ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
double DbarzeroDeltazero_DminusDeltaplus(double sqrts)
D̄⁰Δ⁰ -> D⁻Δ⁺ cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double kplusp_elastic_background(double mandelstam_s)
K+ p elastic background cross section parametrization.
std::optional< double > DbarzeroDeltaplus_elastic(double sqrts)
D̄⁰Δ⁺ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
double piplusp_elastic_AQM(double mandelstam_s, double m1, double m2)
pi+p elactic cross section parametrization.
double piplusp_elastic(double mandelstam_s)
pi+p elastic cross section parametrization, PDG data.
std::optional< double > Dstarzeropizero_elastic(double sqrts)
D*(2007)⁰π⁰ elastic cross section (closest reference Song:2015sfa , data provided by Juan Torres-Rinc...
std::optional< double > DzeroDeltazero_elastic(double sqrts)
D⁰Δ⁰ elastic cross section (Tolos:2013kva ), data provided by Juan Torres-Rincon.
std::optional< double > Dzeroeta_elastic(double sqrts)
D⁰η elastic cross section (Tolos:2013kva , data provided by Juan Torres-Rincon).
constexpr double fm2_mb
mb <-> fm^2 conversion factor.
Definition: constants.h:32
const std::pair< double, double > sqrts_range_Npi
Transition range in N collisions.
const double pipi_offset
Constant offset as to where to turn on the strings and elastic processes for reactions (this is an e...
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...
const double KN_offset
Constant offset as to where to shift from 2to2 to string processes (in GeV) in the case of KN reactio...
const std::pair< double, double > sqrts_range_NN
Transition range in NN collisions.
const double sqrts_range_width
Constant for the range of transition region, in the case of AQM this is added to the sum of masses + ...