32 double spin_factor = (c.
spin() + 1) * (d.
spin() + 1);
33 spin_factor /= (a.
spin() + 1) * (b.
spin() + 1);
34 double symmetry_factor = (1 + (a == b));
35 symmetry_factor /= (1 + (c == d));
38 return spin_factor * symmetry_factor * momentum_factor;
54 double spin_factor = (c.
spin() + 1) * (d.
spin() + 1);
55 spin_factor /= (a.
spin() + 1) * (b.
spin() + 1);
56 double symmetry_factor = (1 + (a == b));
57 symmetry_factor /= (1 + (c == d));
58 const double momentum_factor =
61 return spin_factor * symmetry_factor * momentum_factor;
77 double spin_factor = (c.
spin() + 1) * (d.
spin() + 1);
78 spin_factor /= (a.
spin() + 1) * (b.
spin() + 1);
79 double symmetry_factor = (1 + (a == b));
80 symmetry_factor /= (1 + (c == d));
81 const double momentum_factor =
84 return spin_factor * symmetry_factor * momentum_factor;
92 CollisionBranchList in_list,
double weight = 1.) {
93 main_list.reserve(main_list.size() + in_list.size());
94 for (
auto& proc : in_list) {
95 proc->set_weight(proc->weight() * weight);
96 main_list.emplace_back(std::move(proc));
106 for (
auto& proc : list) {
107 xs_sum += proc->weight();
114 const std::pair<FourVector, FourVector> potentials)
115 : incoming_particles_(incoming_particles),
117 potentials_(potentials),
118 is_BBbar_pair_(incoming_particles_[0].type().is_baryon() &&
119 incoming_particles_[1].type().is_baryon() &&
120 incoming_particles_[0].type().antiparticle_sign() ==
121 -incoming_particles_[1].type().antiparticle_sign()) {}
124 double elastic_parameter,
bool two_to_one_switch,
125 ReactionsBitSet included_2to2,
double low_snn_cut,
bool strings_switch,
126 bool use_AQM,
bool strings_with_probability,
NNbarTreatment nnbar_treatment,
128 CollisionBranchList process_list;
132 double p_pythia = 0.;
133 if (strings_with_probability) {
141 const bool reject_by_nucleon_elastic_cutoff =
145 if (incl_elastic && !reject_by_nucleon_elastic_cutoff) {
146 process_list.emplace_back(
elastic(elastic_parameter, use_AQM));
152 const double sig_current =
sum_xs_of(process_list);
153 const double sig_string = std::max(0.,
high_energy() - sig_current);
160 if (two_to_one_switch) {
164 if (included_2to2.any()) {
187 bool use_AQM)
const {
188 double elastic_xs = 0.;
189 if (elast_par >= 0.) {
191 elastic_xs = elast_par;
202 CollisionBranchList process_list;
205 const auto& pdg_a = data_a.
pdgcode();
206 const auto& pdg_b = data_b.
pdgcode();
207 if ((pdg_a.is_nucleon() && pdg_b.is_pion()) ||
208 (pdg_b.is_nucleon() && pdg_a.is_pion())) {
217 double elastic_xs = 0.0;
225 elastic_xs =
nk_el();
229 elastic_xs =
nn_el();
237 const bool is_deuteron =
239 if (is_deuteron && pdg_other.
is_pion()) {
242 }
else if (is_deuteron && pdg_other.
is_nucleon()) {
246 }
else if (use_AQM) {
251 elastic_xs =
nn_el();
296 std::stringstream ss;
299 ss <<
"problem in CrossSections::elastic: a=" << name_a <<
" b=" << name_b
300 <<
" j_a=" << pdg_a.
spin() <<
" j_b=" << pdg_b.
spin()
301 <<
" sigma=" << sig_el <<
" s=" << s;
302 throw std::runtime_error(ss.str());
312 assert(pion != nucleon);
317 switch (nucleon.
code()) {
319 switch (pion.
code()) {
332 switch (pion.
code()) {
345 switch (pion.
code()) {
358 switch (pion.
code()) {
371 throw std::runtime_error(
372 "only the elastic cross section for proton-pion " 379 std::stringstream ss;
382 ss <<
"problem in CrossSections::elastic: a=" << name_a <<
" b=" << name_b
383 <<
" j_a=" << pdg_a.
spin() <<
" j_b=" << pdg_b.
spin()
384 <<
" sigma=" << sig_el <<
" s=" << s;
385 throw std::runtime_error(ss.str());
395 const auto pdg_nucleon = type_nucleon.
pdgcode().
code();
404 CollisionBranchList process_list;
405 switch (pdg_nucleon) {
413 type_K_p, type_Sigma_p);
424 type_K_p, type_Sigma_m);
426 sqrt_s_, type_K_z, type_Sigma_z);
428 sqrt_s_, type_K_z, type_Lambda);
443 sqrt_s_, type_K_p, type_Sigma_z);
445 sqrt_s_, type_K_z, type_Sigma_p);
448 type_K_p, type_Lambda);
464 type_K_z, type_Sigma_p);
466 sqrt_s_, type_K_p, type_Sigma_z);
468 sqrt_s_, type_K_p, type_Lambda);
476 type_K_z, type_Sigma_m);
491 sqrt_s_, type_K_z, type_Sigma_z);
493 sqrt_s_, type_K_p, type_Sigma_m);
496 type_K_z, type_Lambda);
512 type_K_m, type_Sigma_m_bar);
514 sqrt_s_, type_Kbar_z, type_Sigma_z_bar);
516 sqrt_s_, type_Kbar_z, type_Lambda_bar);
524 type_K_m, type_Sigma_p_bar);
539 sqrt_s_, type_K_m, type_Sigma_z_bar);
541 sqrt_s_, type_Kbar_z, type_Sigma_p_bar);
544 type_K_m, type_Lambda_bar);
557 type_Kbar_z, type_Sigma_m_bar);
568 type_Kbar_z, type_Sigma_p_bar);
570 sqrt_s_, type_K_m, type_Sigma_z_bar);
572 sqrt_s_, type_K_m, type_Lambda_bar);
587 sqrt_s_, type_Kbar_z, type_Sigma_z_bar);
589 sqrt_s_, type_K_m, type_Sigma_m_bar);
592 type_Kbar_z, type_Lambda_bar);
609 assert(kaon != nucleon);
614 switch (nucleon.
code()) {
616 switch (kaon.
code()) {
632 switch (kaon.
code()) {
648 switch (kaon.
code()) {
664 switch (kaon.
code()) {
680 throw std::runtime_error(
681 "elastic cross section for antinucleon-kaon " 688 std::stringstream ss;
691 ss <<
"problem in CrossSections::elastic: a=" << name_a <<
" b=" << name_b
692 <<
" j_a=" << pdg_a.
spin() <<
" j_b=" << pdg_b.
spin()
693 <<
" sigma=" << sig_el <<
" s=" << s;
694 throw std::runtime_error(ss.str());
699 const auto& log = logger<LogArea::CrossSections>();
700 CollisionBranchList resonance_process_list;
711 if (type_resonance.is_stable()) {
717 type_resonance.pdgcode() == type_particle_a.
pdgcode()) ||
719 type_resonance.pdgcode() == type_particle_b.
pdgcode())) {
723 double resonance_xsection =
formation(type_resonance, p_cm_sqr);
727 resonance_process_list.push_back(make_unique<CollisionBranch>(
729 log.debug(
"Found resonance: ", type_resonance);
730 log.debug(type_particle_a.
name(), type_particle_b.
name(),
"->",
731 type_resonance.name(),
" at sqrt(s)[GeV] = ",
sqrt_s_,
732 " with xs[mb] = ", resonance_xsection);
735 return resonance_process_list;
739 double cm_momentum_sqr)
const {
743 if (type_resonance.
charge() !=
757 if (partial_width <= 0.) {
762 const double spinfactor =
763 static_cast<double>(type_resonance.
spin() + 1) /
764 ((type_particle_a.
spin() + 1) * (type_particle_b.
spin() + 1));
765 const int sym_factor =
770 return spinfactor * sym_factor * 2. * M_PI * M_PI / cm_momentum_sqr *
777 CollisionBranchList process_list;
782 const auto& pdg_a = data_a.
pdgcode();
783 const auto& pdg_b = data_b.
pdgcode();
785 if (pdg_a.is_nucleon() && pdg_b.is_nucleon() &&
786 pdg_a.antiparticle_sign() == pdg_b.antiparticle_sign()) {
788 process_list =
nn_xx(included_2to2);
796 if ((pdg_a.is_nucleon() && pdg_b.is_kaon()) ||
797 (pdg_b.is_nucleon() && pdg_a.is_kaon())) {
799 process_list =
nk_xx(included_2to2);
800 }
else if ((pdg_a.is_hyperon() && pdg_b.is_pion()) ||
801 (pdg_b.is_hyperon() && pdg_a.is_pion())) {
803 process_list =
ypi_xx(included_2to2);
804 }
else if ((pdg_a.is_Delta() && pdg_b.is_kaon()) ||
805 (pdg_b.is_Delta() && pdg_a.is_kaon())) {
813 process_list =
dn_xx(included_2to2);
819 process_list =
dpi_xx(included_2to2);
827 CollisionBranchList process_list;
833 if (!same_sign && !any_nucleus) {
853 CollisionBranchList process_list, channel_list;
859 bool both_antinucleons =
862 const ParticleTypePtrList& nuc_or_anti_nuc =
865 const ParticleTypePtrList& delta_or_anti_delta =
875 process_list.reserve(process_list.size() + channel_list.size());
876 std::move(channel_list.begin(), channel_list.end(),
877 std::inserter(process_list, process_list.end()));
878 channel_list.clear();
889 process_list.reserve(process_list.size() + channel_list.size());
890 std::move(channel_list.begin(), channel_list.end(),
891 std::inserter(process_list, process_list.end()));
892 channel_list.clear();
904 if (deutron && antideutron && pim && pi0 && pip) {
905 const ParticleTypePtrList deutron_list = {deutron};
906 const ParticleTypePtrList antideutron_list = {antideutron};
907 const ParticleTypePtrList pion_list = {pim, pi0, pip};
909 (both_antinucleons ? antideutron_list : deutron_list), pion_list,
912 return pCM(sqrts, type_res_1.
mass(), type_res_2.
mass());
914 process_list.reserve(process_list.size() + channel_list.size());
915 std::move(channel_list.begin(), channel_list.end(),
916 std::inserter(process_list, process_list.end()));
917 channel_list.clear();
929 const auto pdg_nucleon = type_nucleon.
pdgcode().
code();
947 bool incl_KN_to_KDelta =
950 bool incl_Strangeness_exchange =
953 CollisionBranchList process_list;
958 switch (pdg_nucleon) {
960 if (incl_Strangeness_exchange) {
970 sqrt_s_, type_pi_m, type_Sigma_p);
973 sqrt_s_, type_pi_p, type_Sigma_m);
976 type_pi_z, type_Sigma_z);
979 type_pi_z, type_Lambda);
990 if (incl_Strangeness_exchange) {
998 type_pi_m, type_Sigma_z);
1001 type_pi_z, type_Sigma_m);
1004 type_pi_m, type_Lambda);
1009 if (incl_KN_to_KDelta) {
1017 type_nucleon, type_kaon,
1021 sqrt_s_, type_Kbar_z, type_Delta_pp_bar);
1025 type_nucleon, type_kaon,
1026 type_K_m, type_Delta_p_bar);
1028 sqrt_s_, type_K_m, type_Delta_p_bar);
1033 if (incl_KN_to_KDelta) {
1041 type_nucleon, type_kaon,
1045 sqrt_s_, type_Kbar_z, type_Delta_p_bar);
1049 type_nucleon, type_kaon,
1050 type_K_m, type_Delta_z_bar);
1052 sqrt_s_, type_K_m, type_Delta_z_bar);
1054 if (incl_KN_to_KN) {
1058 type_Kbar_z, type_p_bar);
1068 switch (pdg_nucleon) {
1070 if (incl_KN_to_KDelta) {
1078 type_nucleon, type_kaon,
1079 type_K_z, type_Delta_pp);
1081 sqrt_s_, type_K_z, type_Delta_pp);
1085 type_nucleon, type_kaon,
1086 type_K_p, type_Delta_p);
1088 sqrt_s_, type_K_p, type_Delta_p);
1093 if (incl_KN_to_KDelta) {
1101 type_nucleon, type_kaon,
1102 type_K_z, type_Delta_p);
1104 sqrt_s_, type_K_z, type_Delta_p);
1108 type_nucleon, type_kaon,
1109 type_K_p, type_Delta_z);
1111 sqrt_s_, type_K_p, type_Delta_z);
1113 if (incl_KN_to_KN) {
1122 if (incl_Strangeness_exchange) {
1132 sqrt_s_, type_pi_p, type_Sigma_p_bar);
1135 sqrt_s_, type_pi_m, type_Sigma_m_bar);
1138 type_pi_z, type_Sigma_z_bar);
1141 type_pi_z, type_Lambda_bar);
1143 if (incl_KN_to_KN) {
1147 sqrt_s_, type_K_z, type_n_bar);
1152 if (incl_Strangeness_exchange) {
1160 type_pi_p, type_Sigma_z_bar);
1163 type_pi_z, type_Sigma_m_bar);
1166 type_pi_p, type_Lambda_bar);
1179 switch (pdg_nucleon) {
1181 if (incl_KN_to_KDelta) {
1189 type_nucleon, type_kaon,
1190 type_K_z, type_Delta_p);
1192 sqrt_s_, type_K_z, type_Delta_p);
1196 type_nucleon, type_kaon,
1197 type_K_p, type_Delta_z);
1199 sqrt_s_, type_K_p, type_Delta_z);
1201 if (incl_KN_to_KN) {
1215 if (incl_KN_to_KDelta) {
1223 type_nucleon, type_kaon,
1224 type_K_z, type_Delta_z);
1226 sqrt_s_, type_K_z, type_Delta_z);
1230 type_nucleon, type_kaon,
1231 type_K_p, type_Delta_m);
1233 sqrt_s_, type_K_p, type_Delta_m);
1238 if (incl_Strangeness_exchange) {
1246 type_pi_m, type_Sigma_z_bar);
1249 type_pi_z, type_Sigma_p_bar);
1252 type_pi_m, type_Lambda_bar);
1257 if (incl_Strangeness_exchange) {
1267 sqrt_s_, type_pi_m, type_Sigma_m_bar);
1270 sqrt_s_, type_pi_p, type_Sigma_p_bar);
1273 type_pi_z, type_Sigma_z_bar);
1276 type_pi_z, type_Lambda_bar);
1278 if (incl_KN_to_KN) {
1282 sqrt_s_, type_K_p, type_p_bar);
1290 switch (pdg_nucleon) {
1292 if (incl_Strangeness_exchange) {
1300 type_pi_z, type_Sigma_p);
1303 type_pi_p, type_Sigma_z);
1306 type_pi_p, type_Lambda);
1311 if (incl_Strangeness_exchange) {
1321 sqrt_s_, type_pi_p, type_Sigma_m);
1324 sqrt_s_, type_pi_m, type_Sigma_p);
1327 type_pi_z, type_Sigma_z);
1330 type_pi_z, type_Lambda);
1332 if (incl_KN_to_KN) {
1341 if (incl_KN_to_KDelta) {
1343 const auto& type_Kbar_z = type_kaon;
1349 type_nucleon, type_kaon,
1353 sqrt_s_, type_Kbar_z, type_Delta_bar_m);
1357 type_nucleon, type_kaon,
1358 type_K_m, type_Delta_bar_z);
1360 sqrt_s_, type_K_m, type_Delta_bar_z);
1362 if (incl_KN_to_KN) {
1371 sqrt_s_, type_K_m, type_n_bar);
1376 if (incl_KN_to_KDelta) {
1384 type_nucleon, type_kaon,
1388 sqrt_s_, type_Kbar_z, type_Delta_z_bar);
1392 type_nucleon, type_kaon,
1393 type_K_m, type_Delta_m_bar);
1395 sqrt_s_, type_K_m, type_Delta_m_bar);
1403 return process_list;
1408 CollisionBranchList process_list;
1410 return process_list;
1417 const auto pdg_delta = type_delta.
pdgcode().
code();
1425 switch (
pack(pdg_delta, pdg_kaon)) {
1436 type_p, type_K_p, type_kaon, type_delta) *
1449 type_kaon, type_p_bar,
1452 type_p_bar, type_K_m, type_kaon, type_delta) *
1455 sqrt_s_, type_p_bar, type_K_m);
1470 type_n, type_K_p, type_kaon, type_delta) *
1481 type_p, type_K_z, type_kaon, type_delta) *
1496 type_kaon, type_n_bar,
1499 type_n_bar, type_K_m, type_kaon, type_delta) *
1502 sqrt_s_, type_n_bar, type_K_m);
1507 type_kaon, type_p_bar,
1510 type_p_bar, type_Kbar_z, type_kaon, type_delta) *
1513 sqrt_s_, type_p_bar, type_Kbar_z);
1526 type_n, type_K_z, type_kaon, type_delta) *
1539 type_kaon, type_n_bar,
1542 type_n_bar, type_Kbar_z, type_kaon, type_delta) *
1545 sqrt_s_, type_n_bar, type_Kbar_z);
1552 return process_list;
1556 CollisionBranchList process_list;
1558 return process_list;
1565 const auto pdg_hyperon = type_hyperon.
pdgcode().
code();
1570 switch (
pack(pdg_hyperon, pdg_pion)) {
1577 s, type_hyperon, type_pion, type_n, type_K_m) *
1593 sqrt_s_, type_p, type_Kbar_z);
1602 type_pion, type_n_bar,
1606 sqrt_s_, type_n_bar, type_K_p);
1615 type_pion, type_p_bar,
1619 sqrt_s_, type_p_bar, type_K_z);
1628 s, type_hyperon, type_pion, type_n, type_K_m) *
1644 sqrt_s_, type_p, type_Kbar_z);
1653 type_pion, type_n_bar,
1657 sqrt_s_, type_n_bar, type_K_p);
1666 type_pion, type_p_bar,
1670 sqrt_s_, type_p_bar, type_K_z);
1679 s, type_hyperon, type_pion, type_n, type_K_m) *
1695 sqrt_s_, type_p, type_Kbar_z);
1704 type_pion, type_n_bar,
1708 sqrt_s_, type_n_bar, type_K_p);
1717 type_pion, type_p_bar,
1721 sqrt_s_, type_p_bar, type_K_z);
1732 s, type_hyperon, type_pion, type_p, type_K_m) *
1743 sqrt_s_, type_n, type_Kbar_z);
1754 type_pion, type_p_bar,
1758 sqrt_s_, type_p_bar, type_K_p);
1762 type_pion, type_n_bar,
1766 sqrt_s_, type_n_bar, type_K_z);
1777 s, type_hyperon, type_pion, type_p, type_K_m) *
1788 sqrt_s_, type_n, type_Kbar_z);
1799 type_pion, type_p_bar,
1803 sqrt_s_, type_p_bar, type_K_p);
1807 type_pion, type_n_bar,
1811 sqrt_s_, type_n_bar, type_K_z);
1822 s, type_hyperon, type_pion, type_p, type_K_m) *
1833 sqrt_s_, type_n, type_Kbar_z);
1844 type_pion, type_p_bar,
1848 sqrt_s_, type_p_bar, type_K_p);
1852 type_pion, type_n_bar,
1856 sqrt_s_, type_n_bar, type_K_z);
1867 s, type_hyperon, type_pion, type_p, type_K_m) *
1878 sqrt_s_, type_n, type_Kbar_z);
1889 type_pion, type_p_bar,
1893 sqrt_s_, type_p_bar, type_K_p);
1897 type_pion, type_n_bar,
1901 sqrt_s_, type_n_bar, type_K_z);
1908 return process_list;
1913 const auto& log = logger<LogArea::ScatterAction>();
1914 CollisionBranchList process_list;
1923 ParticleTypePtrList nuc = (baryon_number > 0)
1930 nuc_a->charge() + nuc_b->charge()) {
1934 for (
const int twoI :
I_tot_range(*nuc_a, *nuc_b)) {
1936 type_a, type_b, *nuc_a, *nuc_b, twoI);
1943 const double matrix_element =
1945 if (matrix_element <= 0.) {
1949 const double spin_factor = (nuc_a->spin() + 1) * (nuc_b->spin() + 1);
1950 const int sym_fac_in =
1952 const int sym_fac_out =
1953 (nuc_a->iso_multiplet() == nuc_b->iso_multiplet()) ? 2 : 1;
1954 double p_cm_final =
pCM_from_s(s, nuc_a->mass(), nuc_b->mass());
1955 const double xsection = isospin_factor * spin_factor * sym_fac_in /
1956 sym_fac_out * p_cm_final * matrix_element /
1960 process_list.push_back(make_unique<CollisionBranch>(
1962 log.debug(type_a.
name(), type_b.
name(),
"->", nuc_a->name(),
1963 nuc_b->name(),
" at sqrts [GeV] = ", sqrts,
1964 " with cs[mb] = ", xsection);
1982 if (produced_nucleus == &type_nucleus ||
1984 produced_nucleus->baryon_number() != type_nucleus.
baryon_number()) {
1991 const double matrix_element =
1992 295.5 + 2.862 / (0.00283735 +
pow_int(sqrts - 2.181, 2)) +
1993 0.0672 /
pow_int(tmp, 2) - 6.61753 / tmp;
1994 const double spin_factor =
1995 (produced_nucleus->spin() + 1) * (type_pi.
spin() + 1);
2000 double xsection = matrix_element * spin_factor / (s *
cm_momentum());
2001 if (produced_nucleus->is_stable()) {
2003 xsection *=
pCM_from_s(s, type_pi.
mass(), produced_nucleus->mass());
2006 const double resonance_integral =
2007 produced_nucleus->iso_multiplet()->get_integral_piR(sqrts);
2008 xsection *= resonance_integral;
2009 log.debug(
"Resonance integral ", resonance_integral,
2010 ", matrix element: ", matrix_element,
2013 process_list.push_back(make_unique<CollisionBranch>(
2015 log.debug(type_pi.
name(), type_nucleus.
name(),
"→ ", type_pi.
name(),
2016 produced_nucleus->name(),
" at ", sqrts,
2017 " GeV, xs[mb] = ", xsection);
2020 return process_list;
2029 CollisionBranchList process_list;
2036 if (produced_nucleus == &type_nucleus ||
2038 produced_nucleus->baryon_number() != type_nucleus.
baryon_number()) {
2041 double matrix_element = 0.0;
2049 matrix_element = 79.0474 / std::pow(tmp, 0.7897) + 654.596 * tmp;
2054 matrix_element = 342.572 / std::pow(tmp, 0.6);
2056 const double spin_factor =
2057 (produced_nucleus->spin() + 1) * (type_N.
spin() + 1);
2061 double xsection = matrix_element * spin_factor / (s *
cm_momentum());
2062 if (produced_nucleus->is_stable()) {
2064 xsection *=
pCM_from_s(s, type_N.
mass(), produced_nucleus->mass());
2067 const double resonance_integral =
2068 produced_nucleus->iso_multiplet()->get_integral_NR(sqrts);
2069 xsection *= resonance_integral;
2071 process_list.push_back(make_unique<CollisionBranch>(
2073 const auto& log = logger<LogArea::ScatterAction>();
2074 log.debug(type_N.
name(), type_nucleus.
name(),
"→ ", type_N.
name(),
2075 produced_nucleus->name(),
" at ", sqrts,
2076 " GeV, xs[mb] = ", xsection);
2078 return process_list;
2082 double total_string_xs,
StringProcess* string_process,
bool use_AQM)
const {
2083 const auto& log = logger<LogArea::CrossSections>();
2085 if (!string_process) {
2086 throw std::runtime_error(
"string_process should be initialized.");
2089 CollisionBranchList channel_list;
2090 if (total_string_xs <= 0.) {
2091 return channel_list;
2100 std::array<int, 2> pdgid;
2101 double AQM_factor = 1.;
2102 for (
int i = 0; i < 2; i++) {
2110 bool can_annihilate =
false;
2113 for (
int iq = 1; iq <= n_q_types; iq++) {
2114 std::array<int, 2> nquark;
2115 for (
int i = 0; i < 2; i++) {
2119 if (nquark[0] != 0 && nquark[1] != 0) {
2120 can_annihilate =
true;
2136 std::array<double, 3> xs =
2139 for (
int ip = 0; ip < 3; ip++) {
2140 xs[ip] *= AQM_factor;
2143 double single_diffr_AX = xs[0], single_diffr_XB = xs[1], double_diffr = xs[2];
2144 double single_diffr = single_diffr_AX + single_diffr_XB;
2145 double diffractive = single_diffr + double_diffr;
2151 double sig_annihilation = 0.0;
2152 if (can_annihilate) {
2158 xs_param *= AQM_factor;
2160 sig_annihilation = std::min(total_string_xs, xs_param);
2163 const double nondiffractive_all =
2164 std::max(0., total_string_xs - sig_annihilation - diffractive);
2165 diffractive = total_string_xs - sig_annihilation - nondiffractive_all;
2166 double_diffr = std::max(0., diffractive - single_diffr);
2167 const double a = (diffractive - double_diffr) / single_diffr;
2168 single_diffr_AX *= a;
2169 single_diffr_XB *= a;
2170 assert(std::abs(single_diffr_AX + single_diffr_XB + double_diffr +
2171 sig_annihilation + nondiffractive_all - total_string_xs) <
2174 double nondiffractive_soft = 0.;
2175 double nondiffractive_hard = 0.;
2176 if (nondiffractive_all > 0.) {
2181 nondiffractive_soft =
2182 nondiffractive_all * std::exp(-hard_xsec / nondiffractive_all);
2183 nondiffractive_hard = nondiffractive_all - nondiffractive_soft;
2185 log.debug(
"String cross sections [mb] are");
2186 log.debug(
"Single-diffractive AB->AX: ", single_diffr_AX);
2187 log.debug(
"Single-diffractive AB->XB: ", single_diffr_XB);
2188 log.debug(
"Double-diffractive AB->XX: ", double_diffr);
2189 log.debug(
"Soft non-diffractive: ", nondiffractive_soft);
2190 log.debug(
"Hard non-diffractive: ", nondiffractive_hard);
2191 log.debug(
"B-Bbar annihilation: ", sig_annihilation);
2194 const double sig_string_soft = total_string_xs - nondiffractive_hard;
2197 if (sig_string_soft > 0.) {
2198 channel_list.push_back(make_unique<CollisionBranch>(
2200 channel_list.push_back(make_unique<CollisionBranch>(
2202 channel_list.push_back(make_unique<CollisionBranch>(
2204 channel_list.push_back(make_unique<CollisionBranch>(
2206 if (can_annihilate) {
2207 channel_list.push_back(make_unique<CollisionBranch>(
2211 if (nondiffractive_hard > 0.) {
2212 channel_list.push_back(make_unique<CollisionBranch>(
2215 return channel_list;
2227 if (pdg_a == pdg_b) {
2247 xs = xs_l * (1. - prob_high) + xs_h * prob_high;
2262 }
else if ((pdg_a.is_meson() && pdg_b.
is_baryon()) ||
2263 (pdg_b.
is_meson() && pdg_a.is_baryon())) {
2269 if (pdg_a.is_meson() && pdg_b.
is_meson()) {
2276 xs *= (1. - 0.4 * pdg_a.frac_strange()) * (1. - 0.4 * pdg_b.
frac_strange());
2282 double cross_sec = 0.;
2306 const double current_xs)
const {
2307 const auto& log = logger<LogArea::CrossSections>();
2311 double nnbar_xsec = std::max(0.,
ppbar_total(s) - current_xs);
2312 log.debug(
"NNbar cross section is: ", nnbar_xsec);
2320 const auto& log = logger<LogArea::CrossSections>();
2321 CollisionBranchList channel_list;
2331 if (
sqrt_s_ - 2 * type_N.mass() < 0) {
2332 return channel_list;
2339 log.debug(
"NNbar reverse cross section is: ", xsection);
2340 channel_list.push_back(make_unique<CollisionBranch>(
2342 channel_list.push_back(make_unique<CollisionBranch>(
2345 return channel_list;
2349 const bool is_anti_particles)
const {
2352 CollisionBranchList process_list;
2358 ParticleTypePtrList nuc_or_anti_nuc;
2359 if (is_anti_particles) {
2370 nuc_a->charge() + nuc_b->charge()) {
2374 for (
const int twoI :
I_tot_range(*nuc_a, *nuc_b)) {
2376 type_a, type_b, *nuc_a, *nuc_b, twoI);
2383 const double matrix_element =
2385 if (matrix_element <= 0.) {
2393 const double spin_factor = (nuc_a->spin() + 1) * (nuc_b->spin() + 1);
2394 const int sym_fac_in =
2396 const int sym_fac_out =
2397 (nuc_a->iso_multiplet() == nuc_b->iso_multiplet()) ? 2 : 1;
2398 const double xsection = isospin_factor * spin_factor * sym_fac_in /
2399 sym_fac_out * p_cm_final * matrix_element /
2403 process_list.push_back(make_unique<CollisionBranch>(
2405 const auto& log = logger<LogArea::CrossSections>();
2406 log.debug(
"2->2 absorption with original particles: ", type_a,
2412 return process_list;
2419 const double m_a = type_a.
mass();
2420 const double m_b = type_b.
mass();
2421 const double msqr = 2. * (m_a * m_a + m_b * m_b);
2427 const double uplmt = m_a + m_b + 3.0 * (w_a + w_b) + 3.0;
2428 if (sqrts > uplmt) {
2435 return 68. / std::pow(sqrts - 1.104, 1.951);
2444 }
else if (twoI == 0) {
2445 const double parametrization = 14. / msqr;
2452 return 6.5 * parametrization;
2454 return parametrization;
2467 }
else if (twoI == 0) {
2481 }
else if (twoI == 0) {
2491 (1.0 - std::exp(-(sqrts - 2.0) * 20.0));
2498 template <
class IntegrationMethod>
2500 const ParticleTypePtrList& list_res_1,
2501 const ParticleTypePtrList& list_res_2,
2502 const IntegrationMethod integrator)
const {
2506 const auto& log = logger<LogArea::CrossSections>();
2507 CollisionBranchList channel_list;
2515 if (type_res_1->charge() + type_res_2->charge() !=
2521 for (
const int twoI :
I_tot_range(type_particle_a, type_particle_b)) {
2523 type_particle_a, type_particle_b, *type_res_1, *type_res_2, twoI);
2530 const double lower_limit = type_res_1->min_mass_kinematic();
2531 const double upper_limit =
sqrt_s_ - type_res_2->mass();
2535 if (upper_limit - lower_limit < 1E-3) {
2541 sqrt_s_, *type_res_1, *type_res_2, twoI);
2542 if (matrix_element <= 0.) {
2549 const double resonance_integral = integrator(*type_res_1, *type_res_2);
2554 const double spin_factor =
2555 (type_res_1->spin() + 1) * (type_res_2->spin() + 1);
2556 const double xsection = isospin_factor * spin_factor * matrix_element *
2560 channel_list.push_back(make_unique<CollisionBranch>(
2562 log.debug(
"Found 2->2 creation process for resonance ", type_res_1,
2564 log.debug(
"2->2 with original particles: ", type_particle_a,
2570 return channel_list;
2574 bool use_transition_probability,
2576 bool treat_BBbar_with_strings)
const {
2579 if (!strings_switch) {
2586 const bool is_NN_scattering =
2589 const bool is_BBbar_scattering =
2597 const bool is_AQM_scattering =
2603 const double mass_sum =
2606 if (!is_NN_scattering && !is_BBbar_scattering && !is_Npi_scattering &&
2607 !is_AQM_scattering) {
2609 }
else if (is_BBbar_scattering) {
2616 const bool is_KplusP =
2628 }
else if (pdg1.
is_pion() && pdg2.is_pion()) {
2633 if (!use_transition_probability) {
2634 return static_cast<double>(
sqrt_s_ > mass_sum + aqm_offset);
2638 double region_lower, region_upper;
2639 if (is_Npi_scattering) {
2642 }
else if (is_NN_scattering) {
2648 region_lower = mass_sum + aqm_offset;
2654 }
else if (
sqrt_s_ < region_lower) {
2664 const double region_lower,
const double region_upper)
const {
2673 double x = (
sqrt_s_ - 0.5 * (region_lower + region_upper)) /
2674 (region_upper - region_lower);
2675 assert(x >= -0.5 && x <= 0.5);
2676 double prob = 0.5 * (std::sin(M_PI * x) + 1.0);
2677 assert(prob >= 0. && prob <= 1.);
static PdgCode from_decimal(const int pdgcode_decimal)
Construct PDG code from decimal number.
unsigned int spin() const
CollisionBranchPtr NNbar_annihilation(const double current_xs) const
Determine the cross section for NNbar annihilation, which is given by the difference between the para...
double piminusp_elastic(double mandelstam_s)
pi-p elastic cross section parametrization Source: GiBUU:parametrizationBarMes_HighEnergy.f90
CollisionBranchList bb_xx_except_nn(ReactionsBitSet included_2to2) const
Find all inelastic 2->2 processes for Baryon-Baryon (BB) Scattering except the more specific Nucleon-...
double string_probability(bool strings_switch, bool use_transition_probability, bool use_AQM, bool treat_nnbar_with_strings) const
double cm_momentum() const
Determine the momenta of the incoming particles in the center-of-mass system.
CollisionBranchList deltak_xx(ReactionsBitSet included_2to2) const
Find all inelastic 2->2 processes for Delta-Kaon (DeltaK) Scattering.
double high_energy() const
Determine the parametrized total cross section at high energies for the given collision, which is non-zero for Baryon-Baryon and Nucleon-Pion scatterings currently.
const bool is_BBbar_pair_
Whether incoming particles are a baryon-antibaryon pair.
bool is_Deltastar() const
double kminusp_kbar0n(double mandelstam_s)
K- p <-> Kbar0 n cross section parametrization.
const double sqrt_s_
Total energy in the center-of-mass frame.
double kminusp_elastic_background(double mandelstam_s)
K- p elastic background cross section parametrization Source: Buss:2011mx, B.3.9. ...
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().
static ParticleTypePtrList & list_baryon_resonances()
int antiparticle_sign() const
IsoParticleType * iso_multiplet() const
static double detailed_balance_factor_RR(double sqrts, double pcm, const ParticleType &a, const ParticleType &b, const ParticleType &c, const ParticleType &d)
Helper function: Calculate the detailed balance factor R such that where and are unstable...
static ParticleTypePtrList & list_Deltas()
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.
double diffractive. Two strings are formed, one from A and one from B.
constexpr double really_small
Numerical error tolerance.
double deuteron_nucleon_elastic(double mandelstam_s)
Deuteron nucleon elastic cross-section [mb] parametrized by Oh:2009gx.
double piplusp_high_energy(double mandelstam_s)
pi+p total cross section at high energies
double piminusp_lambdak0_pdg(double mandelstam_s)
pi- p -> Lambda K0 cross section parametrization, PDG data.
double ppbar_elastic(double mandelstam_s)
ppbar elastic cross section parametrization Source: Bass:1998ca
double piminusp_sigma0k0_res(double mandelstam_s)
pi- p -> Sigma0 K0 cross section parametrization, resonance contribution.
String excitation processes used in SMASH.
double piminusp_sigmaminuskplus_pdg(double mandelstam_s)
pi- p -> Sigma- K+ cross section parametrization, PDG data.
CollisionBranchList dn_xx(ReactionsBitSet included_2to2) const
Find all inelastic 2->2 processes involving Nucleon and (anti-) Deuteron (dN), specifically Nd → Nd'...
CollisionBranchList two_to_two(ReactionsBitSet included_2to2) const
Find all inelastic 2->2 processes for the given scattering.
double kminusp_piminussigmaplus(double sqrts)
K- p <-> pi- Sigma+ cross section parametrization Taken from UrQMD (Graef:2014mra).
constexpr uint64_t pack(int32_t x, int32_t y)
Pack two int32_t into an uint64_t.
double kplusp_elastic_background(double mandelstam_s)
K+ p elastic background cross section parametrization.
static ParticleTypePtrList & list_light_nuclei()
CollisionBranchList rare_two_to_two() const
Find all 2->2 processes which are suppressed at high energies when strings are turned on with probabi...
CollisionBranchList two_to_one() const
Find all resonances that can be produced in a 2->1 collision of the two input particles and the produ...
Collection of useful constants that are known at compile time.
double width_at_pole() const
static ParticleTypePtrList & list_anti_nucleons()
double pp_high_energy(double mandelstam_s)
pp total cross section at high energies
const std::array< double, 2 > sqrts_range_NN
transition range in N-N collisions
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 ...
double kminusn_piminussigma0(double sqrts)
K- n <-> pi- Sigma0 cross section parametrization Follow from the parametrization with the same stran...
double k0n_elastic_background(double mandelstam_s)
K0 n elastic background cross section parametrization Source: Buss:2011mx, B.3.9. ...
a special case of baryon-antibaryon annihilation.
static int pdg_map_for_pythia(PdgCode &pdg)
Take pdg code and map onto particle specie which can be handled by PYTHIA.
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept
T pCM_from_s(const T s, const T mass_a, const T mass_b) noexcept
CollisionBranchList NNbar_creation() const
Determine the cross section for NNbar creation, which is given by detailed balance from the reverse r...
constexpr T pow_int(const T base, unsigned const exponent)
Efficient template for calculating integer powers using squaring.
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...
CollisionBranchList ypi_xx(ReactionsBitSet included_2to2) const
Find all inelastic 2->2 processes for Hyperon-Pion (Ypi) Scattering.
constexpr double fm2_mb
mb <-> fm^2 conversion factor.
2->2 inelastic scattering
double get_integral_RR(const ParticleType &type_res_2, double sqrts)
Look up the tabulated resonance integral for the XX -> RR cross section.
CollisionBranchList dpi_xx(ReactionsBitSet included_2to2) const
Find all inelastic 2->2 processes involving Pion and (anti-) Deuteron (dpi), specifically dπ→ NN...
double NN_string_hard(double mandelstam_s)
nucleon-nucleon hard scattering cross section (with partonic scattering)
double ppbar_high_energy(double mandelstam_s)
ppbar total cross section at high energies
constexpr double nucleon_mass
Nucleon mass in GeV.
double k0p_elastic_background(double mandelstam_s)
K0 p elastic background cross section parametrization Source: Buss:2011mx, B.3.9. ...
T pCM_sqr_from_s(const T s, const T mass_a, const T mass_b) noexcept
CollisionBranchPtr elastic(double elast_par, bool use_AQM) const
Determine the elastic cross section for this collision.
CollisionBranchList string_excitation(double total_string_xs, StringProcess *string_process, bool use_AQM) const
Determine the cross section for string excitations, which is given by the difference between the para...
CollisionBranchList nn_xx(ReactionsBitSet included_2to2) const
Find all inelastic 2->2 processes for Nucelon-Nucelon Scattering.
double spectral_function(double m) const
Full spectral function of the resonance (relativistic Breit-Wigner distribution with mass-dependent ...
constexpr double hbarc
GeV <-> fm conversion factor.
static double sum_xs_of(CollisionBranchList &list)
Helper function: Sum all cross sections of the given process list.
NNbarTreatment
Treatment of N Nbar Annihilation.
double get_integral_NR(double sqrts)
Look up the tabulated resonance integral for the XX -> NR cross section.
double string_hard_cross_section() const
Determine the (parametrized) hard non-diffractive string cross section for this collision.
double get_integral_RK(double sqrts)
Look up the tabulated resonance integral for the XX -> RK cross section.
double nk_el() const
Determine the elastic cross section for a nucleon-kaon (NK) collision.
double kplusn_k0p(double mandelstam_s)
K+ n charge exchange cross section parametrization.
KaonNucleonRatios kaon_nucleon_ratios
double probability_transit_high(const double region_lower, const double region_upper) const
constexpr double pion_mass
Pion mass in GeV.
Range of total isospin for reaction of particle a with particle b.
double kbar0n_elastic_background(double mandelstam_s)
Kbar0 n elastic background cross section parametrization Source: Buss:2011mx, B.3.9.
double nn_el() const
Determine the (parametrized) elastic cross section for a nucleon-nucleon (NN) collision.
CrossSections(const ParticleList &incoming_particles, const double sqrt_s, const std::pair< FourVector, FourVector > potentials)
Construct CrossSections instance.
static const ParticleType & find(PdgCode pdgcode)
Returns the ParticleType object for the given pdgcode.
constexpr int h1
h₁(1170).
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...
double kminusn_piminuslambda(double sqrts)
K- n <-> pi- Lambda cross section parametrization Follow from the parametrization with the same stran...
static const ParticleTypeList & list_all()
elastic scattering: particles remain the same, only momenta change
int charge() const
The charge of the particle.
double pipi_string_hard(double mandelstam_s)
pion-pion hard scattering cross section (with partonic scattering)
double kminusp_pi0sigma0(double sqrts)
K- p <-> pi0 Sigma0 cross section parametrization Fit to Landolt-Börnstein instead of UrQMD values...
int baryon_number() const
constexpr int decimal_d
Deuteron in decimal digits.
constexpr int Delta_pp
Δ⁺⁺.
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.
Particle type contains the static properties of a particle species.
constexpr double minimum_sqrts_pythia_can_handle
Energy in GeV, below which hard reactions via pythia are impossible.
constexpr int decimal_antid
Anti-deuteron in decimal digits.
double piminusp_high_energy(double mandelstam_s)
pi-p total cross section at high energies
double np_elastic(double mandelstam_s)
np elastic cross section parametrization Source: Weil:2013mya, eq.
double kplusn_elastic_background(double mandelstam_s)
K+ n elastic background cross section parametrization sigma(K+n->K+n) = sigma(K+n->K0p) = 0...
CollisionBranchList bar_bar_to_nuc_nuc(const bool is_anti_particles) const
Calculate cross sections for resonance absorption (i.e.
hard string process involving 2->2 QCD process by PYTHIA.
const double KN_offset
Constant offset as to where to shift from 2to2 to string processes (in GeV) in the case of KN reactio...
constexpr double deuteron_mass
Deuteron mass in GeV.
CollisionBranchList generate_collision_list(double elastic_parameter, bool two_to_one_switch, ReactionsBitSet included_2to2, double low_snn_cut, bool strings_switch, bool use_AQM, bool strings_with_probability, NNbarTreatment nnbar_treatment, StringProcess *string_process) const
Generate a list of all possible collisions between the incoming particles with the given c...
PdgCode stores a Particle Data Group Particle Numbering Scheme particle type number.
const double sqrts_range
constant for the range of transition region in the case of AQM this is added to the sum of masses + s...
double piplusp_elastic(double mandelstam_s)
pi+p elastic cross section parametrization, PDG data.
static double detailed_balance_factor_RK(double sqrts, double pcm, const ParticleType &a, const ParticleType &b, const ParticleType &c, const ParticleType &d)
Helper function: Calculate the detailed balance factor R such that where is unstable, is a kaon and are stable.
std::bitset< 6 > ReactionsBitSet
Container for the 2 to 2 reactions in the code.
ParticleTypePtr get_antiparticle() const
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 ParticleList incoming_particles_
List with data of scattering particles.
double pp_elastic(double mandelstam_s)
pp elastic cross section parametrization Source: Weil:2013mya, eq.
double kminusp_piplussigmaminus(double sqrts)
K- p <-> pi+ Sigma- cross section parametrization Taken from UrQMD (Graef:2014mra).
unsigned int spin() const
static double detailed_balance_factor_stable(double s, const ParticleType &a, const ParticleType &b, const ParticleType &c, const ParticleType &d)
Helper function: Calculate the detailed balance factor R such that where are stable.
double kminusn_elastic_background(double mandelstam_s)
K- n elastic background cross section parametrization Source: Buss:2011mx, B.3.9. ...
bool is_antiparticle_of(const PdgCode rhs) const
Use string fragmentation.
Use intermediate Resonances.
const ParticleType & type() const
Get the type of the particle.
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...
double kplusn_inelastic_background(double mandelstam_s)
K+ n inelastic background cross section parametrization Source: Buss:2011mx, B.3.8.
const double pipi_offset
Constant offset as to where to turn on the strings and elastic processes for pi pi reactions (this is...
PdgCode pdgcode() const
Get the pdgcode of the particle.
static const ParticleTypePtr try_find(PdgCode pdgcode)
Returns the ParticleTypePtr for the given pdgcode.
resonance formation (2->1)
double formation(const ParticleType &type_resonance, double cm_momentum_sqr) const
Return the 2-to-1 resonance production cross section for a given resonance.
double ppbar_total(double mandelstam_s)
ppbar total cross section parametrization Source: Bass:1998ca
single diffractive AB->XB.
double kminusp_pi0lambda(double sqrts)
K- p <-> pi0 Lambda cross section parametrization Fit to Landolt-Börnstein instead of UrQMD values...
double npi_el() const
Determine the elastic cross section for a nucleon-pion (Npi) collision.
double piplusp_sigmapluskplus_pdg(double mandelstam_s)
pi+ p to Sigma+ K+ cross section parametrization, PDG data.
A pointer-like interface to global references to ParticleType objects.
double npbar_high_energy(double mandelstam_s)
npbar total cross section at high energies
double piplusp_elastic_high_energy(double mandelstam_s, double m1, double m2)
pi+p elactic cross section parametrization.
double deuteron_pion_elastic(double mandelstam_s)
Deuteron pion elastic cross-section [mb] parametrized to fit pi-d elastic scattering data (the data c...
int antiparticle_sign() const
double frac_strange() const
static ParticleTypePtrList & list_nucleons()
double kbar0p_elastic_background(double mandelstam_s)
Kbar0 p elastic background cross section parametrization Source: Buss:2011mx, B.3.9.
double elastic_parametrization(bool use_AQM) const
Choose the appropriate parametrizations for given incoming particles and return the (parametrized) el...
non-diffractive. Two strings are formed both have ends in A and B.
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
CollisionBranchList npi_yk() const
Find all processes for Nucleon-Pion to Hyperon-Kaon Scattering.
const std::array< double, 2 > sqrts_range_Npi
transition range in N-pi collisions
bool is_Nstar1535() const
CollisionBranchList nk_xx(ReactionsBitSet included_2to2) const
Find all inelastic 2->2 background processes for Nucleon-Kaon (NK) Scattering.
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.
static double nn_to_resonance_matrix_element(double sqrts, const ParticleType &type_a, const ParticleType &type_b, const int twoI)
Scattering matrix amplitude squared (divided by 16π) for resonance production processes like NN → NR...
double np_high_energy(double mandelstam_s)
np total cross section at high energies
ParticleData contains the dynamic information of a certain particle.
(41-45) soft string excitations.
static ParticleTypePtrList & list_anti_Deltas()
std::int32_t code() const
double piplusp_elastic_AQM(double mandelstam_s, double m1, double m2)
An overload of piplusp_elastic_high_energy in which the very low part is replaced by a flat 5 mb cros...
double kplusp_inelastic_background(double mandelstam_s)
K+ p inelastic background cross section parametrization Source: Buss:2011mx, B.3.8.
double Npi_string_hard(double mandelstam_s)
nucleon-pion hard scattering cross section (with partonic scattering)
double xs_ppbar_annihilation(double mandelstam_s)
parametrized cross-section for proton-antiproton annihilation used in the UrQMD model ...
const std::string & name() const