CbmRoot
ThermalModelEVMF.cxx
Go to the documentation of this file.
1 #include "ThermalModelEVMF.h"
2 #include "xMath.h"
3 #include <cmath>
4 #include <cstdlib>
5 #include <iostream>
6 #include <vector>
7 
8 using namespace std;
9 
10 namespace ThermalModelEVMFNamespace {
11  ThermalModelEVMF* gThM;
12 
13  void broyden22(vector<double>& xin,
14  vector<double> (*func)(const vector<double>&)) {
15  const double TOLF = 1.0e-8, EPS = 1.0e-8;
16  const int MAXITS = 200;
17  vector<double> h = xin, xh1 = xin, xh2 = xin;
18  for (unsigned int i = 0; i < xin.size(); ++i) {
19  h[i] = EPS * fabs(h[i]);
20  if (h[i] == 0.0) h[i] = EPS;
21  if (i == 0)
22  xh1[i] = xin[i] + h[i];
23  else
24  xh2[i] = xin[i] + h[i];
25  }
26  vector<double> r1 = func(xin), r21 = func(xh1), r22 = func(xh2);
27  if (fabs(r1[0]) < TOLF && fabs(r1[1]) < TOLF) return;
28  double J[2][2];
29  J[0][0] = (r21[0] - r1[0]) / h[0];
30  J[0][1] = (r22[0] - r1[0]) / h[1];
31  J[1][0] = (r21[1] - r1[1]) / h[0];
32  J[1][1] = (r22[1] - r1[1]) / h[1];
33  double Jinv[2][2];
34  double det = J[0][0] * J[1][1] - J[0][1] * J[1][0];
35  if (det == 0.0) throw("singular Jacobian in broyden22");
36  Jinv[0][0] = J[1][1] / det;
37  Jinv[0][1] = -J[0][1] / det;
38  Jinv[1][0] = -J[1][0] / det;
39  Jinv[1][1] = J[0][0] / det;
40  vector<double> xold = xin;
41  vector<double> rold = r1;
42  vector<double> rprevten = r1;
43  for (int iter = 1; iter <= MAXITS; ++iter) {
44  xin[0] = xold[0] - Jinv[0][0] * rold[0] - Jinv[0][1] * rold[1];
45  xin[1] = xold[1] - Jinv[1][0] * rold[0] - Jinv[1][1] * rold[1];
46  r1 = func(xin);
47  if (fabs(r1[0]) < TOLF && fabs(r1[1]) < TOLF) return;
48  if (iter % 10 == 0) {
49  if (fabs(r1[0] / rprevten[0]) > 1e-2
50  && fabs(r1[1] / rprevten[1]) > 1e-2)
51  throw("broyden is stuck");
52  rprevten = r1;
53  }
54  double JF[2];
55  double DF[2];
56  double dx[2];
57  DF[0] = (r1[0] - rold[0]);
58  DF[1] = (r1[1] - rold[1]);
59  JF[0] = Jinv[0][0] * DF[0] + Jinv[0][1] * DF[1];
60  JF[1] = Jinv[1][0] * DF[0] + Jinv[1][1] * DF[1];
61  dx[0] = xin[0] - xold[0];
62  dx[1] = xin[1] - xold[1];
63  double znam = dx[0] * JF[0] + dx[1] * JF[1];
64  if (znam == 0.0) throw("singular Jacobian in broyden22");
65  double xJ[2];
66  JF[0] = dx[0] - JF[0];
67  JF[1] = dx[1] - JF[1];
68  xJ[0] = dx[0] * Jinv[0][0] + dx[1] * Jinv[1][0];
69  xJ[1] = dx[0] * Jinv[0][1] + dx[1] * Jinv[1][1];
70  Jinv[0][0] = Jinv[0][0] + JF[0] * xJ[0] / znam;
71  Jinv[0][1] = Jinv[0][1] + JF[0] * xJ[1] / znam;
72  Jinv[1][0] = Jinv[1][0] + JF[1] * xJ[0] / znam;
73  Jinv[1][1] = Jinv[1][1] + JF[1] * xJ[1] / znam;
74  xold = xin;
75  rold = r1;
76  }
77  throw("exceed MAXITS in broyden");
78  }
79 
80  vector<double> function22(const vector<double>& xin) {
81  vector<double> ret(2);
82  gThM->Parameters.muQ = xin[0];
83  gThM->Parameters.muS = xin[1];
84  gThM->CalculateDensities();
85  double fBd = gThM->CalculateBaryonDensity();
86  double fCd = gThM->CalculateChargeDensity();
87  double fSd = gThM->CalculateStrangenessDensity();
88  double fASd = gThM->CalculateAbsoluteStrangenessDensity();
89  ret[0] = (fCd / fBd - gThM->QBgoal) / gThM->QBgoal;
90  ret[1] = fSd / fASd;
91  return ret;
92  }
93 }; // namespace ThermalModelEVMFNamespace
94 
95 using namespace ThermalModelEVMFNamespace;
96 
97 
98 ThermalModelEVMF::~ThermalModelEVMF(void) {}
99 
100 void ThermalModelEVMF::FixParameters() {
101  if (fabs(Parameters.muB) < 1e-6) {
102  Parameters.muS = Parameters.muQ = 0.;
103  return;
104  }
105  vector<double> x22(2);
106  x22[0] = Parameters.muQ;
107  x22[1] = Parameters.muS;
108  vector<double> x2(2), xinit(2);
109  xinit[0] = x2[0] = Parameters.muQ;
110  xinit[1] = x2[1] = Parameters.muS;
111  ThermalModelEVMFNamespace::gThM = this;
112  int iter = 0, iterMAX = 10;
113  while (iter < iterMAX) {
114  try {
115  broyden22(x22, ThermalModelEVMFNamespace::function22);
116  } catch (...) {}
117  if (fabs(Parameters.muB / Parameters.muS) > 1.
118  && fabs(Parameters.muB / Parameters.muQ) > 1.)
119  break;
120  xinit[0] = x2[0] = xinit[0] / 10.;
121  xinit[1] = x2[1] = xinit[1] / 10.;
122  iter++;
123  }
124  ThermalModelEVMFNamespace::gThM = NULL;
125 }
126 
127 
128 void ThermalModelEVMF::FixParameters(double QB) {
129  QBgoal = QB;
130  FixParameters();
131 }
132 
133 double ThermalModelEVMF::Density(double n) {
134  double ret = 0.;
135 
136  for (unsigned int i = 0; i < TPS->fParticles.size(); ++i) {
137  double vo = 0.;
138  vo = 4. * 4. * xMath::Pi() / 3. * RHad * RHad * RHad;
139  double dMu = -UVdW(n, Parameters.T, vo);
140  ret += TPS->fParticles[i].CalculateDensity(Parameters.T,
141  Parameters.muB,
142  Parameters.muS,
143  Parameters.muQ,
144  Parameters.gammaS,
145  0,
146  fUseWidth,
147  dMu);
148  }
149 
150  if (fUseHagedorn) {
151  double vo = 0.;
152  vo = 4. * 4. * xMath::Pi() / 3. * RHad * RHad * RHad;
153  double dMu = -UVdW(n, Parameters.T, vo);
154  ret += fHag.CalculateDensity(Parameters.T, Parameters.muB, 0, dMu);
155  }
156  return ret;
157 }
158 
159 void ThermalModelEVMF::CalculateDensities() {
160  fDensity = 0.;
161  {
162  double vo = 4. * 4. * xMath::Pi() / 3. * RHad * RHad * RHad;
163  if (fDensity * vo > 1.) fDensity = 1. / vo / 2.;
164  const double TOLF = 1.0e-8, EPS = 1.0e-8;
165  const int MAXITS = 200;
166  double x = fDensity;
167  double h = x;
168  h = EPS * fabs(h);
169  if (h == 0.0) h = EPS;
170  double xh = x + h;
171  double r1 = x - Density(x);
172  double r2 = xh - Density(xh);
173  double Jinv = h / (r1 - r2);
174  double xold = x, rold = r1;
175  for (int iter = 1; iter <= MAXITS; ++iter) {
176  x = xold - Jinv * rold;
177  if (fMode == 0 && x * vo > 1.)
178  x = 1. / vo * (rand() / static_cast<double>(RAND_MAX));
179  if (fMode == 1 && x * vo > 0.5)
180  x = 0.5 / vo * (rand() / static_cast<double>(RAND_MAX));
181  r1 = x - Density(x);
182  Jinv = (x - xold) / (r1 - rold);
183  if (fabs(r1) < TOLF) break;
184  xold = x;
185  rold = r1;
186  }
187  fDensity = x;
188  }
189  for (unsigned int i = 0; i < TPS->fParticles.size(); ++i) {
190  double vo = 0.;
191  vo = 4. * 4. * xMath::Pi() / 3. * RHad * RHad * RHad;
192  double dMu = -UVdW(fDensity, Parameters.T, vo);
193  densitiesid[i] = TPS->fParticles[i].CalculateDensity(Parameters.T,
194  Parameters.muB,
195  Parameters.muS,
196  Parameters.muQ,
197  Parameters.gammaS,
198  0,
199  fUseWidth,
200  dMu);
201  }
202 
203  if (fUseHagedorn) {
204  double vo = 0.;
205  vo = 4. * 4. * xMath::Pi() / 3. * RHad * RHad * RHad;
206  double dMu = -UVdW(fDensity, Parameters.T, vo);
207  fHagedornDensity =
208  fHag.CalculateDensity(Parameters.T, Parameters.muB, 0, dMu);
209  }
210 
211  for (unsigned int i = 0; i < TPS->fParticles.size(); ++i)
212  densities[i] = densitiesid[i]; // / fSuppression;
213  for (int i = 0; i < static_cast<int>(TPS->fParticles.size()); ++i) {
214  {
215  densitiestotal[i] = densities[i];
216  for (int j = 0;
217  j < static_cast<int>(TPS->fParticles[i].ResonanceBR.size());
218  ++j)
219  if (i != TPS->fParticles[i].ResonanceBR[j].second)
220  densitiestotal[i] +=
221  TPS->fParticles[i].ResonanceBR[j].first
222  * densities[TPS->fParticles[i].ResonanceBR[j].second];
223  }
224  }
225  fCalculated = true;
226 }
227 
228 double ThermalModelEVMF::CalculateShearViscosity() {
229  if (!fCalculated) CalculateDensities();
230  double ret = 0.;
231  double TotalDensity = 0.;
232  for (unsigned int i = 0; i < TPS->fParticles.size(); ++i) {
233  TotalDensity += densities[i];
234  ret += sqrt(TPS->fParticles[i].fMass)
235  * xMath::BesselK(5. / 2., TPS->fParticles[i].fMass / Parameters.T)
236  / xMath::BesselK(2., TPS->fParticles[i].fMass / Parameters.T)
237  * densities[i];
238  }
239 
240  return ret * 5. * sqrt(Parameters.T) / 64. / RHad / RHad / sqrt(xMath::Pi())
241  * 5.06 / TotalDensity;
242 }
ThermalModelParameters::muS
double muS
Definition: ThermalModelBase.h:7
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
ThermalModelEVMF::CalculateShearViscosity
virtual double CalculateShearViscosity()
Definition: ThermalModelEVMF.cxx:228
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
ThermalModelBase::QBgoal
double QBgoal
Definition: ThermalModelBase.h:27
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
ThermalModelEVMFNamespace
Definition: ThermalModelEVMF.cxx:10
ThermalModelEVMFNamespace::function22
vector< double > function22(const vector< double > &xin)
Definition: ThermalModelEVMF.cxx:80
ThermalModelEVMF::Density
double Density(double n)
Definition: ThermalModelEVMF.cxx:133
ThermalModelEVMF::~ThermalModelEVMF
virtual ~ThermalModelEVMF(void)
Definition: ThermalModelEVMF.cxx:98
ThermalModelEVMF::CalculateAbsoluteStrangenessDensity
virtual double CalculateAbsoluteStrangenessDensity()
Definition: ThermalModelEVMF.h:216
ThermalModelEVMF::CalculateChargeDensity
virtual double CalculateChargeDensity()
Definition: ThermalModelEVMF.h:195
ThermalModelEVMF::CalculateBaryonDensity
virtual double CalculateBaryonDensity()
Definition: ThermalModelEVMF.h:188
xMath::Pi
double Pi()
Definition: xMath.h:5
ThermalModelEVMF::FixParameters
virtual void FixParameters()
Definition: ThermalModelEVMF.cxx:100
xMath::BesselK
double BesselK(int n, double x)
Definition: xMath.cxx:180
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
ThermalModelEVMF.h
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
ThermalModelEVMFNamespace::gThM
ThermalModelEVMF * gThM
Definition: ThermalModelEVMF.cxx:11
xMath.h
ThermalModelEVMF::CalculateDensities
virtual void CalculateDensities()
Definition: ThermalModelEVMF.cxx:159
ThermalModelEVMFNamespace::broyden22
void broyden22(vector< double > &xin, vector< double >(*func)(const vector< double > &))
Definition: ThermalModelEVMF.cxx:13
ThermalModelParameters::muQ
double muQ
Definition: ThermalModelBase.h:7
ThermalModelBase::Parameters
ThermalModelParameters Parameters
Definition: ThermalModelBase.h:21
ThermalModelEVMF
Definition: ThermalModelEVMF.h:12
ThermalModelEVMF::CalculateStrangenessDensity
virtual double CalculateStrangenessDensity()
Definition: ThermalModelEVMF.h:202