CbmRoot
ThermalParticle.cxx
Go to the documentation of this file.
1 #include "ThermalParticle.h"
2 #include "xMath.h"
3 #include <cmath>
4 #include <fstream>
5 #include <iostream>
6 #include <sstream>
7 
8 
9 namespace ThermalParticleNamespace {
10 #include "NumericalIntegration.h"
11 }
12 
13 using namespace std;
14 using namespace ThermalParticleNamespace;
15 
16 const double ThermalParticle::GeVtoifm = 1. / 0.197; //5.06773;
17 
18 ThermalParticle::~ThermalParticle(void) {}
19 
20 void ThermalParticle::ReadDecays(string filename) {
21  fDecays.resize(0);
22  ifstream fin(filename.c_str());
23  if (fin.is_open()) {
24  char cc[400];
25  double tmpbr;
26  while (fin >> tmpbr) {
27  ParticleDecay decay;
28  decay.fBratio = tmpbr / 100.;
29  fin.getline(cc, 350);
30  stringstream ss;
31  ss << cc;
32  int tmpid;
33  while (ss >> tmpid) {
34  decay.fDaughters.push_back(tmpid);
35  }
36  fDecays.push_back(decay);
37  }
38  }
39 }
40 
41 void ThermalParticle::NormalizeBranchingRatios() {
42  double sum = 0.;
43  for (unsigned int i = 0; i < fDecays.size(); ++i) {
44  sum += fDecays[i].fBratio;
45  }
46  for (unsigned int i = 0; i < fDecays.size(); ++i) {
47  fDecays[i].fBratio *= 1. / sum;
48  }
49 }
50 
51 double ThermalParticle::CalculateParticleDensity(double T,
52  double muB,
53  double muS,
54  double muQ,
55  double gammaS,
56  double mass,
57  double dMu) const {
58  if (fStatistics == 0) {
59  if (fabs(gammaS - 1.) < 1e-5)
60  return fSpinDegeneracy * xMath::BesselKexp(2, mass / T)
61  * exp((fB * muB + fS * muS + fC * muQ + dMu - mass) / T) / 2.
62  / xMath::Pi() / xMath::Pi() * mass * mass * T * GeVtoifm * GeVtoifm
63  * GeVtoifm;
64  else
65  return fSpinDegeneracy * xMath::BesselKexp(2, mass / T)
66  * pow(gammaS, fAbsS)
67  * exp((fB * muB + fS * muS + fC * muQ + dMu - mass) / T) / 2.
68  / xMath::Pi() / xMath::Pi() * mass * mass * T * GeVtoifm * GeVtoifm
69  * GeVtoifm;
70  }
71  double eps = 1e-2;
72  double ret = 0.;
73  double cur = 0.; //prev = 0., cur = 0.;
74  if (mass < 0.) mass = fMass;
75  muB /= T;
76  muS /= T;
77  muQ /= T;
78  dMu /= T;
79  double tMu = fB * muB + fS * muS + fC * muQ + dMu;
80  double tMass = mass / T;
81  double Sign = 1.;
82  if (tMu > tMass) return 0.;
83  ret = 0.;
84  if (tMass > 70.) return 0.;
85  for (int i = 1; i < 10; ++i) {
86  cur = xMath::BesselK(2, i * tMass)
87  * exp(tMu * i - log(pow(gammaS, -fAbsS)) * i) / i;
88  ret += Sign * cur;
89  if (i != 1 && cur / ret < eps) break;
90  //prev = cur;
91  if (fStatistics == 1) Sign *= -1.;
92  if (fStatistics == 0) break;
93  }
94  return ret * fSpinDegeneracy / 2. / xMath::Pi() / xMath::Pi() * mass * mass
95  * T * GeVtoifm * GeVtoifm * GeVtoifm;
96 }
97 
98 double ThermalParticle::CalculateDensity(double T,
99  double muB,
100  double muS,
101  double muQ,
102  double gammaS,
103  int type,
104  bool useWidth,
105  double dMu) const {
106  if (!useWidth || fWidth / fMass < 1e-2) {
107  if (type == 0)
108  return CalculateParticleDensity(T, muB, muS, muQ, gammaS, fMass, dMu);
109  if (type == 1)
110  return CalculateEnergyDensity(T, muB, muS, muQ, gammaS, fMass, dMu);
111  if (type == 2)
112  return CalculateEntropyDensity(T, muB, muS, muQ, gammaS, fMass, dMu);
113  if (type == 3)
114  return CalculatePressure(T, muB, muS, muQ, gammaS, fMass, dMu);
115  return 0;
116  }
117  vector<double> x, w;
118  double a = max(fThreshold, fMass - 2. * fWidth), b = fMass + 2. * fWidth;
119  int ind = 5;
120  if (fWidth / fMass < 1e-1) {
121  ind = 10;
122  GetCoefsIntegrateLegendre10(a, b, x, w);
123  } else {
124  ind = 32;
125  GetCoefsIntegrateLegendre32(a, b, x, w);
126  }
127  double ret1 = 0., ret2 = 0., tmp = 0.;
128 
129  for (int i = 0; i < ind; i++) {
130  tmp = w[i] * fMass * fWidth * x[i]
131  / ((x[i] * x[i] - fMass * fMass) * (x[i] * x[i] - fMass * fMass)
132  + fMass * fMass * fWidth * fWidth);
133  double dens = 0.;
134  if (type == 0)
135  dens = CalculateParticleDensity(T, muB, muS, muQ, gammaS, x[i], dMu);
136  if (type == 1)
137  dens = CalculateEnergyDensity(T, muB, muS, muQ, gammaS, x[i], dMu);
138  if (type == 2)
139  dens = CalculateEntropyDensity(T, muB, muS, muQ, gammaS, x[i], dMu);
140  if (type == 3)
141  dens = CalculatePressure(T, muB, muS, muQ, gammaS, x[i], dMu);
142  ret1 += tmp * dens;
143  ret2 += tmp;
144  }
145  return ret1 / ret2;
146 }
147 
148 double ThermalParticle::CalculateParticleDensityLaguerre(double T,
149  double muB,
150  double muS,
151  double muQ,
152  double gammaS,
153  double mass,
154  double dMu) const {
155  if (fStatistics == 0)
156  return CalculateParticleDensity(T, muB, muS, muQ, gammaS, mass);
157  double ret = 0.;
158  if (mass < 0.) mass = fMass;
159  muB /= T;
160  muS /= T;
161  muQ /= T;
162  dMu /= T;
163  double tMu = fB * muB + fS * muS + fC * muQ + dMu;
164  double tMass = mass / T;
165  ret = 0.;
166 
167  vector<double> x, w;
168  GetCoefsIntegrateLaguerre32(x, w);
169 
170  if (tMu > tMass) return 0.;
171 
172  for (int i = 0; i < 32; i++) {
173  ret += w[i] * x[i] * x[i]
174  / (exp(sqrt(x[i] * x[i] + tMass * tMass) - tMu) * pow(gammaS, -fAbsS)
175  + fStatistics);
176  }
177 
178  return ret * fSpinDegeneracy / 2. / xMath::Pi() / xMath::Pi() * T * T * T
179  * GeVtoifm * GeVtoifm * GeVtoifm;
180 }
181 
182 double ThermalParticle::CalculatePressure(double T,
183  double muB,
184  double muS,
185  double muQ,
186  double gammaS,
187  double mass,
188  double dMu) const {
189  if (fStatistics == 0)
190  return T * CalculateParticleDensity(T, muB, muS, muQ, mass, dMu);
191  double ret = 0.;
192  if (mass < 0.) mass = fMass;
193  muB /= T;
194  muS /= T;
195  muQ /= T;
196  dMu /= T;
197  double tMu = fB * muB + fS * muS + fC * muQ + dMu;
198  double tMass = mass / T;
199  //double Sign = 1.;
200  ret = 0.;
201 
202  vector<double> x, w;
203  GetCoefsIntegrateLaguerre32(x, w);
204  //double ret = 0.;
205 
206  if (tMu > tMass) return 0.; //cout << "AAAAAAA";
207 
208  for (int i = 0; i < 32; i++) {
209  double x2 = x[i] * x[i];
210  double E = sqrt(x[i] * x[i] + tMass * tMass);
211  ret +=
212  w[i] * x2 * x2 / (exp(E - tMu) * pow(gammaS, -fAbsS) + fStatistics) / E;
213  }
214 
215  return ret * fSpinDegeneracy / 6. / xMath::Pi() / xMath::Pi() * T * T * T * T
216  * GeVtoifm * GeVtoifm * GeVtoifm;
217 }
218 
219 double ThermalParticle::CalculateEntropyDensity(double T,
220  double muB,
221  double muS,
222  double muQ,
223  double gammaS,
224  double mass,
225  double dMu) const {
226  if (fStatistics == 0)
227  return (CalculateEnergyDensity(T, muB, muS, muQ, gammaS, mass, dMu)
228  + CalculatePressure(T, muB, muS, muQ, gammaS, mass, dMu)
229  - (fB * muB + fS * muS + fC * muQ + dMu)
230  * CalculateParticleDensity(T, muB, muS, muQ, gammaS, mass, dMu))
231  / T;
232  double ret = 0.;
233  if (mass < 0.) mass = fMass;
234  muB /= T;
235  muS /= T;
236  muQ /= T;
237  dMu /= T;
238  double tMu = fB * muB + fS * muS + fC * muQ + dMu;
239  double tMass = mass / T;
240  //double Sign = 1.;
241  ret = 0.;
242 
243  vector<double> x, w;
244  GetCoefsIntegrateLaguerre32(x, w);
245  //double ret = 0.;
246 
247  if (tMu > tMass) return 0.; //cout << "AAAAAAA";
248 
249  for (int i = 0; i < 32; i++) {
250  double x2 = x[i] * x[i];
251  double E = sqrt(x[i] * x[i] + tMass * tMass);
252  ret += w[i] * x2 / (exp(E - tMu) * pow(gammaS, -fAbsS) + fStatistics)
253  * (x2 / 3. / E + E - tMu);
254  }
255 
256  return ret * fSpinDegeneracy / 2. / xMath::Pi() / xMath::Pi() * T * T * T
257  * GeVtoifm * GeVtoifm * GeVtoifm;
258 }
259 
260 double ThermalParticle::CalculateEnergyDensity(double T,
261  double muB,
262  double muS,
263  double muQ,
264  double gammaS,
265  double mass,
266  double dMu) const {
267  if (fStatistics == 0)
268  return (3 * T
269  + mass * xMath::BesselK1exp(mass / T)
270  / xMath::BesselKexp(2, mass / T))
271  * CalculateParticleDensity(T, muB, muS, muQ, mass, dMu);
272  double ret = 0.;
273  if (mass < 0.) mass = fMass;
274  muB /= T;
275  muS /= T;
276  muQ /= T;
277  dMu /= T;
278  double tMu = fB * muB + fS * muS + fC * muQ + dMu;
279  double tMass = mass / T;
280  ret = 0.;
281 
282  vector<double> x, w;
283  GetCoefsIntegrateLaguerre32(x, w);
284 
285  if (tMu > tMass) return 0.;
286 
287  for (int i = 0; i < 32; i++) {
288  //double x2 = x[i] * x[i];
289  double E = sqrt(x[i] * x[i] + tMass * tMass);
290  ret += w[i] * x[i] * x[i] * E
291  / (exp(E - tMu) * pow(gammaS, -fAbsS) + fStatistics);
292  }
293 
294  return ret * fSpinDegeneracy / 2. / xMath::Pi() / xMath::Pi() * T * T * T * T
295  * GeVtoifm * GeVtoifm * GeVtoifm;
296 }
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:134
ThermalParticle::CalculateEntropyDensity
double CalculateEntropyDensity(double T, double muB, double muS, double muQ, double gammaS, double mass=-1., double dMu=0.) const
Definition: ThermalParticle.cxx:219
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
ThermalParticle::CalculateEnergyDensity
double CalculateEnergyDensity(double T, double muB, double muS, double muQ, double gammaS, double mass=-1., double dMu=0.) const
Definition: ThermalParticle.cxx:260
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
GetCoefsIntegrateLaguerre32
void GetCoefsIntegrateLaguerre32(std::vector< double > &x, std::vector< double > &w)
Definition: HRGModel/NumericalIntegration.h:762
ThermalParticle::CalculateDensity
double CalculateDensity(double T, double muB, double muS, double muQ, double gammaS, int type=0, bool useWidth=0, double dMu=0.) const
Definition: ThermalParticle.cxx:98
Cbm::Sign
int Sign(const T &x)
Definition: CbmUtils.h:36
ThermalParticle::NormalizeBranchingRatios
void NormalizeBranchingRatios()
Definition: ThermalParticle.cxx:41
ParticleDecay::fBratio
double fBratio
Definition: ThermalParticle.h:8
GetCoefsIntegrateLegendre32
void GetCoefsIntegrateLegendre32(double a, double b, std::vector< double > &x, std::vector< double > &w)
Definition: HRGModel/NumericalIntegration.h:477
GetCoefsIntegrateLegendre10
void GetCoefsIntegrateLegendre10(double a, double b, std::vector< double > &x, std::vector< double > &w)
Definition: HRGModel/NumericalIntegration.h:570
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
ThermalParticleNamespace
Definition: ThermalParticle.cxx:9
xMath::Pi
double Pi()
Definition: xMath.h:5
ParticleDecay
Definition: ThermalParticle.h:7
ThermalModelNoFlowNamespace::GeVtoifm
const Double_t GeVtoifm
Definition: CbmThermalModelNoFlow.cxx:74
ParticleDecay::fDaughters
std::vector< int > fDaughters
Definition: ThermalParticle.h:9
xMath::BesselK1exp
double BesselK1exp(double x)
Definition: xMath.cxx:704
ThermalParticle::GeVtoifm
static const double GeVtoifm
Definition: ThermalParticle.h:17
xMath::BesselK
double BesselK(int n, double x)
Definition: xMath.cxx:180
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
xMath::BesselKexp
double BesselKexp(int n, double x)
Definition: xMath.cxx:747
ThermalParticle::ReadDecays
void ReadDecays(std::string filename="")
Definition: ThermalParticle.cxx:20
ThermalParticle::CalculatePressure
double CalculatePressure(double T, double muB, double muS, double muQ, double gammaS, double mass=-1., double dMu=0.) const
Definition: ThermalParticle.cxx:182
ThermalParticle::CalculateParticleDensity
double CalculateParticleDensity(double T, double muB, double muS, double muQ, double gammaS, double mass=-1., double dMu=0.) const
Definition: ThermalParticle.cxx:51
xMath.h
ThermalParticle.h
NumericalIntegration.h
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:36
ThermalParticle::CalculateParticleDensityLaguerre
double CalculateParticleDensityLaguerre(double T, double muB, double muS, double muQ, double gammaS, double mass=-1., double dMu=0.) const
Definition: ThermalParticle.cxx:148
ThermalParticle::~ThermalParticle
~ThermalParticle(void)
Definition: ThermalParticle.cxx:18