Thermal-FIST 1.6
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 if (m_Gamma < 1e-7)
422 return;
423 double eps = 1e-8;
424 double l = m_Mthr, r = m_M + 2.*m_Gamma;
425 double m1 = l + (r - l) / 3.;
426 double m2 = r - (r - l) / 3.;
427 int MAXITERS = 200;
428 int iter = 0;
429 while (fabs(m2 - m1) < eps && iter < MAXITERS) {
430 if (f(m1) < f(m2)) {
431 l = m1;
432 }
433 else {
434 r = m2;
435 }
436 m1 = l + (r - l) / 3.;
437 m2 = r - (r - l) / 3.;
438 iter++;
439 }
440 m_Max = f((m1 + m2) / 2.);
441 }
442
443 double BreitWignerGenerator::GetRandom() const {
444 //bool fl = true;
445 if (m_Gamma < 1e-7) return m_M;
446 while (1) {
447 double x0 = m_Mthr + (m_M + 2.*m_Gamma - m_Mthr) * randgenMT.rand();
448 double y0 = m_Max * randgenMT.rand();
449 if (y0 < f(x0)) return x0;
450 }
451 return 0.;
452 }
453
454 void BreitWignerGenerator::SetParameters(double M, double gamma, double mthr) {
455 m_M = M;
456 m_Gamma = gamma;
457 m_Mthr = mthr;
458 FixParameters();
459 }
460
461 void ThermalBreitWignerGenerator::SetParameters(ThermalParticle *part, double T, double Mu)
462 {
463 m_part = part;
464 m_T = T;
465 m_Mu = Mu;
466 FixParameters();
467 }
468
469 void ThermalBreitWignerGenerator::FixParameters()
470 {
471 if (m_part->ZeroWidthEnforced()) {
472 m_Xmin = m_part->Mass();
473 m_Xmax = m_part->Mass();
474 m_Max = 1.;
475 return;
476 }
477 double Threshold = m_part->DecayThresholdMass();
478 double Width = m_part->ResonanceWidth();
479 double Mass = m_part->Mass();
480
481 double a = std::max(Threshold, Mass - 2.*Width);
482 double b = Mass + 2.*Width;
483
484 m_Xmin = a;
485 m_Xmax = b + (b - a)*0.2;
486
487 m_Max = 0.;
488
489 int iters = 1000;
490 double dM = (m_Xmax - m_Xmin) / (iters - 1.);
491 for (int i = 0; i < iters; ++i) {
492 double tM = m_Xmin + dM * i;
493 m_Max = std::max(m_Max, f(tM));
494 }
495 m_Max *= 1.2;
496 }
497
498 double ThermalBreitWignerGenerator::f(double M) const
499 {
500 return m_part->ThermalMassDistribution(M, m_T, m_Mu, m_part->ResonanceWidth());
501 }
502
503 double ThermalBreitWignerGenerator::GetRandom() const
504 {
505 if (m_part->ZeroWidthEnforced())
506 return m_part->Mass();
507 while (true) {
508 double x0 = m_Xmin + (m_Xmax - m_Xmin) * randgenMT.rand();
509 double y0 = m_Max * randgenMT.rand();
510 if (y0 < f(x0)) return x0;
511 }
512 return 0.;
513 }
514
515 void ThermalEnergyBreitWignerGenerator::FixParameters()
516 {
517 if (m_part->ZeroWidthEnforced()) {
518 m_Xmin = m_part->Mass();
519 m_Xmax = m_part->Mass();
520 m_Max = 1.;
521 return;
522 }
523 double Threshold = m_part->DecayThresholdMass();
524 double Width = m_part->ResonanceWidth();
525 double Mass = m_part->Mass();
526
527 double a = std::max(Threshold, Mass - 2.*Width);
528 double b = Mass + 2.*Width;
529
530 if (m_part->Decays().size() == 0)
531 a = Mass - 2.*Width + 1.e-6;
532 else
533 a = m_part->DecayThresholdMassDynamical();
534
535 b = Mass + 2.*Width;
536
537 m_Xmin = a;
538 m_Xmax = b + 0.2 * (b - a);
539
540 //if (m_part->PdgId() == 32214)
541 // printf("%lf %lf\n", m_Xmin, m_Xmax);
542
543 m_Max = 0.;
544
545 int iters = 1000;
546 double dM = (m_Xmax - m_Xmin) / (iters - 1.);
547 for (int i = 0; i < iters; ++i) {
548 double tM = m_Xmin + dM * i;
549 m_Max = std::max(m_Max, f(tM));
550 }
551 m_Max *= 1.2;
552 }
553
554 double ThermalEnergyBreitWignerGenerator::f(double M) const
555 {
556 return m_part->ThermalMassDistribution(M, m_T, m_Mu, m_part->TotalWidtheBW(M));
557 }
558
559 //double BesselDistributionGenerator::pn(int n, double a, int nu)
560 //{
561 // if (nu < 0 || n < 0) return 0.0;
562 // double nfact = 1., nnufact = 1.;
563 // for (int i = 1; i <= n; ++i)
564 // nfact *= i;
565
566 // nnufact = nfact;
567 // for (int i = n + 1; i <= n + nu; ++i)
568 // nnufact *= i;
569
570 // return pow(a / 2., 2 * n + nu) / xMath::BesselI(nu, a) / nfact / nnufact;
571 //}
572
573 double BesselDistributionGenerator::pn(int n, double a, int nu)
574 {
575 if (nu < 0 || n < 0) return 0.0;
576
577 double logret = (2. * n + nu) * log(a/2.) - a;
578 for (int i = 1; i <= n; ++i)
579 logret -= 2. * log(i);
580 for (int i = n + 1; i <= n + nu; ++i)
581 logret -= log(i);
582
583 return exp(logret) / xMath::BesselIexp(nu, a);
584 }
585
586 double BesselDistributionGenerator::R(double x, int nu)
587 {
588 double hn2 = 1., hn1 = 0., hn = 0.;
589 double kn2 = 0., kn1 = 1., kn = 0.;
590 int nmax = 20;
591 for (int n = 1; n <= nmax; ++n) {
592 double an = 2. * (nu + n) / x;
593 hn = an * hn1 + hn2;
594 kn = an * kn1 + kn2;
595
596 hn2 = hn1;
597 hn1 = hn;
598
599 kn2 = kn1;
600 kn1 = kn;
601
602 if (n == nmax && nmax <= 1000 && abs(hn2 / kn2 - hn1 / kn1) > 1.e-9)
603 nmax *= 2;
604
605 if (n == 1000) {
606 std::cerr << "**WARNING** BesselDistributionGenerator::R(x,nu): Reached maximum iterations..." << std::endl;
607 }
608 }
609 return hn / kn;
610 }
611
612 double BesselDistributionGenerator::chi2(double a, int nu)
613 {
614 return mu(a, nu) + (1. / 4.) * a * a * R(a, nu) * (R(a,nu+1) - R(a,nu));
615 }
616
617 double BesselDistributionGenerator::sig2(double a, int nu)
618 {
619 double tmp = m(a, nu) - mu(a, nu);
620 return chi2(a, nu) + tmp * tmp;
621 }
622
623 double BesselDistributionGenerator::Q2(double a, int nu)
624 {
625 double A = sqrt(a*a + nu*nu);
626 double B = sqrt(a*a + (nu+1.)*(nu+1.));
627 double tmp = 1. + a * a * (1. + B - A) / 2. / (nu + A) / (nu + 1. + B);
628 return a*a / 2. / (nu + A) + tmp * tmp;
629 }
630
631 int BesselDistributionGenerator::RandomBesselPoisson(double a, int nu, MTRand & rangen)
632 {
633 if (nu < 0 || a <= 0.) return 0;
634
635 while (true) {
636 int n1 = RandomPoisson(a / 2., rangen);
637 int n2 = RandomPoisson(a / 2., rangen);
638 if (n1 - n2 == nu)
639 return n2;
640 }
641
642 return 0;
643 }
644
645 //double BesselDistributionGenerator::pmXmOverpm(int X, int tm, double a, int nu) {
646 // double tf1 = 1., tf2 = 1.;
647 // if (X > 0) {
648 // for (int i = 1; i <= X; ++i) {
649 // tf1 *= tm + i;
650 // tf2 *= tm + nu + i;
651 // }
652 // }
653 // else if (X < 0) {
654 // for (int i = 1; i <= -X; ++i) {
655 // tf1 *= tm + X + i;
656 // tf2 *= tm + nu + X + i;
657 // }
658 // tf1 = 1. / tf1;
659 // tf2 = 1. / tf2;
660 // }
661
662 // return pow(a / 2., 2 * X) / tf1 / tf2;
663 //}
664
665 double BesselDistributionGenerator::pmXmOverpm(int X, int tm, double a, int nu) {
666 double ret = 1.;
667 if (X > 0) {
668 for (int i = 1; i <= X; ++i) {
669 ret *= (a / 2.)/(tm + i);
670 ret *= (a / 2.)/(tm + nu + i);
671 }
672 }
673 else if (X < 0) {
674 for (int i = 1; i <= -X; ++i) {
675 ret *= (tm + X + i) / (a / 2.);
676 ret *= (tm + nu + X + i) / (a / 2.);
677 }
678 }
679
680 return ret;
681 }
682
683 int BesselDistributionGenerator::RandomBesselDevroye1(double a, int nu, MTRand & rangen)
684 {
685 int tm = m(a, nu);
686 double pm = pn(tm, a, nu);
687 double w = 1. + pm / 2.;
688 while (true) {
689 double U = rangen.rand();
690 double W = rangen.rand();
691 int S = 1;
692 if (rangen.rand() < 0.5)
693 S = -1;
694
695 double Y = 0.;
696 if (U <= w / (1. + w)) {
697 double V = rangen.rand();
698 Y = V * w / pm;
699 }
700 else {
701 double E = -log(rangen.randDblExc());
702 Y = (w + E) / pm;
703 }
704 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
705 int X = S * std::lround(Y);
706 #else
707 int X = S * ((Y > 0.0) ? std::floor(Y + 0.5) : std::ceil(Y - 0.5));
708 #endif
709
710 if (tm + X < 0)
711 continue;
712
713 double pratio = pmXmOverpm(X, tm, a, nu);
714
715 if (pratio != pratio) {
716 std::cerr << "**WARNING** BesselDistributionGenerator::RandomBesselDevroye1: Float problem!" << std::endl;
717 continue;
718 }
719
720 if (W * std::min(1., exp(w - pm * Y)) <= pratio)
721 return tm + X;
722 }
723 return 0;
724 }
725
726 int BesselDistributionGenerator::RandomBesselDevroye2(double a, int nu, MTRand & rangen)
727 {
728 int tm = m(a, nu);
729 double tsig = sqrt(sig2(a, nu));
730 double q = std::min(1. / tsig / sqrt(648.), 1. / 3.);
731 while (true) {
732 double U = rangen.rand();
733 double W = rangen.rand();
734 int S = 1;
735 if (rangen.rand() < 0.5)
736 S = -1;
737
738 double Y = 0.;
739 if (U <= (1. + 2. / q) / (1. + 4. / q)) {
740 double V = rangen.rand();
741 Y = V * (1. / 2. + 1. / q);
742 }
743 else {
744 double E = -log(rangen.randDblExc());
745 Y = 1. / 2. + 1. / q + E / q;
746 }
747 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
748 int X = S * std::lround(Y);
749 #else
750 int X = S * ((Y > 0.0) ? std::floor(Y + 0.5) : std::ceil(Y - 0.5));
751 #endif
752
753 if (tm + X < 0)
754 continue;
755
756 double pratio = pmXmOverpm(X, tm, a, nu);
757
758 if (pratio != pratio) {
759 std::cerr << "**WARNING** BesselDistributionGenerator::RandomBesselDevroye2: Float problem!" << std::endl;
760 continue;
761 }
762
763 if (W * std::min(1., exp(1. + q/2. - q*Y)) <= pratio)
764 return tm + X;
765 }
766 return 0;
767 }
768
769 int BesselDistributionGenerator::RandomBesselDevroye3(double a, int nu, MTRand & rangen)
770 {
771 int tm = m(a, nu);
772 double tQ = sqrt(Q2(a, nu));
773 double q = std::min(1. / tQ / sqrt(648.), 1. / 3.);
774 while (true) {
775 double U = rangen.rand();
776 double W = rangen.rand();
777 int S = 1;
778 if (rangen.rand() < 0.5)
779 S = -1;
780
781 double Y = 0.;
782 if (U <= (1. + 2. / q) / (1. + 4. / q)) {
783 double V = rangen.rand();
784 Y = V * (1. / 2. + 1. / q);
785 }
786 else {
787 double E = -log(rangen.randDblExc());
788 Y = 1. / 2. + 1. / q + E / q;
789 }
790 #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900)
791 int X = S * std::lround(Y);
792 #else
793 int X = S * ((Y > 0.0) ? std::floor(Y + 0.5) : std::ceil(Y - 0.5));
794 #endif
795
796 if (tm + X < 0)
797 continue;
798
799 double pratio = pmXmOverpm(X, tm, a, nu);
800
801 if (pratio != pratio) {
802 std::cerr << "**WARNING** BesselDistributionGenerator::RandomBesselDevroye3: Float problem!" << std::endl;
803 continue;
804 }
805
806
807 if (W * std::min(1., exp(1. + q/2. - q*Y)) <= pratio)
808 return tm + X;
809 }
810 return 0;
811 }
812
813 int BesselDistributionGenerator::RandomBesselNormal(double a, int nu, MTRand & rangen)
814 {
815 double ret = -1.;
816 while (ret <= -0.5) {
817 ret = rangen.randNorm(mu(a,nu), sqrt(chi2(a,nu)));
818 }
819 return static_cast<int>(ret + 0.5);
820 }
821
822 int BesselDistributionGenerator::RandomBesselCombined(double a, int nu, MTRand & rangen)
823 {
824 int tm = m(a, nu);
825 if (tm < 6) {
826 if (nu <= a)
827 return RandomBesselPoisson(a, nu, rangen);
828 else
829 return RandomBesselDevroye3(a, nu, rangen);
830 }
831 else
832 return RandomBesselNormal(a, nu, rangen);
833 }
834
835 std::vector<double> SiemensRasmussenMomentumGeneratorGeneralized::GetMomentum(double mass) const
836 {
837 if (mass < 0.)
838 mass = GetMass();
839
840 double ph = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand();
841 double costh = 2. * RandomGenerators::randgenMT.rand() - 1.;
842 double sinth = sqrt(1. - costh * costh);
843
844 double vx = GetBeta() * sinth * cos(ph);
845 double vy = GetBeta() * sinth * sin(ph);
846 double vz = GetBeta() * costh;
847
848 SimpleParticle part(0., 0., 0., mass, 0);
849
850 double tp = m_Generator.GetP(mass);
851 double tphi = 2. * xMath::Pi() * RandomGenerators::randgenMT.rand();
852 double cthe = 2. * RandomGenerators::randgenMT.rand() - 1.;
853 double sthe = sqrt(1. - cthe * cthe);
854 part.px = tp * cos(tphi) * sthe;
855 part.py = tp * sin(tphi) * sthe;
856 part.pz = tp * cthe;
857 part.p0 = sqrt(mass * mass + tp * tp);
858
859 if (GetBeta() != 0.0)
860 part = ParticleDecaysMC::LorentzBoostMomentumOnly(part, -vx, -vy, -vz);
861
862 std::vector<double> ret(7);
863 ret[0] = part.px;
864 ret[1] = part.py;
865 ret[2] = part.pz;
866
867 // Assume a sphere at t = 0
868 ret[3] = 0.;
869 ret[4] = GetR() * sinth * cos(ph);
870 ret[5] = GetR() * sinth * sin(ph);
871 ret[6] = GetR() * costh;
872
873 return ret;
874 }
875
876}
877
878} // namespace thermalfist
879
880
882// Deprecated code below
883
884namespace thermalfist {
885
886 namespace RandomGenerators {
887
888 void SSHMomentumGenerator::FixParameters() {
889 {
890 double eps = 1e-8;
891 double l = 0., r = 1.;
892 double m1 = l + (r - l) / 3.;
893 double m2 = r - (r - l) / 3.;
894 int MAXITERS = 200;
895 int iter = 0;
896 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
897 if (g(m1) < g(m2)) {
898 l = m1;
899 }
900 else {
901 r = m2;
902 }
903 m1 = l + (r - l) / 3.;
904 m2 = r - (r - l) / 3.;
905 iter++;
906 }
907 m_MaxPt = g((m1 + m2) / 2.);
908
909 m_dndpt.clear();
910 double dx = m_dPt;
911 for (double x = 0.5 * dx; x <= 1.; x += dx) {
912 m_dndpt.add_val(x, g(x));
913 }
914 }
915
916 {
917 m_dndy.resize(0);
918 m_MaxYs.resize(0);
919 double dx = m_dPt;
920 for (double x = 0.5 * dx; x <= 1.; x += dx) {
921
922 double pt = -log(x);
923
924 double eps = 1e-8;
925 double l = -4. - m_EtaMax, r = 4. + m_EtaMax;
926 double m1 = l + (r - l) / 3.;
927 double m2 = r - (r - l) / 3.;
928 int MAXITERS = 200;
929 int iter = 0;
930 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
931 if (m_distr.dndy(m1, pt) < m_distr.dndy(m2, pt)) {
932 l = m1;
933 }
934 else {
935 r = m2;
936 }
937 m1 = l + (r - l) / 3.;
938 m2 = r - (r - l) / 3.;
939 iter++;
940 }
941 m_MaxYs.push_back(m_distr.dndy((m1 + m2) / 2., pt));
942
943 m_dndy.push_back(SplineFunction());
944 double dy = m_dy;
945 for (double ty = -4. - m_EtaMax + 0.5 * dy; ty <= 4. + m_EtaMax; ty += dy) {
946 m_dndy[m_dndy.size() - 1].add_val(ty, m_distr.dndy(ty, pt));
947 }
948 }
949 }
950 }
951
952 void SSHMomentumGenerator::FixParameters2() {
953 {
954 double eps = 1e-8;
955 double l = 0., r = 1.;
956 double m1 = l + (r - l) / 3.;
957 double m2 = r - (r - l) / 3.;
958 int MAXITERS = 200;
959 int iter = 0;
960 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
961 if (g(m1) < g(m2)) {
962 l = m1;
963 }
964 else {
965 r = m2;
966 }
967 m1 = l + (r - l) / 3.;
968 m2 = r - (r - l) / 3.;
969 iter++;
970 }
971 m_MaxPt = g((m1 + m2) / 2.);
972
973 m_dndpt.clear();
974 double dx = m_dPt;
975 for (double x = 0.5 * dx; x <= 1.; x += dx) {
976 m_dndpt.add_val(x, g(x));
977 }
978 }
979
980 {
981 m_dndy.resize(0);
982 m_MaxYs.resize(0);
983 double dx = m_dPt;
984 for (double x = 0.5 * dx; x <= 1.; x += dx) {
985
986 double pt = -log(x);
987
988 double eps = 1e-8;
989 double l = -4. - m_EtaMax, r = 4. + m_EtaMax;
990 double m1 = l + (r - l) / 3.;
991 double m2 = r - (r - l) / 3.;
992 int MAXITERS = 200;
993 int iter = 0;
994 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
995 if (m_distr.dndysingle(m1, pt) < m_distr.dndysingle(m2, pt)) {
996 l = m1;
997 }
998 else {
999 r = m2;
1000 }
1001 m1 = l + (r - l) / 3.;
1002 m2 = r - (r - l) / 3.;
1003 iter++;
1004 }
1005 m_MaxYs.push_back(m_distr.dndysingle((m1 + m2) / 2., pt));
1006
1007 m_dndy.push_back(SplineFunction());
1008 double dy = m_dy;
1009 for (double ty = -4. + 0.5 * dy; ty <= 4.; ty += dy) {
1010 m_dndy[m_dndy.size() - 1].add_val(ty, m_distr.dndysingle(ty, pt));
1011 }
1012 }
1013 }
1014
1015 }
1016
1017 void SSHMomentumGenerator::FindMaximumPt() {
1018
1019 double eps = 1e-8;
1020 double l = 0., r = 1.;
1021 double m1 = l + (r - l) / 3.;
1022 double m2 = r - (r - l) / 3.;
1023 int MAXITERS = 200;
1024 int iter = 0;
1025 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
1026 if (g(m1) < g(m2)) {
1027 l = m1;
1028 }
1029 else {
1030 r = m2;
1031 }
1032 m1 = l + (r - l) / 3.;
1033 m2 = r - (r - l) / 3.;
1034 iter++;
1035 }
1036 m_MaxPt = g((m1 + m2) / 2.);
1037
1038 m_dndpt.clearall();
1039 double dx = 0.05;
1040 for (double x = 0.5 * dx; x <= 1.; x += dx) {
1041 m_dndpt.add_val(x, g(x));
1042 }
1043 }
1044
1045 void SSHMomentumGenerator::FindMaximumY(double pt) {
1046 double eps = 1e-8;
1047 double l = -4. - m_EtaMax, r = 4. + m_EtaMax;
1048 double m1 = l + (r - l) / 3.;
1049 double m2 = r - (r - l) / 3.;
1050 int MAXITERS = 200;
1051 int iter = 0;
1052 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
1053 if (m_distr.dndy(m1, pt) < m_distr.dndy(m2, pt)) {
1054 l = m1;
1055 }
1056 else {
1057 r = m2;
1058 }
1059 m1 = l + (r - l) / 3.;
1060 m2 = r - (r - l) / 3.;
1061 iter++;
1062 }
1063 m_MaxY = m_distr.dndy((m1 + m2) / 2., pt);
1064 }
1065
1066 void SSHMomentumGenerator::FindMaximumY2(double pt) {
1067 double eps = 1e-8;
1068 double l = -4. - m_EtaMax, r = 4. + m_EtaMax;
1069 double m1 = l + (r - l) / 3.;
1070 double m2 = r - (r - l) / 3.;
1071 int MAXITERS = 200;
1072 int iter = 0;
1073 while (fabs(m2 - m1) > eps && iter < MAXITERS) {
1074 if (m_distr.dndysingle(m1, pt) < m_distr.dndysingle(m2, pt)) {
1075 l = m1;
1076 }
1077 else {
1078 r = m2;
1079 }
1080 m1 = l + (r - l) / 3.;
1081 m2 = r - (r - l) / 3.;
1082 iter++;
1083 }
1084 m_MaxY = m_distr.dndysingle((m1 + m2) / 2., pt);
1085 }
1086
1087 std::pair<double, double> SSHMomentumGenerator::GetRandom(double mass) {
1088 double tpt = 0., ty = 0.;
1089 while (1) {
1090 double x0 = randgenMT.randDblExc();
1091 double y0 = m_MaxPt * randgenMT.randDblExc();
1092 if (y0 < g2(x0)) {
1093 tpt = -log(x0);
1094 break;
1095 }
1096 }
1097 while (1) {
1098 int ind = (int)(exp(-tpt) / m_dPt);
1099 if (ind < 0) ind = 0;
1100 if (ind >= static_cast<int>(m_dndy.size())) ind = m_dndy.size() - 1;
1101 double x0 = -4. - m_EtaMax + (8. + 2. * m_EtaMax) * randgenMT.randDblExc();
1102 double y0 = m_MaxYs[ind] * randgenMT.randDblExc();
1103 if (y0 < m_dndy[ind].f(x0)) {
1104 ty = x0;
1105 break;
1106 }
1107 }
1108 return std::make_pair(tpt, ty);
1109 }
1110
1111 std::pair<double, double> SSHMomentumGenerator::GetRandom2(double mass) const {
1112 double tpt = 0., ty = 0., teta = 0.;
1113 while (1) {
1114 double x0 = randgenMT.randDblExc();
1115 double y0 = m_MaxPt * randgenMT.randDblExc();
1116 if (y0 < g2(x0)) {
1117 tpt = -log(x0);
1118 break;
1119 }
1120 }
1121 while (1) {
1122 int ind = (int)(exp(-tpt) / m_dPt);
1123 if (ind < 0) ind = 0;
1124 if (ind >= static_cast<int>(m_dndy.size())) ind = m_dndy.size() - 1;
1125 double x0 = -4. + (8.) * randgenMT.randDblExc();
1126 double y0 = m_MaxYs[ind] * randgenMT.randDblExc();
1127
1128 if (y0 < m_dndy[ind].f(x0)) {
1129 ty = x0;
1130 teta = -m_EtaMax + 2. * m_EtaMax * randgenMT.randDblExc();
1131 break;
1132 }
1133 }
1134 return std::make_pair(tpt, ty - teta);
1135 }
1136
1137 std::vector<double> SSHMomentumGenerator::GetMomentum(double mass) const {
1138 std::vector<double> ret(0);
1139 std::pair<double, double> pty = GetRandom2(mass);
1140 double tpt = pty.first;
1141 double ty = pty.second;
1142 double tphi = 2. * xMath::Pi() * randgenMT.rand();
1143 ret.push_back(tpt * cos(tphi)); //px
1144 ret.push_back(tpt * sin(tphi)); //py
1145 ret.push_back(sqrt(tpt * tpt + m_Mass * m_Mass) * sinh(ty)); //pz
1146 return ret;
1147 }
1148
1149 }
1150
1151} // 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:454
thermalfist::RandomGenerators::BreitWignerGenerator::GetRandom
double GetRandom() const
Samples the resonance mass from a relativistic Breit-Wigner distribution with a constant width.
Definition RandomGenerators.cpp:443
thermalfist::RandomGenerators::SSHMomentumGenerator::GetMomentum
std::vector< double > GetMomentum(double mass=-1.) const
Definition RandomGenerators.cpp:1137
thermalfist::RandomGenerators::SiemensRasmussenMomentumGeneratorGeneralized::GetMomentum
virtual std::vector< double > GetMomentum(double mass=-1.) const
Definition RandomGenerators.cpp:835
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:469
thermalfist::RandomGenerators::ThermalBreitWignerGenerator::SetParameters
void SetParameters(ThermalParticle *part, double T, double Mu)
Sets the parameters of the distribution.
Definition RandomGenerators.cpp:461
thermalfist::RandomGenerators::ThermalBreitWignerGenerator::f
virtual double f(double M) const
Unnormalized resonance mass probability density.
Definition RandomGenerators.cpp:498
thermalfist::RandomGenerators::ThermalEnergyBreitWignerGenerator::f
virtual double f(double M) const
Unnormalized resonance mass probability density.
Definition RandomGenerators.cpp:554
thermalfist::RandomGenerators::ThermalEnergyBreitWignerGenerator::FixParameters
virtual void FixParameters()
Computes some auxiliary stuff needed for sampling.
Definition RandomGenerators.cpp:515
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:813
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye1
static int RandomBesselDevroye1(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:683
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye3
static int RandomBesselDevroye3(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:769
thermalfist::RandomGenerators::BesselDistributionGenerator::pmXmOverpm
static double pmXmOverpm(int X, int tm, double a, int nu)
Definition RandomGenerators.cpp:665
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselCombined
static int RandomBesselCombined(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:822
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselPoisson
static int RandomBesselPoisson(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:631
thermalfist::RandomGenerators::BesselDistributionGenerator::RandomBesselDevroye2
static int RandomBesselDevroye2(double a, int nu, MTRand &rangen)
Definition RandomGenerators.cpp:726
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