11 #include <gsl/gsl_sf_bessel.h>
24 double breit_wigner(
const double m,
const double pole,
const double width) {
25 const double msqr = m * m;
26 const double dmsqr = msqr - pole * pole;
27 return 2. * msqr * width / (M_PI * (dmsqr * dmsqr + msqr * width * width));
32 return cauchy(m, pole, width / 2.);
36 double cauchy(
double x,
double pole,
double width) {
37 const double dm = x - pole;
38 return width / (M_PI * (dm * dm + width * width));
43 const double temperature) {
44 return 4.0 * M_PI * momentum_sqr * std::exp(-energy / temperature);
48 const double temperature) {
49 return momentum_sqr * std::sqrt(momentum_sqr) *
50 std::exp(-energy / temperature);
54 const double temperature) {
55 return ((3.0 / 20.0) * (momentum_sqr / (temperature * temperature)) -
56 (6.0 / 5.0) * (energy / temperature) + (14.0 / 5.0)) *
57 std::exp(-energy / temperature) * momentum_sqr;
61 const double temperature) {
62 return (0.75 + 0.125 * (momentum_sqr / (temperature * temperature)) -
63 (1.0 / 30.0) * (momentum_sqr * energy /
64 (temperature * temperature * temperature)) +
66 (momentum_sqr * momentum_sqr /
67 (temperature * temperature * temperature * temperature))) *
68 std::exp(-energy / temperature) * momentum_sqr;
81 double effective_chemical_potential,
83 const double Boltzmann_factor =
84 std::exp(-(std::sqrt(momentum_radial * momentum_radial + mass * mass) -
85 effective_chemical_potential) /
88 return Boltzmann_factor / (1 + statistics * Boltzmann_factor);
97 const double mom_min = 0.0;
98 const double mom_max =
99 std::sqrt(50. * 50. * temperature * temperature - mass * mass);
102 const double p_non_eq_sq =
103 0.5 * (9 * temperature * temperature +
105 std::sqrt(81 * temperature * temperature + 36 * mass * mass));
106 const double e_non_eq = std::sqrt(p_non_eq_sq + mass * mass);
107 const double probability_max =
112 double energy, momentum_radial, probability;
117 energy = std::sqrt(momentum_radial * momentum_radial + mass * mass);
119 energy, momentum_radial * momentum_radial, temperature);
122 return momentum_radial;
129 double energy_average;
132 const double m_over_T = mass / temperature;
133 energy_average = 3 * temperature + mass * gsl_sf_bessel_K1(m_over_T) /
134 gsl_sf_bessel_Kn(2, m_over_T);
137 energy_average = 3 * temperature;
139 const double momentum_average_sqr =
140 (energy_average - mass) * (energy_average + mass);
141 const double energy_min = mass;
142 const double energy_max = 50. * temperature;
145 const double probability_max =
151 double momentum_radial_sqr, probability;
154 momentum_radial_sqr = (energy - mass) * (energy + mass);
159 return std::sqrt(momentum_radial_sqr);
166 double energy_average;
169 const double m_over_T = mass / temperature;
170 energy_average = 3 * temperature + mass * gsl_sf_bessel_K1(m_over_T) /
171 gsl_sf_bessel_Kn(2, m_over_T);
174 energy_average = 3 * temperature;
176 const double momentum_average_sqr =
177 (energy_average - mass) * (energy_average + mass);
178 const double energy_min = mass;
179 const double energy_max = 50. * temperature;
182 const double probability_max =
188 double momentum_radial_sqr, probability;
191 momentum_radial_sqr = (energy - mass) * (energy + mass);
196 return std::sqrt(momentum_radial_sqr);
203 double momentum_radial, energy;
205 if (temperature > 0.6 * mass) {
210 momentum_radial = temperature * (a + b + c);
211 energy = std::sqrt(momentum_radial * momentum_radial + mass * mass);
213 std::exp((momentum_radial - energy) / temperature)) {
220 const double I1 = mass * mass;
221 const double I2 = 2.0 * mass * temperature;
222 const double I3 = 2.0 * temperature * temperature;
223 const double Itot = I1 + I2 + I3;
225 if (r0 < I1 / Itot) {
227 K = -temperature * std::log(r1);
228 }
else if (r0 < (I1 + I2) / Itot) {
231 K = -temperature * std::log(r1 * r2);
236 K = -temperature * std::log(r1 * r2 * r3);
239 momentum_radial = std::sqrt((energy + mass) * (energy - mass));
245 return momentum_radial;
249 double momentum_radial;
254 momentum_radial = (3.0 / 4.0) * temperature * (a + b + c + d);
256 return momentum_radial;