Skip to content

Commit 03bbc53

Browse files
vlvovchvlvovch
authored
Merge pull request #51 from aohlson/master
Fixes for: - infinite loop in event generator at low multiplicity (no particles) - std::invalid_argument errors and other issues due to Breit-Wigner enforcement in event generator Enhancements: - Implemented m->0 limits for ideal gas functions utilizing cluster expansion method
2 parents 8620e34 + b159260 commit 03bbc53

File tree

5 files changed

+149
-38
lines changed

5 files changed

+149
-38
lines changed

‎include/HRGEventGenerator/RandomGenerators.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
#include "HRGEventGenerator/MomentumDistribution.h"
1515
#include "HRGBase/ThermalParticle.h"
1616
#include <random>
17+
#include <utility>
1718

1819
namespace thermalfist {
1920

‎src/library/HRGBase/IdealGasFunctions.cpp‎

Lines changed: 119 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -199,12 +199,21 @@ namespace thermalfist {
199199
double sign = 1.;
200200
double ret = 0.;
201201
for (int i = 1; i <= order; ++i) {
202-
ret += sign * xMath::BesselKexp(2, i*moverT) * cfug / static_cast<double>(i);
202+
if (m == 0.)
203+
ret += sign * cfug / static_cast<double>(i) / static_cast<double>(i) / static_cast<double>(i);
204+
else
205+
ret += sign * xMath::BesselKexp(2, i*moverT) * cfug / static_cast<double>(i);
203206
cfug *= tfug;
204207
if (signchange) sign = -sign;
205208
}
206-
ret *= deg * m * m * T / 2. / xMath::Pi() / xMath::Pi() * xMath::GeVtoifm3();
207-
return ret;
209+
210+
double prefactor = 1.;
211+
if (m == 0.)
212+
prefactor = deg * T * T * T / 2. / xMath::Pi() / xMath::Pi() * 2.;
213+
else
214+
prefactor = deg * m * m * T / 2. / xMath::Pi() / xMath::Pi();
215+
216+
return ret * prefactor * xMath::GeVtoifm3();
208217
}
209218

210219
double QuantumClusterExpansionPressure(int statistics, double T, double mu, double m, double deg, int order,
@@ -254,12 +263,21 @@ namespace thermalfist {
254263
double sign = 1.;
255264
double ret = 0.;
256265
for (int i = 1; i <= order; ++i) {
257-
ret += sign * xMath::BesselKexp(2, i*moverT) * cfug / static_cast<double>(i) / static_cast<double>(i);
266+
if (m == 0.)
267+
ret += sign * cfug / static_cast<double>(i) / static_cast<double>(i) / static_cast<double>(i) / static_cast<double>(i);
268+
else
269+
ret += sign * xMath::BesselKexp(2, i*moverT) * cfug / static_cast<double>(i) / static_cast<double>(i);
258270
cfug *= tfug;
259271
if (signchange) sign = -sign;
260272
}
261-
ret *= deg * m * m * T * T / 2. / xMath::Pi() / xMath::Pi() * xMath::GeVtoifm3();
262-
return ret;
273+
274+
double prefactor = 1.;
275+
if (m == 0.)
276+
prefactor = deg * T * T * T * T / 2. / xMath::Pi() / xMath::Pi() * 2.;
277+
else
278+
prefactor = deg * m * m * T * T / 2. / xMath::Pi() / xMath::Pi();
279+
280+
return ret * prefactor * xMath::GeVtoifm3();
263281
}
264282

265283
double QuantumClusterExpansionEnergyDensity(int statistics, double T, double mu, double m, double deg, int order,
@@ -311,12 +329,21 @@ namespace thermalfist {
311329
double sign = 1.;
312330
double ret = 0.;
313331
for (int i = 1; i <= order; ++i) {
314-
ret += sign * (xMath::BesselKexp(1, i*moverT) + 3. * xMath::BesselKexp(2, i*moverT) / moverT / static_cast<double>(i)) * cfug / static_cast<double>(i);
332+
if (m == 0.)
333+
ret += sign * cfug / static_cast<double>(i) / static_cast<double>(i) / static_cast<double>(i) / static_cast<double>(i);
334+
else
335+
ret += sign * (xMath::BesselKexp(1, i*moverT) + 3. * xMath::BesselKexp(2, i*moverT) / moverT / static_cast<double>(i)) * cfug / static_cast<double>(i);
315336
cfug *= tfug;
316337
if (signchange) sign = -sign;
317338
}
318-
ret *= deg * m * m * m * T / 2. / xMath::Pi() / xMath::Pi() * xMath::GeVtoifm3();
319-
return ret;
339+
340+
double prefactor = 1.;
341+
if (m == 0.)
342+
prefactor = deg * T * T * T * T / 2. / xMath::Pi() / xMath::Pi() * 6.;
343+
else
344+
prefactor = deg * m * m * m * T / 2. / xMath::Pi() / xMath::Pi();
345+
346+
return ret * prefactor * xMath::GeVtoifm3();
320347
}
321348

322349
double QuantumClusterExpansionEntropyDensity(int statistics, double T, double mu, double m, double deg, int order,
@@ -378,6 +405,10 @@ namespace thermalfist {
378405
double moverT = m / T;
379406
double sign = 1.;
380407
double ret = 0.;
408+
409+
if (m == 0.)
410+
return ret;
411+
381412
for (int i = 1; i <= order; ++i) {
382413
ret += sign * xMath::BesselKexp(1, i*moverT) * cfug / static_cast<double>(i);
383414
cfug *= tfug;
@@ -433,12 +464,21 @@ namespace thermalfist {
433464
double sign = 1.;
434465
double ret = 0.;
435466
for (int i = 1; i <= order; ++i) {
436-
ret += sign * xMath::BesselKexp(2, i*moverT) * cfug * pow(static_cast<double>(i), N - 1);
467+
if (m == 0.)
468+
ret += sign * cfug * pow(static_cast<double>(i), N - 3);
469+
else
470+
ret += sign * xMath::BesselKexp(2, i*moverT) * cfug * pow(static_cast<double>(i), N - 1);
437471
cfug *= tfug;
438472
if (signchange) sign = -sign;
439-
}
440-
ret *= deg * m * m * T / 2. / xMath::Pi() / xMath::Pi() * xMath::GeVtoifm3();
441-
return ret;
473+
}
474+
475+
double prefactor = 1.;
476+
if (m == 0.)
477+
prefactor = deg * T * T * T / 2. / xMath::Pi() / xMath::Pi() * 2.;
478+
else
479+
prefactor = deg * m * m * T / 2. / xMath::Pi() / xMath::Pi();
480+
481+
return ret * prefactor * xMath::GeVtoifm3();
442482
}
443483

444484
double QuantumClusterExpansionChiN(int N, int statistics, double T, double mu, double m, double deg, int order,
@@ -1607,6 +1647,8 @@ namespace thermalfist {
16071647
assert(extraConfig.MagneticField.B == 0);
16081648

16091649
// No magnetic field
1650+
if (m == 0.)
1651+
return -deg * mu / T / T / xMath::Pi() / xMath::Pi() * exp(mu/ T);
16101652
return deg * m * m / 2. / T / T / T / T / xMath::Pi() / xMath::Pi()
16111653
* (m * xMath::BesselKexp(1, m / T) - mu * xMath::BesselKexp(2, m / T))
16121654
* exp((mu - m) / T);
@@ -1632,15 +1674,27 @@ namespace thermalfist {
16321674
double sign = 1.;
16331675
double ret = 0.;
16341676
for (int i = 1; i <= order; ++i) {
1635-
ret += sign * (
1677+
if (m == 0.) {
1678+
ret += sign * (
1679+
- (mu - 3. * T/ static_cast<double>(i)) / static_cast<double>(i) / static_cast<double>(i)
1680+
) * cfug;
1681+
}
1682+
else {
1683+
ret += sign * (
16361684
m * xMath::BesselKexp(1, i*moverT)
16371685
- (mu - 3. * T/ static_cast<double>(i)) * xMath::BesselKexp(2, i*moverT)
16381686
) * cfug;
1687+
}
16391688
cfug *= tfug;
16401689
if (signchange) sign = -sign;
16411690
}
1642-
ret *= deg * m * m / T / 2. / xMath::Pi() / xMath::Pi() * xMath::GeVtoifm3();
1643-
return ret;
1691+
double prefactor = 1.;
1692+
if (m == 0.) {
1693+
prefactor = deg * T / xMath::Pi() / xMath::Pi();
1694+
} else {
1695+
prefactor = deg * m * m / T / 2. / xMath::Pi() / xMath::Pi();
1696+
}
1697+
return ret * prefactor * xMath::GeVtoifm3();
16441698
}
16451699

16461700
double QuantumClusterExpansiond2ndT2(int statistics, double T, double mu, double m, double deg, int order,
@@ -1664,15 +1718,24 @@ namespace thermalfist {
16641718
double ret = 0.;
16651719
for (int i = 1; i <= order; ++i) {
16661720
double k = static_cast<double>(i);
1667-
ret += sign * k * (
1721+
if (m == 0.) {
1722+
ret += sign * 2. * T * (mu * mu - 4. * mu * T / k + 6. * T / k * T / k) / k * cfug;
1723+
} else {
1724+
ret += sign * k * (
16681725
(m * (m * m + mu * mu - 4. * mu * T / k + 6. * T * T / k / k) * xMath::BesselKexp(0, i*moverT)
16691726
+ (m * m * (3. * T / k - 2. * mu) + 2. * T / k * (mu * mu - 4. * mu * T / k + 6. * T / k * T / k)) * xMath::BesselKexp(1, i*moverT))
16701727
) * cfug;
1728+
}
16711729
cfug *= tfug;
16721730
if (signchange) sign = -sign;
16731731
}
1674-
ret *= deg * m / T / T / T / 2. / xMath::Pi() / xMath::Pi() * xMath::GeVtoifm3();
1675-
return ret;
1732+
double prefactor = 1.;
1733+
if (m == 0.) {
1734+
prefactor = deg / T / T / 2. / xMath::Pi() / xMath::Pi();
1735+
} else {
1736+
prefactor = deg * m / T / T / T / 2. / xMath::Pi() / xMath::Pi();
1737+
}
1738+
return ret * prefactor * xMath::GeVtoifm3();
16761739
}
16771740

16781741
double QuantumClusterExpansiondedT(int statistics, double T, double mu, double m, double deg, int order,
@@ -1696,15 +1759,24 @@ namespace thermalfist {
16961759
double ret = 0.;
16971760
for (int i = 1; i <= order; ++i) {
16981761
double k = static_cast<double>(i);
1699-
ret += sign * (
1762+
if (m == 0.) {
1763+
ret += sign * 6. * T / k * T / k * (4. * T / k - mu) / k * cfug;
1764+
} else {
1765+
ret += sign * (
17001766
m * (m*m + 3*T/k*(4.*T/k-mu))*xMath::BesselKexp(0, i*moverT)
17011767
+ (6.*T/k*T/k*(4.*T/k-mu) - m*m*(mu-5.*T/k) )*xMath::BesselKexp(1, i*moverT)
17021768
) * cfug;
1769+
}
17031770
cfug *= tfug;
17041771
if (signchange) sign = -sign;
17051772
}
1706-
ret *= deg * m / T / 2. / xMath::Pi() / xMath::Pi() * xMath::GeVtoifm3();
1707-
return ret;
1773+
double prefactor = 1.;
1774+
if (m == 0.) {
1775+
prefactor = deg / 2. / xMath::Pi() / xMath::Pi();
1776+
} else {
1777+
prefactor = deg * m / T / 2. / xMath::Pi() / xMath::Pi();
1778+
}
1779+
return ret * prefactor * xMath::GeVtoifm3();
17081780
}
17091781

17101782
double QuantumClusterExpansiondedmu(int statistics, double T, double mu, double m, double deg, int order,
@@ -1727,12 +1799,21 @@ namespace thermalfist {
17271799
double sign = 1.;
17281800
double ret = 0.;
17291801
for (int i = 1; i <= order; ++i) {
1730-
ret += sign * (xMath::BesselKexp(1, i*moverT) + 3. * xMath::BesselKexp(2, i*moverT) / moverT / static_cast<double>(i)) * cfug / T;
1802+
if (m == 0.) {
1803+
ret += sign * 3. / static_cast<double>(i) / static_cast<double>(i) / static_cast<double>(i) * cfug / T;
1804+
} else {
1805+
ret += sign * (xMath::BesselKexp(1, i*moverT) + 3. * xMath::BesselKexp(2, i*moverT) / moverT / static_cast<double>(i)) * cfug / T;
1806+
}
17311807
cfug *= tfug;
17321808
if (signchange) sign = -sign;
17331809
}
1734-
ret *= deg * m * m * m * T / 2. / xMath::Pi() / xMath::Pi() * xMath::GeVtoifm3();
1735-
return ret;
1810+
double prefactor = 1.;
1811+
if (m == 0.) {
1812+
prefactor = deg * T * T * T * T / xMath::Pi() / xMath::Pi();
1813+
} else {
1814+
prefactor = deg * m * m * m * T / 2. / xMath::Pi() / xMath::Pi();
1815+
}
1816+
return ret * prefactor * xMath::GeVtoifm3();
17361817
}
17371818

17381819

@@ -1757,15 +1838,25 @@ namespace thermalfist {
17571838
double ret = 0.;
17581839
for (int i = 1; i <= order; ++i) {
17591840
double k = static_cast<double>(i);
1760-
ret += sign * k * (
1841+
if (m == 0.) {
1842+
ret += sign * k * (-mu / k / k) * cfug;
1843+
}
1844+
else {
1845+
ret += sign * k * (
17611846
m * xMath::BesselKexp(1, i*moverT)
17621847
- mu * xMath::BesselKexp(2, i*moverT)
17631848
) * cfug;
1849+
}
17641850
cfug *= tfug;
17651851
if (signchange) sign = -sign;
17661852
}
1767-
ret *= deg * m * m / 2. / T / T / T / T / xMath::Pi() / xMath::Pi();
1768-
return ret;
1853+
double prefactor = 1.;
1854+
if (m == 0.) {
1855+
prefactor = deg / T / T / xMath::Pi() / xMath::Pi();
1856+
} else {
1857+
prefactor = deg * m * m / 2. / T / T / T / T / xMath::Pi() / xMath::Pi();
1858+
}
1859+
return ret * prefactor;
17691860
}
17701861

17711862
double QuantumNumericalIntegrationdndT(int statistics, double T, double mu, double m, double deg,

‎src/library/HRGBase/ThermalParticle.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -58,8 +58,8 @@ namespace thermalfist {
5858

5959
bool ThermalParticle::ZeroWidthEnforced() const
6060
{
61-
//if (PdgId() == 223) // omega(782)
62-
// return true;
61+
if (Mass() == 0.)
62+
return true;
6363
return ((ResonanceWidth() / Mass()) < 0.01);
6464
}
6565

@@ -162,7 +162,8 @@ namespace thermalfist {
162162

163163
void ThermalParticle::CalculateThermalBranchingRatios(const ThermalModelParameters & params, bool useWidth, double mu)
164164
{
165-
if (!useWidth || m_Width == 0.0 || m_Width / m_Mass < 1.e-2 || m_ResonanceWidthIntegrationType != eBW) {
165+
//if (!useWidth || m_Width == 0.0 || m_Width / m_Mass < 1.e-2 || m_ResonanceWidthIntegrationType != eBW) {
166+
if (!useWidth || m_Width == 0.0 || ZeroWidthEnforced() || m_ResonanceWidthIntegrationType != eBW) {
166167
for (size_t j = 0; j < m_Decays.size(); ++j) {
167168
m_Decays[j].mBratioAverage = m_Decays[j].mBratio;
168169
}
@@ -498,7 +499,7 @@ namespace thermalfist {
498499
{
499500
//if (m_ResonanceWidthIntegrationType != eBW)
500501
// return m_Width;
501-
if (m_Width / m_Mass < 0.01) {
502+
if (ZeroWidthEnforced()) {
502503
return m_Width;
503504
}
504505

@@ -522,7 +523,7 @@ namespace thermalfist {
522523
{
523524
std::vector<double> ret(m_Decays.size(), 0.);
524525

525-
if (!eBW || m_Width / m_Mass < 0.01) {
526+
if (!eBW || ZeroWidthEnforced()) {
526527
for (size_t i = 0; i < m_Decays.size(); ++i)
527528
ret[i] = m_Decays[i].mBratio;
528529

@@ -650,7 +651,6 @@ namespace thermalfist {
650651
if (!(params.gammaS == 1. || m_AbsS == 0.)) mu += log(params.gammaS) * m_AbsS * params.T;
651652
if (!(params.gammaC == 1. || m_AbsC == 0.)) mu += log(params.gammaC) * m_AbsC * params.T;
652653

653-
//if (!useWidth || m_Width / m_Mass < 1.e-2 || m_ResonanceWidthIntegrationType == ZeroWidth) {
654654
if (!useWidth || ZeroWidthEnforced() || m_ResonanceWidthIntegrationType == ZeroWidth) {
655655
return mn * IdealGasFunctions::IdealGasQuantity(type, m_QuantumStatisticsCalculationType, 0, params.T / static_cast<double>(n), mu, m_Mass, m_Degeneracy, 1, m_IGFExtraConfig);
656656
}

‎src/library/HRGEventGenerator/EventGeneratorBase.cpp‎

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1251,6 +1251,11 @@ namespace thermalfist {
12511251
bool flOverlap = true;
12521252
while (flOverlap) {
12531253
int sampled = 0;
1254+
1255+
// Check if the event has no particles at all
1256+
if(ids.size() == 0)
1257+
break;
1258+
12541259
while (sampled < ids.size()) {
12551260
flOverlap = false;
12561261
int i = ids[sampled];
@@ -1517,7 +1522,7 @@ namespace thermalfist {
15171522
}
15181523
}
15191524

1520-
// Sample lifetime in rest frame
1525+
// Sample the lifetime in the rest frame
15211526
ct = -ct * log(1. - RandomGenerators::randgenMT.randDblExc());
15221527

15231528
// Lorentz time delay
@@ -1564,7 +1569,7 @@ namespace thermalfist {
15641569
primParticles[i][j].processed = true;
15651570
}
15661571
else {
1567-
// Decay through unknown branching ratio, presumably radiative, no hadrons, just ignore decay products
1572+
// Decay through unknown branching ratio, presumably radiative, no hadrons, just ignore the decay products
15681573
primParticles[i][j].processed = true;
15691574
}
15701575
}
@@ -2032,4 +2037,4 @@ namespace thermalfist {
20322037
fUseGCEConservedCharges = false;
20332038
}
20342039

2035-
} // namespace thermalfist
2040+
} // namespace thermalfist

‎src/library/HRGEventGenerator/RandomGenerators.cpp‎

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -418,6 +418,8 @@ namespace thermalfist {
418418
}
419419

420420
void BreitWignerGenerator::FixParameters() {
421+
if (m_Gamma < 1e-7)
422+
return;
421423
double eps = 1e-8;
422424
double l = m_Mthr, r = m_M + 2.*m_Gamma;
423425
double m1 = l + (r - l) / 3.;
@@ -466,6 +468,12 @@ namespace thermalfist {
466468

467469
void ThermalBreitWignerGenerator::FixParameters()
468470
{
471+
if (m_part->ZeroWidthEnforced()) {
472+
m_Xmin = m_part->Mass();
473+
m_Xmax = m_part->Mass();
474+
m_Max = 1.;
475+
return;
476+
}
469477
double Threshold = m_part->DecayThresholdMass();
470478
double Width = m_part->ResonanceWidth();
471479
double Mass = m_part->Mass();
@@ -494,7 +502,7 @@ namespace thermalfist {
494502

495503
double ThermalBreitWignerGenerator::GetRandom() const
496504
{
497-
if (m_part->ResonanceWidth() / m_part->Mass() < 1.e-2)
505+
if (m_part->ZeroWidthEnforced())
498506
return m_part->Mass();
499507
while (true) {
500508
double x0 = m_Xmin + (m_Xmax - m_Xmin) * randgenMT.rand();
@@ -506,6 +514,12 @@ namespace thermalfist {
506514

507515
void ThermalEnergyBreitWignerGenerator::FixParameters()
508516
{
517+
if (m_part->ZeroWidthEnforced()) {
518+
m_Xmin = m_part->Mass();
519+
m_Xmax = m_part->Mass();
520+
m_Max = 1.;
521+
return;
522+
}
509523
double Threshold = m_part->DecayThresholdMass();
510524
double Width = m_part->ResonanceWidth();
511525
double Mass = m_part->Mass();

0 commit comments

Comments
 (0)