Thermal-FIST 1.5
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
RandomGenerators.cpp
Go to the documentation of this file.
1/*
2 * Thermal-FIST package
3 *
4 * Copyright (c) 2015-2019 Volodymyr Vovchenko
5 *
6 * GNU General Public License (GPLv3 or later)
7 */
8#include "HRGEventGenerator/RandomGenerators.h"
9
10#include "HRGBase/xMath.h"
11#include "HRGEventGenerator/SimpleParticle.h"
12#include "HRGEventGenerator/ParticleDecaysMC.h"
13#include "HRGEventGenerator/EventGeneratorBase.h"
14
15namespace thermalfist {
16
17 namespace RandomGenerators {
18
19 MTRand randgenMT;
20
21 std::mt19937 rng_std(std::time(nullptr));
22
23 void SetSeed(const unsigned int seed) {
24 randgenMT.seed(seed);
25 }
26
27 int RandomPoisson(double mean) {
28 int n;
29 if (mean <= 0) return 0;
30 if (mean < 25) {
31 double expmean = exp(-mean);
32 double pir = 1;
33 n = -1;
34 while (1) {
35 n++;
36 pir *= randgenMT.rand();
37 if (pir <= expmean) break;
38 }
39 return n;
40 }
41 // for large value we use inversion method
42 else {//if (mean < 1E9) {
43 double em, t, y;
44 double sq, alxm, g;
45 double pi = xMath::Pi();
46
47 sq = sqrt(2.0*mean);
48 alxm = log(mean);
49 g = mean * alxm - xMath::LogGamma(mean + 1.0);
50
51 do {
52 do {
53 y = tan(pi*randgenMT.rand());
54 em = sq * y + mean;
55 } while (em < 0.0);
56
57 em = floor(em);
58 t = 0.9*(1.0 + y * y)* exp(em*alxm - xMath::LogGamma(em + 1.0) - g);
59 } while (randgenMT.rand() > t);
60
61 return static_cast<int> (em);
62
63 }
64 //else {
65 // // use Gaussian approximation vor very large values
66 // n = Int_t(Gaus(0,1)*TMath::Sqrt(mean) + mean + 0.5);
67 // return n;
68 //}
69 }
70
71
72 int RandomPoisson(double mean, MTRand &rangen) {
73 int n;
74 if (mean <= 0) return 0;
75 if (mean < 25) {
76 double expmean = exp(-mean);
77 double pir = 1;
78 n = -1;
79 while (1) {
80 n++;
81 pir *= rangen.rand();
82 if (pir <= expmean) break;
83 }
84 return n;
85 }
86 // for large value we use inversion method
87 else {//if (mean < 1E9) {
88 double em, t, y;
89 double sq, alxm, g;
90 double pi = xMath::Pi();
91
92 sq = sqrt(2.0*mean);
93 alxm = log(mean);
94 g = mean * alxm - xMath::LogGamma(mean + 1.0);
95
96 do {
97 do {
98 y = tan(pi*rangen.rand());
99 em = sq * y + mean;
100 } while (em < 0.0);
101
102 em = floor(em);
103 t = 0.9*(1.0 + y * y)* exp(em*alxm - xMath::LogGamma(em + 1.0) - g);
104 } while (rangen.rand() > t);
105
106 return static_cast<int> (em);
107
108 }
109 //else {
110 // // use Gaussian approximation vor very large values
111 // n = Int_t(Gaus(0,1)*TMath::Sqrt(mean) + mean +0.5);
112 // return n;
113 //}
114 }
115
116 double SkellamProbability(int k, double mu1, double mu2)
117 {
118 return exp(-(mu1 + mu2)) * pow(sqrt(mu1 / mu2), k) * xMath::BesselI(k, 2. * sqrt(mu1 * mu2));
119 }
120
121
122 double SiemensRasmussenMomentumGenerator::g(double x, double mass) const {
123 if (mass < 0.)
124 mass = m_Mass;
125 double tp = -log(x);
126 double talpha = alpha(tp);
127 //double en = w(tp);
128 double en = sqrt(tp * tp + mass * mass);
129 double sh = sinh(talpha);
130 double shtalpha = 1.;
131 if (talpha != 0.0)
132 shtalpha = sh / talpha;
133 double ch = sqrt(1. + sh * sh);
134 return tp * tp * exp(-m_Gamma * en / m_T) * ((1. + m_T / m_Gamma / en)*shtalpha - m_T / m_Gamma / en * ch) / x;
135 }
136
137 double SiemensRasmussenMomentumGenerator::g2(double x, double tp, double mass) const {
138 if (mass < 0.)
139 mass = m_Mass;
140 double talpha = alpha(tp);
141 //double en = w(tp);
142 double en = sqrt(tp * tp + mass * mass);
143 double sh = sinh(talpha);
144 double shtalpha = 1.;
145 if (talpha != 0.0)
146 shtalpha = sh / talpha;
147 double ch = sqrt(1. + sh * sh);
148 return tp * tp * exp(-m_Gamma * en / m_T) * ((1. + m_T / m_Gamma / en)*shtalpha - m_T / m_Gamma / en * ch) / x;
149 }
150
151 void SiemensRasmussenMomentumGenerator::FixParameters() {
152 double eps = 1e-8;
153 double l = 0., r = 1.;
154 double m1 = l + (r - l) / 3.;
155 double m2 = r - (r - l) / 3.;
156 int MAXITERS = 200;
157 int iter = 0;
158 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
159 if (g(m1) < g(m2)) {
160 l = m1;
161 }
162 else {
163 r = m2;
164 }
165 m1 = l + (r - l) / 3.;
166 m2 = r - (r - l) / 3.;
167 iter++;
168 }
169 m_Max = g((m1 + m2) / 2.);
170 }
171
172 double SiemensRasmussenMomentumGenerator::GetRandom(double mass) const {
173 if (mass < 0.)
174 mass = m_Mass;
175 while (1) {
176 double x0 = randgenMT.randDblExc();
177 double y0 = m_Max * randgenMT.randDblExc();
178 double mn = 1.;
179 if (mass != m_Mass)
180 mn = 10.;
181 if (y0 < mn * g(x0)) return -log(x0);
182 }
183 return 0.;
184 }
185
186 std::vector<double> SiemensRasmussenMomentumGenerator::GetMomentum(double mass) const {
187 std::vector<double> ret(0);
188 double tp = GetRandom(mass);
189 double tphi = 2. * xMath::Pi() * randgenMT.rand();
190 double cthe = 2. * randgenMT.rand() - 1.;
191 double sthe = sqrt(1. - cthe * cthe);
192 ret.push_back(tp*cos(tphi)*sthe); //px
193 ret.push_back(tp*sin(tphi)*sthe); //py
194 ret.push_back(tp*cthe); //pz
195 // TODO: proper Cartesian coordinates
196 ret.push_back(0.); //r0
197 ret.push_back(0.); //rx
198 ret.push_back(0.); //ry
199 ret.push_back(0.); //rz
200 return ret;
201 }
202
203
204 double ThermalMomentumGenerator::g(double x, double mass) const
205 {
206 if (mass < 0.)
207 mass = m_Mass;
208 double tp = -log(x);
209 double texp = exp((sqrt(mass * mass + tp * tp) - m_Mu) / m_T);
210 return tp * tp / (texp + m_Statistics) / x;
211 }
212
213 //void ThermalMomentumGenerator::FixParameters()
214 //{
215 // //double eps = 1e-8;
216 // //double l = 0., r = 1.;
217 // //double m1 = l + (r - l) / 3.;
218 // //double m2 = r - (r - l) / 3.;
219 // //int MAXITERS = 200;
220 // //int iter = 0;
221 // //while (fabs(m2 - m1) > eps && iter < MAXITERS) {
222 // // if (g(m1) < g(m2)) {
223 // // l = m1;
224 // // }
225 // // else {
226 // // r = m2;
227 // // }
228 // // m1 = l + (r - l) / 3.;
229 // // m2 = r - (r - l) / 3.;
230 // // iter++;
231 // //}
232 // //m_Max = g((m1 + m2) / 2.);
233 // m_Max = ComputeMaximum(m_Mass);
234 //}
235
236 double ThermalMomentumGenerator::ComputeMaximum(double mass) const
237 {
238 double eps = 1e-8;
239 double l = 0., r = 1.;
240 double m1 = l + (r - l) / 3.;
241 double m2 = r - (r - l) / 3.;
242 int MAXITERS = 200;
243 int iter = 0;
244 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
245 if (g(m1, mass) < g(m2, mass)) {
246 l = m1;
247 }
248 else {
249 r = m2;
250 }
251 m1 = l + (r - l) / 3.;
252 m2 = r - (r - l) / 3.;
253 iter++;
254 }
255 return g((m1 + m2) / 2., mass);
256 }
257
258 double ThermalMomentumGenerator::GetP(double mass) const
259 {
260 if (mass < 0.)
261 mass = m_Mass;
262 while (1) {
263 double x0 = randgenMT.randDblExc();
264
265 if (mass < m_Mu && m_Statistics == -1)
266 std::cerr << "**WARNING** ThermalMomentumGenerator::GetP: Bose-condensation mu " << m_Mu << " > mass " << mass << std::endl;
267
268 double M = m_Max;
269 if (mass != m_Mass)
270 M = ComputeMaximum(mass);
271
272 double prob = g(x0, mass) / M;
273
274 if (prob > 1.)
275 std::cerr << "**WARNING** ThermalMomentumGenerator::GetP: Probability exceeds unity by " << prob - 1. << std::endl;
276
277 if (randgenMT.randDblExc() < prob) return -log(x0);
278 }
279 return 0.;
280 }
281
282
283 BoostInvariantMomentumGenerator::BoostInvariantMomentumGenerator(BoostInvariantFreezeoutParametrization* FreezeoutModel,
284 double Tkin, double etamax, double mass, int statistics, double mu) :
285 m_FreezeoutModel(FreezeoutModel),
286 m_Tkin(Tkin), m_EtaMax(etamax), m_Mass(mass),
287 m_Generator(mass, statistics, Tkin, mu)
288 {
289 if (m_FreezeoutModel == NULL) {
290 //m_FreezeoutModel = new BoostInvariantFreezeoutParametrization();
291 }
292 }
293
294 BoostInvariantMomentumGenerator::~BoostInvariantMomentumGenerator()
295 {
296 if (m_FreezeoutModel != NULL)
297 delete m_FreezeoutModel;
298 }
299
300 std::vector<double> BoostInvariantMomentumGenerator::GetMomentum(double mass) const
301 {
302 if (mass < 0.)
303 mass = Mass();
304
305
306 double zetacand = GetRandomZeta(RandomGenerators::randgenMT);
307 double eta = -EtaMax() + 2. * EtaMax() * RandomGenerators::randgenMT.rand();
308 double ph = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand();
309
310 double betar = m_FreezeoutModel->tanhetaperp(zetacand);
311 double cosheta = cosh(eta);
312 double sinheta = sinh(eta);
313
314 double cosphi = cos(ph);
315 double sinphi = sin(ph);
316
317 double vx = betar * cosphi / cosheta;
318 double vy = betar * sinphi / cosheta;
319 double vz = tanh(eta);
320
321 double dRdZeta = m_FreezeoutModel->dRdZeta(zetacand);
322 double dtaudZeta = m_FreezeoutModel->dtaudZeta(zetacand);
323
324 double coshetaperp = m_FreezeoutModel->coshetaperp(zetacand);
325 double sinhetaperp = m_FreezeoutModel->sinhetaperp(zetacand);
326
327 std::vector<double> dsigma_lab;
328 dsigma_lab.push_back(dRdZeta * cosheta);
329 dsigma_lab.push_back(dtaudZeta * cosphi);
330 dsigma_lab.push_back(dtaudZeta * sinphi);
331 dsigma_lab.push_back(dRdZeta * sinheta);
332
333 // dsigma^\mu in the local rest frame
334 std::vector<double> dsigma_loc = LorentzBoost(dsigma_lab, vx, vy, vz);
335
336 // Maximum weight for the rejection sampling of the momentum
337 double maxWeight = 1. + std::abs(dsigma_loc[1] / dsigma_loc[0]) + std::abs(dsigma_loc[2] / dsigma_loc[0]) + std::abs(dsigma_loc[3] / dsigma_loc[0]);
338 maxWeight += 1.e-5; // To stabilize models where maxWeight is equal to unity, like Cracow model
339
340 SimpleParticle part(0., 0., 0., mass, 0);
341
342 while (true) {
343
344 double tp = m_Generator.GetP(mass);
345 double tphi = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand();
346 double cthe = 2. * RandomGenerators::randgenMT.rand() - 1.;
347 double sthe = sqrt(1. - cthe * cthe);
348 part.px = tp * cos(tphi) * sthe;
349 part.py = tp * sin(tphi) * sthe;
350 part.pz = tp * cthe;
351 part.p0 = sqrt(mass * mass + tp * tp);
352
353 double p0LRF = part.p0;
354
355 double dsigmamu_pmu_loc = dsigma_loc[0] * part.p0
356 - dsigma_loc[1] * part.px - dsigma_loc[2] * part.py - dsigma_loc[3] * part.pz;
357
358
359 double dsigmamu_umu_loc = dsigma_loc[0];
360
361 double dumu_pmu_loc = p0LRF;
362
363 double Weight = dsigmamu_pmu_loc / dsigmamu_umu_loc / dumu_pmu_loc / maxWeight;
364
365 if (Weight > 1.) {
366 std::cerr << "**WARNING** BoostInvariantHypersurfaceMomentumGenerator::GetMomentum: Weight exceeds unity by " << Weight - 1. << std::endl;
367 }
368
369 if (RandomGenerators::randgenMT.rand() < Weight)
370 break;
371
372 }
373
374 if (betar != 0.0 || eta != 0.0)
375 part = ParticleDecaysMC::LorentzBoostMomentumOnly(part, -vx, -vy, -vz);
376
377
378 std::vector<double> ret(3, 0.);
379 ret[0] = part.px;
380 ret[1] = part.py;
381 ret[2] = part.pz;
382
383 // Space-time coordinates
384 double tau = m_FreezeoutModel->taufunc(zetacand);
385 double r0 = tau * cosheta;
386 double rz = tau * sinheta;
387
388 double Rperp = m_FreezeoutModel->Rfunc(zetacand);
389 double rx = Rperp * cosphi;
390 double ry = Rperp * sinphi;
391
392 ret.push_back(r0);
393 ret.push_back(rx);
394 ret.push_back(ry);
395 ret.push_back(rz);
396 return ret;
397 }
398
399 double BoostInvariantMomentumGenerator::GetRandomZeta(MTRand& rangen) const
400 {
401 if (m_FreezeoutModel->InverseZetaDistributionIsExplicit())
402 return m_FreezeoutModel->InverseZetaDistribution(rangen.rand());
403
404 while (1) {
405 double zetacand = rangen.rand();
406
407 double prob = m_FreezeoutModel->ZetaProbability(zetacand) / m_FreezeoutModel->ProbabilityMaximum();
408
409 if (rangen.rand() < prob)
410 return zetacand;
411 }
412 return 0.;
413 }
414
415
416 double BreitWignerGenerator::f(double x) const {
417 return x / ((x*x - m_M * m_M)*(x*x - m_M * m_M) + m_M * m_M*m_Gamma*m_Gamma);
418 }
419
420 void BreitWignerGenerator::FixParameters() {
421 double eps = 1e-8;
422 double l = m_Mthr, r = m_M + 2.*m_Gamma;
423 double m1 = l + (r - l) / 3.;
424 double m2 = r - (r - l) / 3.;
425 int MAXITERS = 200;
426 int iter = 0;
427 while (fabs(m2 - m1) < eps && iter < MAXITERS) {
428 if (f(m1) < f(m2)) {
429 l = m1;
430 }
431 else {
432 r = m2;
433 }
434 m1 = l + (r - l) / 3.;
435 m2 = r - (r - l) / 3.;
436 iter++;
437 }
438 m_Max = f((m1 + m2) / 2.);
439 }
440
441 double BreitWignerGenerator::GetRandom() const {
442 //bool fl = true;
443 if (m_Gamma < 1e-7) return m_M;
444 while (1) {
445 double x0 = m_Mthr + (m_M + 2.*m_Gamma - m_Mthr) * randgenMT.rand();
446 double y0 = m_Max * randgenMT.rand();
447 if (y0 < f(x0)) return x0;
448 }
449 return 0.;
450 }
451
452 void BreitWignerGenerator::SetParameters(double M, double gamma, double mthr) {
453 m_M = M;
454 m_Gamma = gamma;
455 m_Mthr = mthr;
456 FixParameters();
457 }
458
459 void ThermalBreitWignerGenerator::SetParameters(ThermalParticle *part, double T, double Mu)
460 {
461 m_part = part;
462 m_T = T;
463 m_Mu = Mu;
464 FixParameters();
465 }
466
467 void ThermalBreitWignerGenerator::FixParameters()
468 {
469 double Threshold = m_part->DecayThresholdMass();
470 double Width = m_part->ResonanceWidth();
471 double Mass = m_part->Mass();
472
473 double a = std::max(Threshold, Mass - 2.*Width);
474 double b = Mass + 2.*Width;
475
476 m_Xmin = a;
477 m_Xmax = b + (b - a)*0.2;
478
479 m_Max = 0.;
480
481 int iters = 1000;
482 double dM = (m_Xmax - m_Xmin) / (iters - 1.);
483 for (int i = 0; i < iters; ++i) {
484 double tM = m_Xmin + dM * i;
485 m_Max = std::max(m_Max, f(tM));
486 }
487 m_Max *= 1.2;
488 }
489
490 double ThermalBreitWignerGenerator::f(double M) const
491 {
492 return m_part->ThermalMassDistribution(M, m_T, m_Mu, m_part->ResonanceWidth());
493 }
494
495 double ThermalBreitWignerGenerator::GetRandom() const
496 {
497 if (m_part->ResonanceWidth() / m_part->Mass() < 1.e-2)
498 return m_part->Mass();
499 while (true) {
500 double x0 = m_Xmin + (m_Xmax - m_Xmin) * randgenMT.rand();
501 double y0 = m_Max * randgenMT.rand();
502 if (y0 < f(x0)) return x0;
503 }
504 return 0.;
505 }
506
507 void ThermalEnergyBreitWignerGenerator::FixParameters()
508 {
509 double Threshold = m_part->DecayThresholdMass();
510 double Width = m_part->ResonanceWidth();
511 double Mass = m_part->Mass();
512
513 double a = std::max(Threshold, Mass - 2.*Width);
514 double b = Mass + 2.*Width;
515
516 if (m_part->Decays().size() == 0)
517 a = Mass - 2.*Width + 1.e-6;
518 else
519 a = m_part->DecayThresholdMassDynamical();
520
521 b = Mass + 2.*Width;
522
523 m_Xmin = a;
524 m_Xmax = b + 0.2 * (b - a);
525
526 //if (m_part->PdgId() == 32214)
527 // printf("%lf %lf\n", m_Xmin, m_Xmax);
528
529 m_Max = 0.;
530
531 int iters = 1000;
532 double dM = (m_Xmax - m_Xmin) / (iters - 1.);
533 for (int i = 0; i < iters; ++i) {
534 double tM = m_Xmin + dM * i;
535 m_Max = std::max(m_Max, f(tM));
536 }
537 m_Max *= 1.2;
538 }
539
540 double ThermalEnergyBreitWignerGenerator::f(double M) const
541 {
542 return m_part->ThermalMassDistribution(M, m_T, m_Mu, m_part->TotalWidtheBW(M));
543 }
544
545 //double BesselDistributionGenerator::pn(int n, double a, int nu)
546 //{
547 // if (nu < 0 || n < 0) return 0.0;
548 // double nfact = 1., nnufact = 1.;
549 // for (int i = 1; i <= n; ++i)
550 // nfact *= i;
551
552 // nnufact = nfact;
553 // for (int i = n + 1; i <= n + nu; ++i)
554 // nnufact *= i;
555
556 // return pow(a / 2., 2 * n + nu) / xMath::BesselI(nu, a) / nfact / nnufact;
557 //}
558
559 double BesselDistributionGenerator::pn(int n, double a, int nu)
560 {
561 if (nu < 0 || n < 0) return 0.0;
562
563 double logret = (2. * n + nu) * log(a/2.) - a;
564 for (int i = 1; i <= n; ++i)
565 logret -= 2. * log(i);
566 for (int i = n + 1; i <= n + nu; ++i)
567 logret -= log(i);
568
569 return exp(logret) / xMath::BesselIexp(nu, a);
570 }
571
572 double BesselDistributionGenerator::R(double x, int nu)
573 {
574 double hn2 = 1., hn1 = 0., hn = 0.;
575 double kn2 = 0., kn1 = 1., kn = 0.;
576 int nmax = 20;
577 for (int n = 1; n <= nmax; ++n) {
578 double an = 2. * (nu + n) / x;
579 hn = an * hn1 + hn2;
580 kn = an * kn1 + kn2;
581
582 hn2 = hn1;
583 hn1 = hn;
584
585 kn2 = kn1;
586 kn1 = kn;
587
588 if (n == nmax && nmax <= 1000 && abs(hn2 / kn2 - hn1 / kn1) > 1.e-9)
589 nmax *= 2;
590
591 if (n == 1000) {
592 std::cerr << "**WARNING** BesselDistributionGenerator::R(x,nu): Reached maximum iterations..." << std::endl;
593 }
594 }
595 return hn / kn;
596 }
597
598 double BesselDistributionGenerator::chi2(double a, int nu)
599 {
600 return mu(a, nu) + (1. / 4.) * a * a * R(a, nu) * (R(a,nu+1) - R(a,nu));
601 }
602
603 double BesselDistributionGenerator::sig2(double a, int nu)
604 {
605 double tmp = m(a, nu) - mu(a, nu);
606 return chi2(a, nu) + tmp * tmp;
607 }
608
609 double BesselDistributionGenerator::Q2(double a, int nu)
610 {
611 double A = sqrt(a*a + nu*nu);
612 double B = sqrt(a*a + (nu+1.)*(nu+1.));
613 double tmp = 1. + a * a * (1. + B - A) / 2. / (nu + A) / (nu + 1. + B);
614 return a*a / 2. / (nu + A) + tmp * tmp;
615 }
616
617 int BesselDistributionGenerator::RandomBesselPoisson(double a, int nu, MTRand & rangen)
618 {
619 if (nu < 0 || a <= 0.) return 0;
620
621 while (true) {
622 int n1 = RandomPoisson(a / 2., rangen);
623 int n2 = RandomPoisson(a / 2., rangen);
624 if (n1 - n2 == nu)
625 return n2;
626 }
627
628 return 0;
629 }
630
631 //double BesselDistributionGenerator::pmXmOverpm(int X, int tm, double a, int nu) {
632 // double tf1 = 1., tf2 = 1.;
633 // if (X > 0) {
634 // for (int i = 1; i <= X; ++i) {
635 // tf1 *= tm + i;
636 // tf2 *= tm + nu + i;
637 // }
638 // }
639 // else if (X < 0) {
640 // for (int i = 1; i <= -X; ++i) {
641 // tf1 *= tm + X + i;
642 // tf2 *= tm + nu + X + i;
643 // }
644 // tf1 = 1. / tf1;
645 // tf2 = 1. / tf2;
646 // }
647
648 // return pow(a / 2., 2 * X) / tf1 / tf2;
649 //}
650
651 double BesselDistributionGenerator::pmXmOverpm(int X, int tm, double a, int nu) {
652 double ret = 1.;
653 if (X > 0) {
654 for (int i = 1; i <= X; ++i) {
655 ret *= (a / 2.)/(tm + i);
656 ret *= (a / 2.)/(tm + nu + i);
657 }
658 }
659 else if (X < 0) {
660 for (int i = 1; i <= -X; ++i) {
661 ret *= (tm + X + i) / (a / 2.);
662 ret *= (tm + nu + X + i) / (a / 2.);
663 }
664 }
665
666 return ret;
667 }
668
669 int BesselDistributionGenerator::RandomBesselDevroye1(double a, int nu, MTRand & rangen)
670 {
671 int tm = m(a, nu);
672 double pm = pn(tm, a, nu);
673 double w = 1. + pm / 2.;
674 while (true) {
675 double U = rangen.rand();
676 double W = rangen.rand();
677 int S = 1;
678 if (rangen.rand() < 0.5)
679 S = -1;
680
681 double Y = 0.;
682 if (U <= w / (1. + w)) {
683 double V = rangen.rand();
684 Y = V * w / pm;
685 }
686 else {
687 double E = -log(rangen.randDblExc());
688 Y = (w + E) / pm;
689 }
690 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
691 int X = S * std::lround(Y);
692 #else
693 int X = S * ((Y > 0.0) ? std::floor(Y + 0.5) : std::ceil(Y - 0.5));
694 #endif
695
696 if (tm + X < 0)
697 continue;
698
699 double pratio = pmXmOverpm(X, tm, a, nu);
700
701 if (pratio != pratio) {
702 std::cerr << "**WARNING** BesselDistributionGenerator::RandomBesselDevroye1: Float problem!" << std::endl;
703 continue;
704 }
705
706 if (W * std::min(1., exp(w - pm * Y)) <= pratio)
707 return tm + X;
708 }
709 return 0;
710 }
711
712 int BesselDistributionGenerator::RandomBesselDevroye2(double a, int nu, MTRand & rangen)
713 {
714 int tm = m(a, nu);
715 double tsig = sqrt(sig2(a, nu));
716 double q = std::min(1. / tsig / sqrt(648.), 1. / 3.);
717 while (true) {
718 double U = rangen.rand();
719 double W = rangen.rand();
720 int S = 1;
721 if (rangen.rand() < 0.5)
722 S = -1;
723
724 double Y = 0.;
725 if (U <= (1. + 2. / q) / (1. + 4. / q)) {
726 double V = rangen.rand();
727 Y = V * (1. / 2. + 1. / q);
728 }
729 else {
730 double E = -log(rangen.randDblExc());
731 Y = 1. / 2. + 1. / q + E / q;
732 }
733 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
734 int X = S * std::lround(Y);
735 #else
736 int X = S * ((Y > 0.0) ? std::floor(Y + 0.5) : std::ceil(Y - 0.5));
737 #endif
738
739 if (tm + X < 0)
740 continue;
741
742 double pratio = pmXmOverpm(X, tm, a, nu);
743
744 if (pratio != pratio) {
745 std::cerr << "**WARNING** BesselDistributionGenerator::RandomBesselDevroye2: Float problem!" << std::endl;
746 continue;
747 }
748
749 if (W * std::min(1., exp(1. + q/2. - q*Y)) <= pratio)
750 return tm + X;
751 }
752 return 0;
753 }
754
755 int BesselDistributionGenerator::RandomBesselDevroye3(double a, int nu, MTRand & rangen)
756 {
757 int tm = m(a, nu);
758 double tQ = sqrt(Q2(a, nu));
759 double q = std::min(1. / tQ / sqrt(648.), 1. / 3.);
760 while (true) {
761 double U = rangen.rand();
762 double W = rangen.rand();
763 int S = 1;
764 if (rangen.rand() < 0.5)
765 S = -1;
766
767 double Y = 0.;
768 if (U <= (1. + 2. / q) / (1. + 4. / q)) {
769 double V = rangen.rand();
770 Y = V * (1. / 2. + 1. / q);
771 }
772 else {
773 double E = -log(rangen.randDblExc());
774 Y = 1. / 2. + 1. / q + E / q;
775 }
776 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
777 int X = S * std::lround(Y);
778 #else
779 int X = S * ((Y > 0.0) ? std::floor(Y + 0.5) : std::ceil(Y - 0.5));
780 #endif
781
782 if (tm + X < 0)
783 continue;
784
785 double pratio = pmXmOverpm(X, tm, a, nu);
786
787 if (pratio != pratio) {
788 std::cerr << "**WARNING** BesselDistributionGenerator::RandomBesselDevroye3: Float problem!" << std::endl;
789 continue;
790 }
791
792
793 if (W * std::min(1., exp(1. + q/2. - q*Y)) <= pratio)
794 return tm + X;
795 }
796 return 0;
797 }
798
799 int BesselDistributionGenerator::RandomBesselNormal(double a, int nu, MTRand & rangen)
800 {
801 double ret = -1.;
802 while (ret <= -0.5) {
803 ret = rangen.randNorm(mu(a,nu), sqrt(chi2(a,nu)));
804 }
805 return static_cast<int>(ret + 0.5);
806 }
807
808 int BesselDistributionGenerator::RandomBesselCombined(double a, int nu, MTRand & rangen)
809 {
810 int tm = m(a, nu);
811 if (tm < 6) {
812 if (nu <= a)
813 return RandomBesselPoisson(a, nu, rangen);
814 else
815 return RandomBesselDevroye3(a, nu, rangen);
816 }
817 else
818 return RandomBesselNormal(a, nu, rangen);
819 }
820
821 std::vector<double> SiemensRasmussenMomentumGeneratorGeneralized::GetMomentum(double mass) const
822 {
823 if (mass < 0.)
824 mass = GetMass();
825
826 double ph = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand();
827 double costh = 2. * RandomGenerators::randgenMT.rand() - 1.;
828 double sinth = sqrt(1. - costh * costh);
829
830 double vx = GetBeta() * sinth * cos(ph);
831 double vy = GetBeta() * sinth * sin(ph);
832 double vz = GetBeta() * costh;
833
834 SimpleParticle part(0., 0., 0., mass, 0);
835
836 double tp = m_Generator.GetP(mass);
837 double tphi = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand();
838 double cthe = 2. * RandomGenerators::randgenMT.rand() - 1.;
839 double sthe = sqrt(1. - cthe * cthe);
840 part.px = tp * cos(tphi) * sthe;
841 part.py = tp * sin(tphi) * sthe;
842 part.pz = tp * cthe;
843 part.p0 = sqrt(mass * mass + tp * tp);
844
845 if (GetBeta() != 0.0)
846 part = ParticleDecaysMC::LorentzBoostMomentumOnly(part, -vx, -vy, -vz);
847
848 std::vector<double> ret(7);
849 ret[0] = part.px;
850 ret[1] = part.py;
851 ret[2] = part.pz;
852
853 // Assume a sphere at t = 0
854 ret[3] = 0.;
855 ret[4] = GetR() * sinth * cos(ph);
856 ret[5] = GetR() * sinth * sin(ph);
857 ret[6] = GetR() * costh;
858
859 return ret;
860 }
861
862}
863
864} // namespace thermalfist
865
866
868// Deprecated code below
869
870namespace thermalfist {
871
872 namespace RandomGenerators {
873
874 void SSHMomentumGenerator::FixParameters() {
875 {
876 double eps = 1e-8;
877 double l = 0., r = 1.;
878 double m1 = l + (r - l) / 3.;
879 double m2 = r - (r - l) / 3.;
880 int MAXITERS = 200;
881 int iter = 0;
882 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
883 if (g(m1) < g(m2)) {
884 l = m1;
885 }
886 else {
887 r = m2;
888 }
889 m1 = l + (r - l) / 3.;
890 m2 = r - (r - l) / 3.;
891 iter++;
892 }
893 m_MaxPt = g((m1 + m2) / 2.);
894
895 m_dndpt.clear();
896 double dx = m_dPt;
897 for (double x = 0.5 * dx; x <= 1.; x += dx) {
898 m_dndpt.add_val(x, g(x));
899 }
900 }
901
902 {
903 m_dndy.resize(0);
904 m_MaxYs.resize(0);
905 double dx = m_dPt;
906 for (double x = 0.5 * dx; x <= 1.; x += dx) {
907
908 double pt = -log(x);
909
910 double eps = 1e-8;
911 double l = -4. - m_EtaMax, r = 4. + m_EtaMax;
912 double m1 = l + (r - l) / 3.;
913 double m2 = r - (r - l) / 3.;
914 int MAXITERS = 200;
915 int iter = 0;
916 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
917 if (m_distr.dndy(m1, pt) < m_distr.dndy(m2, pt)) {
918 l = m1;
919 }
920 else {
921 r = m2;
922 }
923 m1 = l + (r - l) / 3.;
924 m2 = r - (r - l) / 3.;
925 iter++;
926 }
927 m_MaxYs.push_back(m_distr.dndy((m1 + m2) / 2., pt));
928
929 m_dndy.push_back(SplineFunction());
930 double dy = m_dy;
931 for (double ty = -4. - m_EtaMax + 0.5 * dy; ty <= 4. + m_EtaMax; ty += dy) {
932 m_dndy[m_dndy.size() - 1].add_val(ty, m_distr.dndy(ty, pt));
933 }
934 }
935 }
936 }
937
938 void SSHMomentumGenerator::FixParameters2() {
939 {
940 double eps = 1e-8;
941 double l = 0., r = 1.;
942 double m1 = l + (r - l) / 3.;
943 double m2 = r - (r - l) / 3.;
944 int MAXITERS = 200;
945 int iter = 0;
946 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
947 if (g(m1) < g(m2)) {
948 l = m1;
949 }
950 else {
951 r = m2;
952 }
953 m1 = l + (r - l) / 3.;
954 m2 = r - (r - l) / 3.;
955 iter++;
956 }
957 m_MaxPt = g((m1 + m2) / 2.);
958
959 m_dndpt.clear();
960 double dx = m_dPt;
961 for (double x = 0.5 * dx; x <= 1.; x += dx) {
962 m_dndpt.add_val(x, g(x));
963 }
964 }
965
966 {
967 m_dndy.resize(0);
968 m_MaxYs.resize(0);
969 double dx = m_dPt;
970 for (double x = 0.5 * dx; x <= 1.; x += dx) {
971
972 double pt = -log(x);
973
974 double eps = 1e-8;
975 double l = -4. - m_EtaMax, r = 4. + m_EtaMax;
976 double m1 = l + (r - l) / 3.;
977 double m2 = r - (r - l) / 3.;
978 int MAXITERS = 200;
979 int iter = 0;
980 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
981 if (m_distr.dndysingle(m1, pt) < m_distr.dndysingle(m2, pt)) {
982 l = m1;
983 }
984 else {
985 r = m2;
986 }
987 m1 = l + (r - l) / 3.;
988 m2 = r - (r - l) / 3.;
989 iter++;
990 }
991 m_MaxYs.push_back(m_distr.dndysingle((m1 + m2) / 2., pt));
992
993 m_dndy.push_back(SplineFunction());
994 double dy = m_dy;
995 for (double ty = -4. + 0.5 * dy; ty <= 4.; ty += dy) {
996 m_dndy[m_dndy.size() - 1].add_val(ty, m_distr.dndysingle(ty, pt));
997 }
998 }
999 }
1000
1001 }
1002
1003 void SSHMomentumGenerator::FindMaximumPt() {
1004
1005 double eps = 1e-8;
1006 double l = 0., r = 1.;
1007 double m1 = l + (r - l) / 3.;
1008 double m2 = r - (r - l) / 3.;
1009 int MAXITERS = 200;
1010 int iter = 0;
1011 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
1012 if (g(m1) < g(m2)) {
1013 l = m1;
1014 }
1015 else {
1016 r = m2;
1017 }
1018 m1 = l + (r - l) / 3.;
1019 m2 = r - (r - l) / 3.;
1020 iter++;
1021 }
1022 m_MaxPt = g((m1 + m2) / 2.);
1023
1024 m_dndpt.clearall();
1025 double dx = 0.05;
1026 for (double x = 0.5 * dx; x <= 1.; x += dx) {
1027 m_dndpt.add_val(x, g(x));
1028 }
1029 }
1030
1031 void SSHMomentumGenerator::FindMaximumY(double pt) {
1032 double eps = 1e-8;
1033 double l = -4. - m_EtaMax, r = 4. + m_EtaMax;
1034 double m1 = l + (r - l) / 3.;
1035 double m2 = r - (r - l) / 3.;
1036 int MAXITERS = 200;
1037 int iter = 0;
1038 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
1039 if (m_distr.dndy(m1, pt) < m_distr.dndy(m2, pt)) {
1040 l = m1;
1041 }
1042 else {
1043 r = m2;
1044 }
1045 m1 = l + (r - l) / 3.;
1046 m2 = r - (r - l) / 3.;
1047 iter++;
1048 }
1049 m_MaxY = m_distr.dndy((m1 + m2) / 2., pt);
1050 }
1051
1052 void SSHMomentumGenerator::FindMaximumY2(double pt) {
1053 double eps = 1e-8;
1054 double l = -4. - m_EtaMax, r = 4. + m_EtaMax;
1055 double m1 = l + (r - l) / 3.;
1056 double m2 = r - (r - l) / 3.;
1057 int MAXITERS = 200;
1058 int iter = 0;
1059 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
1060 if (m_distr.dndysingle(m1, pt) < m_distr.dndysingle(m2, pt)) {
1061 l = m1;
1062 }
1063 else {
1064 r = m2;
1065 }
1066 m1 = l + (r - l) / 3.;
1067 m2 = r - (r - l) / 3.;
1068 iter++;
1069 }
1070 m_MaxY = m_distr.dndysingle((m1 + m2) / 2., pt);
1071 }
1072
1073 std::pair<double, double> SSHMomentumGenerator::GetRandom(double mass) {
1074 double tpt = 0., ty = 0.;
1075 while (1) {
1076 double x0 = randgenMT.randDblExc();
1077 double y0 = m_MaxPt * randgenMT.randDblExc();
1078 if (y0 < g2(x0)) {
1079 tpt = -log(x0);
1080 break;
1081 }
1082 }
1083 while (1) {
1084 int ind = (int)(exp(-tpt) / m_dPt);
1085 if (ind < 0) ind = 0;
1086 if (ind >= static_cast<int>(m_dndy.size())) ind = m_dndy.size() - 1;
1087 double x0 = -4. - m_EtaMax + (8. + 2. * m_EtaMax) * randgenMT.randDblExc();
1088 double y0 = m_MaxYs[ind] * randgenMT.randDblExc();
1089 if (y0 < m_dndy[ind].f(x0)) {
1090 ty = x0;
1091 break;
1092 }
1093 }
1094 return std::make_pair(tpt, ty);
1095 }
1096
1097 std::pair<double, double> SSHMomentumGenerator::GetRandom2(double mass) const {
1098 double tpt = 0., ty = 0., teta = 0.;
1099 while (1) {
1100 double x0 = randgenMT.randDblExc();
1101 double y0 = m_MaxPt * randgenMT.randDblExc();
1102 if (y0 < g2(x0)) {
1103 tpt = -log(x0);
1104 break;
1105 }
1106 }
1107 while (1) {
1108 int ind = (int)(exp(-tpt) / m_dPt);
1109 if (ind < 0) ind = 0;
1110 if (ind >= static_cast<int>(m_dndy.size())) ind = m_dndy.size() - 1;
1111 double x0 = -4. + (8.) * randgenMT.randDblExc();
1112 double y0 = m_MaxYs[ind] * randgenMT.randDblExc();
1113
1114 if (y0 < m_dndy[ind].f(x0)) {
1115 ty = x0;
1116 teta = -m_EtaMax + 2. * m_EtaMax * randgenMT.randDblExc();
1117 break;
1118 }
1119 }
1120 return std::make_pair(tpt, ty - teta);
1121 }
1122
1123 std::vector<double> SSHMomentumGenerator::GetMomentum(double mass) const {
1124 std::vector<double> ret(0);
1125 std::pair<double, double> pty = GetRandom2(mass);
1126 double tpt = pty.first;
1127 double ty = pty.second;
1128 double tphi = 2. * xMath::Pi() * randgenMT.rand();
1129 ret.push_back(tpt * cos(tphi)); //px
1130 ret.push_back(tpt * sin(tphi)); //py
1131 ret.push_back(sqrt(tpt * tpt + m_Mass * m_Mass) * sinh(ty)); //pz
1132 return ret;
1133 }
1134
1135 }
1136
1137} // namespace thermalfist
thermalfist::BoostInvariantFreezeoutParametrization
Base class implementing a longitudinally boost-invariant azimuthally symmetric freeze-out parametriza...
Definition FreezeoutModels.h:28
thermalfist::RandomGenerators::BoostInvariantMomentumGenerator::GetMomentum
virtual std::vector< double > GetMomentum(double mass=-1.) const
Definition RandomGenerators.cpp:300
thermalfist::RandomGenerators::BoostInvariantMomentumGenerator::BoostInvariantMomentumGenerator
BoostInvariantMomentumGenerator(BoostInvariantFreezeoutParametrization *FreezeoutModel=NULL, double Tkin=0.100, double etamax=3.0, double mass=0.938, int statistics=0, double mu=0)
Construct a new BoostInvariantMomentumGenerator object.
Definition RandomGenerators.cpp:283
thermalfist::RandomGenerators::BoostInvariantMomentumGenerator::GetRandomZeta
virtual double GetRandomZeta(MTRand &rangen=RandomGenerators::randgenMT) const
Samples zeta for use in Monte Carlo event generator.
Definition RandomGenerators.cpp:399
thermalfist::RandomGenerators::BoostInvariantMomentumGenerator::~BoostInvariantMomentumGenerator
virtual ~BoostInvariantMomentumGenerator()
BoostInvariantMomentumGenerator desctructor.
Definition RandomGenerators.cpp:294
thermalfist::RandomGenerators::BreitWignerGenerator::SetParameters
void SetParameters(double M, double gamma, double mthr)
Set the Breit-Wigner spectral function parameters.
Definition RandomGenerators.cpp:452
thermalfist::RandomGenerators::BreitWignerGenerator::GetRandom
double GetRandom() const
Samples the resonance mass from a relativistic Breit-Wigner distribution with a constant width.
Definition RandomGenerators.cpp:441
thermalfist::RandomGenerators::SSHMomentumGenerator::GetMomentum
std::vector< double > GetMomentum(double mass=-1.) const
Definition RandomGenerators.cpp:1123
thermalfist::RandomGenerators::SiemensRasmussenMomentumGeneratorGeneralized::GetMomentum
virtual std::vector< double > GetMomentum(double mass=-1.) const
Definition RandomGenerators.cpp:821
thermalfist::RandomGenerators::SiemensRasmussenMomentumGenerator::GetMomentum
virtual std::vector< double > GetMomentum(double mass=-1.) const
Definition RandomGenerators.cpp:186
thermalfist::RandomGenerators::ThermalBreitWignerGenerator::FixParameters
virtual void FixParameters()
Computes some auxiliary stuff needed for sampling.
Definition RandomGenerators.cpp:467
thermalfist::RandomGenerators::ThermalBreitWignerGenerator::SetParameters
void SetParameters(ThermalParticle *part, double T, double Mu)
Sets the parameters of the distribution.
Definition RandomGenerators.cpp:459
thermalfist::RandomGenerators::ThermalBreitWignerGenerator::f
virtual double f(double M) const
Unnormalized resonance mass probability density.
Definition RandomGenerators.cpp:490
thermalfist::RandomGenerators::ThermalEnergyBreitWignerGenerator::f
virtual double f(double M) const
Unnormalized resonance mass probability density.
Definition RandomGenerators.cpp:540
thermalfist::RandomGenerators::ThermalEnergyBreitWignerGenerator::FixParameters
virtual void FixParameters()
Computes some auxiliary stuff needed for sampling.
Definition RandomGenerators.cpp:507
thermalfist::RandomGenerators::ThermalMomentumGenerator::GetP
double GetP(double mass=-1.) const
Samples the momentum of a particle.
Definition RandomGenerators.cpp:258
thermalfist::ThermalParticle
Class containing all information about a particle specie.
Definition ThermalParticle.h:64
thermalfist::ParticleDecaysMC::LorentzBoostMomentumOnly
SimpleParticle LorentzBoostMomentumOnly(const SimpleParticle &part, double vx, double vy, double vz)
Lorentz boost of the 4-momentum of a particle.
Definition ParticleDecaysMC.cpp:64
thermalfist::RandomGenerators
Contains random generator functions used in the Monte Carlo Thermal Event Generator.
Definition HypersurfaceSampler.h:33
thermalfist::RandomGenerators::SkellamProbability
double SkellamProbability(int k, double mu1, double mu2)
Probability of a Skellam distributed random variable with Poisson means mu1 and mu2 to have the value...
Definition RandomGenerators.cpp:116
thermalfist::RandomGenerators::RandomPoisson
int RandomPoisson(double mean)
Generates random integer distributed by Poisson with specified mean Uses randgenMT.
Definition RandomGenerators.cpp:27
thermalfist::RandomGenerators::randgenMT
MTRand randgenMT
The Mersenne Twister random number generator.
Definition RandomGenerators.cpp:19
thermalfist::RandomGenerators::SetSeed
void SetSeed(const unsigned int seed)
Set the seed of the random number generator randgenMT.
Definition RandomGenerators.cpp:23
thermalfist::RandomGenerators::rng_std
std::mt19937 rng_std
The Mersenne Twister random number generator from STD library.
thermalfist::xMath::Pi
constexpr double Pi()
Pi constant.
Definition xMath.h:23
thermalfist::xMath::BesselIexp
double BesselIexp(int n, double x)
Modified Bessel function I_n(x), divided by exponential factor.
Definition xMath.cpp:791
thermalfist::xMath::LogGamma
double LogGamma(double x)
Computes the logarithm of the Gamma function.
Definition xMath.cpp:955
thermalfist::xMath::BesselI
double BesselI(int n, double x)
Integer order modified Bessel function I_n(x)
Definition xMath.cpp:192
thermalfist
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition CosmicEoS.h:9
thermalfist::LorentzBoost
std::vector< double > LorentzBoost(const std::vector< double > &fourvector, double vx, double vy, double vz)
Performs a Lorentz boost on a four-vector.
Definition EventGeneratorBase.cpp:32
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselNormal
static int RandomBesselNormal(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:799
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye1
static int RandomBesselDevroye1(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:669
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye3
static int RandomBesselDevroye3(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:755
thermalfist::RandomGenerators::BesselDistributionGenerator::pmXmOverpm
static double pmXmOverpm(int X, int tm, double a, int nu)
Definition RandomGenerators.cpp:651
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselCombined
static int RandomBesselCombined(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:808
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselPoisson
static int RandomBesselPoisson(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:617
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye2
static int RandomBesselDevroye2(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:712
thermalfist::SimpleParticle
Structure holding information about a single particle in the event generator.
Definition SimpleParticle.h:15
thermalfist::SimpleParticle::p0
double p0
Energy (in GeV)
Definition SimpleParticle.h:18
thermalfist::SimpleParticle::pz
double pz
3-momentum components (in GeV)
Definition SimpleParticle.h:16
xMath.h
Contains some extra mathematical functions used in the code.
Generated by doxygen 1.13.2