Thermal-FIST 1.5
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
ThermalModelRealGas.cpp
Go to the documentation of this file.
1#include "HRGRealGas/ThermalModelRealGas.h"
2
3#include <vector>
4#include <cmath>
5#include <iostream>
6#include <ctime>
7#include <iostream>
8
9#include <Eigen/Dense>
10
11#include "HRGBase/xMath.h"
12#include "HRGEV/ExcludedVolumeHelper.h"
13
14#ifdef USE_OPENMP
15#include <omp.h>
16#include "ThermalModelRealGas.h"
17#endif
18
19
20using namespace Eigen;
21using namespace std;
22
23namespace thermalfist {
24
25 ThermalModelRealGas::ThermalModelRealGas(ThermalParticleSystem* TPS_, const ThermalModelParameters& params) :
26 ThermalModelBase(TPS_, params), m_SearchMultipleSolutions(false), m_ComponentMapCalculated(false)
27 {
28 m_chi.resize(6);
29 for (int i = 0; i < 6; ++i) m_chi[i].resize(3);
30 m_chiarb.resize(6);
31 m_DensitiesId.resize(m_densities.size());
32 m_MuStar.resize(m_densities.size());
33 m_Volume = params.V;
34 m_TAG = "ThermalModelRealGas";
35
36 m_exvolmodideal = new ExcludedVolumeModelMultiBase(m_densities.size());
37 m_mfmodideal = new MeanFieldModelMultiBase(m_densities.size());
38 m_exvolmod = m_exvolmodideal;
39 m_mfmod = m_mfmodideal;
40
41 m_Ensemble = GCE;
42 m_InteractionModel = RealGas;
43 }
44
45 ThermalModelRealGas::~ThermalModelRealGas(void)
46 {
47 if (m_exvolmodideal != NULL) {
48 delete m_exvolmodideal;
49 if (m_exvolmod == m_exvolmodideal)
50 m_exvolmod = NULL;
51 m_exvolmodideal = NULL;
52 }
53 if (m_mfmodideal != NULL) {
54 delete m_mfmodideal;
55 if (m_mfmod == m_mfmodideal)
56 m_mfmod = NULL;
57 m_mfmodideal = NULL;
58 }
59 if (m_exvolmod != NULL) {
60 delete m_exvolmod;
61 m_exvolmod = NULL;
62 }
63 if (m_mfmod != NULL) {
64 delete m_mfmod;
65 m_mfmod = NULL;
66 }
67 }
68
69 void ThermalModelRealGas::FillChemicalPotentials() {
70 ThermalModelBase::FillChemicalPotentials();
71 for (size_t i = 0; i < m_MuStar.size(); ++i)
72 m_MuStar[i] = m_Chem[i];
73 }
74
75 void ThermalModelRealGas::SetChemicalPotentials(const vector<double>& chem)
76 {
77 ThermalModelBase::SetChemicalPotentials(chem);
78 for (size_t i = 0; i < m_MuStar.size(); ++i)
79 m_MuStar[i] = m_Chem[i];
80 }
81
82 void ThermalModelRealGas::CalculateComponentsMap()
83 {
84 map< vector<int>, int> MapComp;
85
86 int NN = m_densities.size();
87
88 m_MapTodMuStar.resize(NN);
89 m_MapFromdMuStar.clear();
90 MapComp.clear();
91 m_dMuStarIndices.clear();
92
93 int tind = 0;
94 for (int i = 0; i < NN; ++i) {
95 vector<int> IntParam(0);
96 if (m_exvolmod == NULL)
97 IntParam.push_back(0);
98 else
99 IntParam.push_back(m_exvolmod->ComponentIndices()[i]);
100 if (m_mfmod == NULL)
101 IntParam.push_back(0);
102 else
103 IntParam.push_back(m_mfmod->ComponentIndices()[i]);
104
105 if (MapComp.count(IntParam) == 0) {
106 MapComp[IntParam] = tind;
107 m_MapTodMuStar[i] = tind;
108 m_MapFromdMuStar.push_back(i);
109 m_dMuStarIndices.push_back(vector<int>(1, i));
110 tind++;
111 }
112 else {
113 m_MapTodMuStar[i] = MapComp[IntParam];
114 m_dMuStarIndices[MapComp[IntParam]].push_back(i);
115 }
116 }
117
118 m_ComponentMapCalculated = true;
119 }
120
121 std::vector<double> ThermalModelRealGas::SearchSingleSolution(const std::vector<double> &muStarInit)
122 {
123 return SearchSingleSolutionUsingComponents(muStarInit);
124 //return SearchSingleSolutionDirect(muStarInit);
125 }
126
127 vector<double> ThermalModelRealGas::SearchSingleSolutionUsingComponents(const vector<double> &muStarInit)
128 {
129 if (!m_ComponentMapCalculated)
130 CalculateComponentsMap();
131
132
133 int NN = m_densities.size();
134 int NNdmu = m_MapFromdMuStar.size();
135
136 vector<double> dmuscur(NNdmu, 0.);
137 for (int i = 0; i < NNdmu; ++i)
138 dmuscur[i] = muStarInit[m_MapFromdMuStar[i]] - m_Chem[m_MapFromdMuStar[i]];
139
140 BroydenEquationsRealGasComponents eqs(this);
141 BroydenJacobianRealGasComponents jac(this);
142 Broyden broydn(&eqs, &jac);
143 BroydenSolutionCriteriumRealGas crit(this);
144
145 dmuscur = broydn.Solve(dmuscur, &crit);
146
147 if (broydn.Iterations() == broydn.MaxIterations())
148 m_LastBroydenSuccessFlag = false;
149 else m_LastBroydenSuccessFlag = true;
150
151 m_MaxDiff = broydn.MaxDifference();
152
153 vector<double> ret(NN);
154 for (int i = 0; i < NN; ++i)
155 ret[i] = m_Chem[i] + dmuscur[m_MapTodMuStar[i]];
156
157 return ret;
158 }
159
160 vector<double> ThermalModelRealGas::SearchSingleSolutionDirect(const vector<double> &muStarInit)
161 {
162 int NN = m_densities.size();
163
164 vector<double> dmuscur(NN, 0.);
165 for (int i = 0; i < NN; ++i)
166 dmuscur[i] = muStarInit[i] - m_Chem[i];
167
168 BroydenEquationsRealGas eqs(this);
169 BroydenJacobianRealGas jac(this);
170 Broyden broydn(&eqs, &jac);
171 BroydenSolutionCriteriumRealGas crit(this);
172
173 dmuscur = broydn.Solve(dmuscur, &crit);
174
175 if (broydn.Iterations() == broydn.MaxIterations())
176 m_LastBroydenSuccessFlag = false;
177 else m_LastBroydenSuccessFlag = true;
178
179 m_MaxDiff = broydn.MaxDifference();
180
181 vector<double> ret(NN);
182 for (int i = 0; i < NN; ++i)
183 ret[i] = m_Chem[i] + dmuscur[i];
184
185 return ret;
186 }
187
188 vector<double> ThermalModelRealGas::SearchMultipleSolutions(int iters) {
189 vector<double> csol(m_densities.size(), 0.);
190 double Psol = 0.;
191 bool solved = false;
192 double muBmin = m_Parameters.muB - 0.5 * xMath::mnucleon();
193 double muBmax = m_Parameters.muB + 0.5 * xMath::mnucleon();
194 double dmu = (muBmax - muBmin) / iters;
195 vector<double> curmust(m_densities.size(), 0.);
196 double maxdif = 0.;
197 for (int isol = 0; isol < iters; ++isol) {
198 double tmu = muBmin + (0.5 + isol) * dmu;
199 for (size_t j = 0; j < curmust.size(); ++j) {
200 curmust[j] = m_Chem[j] + (tmu - m_Parameters.muB) * m_Chem[j] / m_Parameters.muB;
201 if (m_TPS->Particles()[j].Statistics() == -1 && curmust[j] > m_TPS->Particles()[j].Mass())
202 curmust[j] = 0.98 * m_TPS->Particles()[j].Mass();
203 }
204
205 vector<double> sol = SearchSingleSolution(curmust);
206
207 bool fl = true;
208 for (size_t i = 0; i < sol.size(); ++i)
209 if (sol[i] != sol[i]) fl = false;
210 fl &= m_LastBroydenSuccessFlag;
211 if (!fl) continue;
212
213 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
214 m_DensitiesId[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, sol[i]);
215
216 m_densities = m_exvolmod->nsol(m_DensitiesId);
217 m_exvolmod->SetDensities(m_densities);
218 m_mfmod->SetDensities(m_densities);
219
220 double tP = 0.;
221 for (size_t i = 0; i < m_densities.size(); ++i) {
222 double tfsum = 0.;
223 for (size_t j = 0; j < m_densities.size(); ++j) {
224 tfsum += m_densities[j] * m_exvolmod->df(i, j);
225 }
226 tfsum = m_exvolmod->f(i) - tfsum;
227 tP += tfsum * m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, sol[i]);
228 }
229
230 tP += -m_mfmod->v();
231 for (size_t i = 0; i < m_densities.size(); ++i)
232 tP += m_mfmod->dv(i) * m_densities[i];
233
234 if (!solved || tP > Psol) {
235 solved = true;
236 Psol = tP;
237 csol = sol;
238 maxdif = m_MaxDiff;
239 }
240 }
241 m_LastBroydenSuccessFlag = solved;
242 m_MaxDiff = maxdif;
243 return csol;
244 }
245
246 void ThermalModelRealGas::SolveEquations() {
247 if (!m_SearchMultipleSolutions) {
248
249 vector<double> muStarInit = m_MuStar;
250
251 for (size_t i = 0; i < muStarInit.size(); ++i) {
252 if (m_TPS->Particles()[i].Statistics() == -1 && muStarInit[i] > m_TPS->Particles()[i].Mass())
253 muStarInit[i] = 0.98 * m_TPS->Particles()[i].Mass();
254 }
255 m_MuStar = SearchSingleSolution(muStarInit);
256 }
257 else {
258 m_MuStar = SearchMultipleSolutions(100);
259 }
260 }
261
262 void ThermalModelRealGas::CalculatePrimordialDensities() {
263 m_FluctuationsCalculated = false;
264
265 int NN = m_densities.size();
266
267 SolveEquations();
268
269 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
270 m_DensitiesId[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, m_MuStar[i]);
271
272 m_densities = m_exvolmod->nsol(m_DensitiesId);
273 m_exvolmod->SetDensities(m_densities);
274 m_mfmod->SetDensities(m_densities);
275
276 // TODO: Scalar density properly
277 m_scaldens = m_densities;
278
279 m_Calculated = true;
280
281 ValidateCalculation();
282 }
283
284 vector<double> ThermalModelRealGas::CalculateChargeFluctuations(const vector<double>& chgs, int order) {
285 vector<double> ret(order + 1, 0.);
286
287 // chi1
288 for (size_t i = 0; i < m_densities.size(); ++i)
289 ret[0] += chgs[i] * m_densities[i];
290
291 ret[0] /= pow(m_Parameters.T * xMath::GeVtoifm(), 3);
292
293 if (order < 2) return ret;
294 // Preparing matrix for system of linear equations
295 int NN = m_densities.size();
296 int Nevcomp = m_exvolmod->ComponentsNumber();
297 const vector<int>& evinds = m_exvolmod->ComponentIndices();
298 const vector<int>& evindsfrom = m_exvolmod->ComponentIndicesFrom();
299
300 int Nmfcomp = m_mfmod->ComponentsNumber();
301 const vector<int>& mfinds = m_mfmod->ComponentIndices();
302 const vector<int>& mfindsfrom = m_mfmod->ComponentIndicesFrom();
303
304 MatrixXd densMatrix(2 * NN, 2 * NN);
305 VectorXd solVector(2 * NN), xVector(2 * NN);
306
307 vector<double> chi2id(m_densities.size()), Ps(m_densities.size());
308 for (int i = 0; i < NN; ++i) {
309 //chi2id[i] = m_TPS->Particles()[i].chi(2, m_Parameters, m_UseWidth, m_MuStar[i]);
310 Ps[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_MuStar[i]);
311 chi2id[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_MuStar[i]) * pow(xMath::GeVtoifm(), 3);
312 }
313
314 vector<double> evc_chi2id(Nevcomp, 0.), evc_Ps(Nevcomp, 0.), evc_ns(Nevcomp, 0.);
315 for (int i = 0; i < NN; ++i) {
316 evc_Ps[evinds[i]] += Ps[i];
317 evc_ns[evinds[i]] += m_DensitiesId[i];
318 evc_chi2id[evinds[i]] += chi2id[i];
319 }
320
321 vector<vector<double>> pkd2fkij(Nevcomp, vector<double>(Nevcomp, 0.));
322 for (int indi = 0; indi < Nevcomp; ++indi) {
323 //int indi = evinds[i];
324 int i = evindsfrom[indi];
325 for (int indj = 0; indj < Nevcomp; ++indj) {
326 //int indj = evinds[j];
327 int j = evindsfrom[indj];
328 //for (int k = 0; k < NN; ++k) {
329 // pkd2fkij[indi][indj] += Ps[k] * m_exvolmod->d2f(k, i, j);
330 //}
331 for (int k = 0; k < Nevcomp; ++k) {
332 pkd2fkij[indi][indj] += evc_Ps[k] * m_exvolmod->d2f(evindsfrom[k], i, j);
333 }
334 }
335 }
336
337 for (int i = 0; i < NN; ++i)
338 for (int j = 0; j < NN; ++j) {
339 densMatrix(i, j) = -m_exvolmod->df(i, j) * m_DensitiesId[i];
340 if (i == j) densMatrix(i, j) += 1.;
341 }
342
343 for (int i = 0; i < NN; ++i)
344 for (int j = 0; j < NN; ++j)
345 densMatrix(i, NN + j) = 0.;
346
347 for (int i = 0; i < NN; ++i) {
348 densMatrix(i, NN + i) = -m_exvolmod->f(i) * chi2id[i];
349 }
350
351 for (int i = 0; i < NN; ++i)
352 for (int j = 0; j < NN; ++j) {
353 densMatrix(NN + i, j) = m_mfmod->d2v(i, j);
354 //for (int k = 0; k < NN; ++k)
355 // densMatrix(NN + i, j) += -m_exvolmod->d2f(k, i, j) * Ps[k];
356 densMatrix(NN + i, j) += -pkd2fkij[evinds[i]][evinds[j]];
357 }
358
359 for (int i = 0; i < NN; ++i)
360 for (int j = 0; j < NN; ++j) {
361 densMatrix(NN + i, NN + j) = -m_exvolmod->df(j, i) * m_DensitiesId[j];
362 if (i == j) densMatrix(NN + i, NN + j) += 1.;
363 }
364
365 PartialPivLU<MatrixXd> decomp(densMatrix);
366
367 // chi2
368 vector<double> dni(NN, 0.), dmus(NN, 0.);
369
370 for (int i = 0; i < NN; ++i) {
371 xVector[i] = 0.;
372 xVector[NN + i] = chgs[i];
373 }
374
375 solVector = decomp.solve(xVector);
376
377 for (int i = 0; i < NN; ++i) {
378 dni[i] = solVector[i];
379 dmus[i] = solVector[NN + i];
380 }
381
382 for (int i = 0; i < NN; ++i)
383 ret[1] += chgs[i] * dni[i];
384
385 ret[1] /= pow(m_Parameters.T, 2) * pow(xMath::GeVtoifm(), 3);
386
387 if (order < 3) return ret;
388
389 vector<double> evc_dn(Nevcomp, 0.), evc_dmus(Nevcomp, 0.), evc_nsdmus(Nevcomp, 0.);
390 for (int i = 0; i < NN; ++i) {
391 evc_dn[evinds[i]] += dni[i];
392 evc_dmus[evinds[i]] += dmus[i];
393 evc_nsdmus[evinds[i]] += m_DensitiesId[i] * dmus[i];
394 }
395
396 vector<double> mfc_dn(Nmfcomp, 0.);
397 for (int i = 0; i < NN; ++i) {
398 mfc_dn[mfinds[i]] += dni[i];
399 }
400
401 // chi3
402 vector<double> d2ni(NN, 0.), d2mus(NN, 0.);
403
404 vector<double> chi3id(m_densities.size());
405 for (int i = 0; i < NN; ++i)
406 chi3id[i] = m_TPS->Particles()[i].chiDimensionfull(3, m_Parameters, m_UseWidth, m_MuStar[i]) * pow(xMath::GeVtoifm(), 3);
407
408
409 vector<vector<double>> d2fijkdnk(Nevcomp, vector<double>(Nevcomp, 0.));
410 for (int indi = 0; indi < Nevcomp; ++indi) {
411 int i = evindsfrom[indi];
412 for (int indj = 0; indj < Nevcomp; ++indj) {
413 int j = evindsfrom[indj];
414 //for (int k = 0; k < NN; ++k) {
415 // d2fijkdnk[indi][indj] += dni[k] * m_exvolmod->d2f(i, j, k);
416 //}
417 for (int indk = 0; indk < Nevcomp; ++indk) {
418 int k = evindsfrom[indk];
419 d2fijkdnk[indi][indj] += evc_dn[indk] * m_exvolmod->d2f(i, j, k);
420 }
421 }
422 }
423
424 vector<double> dfikdnk(Nevcomp, 0.);
425 for (int indi = 0; indi < Nevcomp; ++indi) {
426 int i = evindsfrom[indi];
427 //for (int k = 0; k < NN; ++k) {
428 // dfikdnk[indi] += dni[k] * m_exvolmod->df(i, k);
429 //}
430 for (int indk = 0; indk < Nevcomp; ++indk) {
431 int k = evindsfrom[indk];
432 dfikdnk[indi] += evc_dn[indk] * m_exvolmod->df(i, k);
433 }
434 }
435
436 vector<vector<double>> d2fkijnskmusk(Nevcomp, vector<double>(Nevcomp, 0.));
437 for (int indi = 0; indi < Nevcomp; ++indi) {
438 int i = evindsfrom[indi];
439 for (int indj = 0; indj < Nevcomp; ++indj) {
440 int j = evindsfrom[indj];
441 //for (int k = 0; k < NN; ++k) {
442 // d2fkijnskmusk[indi][indj] += m_exvolmod->d2f(k, i, j) * m_DensitiesId[k] * dmus[k];
443 //}
444 for (int indk = 0; indk < Nevcomp; ++indk) {
445 int k = evindsfrom[indk];
446 d2fkijnskmusk[indi][indj] += m_exvolmod->d2f(k, i, j) * evc_nsdmus[indk];
447 }
448 }
449 }
450
451 vector<vector<double>> pkd3fkijmdnm(Nevcomp, vector<double>(Nevcomp, 0.));
452 for (int indi = 0; indi < Nevcomp; ++indi) {
453 int i = evindsfrom[indi];
454 for (int indj = 0; indj < Nevcomp; ++indj) {
455 int j = evindsfrom[indj];
456 //for (int k = 0; k < NN; ++k) {
457 // for (int m = 0; m < NN; ++m) {
458 // pkd3fkijmdnm[indi][indj] += Ps[k] * m_exvolmod->d3f(k, i, j, m) * dni[m];
459 // }
460 //}
461 for (int indk = 0; indk < Nevcomp; ++indk) {
462 int k = evindsfrom[indk];
463 for (int indm = 0; indm < Nevcomp; ++indm) {
464 int m = evindsfrom[indm];
465 pkd3fkijmdnm[indi][indj] += evc_Ps[indk] * m_exvolmod->d3f(k, i, j, m) * evc_dn[indm];
466 }
467 }
468 }
469 }
470
471 vector<vector<double>> d3vijkdnk(Nmfcomp, vector<double>(Nmfcomp, 0.));
472 for (int indi = 0; indi < Nmfcomp; ++indi) {
473 int i = mfindsfrom[indi];
474 for (int indj = 0; indj < Nmfcomp; ++indj) {
475 int j = mfindsfrom[indj];
476 //for (int k = 0; k < NN; ++k) {
477 // d3vijkdnk[indi][indj] += dni[k] * m_mfmod->d3v(i, j, k);
478 //}
479 for (int indk = 0; indk < Nmfcomp; ++indk) {
480 int k = mfindsfrom[indk];
481 d3vijkdnk[indi][indj] += mfc_dn[indk] * m_mfmod->d3v(i, j, k);
482 }
483 }
484 }
485
486 vector< vector<double> > daij11, daij12, daij21, daij22;
487
488 daij11.resize(NN);
489 daij12.resize(NN);
490 daij21.resize(NN);
491 daij22.resize(NN);
492 for (int i = 0; i < NN; ++i) {
493 //cout << "chi3 iter: " << i << "\n";
494 daij11[i].resize(NN);
495 daij12[i].resize(NN);
496 daij21[i].resize(NN);
497 daij22[i].resize(NN);
498 for (int j = 0; j < NN; ++j) {
499 daij11[i][j] = 0.;
500 //for (int k = 0; k < NN; ++k)
501 // daij11[i][j] += -m_exvolmod->d2f(i, j, k) * dni[k] * m_DensitiesId[i];
502 daij11[i][j] += -d2fijkdnk[evinds[i]][evinds[j]] * m_DensitiesId[i];
503 daij11[i][j] += -m_exvolmod->df(i, j) * chi2id[i] * dmus[i];
504
505 daij12[i][j] = 0.;
506 if (i == j) {
507 //for (int k = 0; k < NN; ++k)
508 // daij12[i][j] += -m_exvolmod->df(i, k) * chi2id[i] * dni[k];
509 daij12[i][j] += -dfikdnk[evinds[i]] * chi2id[i];
510 daij12[i][j] += -m_exvolmod->f(i) * chi3id[i] * dmus[i];
511 }
512
513
514 daij21[i][j] = 0.;
515 daij21[i][j] += d3vijkdnk[mfinds[i]][mfinds[j]];
516 //for (int k = 0; k < NN; ++k) {
517 // daij21[i][j] += m_mfmod->d3v(i, j, k) * dni[k];
518 // //daij21[i][j] += -m_exvolmod->d2f(k, i, j) * m_DensitiesId[k] * dmus[k];
519 // //for (int m = 0; m < NN; ++m)
520 // // daij21[i][j] += -Ps[k] * m_exvolmod->d3f(k, i, j, m) * dni[m];
521 //}
522 daij21[i][j] += -d2fkijnskmusk[evinds[i]][evinds[j]];
523 daij21[i][j] += -pkd3fkijmdnm[evinds[i]][evinds[j]];
524
525 daij22[i][j] = 0.;
526 daij22[i][j] += -m_exvolmod->df(j, i) * chi2id[j] * dmus[j];
527 //for (int k = 0; k < NN; ++k)
528 // daij22[i][j] += -m_exvolmod->d2f(j, i, k) * m_DensitiesId[j] * dni[k];
529 daij22[i][j] += -m_DensitiesId[j] * d2fijkdnk[evinds[j]][evinds[i]];
530 }
531 }
532
533
534 for (int i = 0; i < NN; ++i) {
535 xVector[i] = 0.;
536
537 for (int j = 0; j < NN; ++j)
538 xVector[i] += -daij11[i][j] * dni[j];
539
540 for (int j = 0; j < NN; ++j)
541 xVector[i] += -daij12[i][j] * dmus[j];
542 }
543 for (int i = 0; i < NN; ++i) {
544 xVector[NN + i] = 0.;
545
546 for (int j = 0; j < NN; ++j)
547 xVector[NN + i] += -daij21[i][j] * dni[j];
548
549 for (int j = 0; j < NN; ++j)
550 xVector[NN + i] += -daij22[i][j] * dmus[j];
551 }
552
553 solVector = decomp.solve(xVector);
554
555 for (int i = 0; i < NN; ++i) {
556 d2ni[i] = solVector[i];
557 d2mus[i] = solVector[NN + i];
558 }
559
560 for (int i = 0; i < NN; ++i)
561 ret[2] += chgs[i] * d2ni[i];
562
563 ret[2] /= m_Parameters.T * pow(xMath::GeVtoifm(), 3);
564
565 if (order < 4) return ret;
566
567 vector<double> evc_d2n(Nevcomp, 0.), evc_d2mus(Nevcomp, 0.), evc_nsd2mus(Nevcomp, 0.), evc_chi2iddmus2(Nevcomp, 0.);
568 for (int i = 0; i < NN; ++i) {
569 evc_d2n[evinds[i]] += d2ni[i];
570 evc_d2mus[evinds[i]] += d2mus[i];
571 evc_nsd2mus[evinds[i]] += m_DensitiesId[i] * d2mus[i];
572 evc_chi2iddmus2[evinds[i]] += chi2id[i] * dmus[i] * dmus[i];
573 }
574
575 vector<double> mfc_d2n(Nmfcomp, 0.);
576 for (int i = 0; i < NN; ++i) {
577 mfc_d2n[mfinds[i]] += d2ni[i];
578 }
579
580 // chi4
581 vector<double> d3ni(NN, 0.), d3mus(NN, 0.);
582
583 vector<double> chi4id(m_densities.size());
584 for (int i = 0; i < NN; ++i)
585 chi4id[i] = m_TPS->Particles()[i].chiDimensionfull(4, m_Parameters, m_UseWidth, m_MuStar[i]) * pow(xMath::GeVtoifm(), 3);
586
587 vector<vector<double>> d2fijkd2nk(Nevcomp, vector<double>(Nevcomp, 0.));
588 for (int indi = 0; indi < Nevcomp; ++indi) {
589 int i = evindsfrom[indi];
590 for (int indj = 0; indj < Nevcomp; ++indj) {
591 int j = evindsfrom[indj];
592 //for (int k = 0; k < NN; ++k) {
593 // d2fijkd2nk[indi][indj] += d2ni[k] * m_exvolmod->d2f(i, j, k);
594 //}
595 for (int indk = 0; indk < Nevcomp; ++indk) {
596 int k = evindsfrom[indk];
597 d2fijkd2nk[indi][indj] += evc_d2n[indk] * m_exvolmod->d2f(i, j, k);
598 }
599 }
600 }
601
602 vector<vector<double>> d3fijkmdnkdnm(Nevcomp, vector<double>(Nevcomp, 0.));
603 for (int indi = 0; indi < Nevcomp; ++indi) {
604 int i = evindsfrom[indi];
605 for (int indj = 0; indj < Nevcomp; ++indj) {
606 int j = evindsfrom[indj];
607 //for (int k = 0; k < NN; ++k) {
608 // for (int m = 0; m < NN; ++m) {
609 // d3fijkmdnkdnm[indi][indj] += m_exvolmod->d3f(i, j, k, m) * dni[k] * dni[m];
610 // }
611 //}
612 for (int indk = 0; indk < Nevcomp; ++indk) {
613 int k = evindsfrom[indk];
614 for (int indm = 0; indm < Nevcomp; ++indm) {
615 int m = evindsfrom[indm];
616 d3fijkmdnkdnm[indi][indj] += m_exvolmod->d3f(i, j, k, m) * evc_dn[indk] * evc_dn[indm];
617 }
618 }
619 }
620 }
621
622 vector<double> dfikd2nk(Nevcomp, 0.);
623 for (int indi = 0; indi < Nevcomp; ++indi) {
624 int i = evindsfrom[indi];
625 //for (int k = 0; k < NN; ++k) {
626 // dfikd2nk[indi] += d2ni[k] * m_exvolmod->df(i, k);
627 //}
628 for (int indk = 0; indk < Nevcomp; ++indk) {
629 int k = evindsfrom[indk];
630 dfikd2nk[indi] += evc_d2n[indk] * m_exvolmod->df(i, k);
631 }
632 }
633
634 vector<double> d2fikmdnkdnm(Nevcomp, 0.);
635 for (int indi = 0; indi < Nevcomp; ++indi) {
636 int i = evindsfrom[indi];
637 //for (int k = 0; k < NN; ++k) {
638 // for (int m = 0; m < NN; ++m) {
639 // d2fikmdnkdnm[indi] += m_exvolmod->d2f(i, k, m) * dni[k] * dni[m];
640 // }
641 //}
642 for (int indk = 0; indk < Nevcomp; ++indk) {
643 int k = evindsfrom[indk];
644 for (int indm = 0; indm < Nevcomp; ++indm) {
645 int m = evindsfrom[indm];
646 d2fikmdnkdnm[indi] += m_exvolmod->d2f(i, k, m) * evc_dn[indk] * evc_dn[indm];
647 }
648 }
649 }
650
651 vector<vector<double>> d2fkijnskd2musk(Nevcomp, vector<double>(Nevcomp, 0.));
652 for (int indi = 0; indi < Nevcomp; ++indi) {
653 int i = evindsfrom[indi];
654 for (int indj = 0; indj < Nevcomp; ++indj) {
655 int j = evindsfrom[indj];
656 //for (int k = 0; k < NN; ++k) {
657 // d2fkijnskd2musk[indi][indj] += m_exvolmod->d2f(k, i, j) * m_DensitiesId[k] * d2mus[k];
658 //}
659 for (int indk = 0; indk < Nevcomp; ++indk) {
660 int k = evindsfrom[indk];
661 d2fkijnskd2musk[indi][indj] += m_exvolmod->d2f(k, i, j) * evc_nsd2mus[indk];
662 }
663 }
664 }
665
666 vector<vector<double>> d2fkijc2kdmuskdmusk(Nevcomp, vector<double>(Nevcomp, 0.));
667 for (int indi = 0; indi < Nevcomp; ++indi) {
668 int i = evindsfrom[indi];
669 for (int indj = 0; indj < Nevcomp; ++indj) {
670 int j = evindsfrom[indj];
671 //for (int k = 0; k < NN; ++k) {
672 // d2fkijc2kdmuskdmusk[indi][indj] += m_exvolmod->d2f(k, i, j) * chi2id[k] * dmus[k] * dmus[k];
673 //}
674 for (int indk = 0; indk < Nevcomp; ++indk) {
675 int k = evindsfrom[indk];
676 d2fkijc2kdmuskdmusk[indi][indj] += m_exvolmod->d2f(k, i, j) * evc_chi2iddmus2[indk];
677 }
678 }
679 }
680
681 vector<vector<double>> nskd3fkijmdmuskdnm(Nevcomp, vector<double>(Nevcomp, 0.));
682 for (int indi = 0; indi < Nevcomp; ++indi) {
683 int i = evindsfrom[indi];
684 for (int indj = 0; indj < Nevcomp; ++indj) {
685 int j = evindsfrom[indj];
686 //for (int k = 0; k < NN; ++k) {
687 // for (int m = 0; m < NN; ++m) {
688 // nskd3fkijmdmuskdnm[indi][indj] += m_DensitiesId[k] * m_exvolmod->d3f(k, i, j, m) * dmus[k] * dni[m];
689 // }
690 //}
691 for (int indk = 0; indk < Nevcomp; ++indk) {
692 int k = evindsfrom[indk];
693 for (int indm = 0; indm < Nevcomp; ++indm) {
694 int m = evindsfrom[indm];
695 nskd3fkijmdmuskdnm[indi][indj] += evc_nsdmus[indk] * m_exvolmod->d3f(k, i, j, m) * evc_dn[indm];
696 }
697 }
698 }
699 }
700
701 vector<vector<double>> pkd3fkijmd2nm(Nevcomp, vector<double>(Nevcomp, 0.));
702 for (int indi = 0; indi < Nevcomp; ++indi) {
703 int i = evindsfrom[indi];
704 for (int indj = 0; indj < Nevcomp; ++indj) {
705 int j = evindsfrom[indj];
706 //for (int k = 0; k < NN; ++k) {
707 // for (int m = 0; m < NN; ++m) {
708 // pkd3fkijmd2nm[indi][indj] += Ps[k] * m_exvolmod->d3f(k, i, j, m) * d2ni[m];
709 // }
710 //}
711 for (int indk = 0; indk < Nevcomp; ++indk) {
712 int k = evindsfrom[indk];
713 for (int indm = 0; indm < Nevcomp; ++indm) {
714 int m = evindsfrom[indm];
715 pkd3fkijmd2nm[indi][indj] += evc_Ps[indk] * m_exvolmod->d3f(k, i, j, m) * evc_d2n[indm];
716 }
717 }
718 }
719 }
720
721 vector<vector<double>> pkd4fkijmldnmdnl(Nevcomp, vector<double>(Nevcomp, 0.));
722 for (int indi = 0; indi < Nevcomp; ++indi) {
723 int i = evindsfrom[indi];
724 for (int indj = 0; indj < Nevcomp; ++indj) {
725 int j = evindsfrom[indj];
726 //for (int k = 0; k < NN; ++k) {
727 // for (int m = 0; m < NN; ++m) {
728 // for (int l = 0; m < NN; ++m) {
729 // pkd4fkijmldnmdnl[indi][indj] += Ps[k] * m_exvolmod->d4f(k, i, j, m, l) * dni[m] * dni[l];
730 // }
731 // }
732 //}
733 for (int indk = 0; indk < Nevcomp; ++indk) {
734 int k = evindsfrom[indk];
735 for (int indm = 0; indm < Nevcomp; ++indm) {
736 int m = evindsfrom[indm];
737 for (int indl = 0; indl < Nevcomp; ++indl) {
738 int l = evindsfrom[indl];
739 pkd4fkijmldnmdnl[indi][indj] += evc_Ps[indk] * m_exvolmod->d4f(k, i, j, m, l) * evc_dn[indm] * evc_dn[indl];
740 }
741 }
742 }
743 }
744 }
745
746 vector<vector<double>> d3vijkd2nk(Nmfcomp, vector<double>(Nmfcomp, 0.));
747 for (int indi = 0; indi < Nmfcomp; ++indi) {
748 int i = mfindsfrom[indi];
749 for (int indj = 0; indj < Nmfcomp; ++indj) {
750 int j = mfindsfrom[indj];
751 //for (int k = 0; k < NN; ++k) {
752 // d3vijkd2nk[indi][indj] += d2ni[k] * m_mfmod->d3v(i, j, k);
753 //}
754 for (int indk = 0; indk < Nmfcomp; ++indk) {
755 int k = mfindsfrom[indk];
756 d3vijkd2nk[indi][indj] += mfc_d2n[indk] * m_mfmod->d3v(i, j, k);
757 }
758 }
759 }
760
761 vector<vector<double>> d4vijkmdnkdnm(Nmfcomp, vector<double>(Nmfcomp, 0.));
762 for (int indi = 0; indi < Nmfcomp; ++indi) {
763 int i = mfindsfrom[indi];
764 for (int indj = 0; indj < Nmfcomp; ++indj) {
765 int j = mfindsfrom[indj];
766 //for (int k = 0; k < NN; ++k) {
767 // for (int m = 0; m < NN; ++m) {
768 // d4vijkmdnkdnm[indi][indj] += dni[k] * dni[m] * m_mfmod->d4v(i, j, k, m);
769 // }
770 //}
771 for (int indk = 0; indk < Nmfcomp; ++indk) {
772 int k = mfindsfrom[indk];
773 for (int indm = 0; indm < Nmfcomp; ++indm) {
774 int m = mfindsfrom[indm];
775 d4vijkmdnkdnm[indi][indj] += mfc_dn[indk] * mfc_dn[indm] * m_mfmod->d4v(i, j, k, m);
776 }
777 }
778 }
779 }
780
781 vector< vector<double> > d2aij11, d2aij12, d2aij21, d2aij22;
782
783 d2aij11.resize(NN);
784 d2aij12.resize(NN);
785 d2aij21.resize(NN);
786 d2aij22.resize(NN);
787 for (int i = 0; i < NN; ++i) {
788 //cout << "chi4 iter: " << i << "\n";
789 d2aij11[i].resize(NN);
790 d2aij12[i].resize(NN);
791 d2aij21[i].resize(NN);
792 d2aij22[i].resize(NN);
793 for (int j = 0; j < NN; ++j) {
794 d2aij11[i][j] = 0.;
795 d2aij11[i][j] += -m_exvolmod->df(i, j) * chi3id[i] * dmus[i] * dmus[i];
796 d2aij11[i][j] += -m_exvolmod->df(i, j) * chi2id[i] * d2mus[i];
797
798 d2aij11[i][j] += -2. * d2fijkdnk[evinds[i]][evinds[j]] * chi2id[i] * dmus[i];
799 d2aij11[i][j] += -d2fijkd2nk[evinds[i]][evinds[j]] * m_DensitiesId[i];
800 d2aij11[i][j] += -d3fijkmdnkdnm[evinds[i]][evinds[j]] * m_DensitiesId[i];
801
802 //for (int k = 0; k < NN; ++k) {
803 // //d2aij11[i][j] += -m_exvolmod->d2f(i, j, k) * chi2id[i] * dmus[i] * dni[k];
804 // //d2aij11[i][j] += -m_exvolmod->d2f(i, j, k) * m_DensitiesId[i] * d2ni[k];
805 // //d2aij11[i][j] += -m_exvolmod->d2f(i, j, k) * chi2id[i] * dmus[i] * dni[k];
806 // //for (int m = 0; m < NN; ++m) {
807 // // d2aij11[i][j] += -m_exvolmod->d3f(i, j, k, m) * m_DensitiesId[i] * dni[k] * dni[m];
808 // //}
809 //}
810
811 d2aij12[i][j] = 0.;
812 if (i == j) {
813 d2aij12[i][j] += -m_exvolmod->f(i) * chi3id[i] * d2mus[i];
814 d2aij12[i][j] += -m_exvolmod->f(i) * chi4id[i] * dmus[i] * dmus[i];
815
816 d2aij12[i][j] += -2. * dfikdnk[evinds[i]] * chi3id[i] * dmus[i];
817 d2aij12[i][j] += -dfikd2nk[evinds[i]] * chi2id[i];
818 d2aij12[i][j] += -d2fikmdnkdnm[evinds[i]] * chi2id[i];
819
820 //for (int k = 0; k < NN; ++k) {
821 // d2aij12[i][j] += -2. * m_exvolmod->df(i, k) * chi3id[i] * dmus[i] * dni[k];
822 // d2aij12[i][j] += -m_exvolmod->df(i, k) * chi2id[i] * d2ni[k];
823
824 // //for (int m = 0; m < NN; ++m) {
825 // // d2aij12[i][j] += -m_exvolmod->d2f(i, k, m) * chi2id[i] * dni[k] * dni[m];
826 // //}
827 //}
828 }
829
830 d2aij21[i][j] = 0.;
831
832 d2aij21[i][j] += -d2fkijnskd2musk[evinds[i]][evinds[j]];
833 d2aij21[i][j] += -d2fkijc2kdmuskdmusk[evinds[i]][evinds[j]];
834 d2aij21[i][j] += -2. * nskd3fkijmdmuskdnm[evinds[i]][evinds[j]];
835 d2aij21[i][j] += -pkd3fkijmd2nm[evinds[i]][evinds[j]];
836 d2aij21[i][j] += -pkd4fkijmldnmdnl[evinds[i]][evinds[j]];
837
838 d2aij21[i][j] += d3vijkd2nk[mfinds[i]][mfinds[j]];
839 d2aij21[i][j] += d4vijkmdnkdnm[mfinds[i]][mfinds[j]];
840
841 //for (int k = 0; k < NN; ++k) {
842 // ////daij21[i][j] += m_mfmod->d3v(i, j, k) * dni[k];
843 // d2aij21[i][j] += m_mfmod->d3v(i, j, k) * d2ni[k];
844 // for (int m = 0; m < NN; ++m)
845 // d2aij21[i][j] += m_mfmod->d4v(i, j, k, m) * dni[k] * dni[m];
846 // ////daij21[i][j] += -m_exvolmod->d2f(k, i, j) * m_DensitiesId[k] * dmus[k];
847 // //d2aij21[i][j] += -m_exvolmod->d2f(k, i, j) * m_DensitiesId[k] * d2mus[k];
848 // //d2aij21[i][j] += -m_exvolmod->d2f(k, i, j) * chi2id[k] * dmus[k] * dmus[k];
849 // //for (int m = 0; m < NN; ++m)
850 // // d2aij21[i][j] += -m_exvolmod->d3f(k, i, j, m) * m_DensitiesId[k] * dmus[k] * dni[m];
851
852 // //for (int m = 0; m < NN; ++m) {
853 // // ////daij21[i][j] += -Ps[k] * m_exvolmod->d3f(k, i, j, m) * dni[m];
854 // // //d2aij21[i][j] += -m_DensitiesId[k] * m_exvolmod->d3f(k, i, j, m) * dni[m] * dmus[k];
855 // // //d2aij21[i][j] += -Ps[k] * m_exvolmod->d3f(k, i, j, m) * d2ni[m];
856 // // //for (int l = 0; l < NN; ++l)
857 // // // d2aij21[i][j] += -Ps[k] * m_exvolmod->d4f(k, i, j, m, l) * dni[m] * dni[l];
858 // //}
859 //}
860
861 d2aij22[i][j] = 0.;
862 //daij22[i][j] += -m_exvolmod->df(j, i) * chi2id[j] * dmus[j];
863 d2aij22[i][j] += -m_exvolmod->df(j, i) * chi3id[j] * dmus[j] * dmus[j];
864 d2aij22[i][j] += -m_exvolmod->df(j, i) * chi2id[j] * d2mus[j];
865
866 d2aij22[i][j] += -2. * d2fijkdnk[evinds[j]][evinds[i]] * chi2id[j] * dmus[j];
867
868 //for (int k = 0; k < NN; ++k)
869 // d2aij22[i][j] += -m_exvolmod->d2f(j, i, k) * chi2id[j] * dmus[j] * dni[k];
870
871 d2aij22[i][j] += -d2fijkd2nk[evinds[j]][evinds[i]] * m_DensitiesId[j];
872 d2aij22[i][j] += -d3fijkmdnkdnm[evinds[j]][evinds[i]] * m_DensitiesId[j];
873
874 //for (int k = 0; k < NN; ++k) {
875 // ////daij22[i][j] += -m_exvolmod->d2f(j, i, k) * m_DensitiesId[j] * dni[k];
876 // //d2aij22[i][j] += -m_exvolmod->d2f(j, i, k) * chi2id[j] * dni[k] * dmus[j];
877 // //d2aij22[i][j] += -m_exvolmod->d2f(j, i, k) * m_DensitiesId[j] * d2ni[k];
878 // //for (int m = 0; m < NN; ++m) {
879 // // d2aij22[i][j] += -m_exvolmod->d3f(j, i, k, m) * m_DensitiesId[j] * dni[k] * dni[m];
880 // //}
881 //}
882 }
883 }
884
885
886 for (int i = 0; i < NN; ++i) {
887 xVector[i] = 0.;
888
889 for (int j = 0; j < NN; ++j)
890 xVector[i] += -2. * daij11[i][j] * d2ni[j] - d2aij11[i][j] * dni[j];
891
892 for (int j = 0; j < NN; ++j)
893 xVector[i] += -2. * daij12[i][j] * d2mus[j] - d2aij12[i][j] * dmus[j];
894 }
895 for (int i = 0; i < NN; ++i) {
896 xVector[NN + i] = 0.;
897
898 for (int j = 0; j < NN; ++j)
899 xVector[NN + i] += -2. * daij21[i][j] * d2ni[j] - d2aij21[i][j] * dni[j];
900
901 for (int j = 0; j < NN; ++j)
902 xVector[NN + i] += -2. * daij22[i][j] * d2mus[j] - d2aij22[i][j] * dmus[j];
903 }
904
905 solVector = decomp.solve(xVector);
906
907 for (int i = 0; i < NN; ++i) {
908 d3ni[i] = solVector[i];
909 d3mus[i] = solVector[NN + i];
910 }
911
912 for (int i = 0; i < NN; ++i)
913 ret[3] += chgs[i] * d3ni[i];
914
915 ret[3] /= pow(xMath::GeVtoifm(), 3);
916
917 return ret;
918 }
919
920 vector<double> ThermalModelRealGas::CalculateChargeFluctuationsOld(const vector<double>& chgs, int order) {
921 vector<double> ret(order + 1, 0.);
922
923 // chi1
924 for (size_t i = 0; i < m_densities.size(); ++i)
925 ret[0] += chgs[i] * m_densities[i];
926
927 ret[0] /= pow(m_Parameters.T * xMath::GeVtoifm(), 3);
928
929 if (order < 2) return ret;
930 // Preparing matrix for system of linear equations
931 int NN = m_densities.size();
932
933 MatrixXd densMatrix(2 * NN, 2 * NN);
934 VectorXd solVector(2 * NN), xVector(2 * NN);
935
936 vector<double> chi2id(m_densities.size()), Ps(m_densities.size());
937 for (int i = 0; i < NN; ++i) {
938 //chi2id[i] = m_TPS->Particles()[i].chi(2, m_Parameters, m_UseWidth, m_MuStar[i]);
939 Ps[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_MuStar[i]);
940 chi2id[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_MuStar[i]) * pow(xMath::GeVtoifm(), 3);
941 }
942
943 for (int i = 0; i < NN; ++i)
944 for (int j = 0; j < NN; ++j) {
945 densMatrix(i, j) = -m_exvolmod->df(i, j) * m_DensitiesId[i];
946 if (i == j) densMatrix(i, j) += 1.;
947 }
948
949 for (int i = 0; i < NN; ++i)
950 for (int j = 0; j < NN; ++j)
951 densMatrix(i, NN + j) = 0.;
952
953 for (int i = 0; i < NN; ++i) {
954 densMatrix(i, NN + i) = -m_exvolmod->f(i) * chi2id[i];
955 }
956
957 for (int i = 0; i < NN; ++i)
958 for (int j = 0; j < NN; ++j) {
959 densMatrix(NN + i, j) = m_mfmod->d2v(i, j);
960 for (int k = 0; k < NN; ++k) {
961 densMatrix(NN + i, j) += -m_exvolmod->d2f(k, i, j) * Ps[k];
962 }
963 }
964
965 for (int i = 0; i < NN; ++i)
966 for (int j = 0; j < NN; ++j) {
967 densMatrix(NN + i, NN + j) = -m_exvolmod->df(j, i) * m_DensitiesId[j];
968 if (i == j) densMatrix(NN + i, NN + j) += 1.;
969 }
970
971 PartialPivLU<MatrixXd> decomp(densMatrix);
972
973 // chi2
974 vector<double> dni(NN, 0.), dmus(NN, 0.);
975
976 for (int i = 0; i < NN; ++i) {
977 xVector[i] = 0.;
978 xVector[NN + i] = chgs[i];
979 }
980
981 solVector = decomp.solve(xVector);
982
983 for (int i = 0; i < NN; ++i) {
984 dni[i] = solVector[i];
985 dmus[i] = solVector[NN + i];
986 }
987
988 for (int i = 0; i < NN; ++i)
989 ret[1] += chgs[i] * dni[i];
990
991 ret[1] /= pow(m_Parameters.T, 2) * pow(xMath::GeVtoifm(), 3);
992
993 if (order < 3) return ret;
994
995 // chi3
996 vector<double> d2ni(NN, 0.), d2mus(NN, 0.);
997
998 vector<double> chi3id(m_densities.size());
999 for (int i = 0; i < NN; ++i)
1000 chi3id[i] = m_TPS->Particles()[i].chiDimensionfull(3, m_Parameters, m_UseWidth, m_MuStar[i]) * pow(xMath::GeVtoifm(), 3);
1001
1002 vector< vector<double> > daij11, daij12, daij21, daij22;
1003
1004 daij11.resize(NN);
1005 daij12.resize(NN);
1006 daij21.resize(NN);
1007 daij22.resize(NN);
1008 for (int i = 0; i < NN; ++i) {
1009 cout << "chi3 iter: " << i << "\n";
1010 daij11[i].resize(NN);
1011 daij12[i].resize(NN);
1012 daij21[i].resize(NN);
1013 daij22[i].resize(NN);
1014 for (int j = 0; j < NN; ++j) {
1015 daij11[i][j] = 0.;
1016 for (int k = 0; k < NN; ++k)
1017 daij11[i][j] += -m_exvolmod->d2f(i, j, k) * dni[k] * m_DensitiesId[i];
1018 daij11[i][j] += -m_exvolmod->df(i, j) * chi2id[i] * dmus[i];
1019
1020 daij12[i][j] = 0.;
1021 if (i == j) {
1022 for (int k = 0; k < NN; ++k)
1023 daij12[i][j] += -m_exvolmod->df(i, k) * chi2id[i] * dni[k];
1024 daij12[i][j] += -m_exvolmod->f(i) * chi3id[i] * dmus[i];
1025 }
1026
1027
1028 daij21[i][j] = 0.;
1029 for (int k = 0; k < NN; ++k) {
1030 daij21[i][j] += m_mfmod->d3v(i, j, k) * dni[k];
1031 daij21[i][j] += -m_exvolmod->d2f(k, i, j) * m_DensitiesId[k] * dmus[k];
1032 for (int m = 0; m < NN; ++m)
1033 daij21[i][j] += -Ps[k] * m_exvolmod->d3f(k, i, j, m) * dni[m];
1034 }
1035
1036 daij22[i][j] = 0.;
1037 daij22[i][j] += -m_exvolmod->df(j, i) * chi2id[j] * dmus[j];
1038 for (int k = 0; k < NN; ++k)
1039 daij22[i][j] += -m_exvolmod->d2f(j, i, k) * m_DensitiesId[j] * dni[k];
1040 }
1041 }
1042
1043
1044 for (int i = 0; i < NN; ++i) {
1045 xVector[i] = 0.;
1046
1047 for (int j = 0; j < NN; ++j)
1048 xVector[i] += -daij11[i][j] * dni[j];
1049
1050 for (int j = 0; j < NN; ++j)
1051 xVector[i] += -daij12[i][j] * dmus[j];
1052 }
1053 for (int i = 0; i < NN; ++i) {
1054 xVector[NN + i] = 0.;
1055
1056 for (int j = 0; j < NN; ++j)
1057 xVector[NN + i] += -daij21[i][j] * dni[j];
1058
1059 for (int j = 0; j < NN; ++j)
1060 xVector[NN + i] += -daij22[i][j] * dmus[j];
1061 }
1062
1063 solVector = decomp.solve(xVector);
1064
1065 for (int i = 0; i < NN; ++i) {
1066 d2ni[i] = solVector[i];
1067 d2mus[i] = solVector[NN + i];
1068 }
1069
1070 for (int i = 0; i < NN; ++i)
1071 ret[2] += chgs[i] * d2ni[i];
1072
1073 ret[2] /= m_Parameters.T * pow(xMath::GeVtoifm(), 3);
1074
1075 if (order < 4) return ret;
1076
1077 // chi4
1078 vector<double> d3ni(NN, 0.), d3mus(NN, 0.);
1079
1080 vector<double> chi4id(m_densities.size());
1081 for (int i = 0; i < NN; ++i)
1082 chi4id[i] = m_TPS->Particles()[i].chiDimensionfull(4, m_Parameters, m_UseWidth, m_MuStar[i]) * pow(xMath::GeVtoifm(), 3);
1083
1084 vector< vector<double> > d2aij11, d2aij12, d2aij21, d2aij22;
1085
1086 d2aij11.resize(NN);
1087 d2aij12.resize(NN);
1088 d2aij21.resize(NN);
1089 d2aij22.resize(NN);
1090 for (int i = 0; i < NN; ++i) {
1091 cout << "chi4 iter: " << i << "\n";
1092 d2aij11[i].resize(NN);
1093 d2aij12[i].resize(NN);
1094 d2aij21[i].resize(NN);
1095 d2aij22[i].resize(NN);
1096 for (int j = 0; j < NN; ++j) {
1097 d2aij11[i][j] = 0.;
1098 d2aij11[i][j] += -m_exvolmod->df(i, j) * chi3id[i] * dmus[i] * dmus[i];
1099 d2aij11[i][j] += -m_exvolmod->df(i, j) * chi2id[i] * d2mus[i];
1100 for (int k = 0; k < NN; ++k) {
1101 d2aij11[i][j] += -m_exvolmod->d2f(i, j, k) * chi2id[i] * dmus[i] * dni[k];
1102 d2aij11[i][j] += -m_exvolmod->d2f(i, j, k) * m_DensitiesId[i] * d2ni[k];
1103 d2aij11[i][j] += -m_exvolmod->d2f(i, j, k) * chi2id[i] * dmus[i] * dni[k];
1104 for (int m = 0; m < NN; ++m) {
1105 d2aij11[i][j] += -m_exvolmod->d3f(i, j, k, m) * m_DensitiesId[i] * dni[k] * dni[m];
1106 }
1107 }
1108
1109 d2aij12[i][j] = 0.;
1110 if (i == j) {
1111 d2aij12[i][j] += -m_exvolmod->f(i) * chi3id[i] * d2mus[i];
1112 d2aij12[i][j] += -m_exvolmod->f(i) * chi4id[i] * dmus[i] * dmus[i];
1113 for (int k = 0; k < NN; ++k) {
1114 d2aij12[i][j] += -2. * m_exvolmod->df(i, k) * chi3id[i] * dmus[i] * dni[k];
1115 d2aij12[i][j] += -m_exvolmod->df(i, k) * chi2id[i] * d2ni[k];
1116 for (int m = 0; m < NN; ++m) {
1117 d2aij12[i][j] += -m_exvolmod->d2f(i, k, m) * chi2id[i] * dni[k] * dni[m];
1118 }
1119 }
1120 }
1121
1122 d2aij21[i][j] = 0.;
1123 for (int k = 0; k < NN; ++k) {
1124 //daij21[i][j] += m_mfmod->d3v(i, j, k) * dni[k];
1125 d2aij21[i][j] += m_mfmod->d3v(i, j, k) * d2ni[k];
1126 for (int m = 0; m < NN; ++m)
1127 d2aij21[i][j] += m_mfmod->d4v(i, j, k, m) * dni[k] * dni[m];
1128 //daij21[i][j] += -m_exvolmod->d2f(k, i, j) * m_DensitiesId[k] * dmus[k];
1129 d2aij21[i][j] += -m_exvolmod->d2f(k, i, j) * m_DensitiesId[k] * d2mus[k];
1130 d2aij21[i][j] += -m_exvolmod->d2f(k, i, j) * chi2id[k] * dmus[k] * dmus[k];
1131 for (int m = 0; m < NN; ++m)
1132 d2aij21[i][j] += -m_exvolmod->d3f(k, i, j, m) * m_DensitiesId[k] * dmus[k] * dni[m];
1133
1134 for (int m = 0; m < NN; ++m) {
1135 //daij21[i][j] += -Ps[k] * m_exvolmod->d3f(k, i, j, m) * dni[m];
1136 d2aij21[i][j] += -m_DensitiesId[k] * m_exvolmod->d3f(k, i, j, m) * dni[m] * dmus[k];
1137 d2aij21[i][j] += -Ps[k] * m_exvolmod->d3f(k, i, j, m) * d2ni[m];
1138 for (int l = 0; l < NN; ++l)
1139 d2aij21[i][j] += -Ps[k] * m_exvolmod->d4f(k, i, j, m, l) * dni[m] * dni[l];
1140 }
1141 }
1142
1143 d2aij22[i][j] = 0.;
1144 //daij22[i][j] += -m_exvolmod->df(j, i) * chi2id[j] * dmus[j];
1145 d2aij22[i][j] += -m_exvolmod->df(j, i) * chi3id[j] * dmus[j] * dmus[j];
1146 d2aij22[i][j] += -m_exvolmod->df(j, i) * chi2id[j] * d2mus[j];
1147 for (int k = 0; k < NN; ++k)
1148 d2aij22[i][j] += -m_exvolmod->d2f(j, i, k) * chi2id[j] * dmus[j] * dni[k];
1149 for (int k = 0; k < NN; ++k) {
1150 //daij22[i][j] += -m_exvolmod->d2f(j, i, k) * m_DensitiesId[j] * dni[k];
1151 d2aij22[i][j] += -m_exvolmod->d2f(j, i, k) * chi2id[j] * dni[k] * dmus[j];
1152 d2aij22[i][j] += -m_exvolmod->d2f(j, i, k) * m_DensitiesId[j] * d2ni[k];
1153 for (int m = 0; m < NN; ++m) {
1154 d2aij22[i][j] += -m_exvolmod->d3f(j, i, k, m) * m_DensitiesId[j] * dni[k] * dni[m];
1155 }
1156 }
1157 }
1158 }
1159
1160
1161 for (int i = 0; i < NN; ++i) {
1162 xVector[i] = 0.;
1163
1164 for (int j = 0; j < NN; ++j)
1165 xVector[i] += -2. * daij11[i][j] * d2ni[j] - d2aij11[i][j] * dni[j];
1166
1167 for (int j = 0; j < NN; ++j)
1168 xVector[i] += -2. * daij12[i][j] * d2mus[j] - d2aij12[i][j] * dmus[j];
1169 }
1170 for (int i = 0; i < NN; ++i) {
1171 xVector[NN + i] = 0.;
1172
1173 for (int j = 0; j < NN; ++j)
1174 xVector[NN + i] += -2. * daij21[i][j] * d2ni[j] - d2aij21[i][j] * dni[j];
1175
1176 for (int j = 0; j < NN; ++j)
1177 xVector[NN + i] += -2. * daij22[i][j] * d2mus[j] - d2aij22[i][j] * dmus[j];
1178 }
1179
1180 solVector = decomp.solve(xVector);
1181
1182 for (int i = 0; i < NN; ++i) {
1183 d3ni[i] = solVector[i];
1184 d3mus[i] = solVector[NN + i];
1185 }
1186
1187 for (int i = 0; i < NN; ++i)
1188 ret[3] += chgs[i] * d3ni[i];
1189
1190 ret[3] /= pow(xMath::GeVtoifm(), 3);
1191
1192 return ret;
1193 }
1194
1195
1196 // TODO include correlations
1197 vector< vector<double> > ThermalModelRealGas::CalculateFluctuations(int order) {
1198 if (order < 1) return m_chi;
1199
1200 vector<double> chgs(m_densities.size());
1201 vector<double> chis;
1202
1203 // Baryon charge
1204 for (size_t i = 0; i < chgs.size(); ++i)
1205 chgs[i] = m_TPS->Particles()[i].BaryonCharge();
1206 chis = CalculateChargeFluctuations(chgs, order);
1207 for (int i = 0; i < order; ++i) m_chi[i][0] = chis[i];
1208
1209 // Electric charge
1210 for (size_t i = 0; i < chgs.size(); ++i)
1211 chgs[i] = m_TPS->Particles()[i].ElectricCharge();
1212 chis = CalculateChargeFluctuations(chgs, order);
1213 for (int i = 0; i < order; ++i) m_chi[i][1] = chis[i];
1214
1215 // Strangeness charge
1216 for (size_t i = 0; i < chgs.size(); ++i)
1217 chgs[i] = m_TPS->Particles()[i].Strangeness();
1218 chis = CalculateChargeFluctuations(chgs, order);
1219 for (int i = 0; i < order; ++i) m_chi[i][2] = chis[i];
1220
1221 // Arbitrary charge
1222 for (size_t i = 0; i < chgs.size(); ++i)
1223 chgs[i] = m_TPS->Particles()[i].ArbitraryCharge();
1224 chis = CalculateChargeFluctuations(chgs, order);
1225 for (int i = 0; i < order; ++i) m_chiarb[i] = chis[i];
1226
1227 return m_chi;
1228 }
1229
1230 void ThermalModelRealGas::CalculateTwoParticleCorrelations()
1231 {
1232 int NN = m_densities.size();
1233 int Nevcomp = m_exvolmod->ComponentsNumber();
1234 const vector<int>& evinds = m_exvolmod->ComponentIndices();
1235 const vector<int>& evindsfrom = m_exvolmod->ComponentIndicesFrom();
1236
1237
1238 m_PrimCorrel.resize(NN);
1239 for (int i = 0; i < NN; ++i)
1240 m_PrimCorrel[i].resize(NN);
1241 m_dmusdmu = m_PrimCorrel;
1242 m_TotalCorrel = m_PrimCorrel;
1243
1244 MatrixXd densMatrix(2 * NN, 2 * NN);
1245 VectorXd solVector(2 * NN), xVector(2 * NN);
1246
1247 vector<double> chi2id(m_densities.size()), Ps(m_densities.size());
1248 for (int i = 0; i < NN; ++i) {
1249 //chi2id[i] = m_TPS->Particles()[i].chi(2, m_Parameters, m_UseWidth, m_MuStar[i]);
1250 Ps[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_MuStar[i]);
1251 chi2id[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_MuStar[i]);
1252 }
1253
1254 vector<double> evc_Ps(Nevcomp, 0.);
1255 for (int i = 0; i < NN; ++i) {
1256 evc_Ps[evinds[i]] += Ps[i];
1257 }
1258
1259 for (int i = 0; i < NN; ++i)
1260 for (int j = 0; j < NN; ++j) {
1261 densMatrix(i, j) = -m_exvolmod->df(i, j) * m_DensitiesId[i];
1262 if (i == j) densMatrix(i, j) += 1.;
1263 }
1264
1265 for (int i = 0; i < NN; ++i)
1266 for (int j = 0; j < NN; ++j)
1267 densMatrix(i, NN + j) = 0.;
1268
1269 for (int i = 0; i < NN; ++i) {
1270 densMatrix(i, NN + i) = -m_exvolmod->f(i) * chi2id[i] * pow(xMath::GeVtoifm(), 3);
1271 }
1272
1273 for (int i = 0; i < NN; ++i)
1274 for (int j = 0; j < NN; ++j) {
1275 densMatrix(NN + i, j) = m_mfmod->d2v(i, j);
1276 //for (int k = 0; k < NN; ++k) {
1277 // densMatrix(NN + i, j) += -m_exvolmod->d2f(k, i, j) * Ps[k];
1278 //}
1279 for (int indk = 0; indk < Nevcomp; ++indk) {
1280 int k = evindsfrom[indk];
1281 densMatrix(NN + i, j) += -m_exvolmod->d2f(k, i, j) * evc_Ps[indk];
1282 }
1283 }
1284
1285 for (int i = 0; i < NN; ++i)
1286 for (int j = 0; j < NN; ++j) {
1287 densMatrix(NN + i, NN + j) = -m_exvolmod->df(j,i) * m_DensitiesId[j];
1288 if (i == j) densMatrix(NN + i, NN + j) += 1.;
1289 }
1290
1291 PartialPivLU<MatrixXd> decomp(densMatrix);
1292
1293 for (int k = 0; k < NN; ++k) {
1294 vector<double> dni(NN, 0.), dmus(NN, 0.);
1295
1296 for (int i = 0; i < NN; ++i) {
1297 xVector[i] = 0.;
1298 xVector[NN + i] = static_cast<int>(i == k);
1299 }
1300
1301 solVector = decomp.solve(xVector);
1302
1303 for (int i = 0; i < NN; ++i) {
1304 dni[i] = solVector[i];
1305 dmus[i] = solVector[NN + i];
1306 m_dmusdmu[i][k] = dmus[i];
1307 }
1308
1309 for (int j = 0; j < NN; ++j) {
1310 m_PrimCorrel[j][k] = dni[j];
1311 }
1312 }
1313
1314 for (int i = 0; i < NN; ++i) {
1315 m_wprim[i] = m_PrimCorrel[i][i];
1316 if (m_densities[i] > 0.) m_wprim[i] *= m_Parameters.T / m_densities[i];
1317 else m_wprim[i] = 1.;
1318 }
1319
1320 }
1321
1322 void ThermalModelRealGas::CalculateTemperatureDerivatives() {
1323 int N = m_TPS->ComponentsNumber();
1324 int Nevcomp = m_exvolmod->ComponentsNumber();
1325 const vector<int>& evinds = m_exvolmod->ComponentIndices();
1326 const vector<int>& evindsfrom = m_exvolmod->ComponentIndicesFrom();
1327 int Nmfcomp = m_mfmod->ComponentsNumber();
1328 const vector<int>& mfinds = m_mfmod->ComponentIndices();
1329 const vector<int>& mfindsfrom = m_mfmod->ComponentIndicesFrom();
1330
1331 double T = m_Parameters.T;
1332
1333 m_dndT = vector<double>(N, 0.);
1334 m_dmusdT = vector<double>(N, 0.);
1335 m_PrimChi2sdT = vector<vector<double>>(N, vector<double>(N, 0.));
1336
1337 vector<double> nids(N, 0.), dniddTs(N, 0.), chi2ids(N, 0.), dchi2idsdT(N, 0.), chi3ids(N, 0.);
1338 for (int i = 0; i < N; ++i) {
1339 nids[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, m_MuStar[i]);
1340 dniddTs[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::dndT, m_UseWidth, m_MuStar[i]);
1341 chi2ids[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_MuStar[i]) * xMath::GeVtoifm3();
1342 // chi2til = chi2 * T^2
1343 // dchi2til/dT = T^2 * dchi2/dT + 2 * T * chi2 = T^2 * dchi2/dT + 2 * chi2til / T
1344 dchi2idsdT[i] = (T * T * m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::dchi2dT, m_UseWidth, m_MuStar[i]) * xMath::GeVtoifm3() + 2. * T * chi2ids[i] / T / T);
1345 chi3ids[i] = m_TPS->Particles()[i].chiDimensionfull(3, m_Parameters, m_UseWidth, m_MuStar[i]) * xMath::GeVtoifm3();
1346 }
1347
1348 int NN = m_densities.size();
1349
1350 MatrixXd densMatrix(2 * NN, 2 * NN);
1351 VectorXd solVector(2 * NN), xVector(2 * NN);
1352
1353 vector<double> chi2id(m_densities.size()), Ps(m_densities.size()), sids(m_densities.size());
1354 for (int i = 0; i < NN; ++i) {
1355 //chi2id[i] = m_TPS->Particles()[i].chi(2, m_Parameters, m_UseWidth, m_MuStar[i]);
1356 Ps[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_MuStar[i]);
1357 sids[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::EntropyDensity, m_UseWidth, m_MuStar[i]);
1358 //chi2id[i] = m_TPS->Particles()[i].chiDimensionfull(2, m_Parameters, m_UseWidth, m_MuStar[i]);
1359 }
1360
1361 vector<double> evc_Ps(Nevcomp, 0.);
1362 for (int i = 0; i < NN; ++i) {
1363 evc_Ps[evinds[i]] += Ps[i];
1364 }
1365
1366 for (int i = 0; i < NN; ++i)
1367 for (int j = 0; j < NN; ++j) {
1368 densMatrix(i, j) = -m_exvolmod->df(i, j) * m_DensitiesId[i];
1369 if (i == j) densMatrix(i, j) += 1.;
1370 }
1371
1372 for (int i = 0; i < NN; ++i)
1373 for (int j = 0; j < NN; ++j)
1374 densMatrix(i, NN + j) = 0.;
1375
1376 for (int i = 0; i < NN; ++i) {
1377 densMatrix(i, NN + i) = -m_exvolmod->f(i) * chi2ids[i];
1378 }
1379
1380 for (int i = 0; i < NN; ++i)
1381 for (int j = 0; j < NN; ++j) {
1382 densMatrix(NN + i, j) = m_mfmod->d2v(i, j);
1383 for (int indk = 0; indk < Nevcomp; ++indk) {
1384 int k = evindsfrom[indk];
1385 densMatrix(NN + i, j) += -m_exvolmod->d2f(k, i, j) * evc_Ps[indk];
1386 }
1387 }
1388
1389 for (int i = 0; i < NN; ++i)
1390 for (int j = 0; j < NN; ++j) {
1391 densMatrix(NN + i, NN + j) = -m_exvolmod->df(j,i) * m_DensitiesId[j];
1392 if (i == j) densMatrix(NN + i, NN + j) += 1.;
1393 }
1394
1395 PartialPivLU<MatrixXd> decomp(densMatrix);
1396
1397 for(int i = 0; i < NN; ++i) {
1398 double fi = m_exvolmod->f(i);
1399 xVector[i] = fi * dniddTs[i];
1400 xVector[NN + i] = 0.;
1401 for(int j = 0; j < NN; ++j) {
1402 xVector[NN + i] += m_exvolmod->df(j, i) * sids[j];
1403 }
1404 }
1405
1406 solVector = decomp.solve(xVector);
1407
1408 for(int i=0;i<NN;++i) {
1409 m_dndT[i] = solVector[i];
1410 m_dmusdT[i] = solVector[NN+i];
1411 }
1412
1413 // dchi2's
1414 if (IsSusceptibilitiesCalculated()) {
1415 vector<double> dnevdT(Nevcomp, 0.);
1416 for(int l = 0; l < NN; ++l) {
1417 dnevdT[evinds[l]] += m_dndT[l];
1418 }
1419
1420 vector<vector<double>> chikjd2fikldnldT(Nevcomp, vector<double>(N, 0.));
1421 for (int indi = 0; indi < Nevcomp; ++indi) {
1422 int i = evindsfrom[indi];
1423 for (int j = 0; j < N; ++j) {
1424 for (int k = 0; k < N; ++k) {
1425 for(int indl = 0; indl < Nevcomp; ++indl) {
1426 int l = evindsfrom[indl];
1427 chikjd2fikldnldT[indi][j] += m_PrimCorrel[k][j] * m_exvolmod->d2f(i, k, l) * dnevdT[indl];
1428 }
1429 }
1430 }
1431 }
1432
1433 vector<vector<double>> chikjdfik(Nevcomp, vector<double>(N, 0.));
1434 for (int indi = 0; indi < Nevcomp; ++indi) {
1435 int i = evindsfrom[indi];
1436 for (int j = 0; j < N; ++j) {
1437 for (int k = 0; k < N; ++k) {
1438 chikjdfik[indi][j] += m_PrimCorrel[k][j] * m_exvolmod->df(i, k);
1439 }
1440 }
1441 }
1442
1443 vector<double> dfikdnkdT(Nevcomp, 0.);
1444 for (int indi = 0; indi < Nevcomp; ++indi) {
1445 int i = evindsfrom[indi];
1446 for (int k = 0; k < N; ++k) {
1447 dfikdnkdT[indi] += m_exvolmod->df(i, k) * m_dndT[k];
1448 }
1449 }
1450
1451 vector<vector<double>> chikjd3vikldnldT(Nmfcomp, vector<double>(N, 0.));
1452 vector<double> dnmfdT(Nmfcomp, 0.);
1453 for(int l = 0; l < NN; ++l) {
1454 dnmfdT[mfinds[l]] += m_dndT[l];
1455 }
1456
1457 for (int indi = 0; indi < Nmfcomp; ++indi) {
1458 int i = mfindsfrom[indi];
1459 for (int j = 0; j < N; ++j) {
1460 for (int k = 0; k < N; ++k) {
1461 for(int indl = 0; indl < Nmfcomp; ++indl) {
1462 int l = mfindsfrom[indl];
1463 chikjd3vikldnldT[indi][j] += m_PrimCorrel[k][j] * m_mfmod->d3v(i, k, l) * dnmfdT[indl];
1464 }
1465 }
1466 }
1467 }
1468
1469 vector<vector<double>> chikjd3fmikldnldTPsm(Nevcomp, vector<double>(N, 0.));
1470
1471 for (int indi = 0; indi < Nevcomp; ++indi) {
1472 int i = evindsfrom[indi];
1473 for (int j = 0; j < N; ++j) {
1474 for (int k = 0; k < N; ++k) {
1475 for(int indl = 0; indl < Nevcomp; ++indl) {
1476 int l = evindsfrom[indl];
1477 for(int indm = 0; indm < Nevcomp; ++indm) {
1478 int m = evindsfrom[indm];
1479 chikjd3fmikldnldTPsm[indi][j] += m_PrimCorrel[k][j] * m_exvolmod->d3f(m, i, k, l) * dnevdT[indl] * evc_Ps[indm];
1480 }
1481 }
1482 }
1483 }
1484 }
1485
1486 vector<double> splusnidsums(Nevcomp, 0.);
1487 for(int m = 0; m < NN; ++m) {
1488 splusnidsums[evinds[m]] += sids[m] + nids[m] * m_dmusdT[m];
1489 }
1490 vector<vector<double>> chikjd2fmiketc(Nevcomp, vector<double>(N, 0.));
1491 for (int indi = 0; indi < Nevcomp; ++indi) {
1492 int i = evindsfrom[indi];
1493 for (int j = 0; j < N; ++j) {
1494 for (int k = 0; k < N; ++k) {
1495 for (int indm = 0; indm < Nevcomp; ++indm) {
1496 int m = evindsfrom[indm];
1497 chikjd2fmiketc[indi][j] += m_PrimCorrel[k][j] * m_exvolmod->d2f(m, i, k) * splusnidsums[indm];
1498 }
1499 }
1500 }
1501 }
1502
1503 vector<vector<double>> dmusmukjd2fkildnldTnidk(Nevcomp, vector<double>(N, 0.));
1504 for (int indi = 0; indi < Nevcomp; ++indi) {
1505 int i = evindsfrom[indi];
1506 for (int j = 0; j < N; ++j) {
1507 for (int k = 0; k < N; ++k) {
1508 for(int indl = 0; indl < Nevcomp; ++indl) {
1509 int l = evindsfrom[indl];
1510 dmusmukjd2fkildnldTnidk[indi][j] += m_dmusdmu[k][j] * m_exvolmod->d2f(k, i, l) * dnevdT[indl] * nids[k];
1511 }
1512 }
1513 }
1514 }
1515
1516 vector<vector<double>> dmusdmukjdfkietc(Nevcomp, vector<double>(N, 0.));
1517 for (int indi = 0; indi < Nevcomp; ++indi) {
1518 int i = evindsfrom[indi];
1519 for (int j = 0; j < N; ++j) {
1520 for (int k = 0; k < N; ++k) {
1521 dmusdmukjdfkietc[indi][j] += m_dmusdmu[k][j] * m_exvolmod->df(k, i) * (dniddTs[k] + chi2ids[k] * m_dmusdT[k]);
1522 }
1523 }
1524 }
1525
1526
1527 for (int j = 0; j < NN; ++j) {
1528 for (int i = 0; i < NN; ++i) {
1529 // a1
1530 double a1 = 0.;
1531
1532 a1 += chikjd2fikldnldT[evinds[i]][j] * nids[i];
1533 a1 += chikjdfik[evinds[i]][j] * (dniddTs[i] + chi2ids[i] * m_dmusdT[i]);
1534 a1 += m_dmusdmu[i][j] * dfikdnkdT[evinds[i]] * chi2ids[i];
1535
1536// for (int k = 0; k < NN; ++k) {
1542// }
1543 a1 += m_dmusdmu[i][j] * m_exvolmod->f(i) * (dchi2idsdT[i] + chi3ids[i] * m_dmusdT[i]);
1544
1545 xVector[i] = a1;
1546
1547
1548 // a2
1549 double a2 = 0.;
1550
1551
1552 a2 += -chikjd3vikldnldT[mfinds[i]][j];
1553
1554 a2 += chikjd3fmikldnldTPsm[evinds[i]][j];
1555
1556 a2 += chikjd2fmiketc[evinds[i]][j];
1557
1558 a2 += dmusmukjd2fkildnldTnidk[evinds[i]][j];
1559
1560 a2 += dmusdmukjdfkietc[evinds[i]][j];
1561
1562// for (int k = 0; k < NN; ++k) {
1570//
1575// }
1576
1577 xVector[i + NN] = a2;
1578 }
1579
1580 solVector = decomp.solve(xVector);
1581
1582 for (int i = 0; i < NN; ++i) {
1583 // chi2 = chi2til / T^2
1584 // dchi2/dT = dchi2til/dT / T^2 - 2 * chi2til / T / T / T
1585 m_PrimChi2sdT[i][j] = solVector[i] / (xMath::GeVtoifm3() * T * T) - 2. * m_PrimCorrel[i][j] / T / T / T / xMath::GeVtoifm3();
1586 }
1587 }
1588
1589 m_TemperatureDerivativesCalculated = true;
1590 CalculateSusceptibilityMatrix();
1591 }
1592
1593 m_TemperatureDerivativesCalculated = true;
1594 }
1595
1596 void ThermalModelRealGas::CalculateFluctuations()
1597 {
1598 CalculateTwoParticleCorrelations();
1599 CalculateSusceptibilityMatrix();
1600 CalculateTwoParticleFluctuationsDecays();
1601 CalculateProxySusceptibilityMatrix();
1602 CalculateParticleChargeCorrelationMatrix();
1603
1604 for (size_t i = 0; i < m_wprim.size(); ++i) {
1605 m_skewprim[i] = 1.;
1606 m_kurtprim[i] = 1.;
1607 m_skewtot[i] = 1.;
1608 m_kurttot[i] = 1.;
1609 }
1610
1611 m_FluctuationsCalculated = true;
1612
1613 if (IsTemperatureDerivativesCalculated()) {
1614 m_TemperatureDerivativesCalculated = false;
1615 CalculateTemperatureDerivatives();
1616 }
1617 }
1618
1619 double ThermalModelRealGas::CalculateEnergyDensityDerivativeT() {
1620 if (!IsTemperatureDerivativesCalculated())
1621 CalculateTemperatureDerivatives();
1622
1623 // Compute dE/dT
1624 double ret = 0.;
1625
1626 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
1627 const ThermalParticle &part = m_TPS->Particles()[i];
1628 // Term 1
1629 for (int k = 0; k < m_TPS->ComponentsNumber(); ++k) {
1630 double dfik = m_exvolmod->df(i, k);
1631 ret += dfik * m_dndT[k] * part.Density(m_Parameters, IdealGasFunctions::EnergyDensity, m_UseWidth, m_MuStar[i]);
1632 }
1633
1634 double fi = m_exvolmod->f(i);
1635 double dvi = m_mfmod->dv(i);
1636 // Term 2
1637 ret += fi * part.Density(m_Parameters, IdealGasFunctions::dedT, m_UseWidth, m_MuStar[i]);
1638
1639 // Term 3
1640 ret += fi * part.Density(m_Parameters, IdealGasFunctions::dedmu, m_UseWidth, m_MuStar[i]) * m_dmusdT[i];
1641
1642 // Term 4
1643 ret += dvi * m_dndT[i];
1644 }
1645
1646 // No T dependence implemented yet
1647 assert(m_mfmod->dvdT() == 0.0);
1648
1649 return ret;
1650 }
1651
1652 double ThermalModelRealGas::CalculateEnergyDensity() {
1653 if (!m_Calculated) CalculateDensities();
1654 double ret = 0.;
1655 for (size_t i = 0; i < m_densities.size(); ++i)
1656 ret += m_exvolmod->f(i) * m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::EnergyDensity, m_UseWidth, m_MuStar[i]);
1657 ret += m_mfmod->v();
1658
1659 if (1 /*&& m_TemperatureDependentAB*/) {
1660 for (size_t i = 0; i < m_densities.size(); ++i) {
1661 double tPid = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_MuStar[i]);
1662 ret += -m_Parameters.T * tPid * m_exvolmod->dfdT(i);
1663 }
1664 ret += m_Parameters.T * m_mfmod->dvdT();
1665 }
1666
1667 return ret;
1668 }
1669
1670 double ThermalModelRealGas::CalculateEntropyDensity() {
1671 if (!m_Calculated) CalculateDensities();
1672 double ret = 0.;
1673 for (size_t i = 0; i < m_densities.size(); ++i)
1674 ret += m_exvolmod->f(i) * m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::EntropyDensity, m_UseWidth, m_MuStar[i]);
1675
1676 if (1 /*&& m_TemperatureDependentAB*/) {
1677 for (size_t i = 0; i < m_densities.size(); ++i) {
1678 double tPid = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_MuStar[i]);
1679 ret += -m_Parameters.T * tPid * m_exvolmod->dfdT(i);
1680 }
1681 ret += m_Parameters.T * m_mfmod->dvdT();
1682 }
1683
1684 return ret;
1685 }
1686
1687 double ThermalModelRealGas::CalculatePressure() {
1688 if (!m_Calculated) CalculateDensities();
1689 double ret = 0.;
1690 for (size_t i = 0; i < m_densities.size(); ++i) {
1691 double tfsum = 0.;
1692 for (size_t j = 0; j < m_densities.size(); ++j) {
1693 tfsum += m_densities[j] * m_exvolmod->df(i, j);
1694 }
1695 tfsum = m_exvolmod->f(i) - tfsum;
1696 ret += tfsum * m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, m_MuStar[i]);
1697 }
1698
1699 ret += -m_mfmod->v();
1700 for (size_t i = 0; i < m_densities.size(); ++i)
1701 ret += m_mfmod->dv(i) * m_densities[i];
1702
1703 return ret;
1704 }
1705
1706 double ThermalModelRealGas::ParticleScalarDensity(int part) {
1707 if (!m_Calculated) CalculateDensities();
1708
1709 return m_scaldens[part];
1710 }
1711
1712 double ThermalModelRealGas::MuShift(int id) const
1713 {
1714 if (id >= 0. && id < static_cast<int>(ComponentsNumber()))
1715 return m_MuStar[id] - m_Chem[id];
1716 else
1717 return 0.0;
1718 }
1719
1720 std::vector<double> ThermalModelRealGas::BroydenEquationsRealGas::Equations(const std::vector<double>& x)
1721 {
1722 int NN = m_THM->Densities().size();
1723 vector<double> Ps(NN, 0.);
1724 for (int i = 0; i < NN; ++i) {
1725 Ps[i] = m_THM->TPS()->Particles()[i].Density(m_THM->Parameters(),
1726 IdealGasFunctions::Pressure,
1727 m_THM->UseWidth(),
1728 m_THM->ChemicalPotential(i) + x[i]
1729 );
1730 }
1731
1732 vector<double> ns(NN, 0.);
1733 for (int i = 0; i < NN; ++i) {
1734 ns[i] = m_THM->TPS()->Particles()[i].Density(m_THM->Parameters(),
1735 IdealGasFunctions::ParticleDensity,
1736 m_THM->UseWidth(),
1737 m_THM->ChemicalPotential(i) + x[i]
1738 );
1739 }
1740
1741 vector<double> np = m_THM->m_exvolmod->nsol(ns);
1742 m_THM->m_exvolmod->SetDensities(np);
1743 m_THM->m_mfmod->SetDensities(np);
1744
1745 vector<double> ret(m_N, 0.);
1746 for (size_t i = 0; i < ret.size(); ++i) {
1747 ret[i] = x[i];
1748 for (int j = 0; j < NN; ++j)
1749 ret[i] += -m_THM->m_exvolmod->df(j, i) * Ps[j];
1750 ret[i] += m_THM->m_mfmod->dv(i);
1751 }
1752 return ret;
1753 }
1754
1755 std::vector<double> ThermalModelRealGas::BroydenJacobianRealGas::Jacobian(const std::vector<double>& x)
1756 {
1757 int NN = m_THM->m_densities.size();
1758
1759 MatrixXd densMatrix(NN, NN);
1760 VectorXd solVector(NN), xVector(NN);
1761
1762 std::vector<double> ret(NN * NN, 0.);
1763 {
1764 vector<double> Ps(NN, 0.);
1765 for (int i = 0; i < NN; ++i)
1766 Ps[i] = m_THM->TPS()->Particles()[i].Density(m_THM->Parameters(),
1767 IdealGasFunctions::Pressure,
1768 m_THM->UseWidth(),
1769 m_THM->ChemicalPotential(i) + x[i]
1770 );
1771
1772 vector<double> ns(NN, 0.);
1773 for (int i = 0; i < NN; ++i)
1774 ns[i] = m_THM->TPS()->Particles()[i].Density(m_THM->Parameters(),
1775 IdealGasFunctions::ParticleDensity,
1776 m_THM->UseWidth(),
1777 m_THM->ChemicalPotential(i) + x[i]
1778 );
1779
1780 vector<double> chi2s(NN, 0.);
1781 for (int i = 0; i < NN; ++i)
1782 chi2s[i] = m_THM->TPS()->Particles()[i].chiDimensionfull(2, m_THM->Parameters(),
1783 m_THM->UseWidth(),
1784 m_THM->ChemicalPotential(i) + x[i]
1785 );
1786
1787 vector<double> np = m_THM->m_exvolmod->nsol(ns);
1788 m_THM->m_exvolmod->SetDensities(np);
1789 m_THM->m_mfmod->SetDensities(np);
1790
1791 for (int i = 0; i < NN; ++i) {
1792 for (int j = 0; j < NN; ++j) {
1793 densMatrix(i, j) = 0.;
1794 if (i == j)
1795 densMatrix(i, j) += 1.;
1796
1797 densMatrix(i, j) += -m_THM->m_exvolmod->df(i, j) * ns[i];
1798 }
1799 }
1800
1801 int Nevcomp = m_THM->m_exvolmod->ComponentsNumber();
1802 const vector<int>& evinds = m_THM->m_exvolmod->ComponentIndices();
1803 const vector<int>& evindsfrom = m_THM->m_exvolmod->ComponentIndicesFrom();
1804
1805 int Nmfcomp = m_THM->m_mfmod->ComponentsNumber();
1806 const vector<int>& mfinds = m_THM->m_mfmod->ComponentIndices();
1807 const vector<int>& mfindsfrom = m_THM->m_mfmod->ComponentIndicesFrom();
1808
1809 vector<double> evc_Ps(Nevcomp, 0.);
1810 for (int i = 0; i < NN; ++i)
1811 evc_Ps[m_THM->m_exvolmod->ComponentIndices()[i]] += Ps[i];
1812
1813 PartialPivLU<MatrixXd> decomp(densMatrix);
1814
1815 for (int kp = 0; kp < NN; ++kp) {
1816
1817 if (1 /*attrfl*/) {
1818 for (int l = 0; l < NN; ++l) {
1819 xVector[l] = 0.;
1820 if (l == kp) {
1821 xVector[l] = chi2s[l] * pow(xMath::GeVtoifm(), 3) * m_THM->m_exvolmod->f(l);
1822 }
1823 }
1824
1825 solVector = decomp.solve(xVector);
1826 for (int i = 0; i < NN; ++i)
1827 if (solVector[i] > 1.) solVector[i] = 1.; // Stabilizer
1828 }
1829
1830 vector<double> dnjdmukp(NN, 0.);
1831 vector<double> evc_dnjdmukp(Nevcomp, 0.);
1832 vector<double> evc_d2fmul(Nevcomp, 0.);
1833 vector<double> mfc_dnjdmukp(Nmfcomp, 0.);
1834 vector<double> mfc_d2vmul(Nmfcomp, 0.);
1835 if (1 /*attrfl*/) {
1836 for (int j = 0; j < NN; ++j) {
1837 dnjdmukp[j] = solVector[j];
1838 evc_dnjdmukp[evinds[j]] += solVector[j];
1839 mfc_dnjdmukp[mfinds[j]] += solVector[j];
1840 }
1841 for (int nn = 0; nn < Nevcomp; ++nn) {
1842 for (int in = 0; in < Nevcomp; ++in) {
1843 for (int inp = 0; inp < Nevcomp; ++inp) {
1844 int n = evindsfrom[in];
1845 int np = evindsfrom[inp];
1846 int k = evindsfrom[nn];
1847 evc_d2fmul[nn] += -evc_Ps[in] * m_THM->m_exvolmod->d2f(n, k, np) * evc_dnjdmukp[inp];
1848 }
1849 }
1850 }
1851
1852 for (int nn = 0; nn < Nmfcomp; ++nn) {
1853 for (int in = 0; in < Nmfcomp; ++in) {
1854 int n = mfindsfrom[in];
1855 int k = mfindsfrom[nn];
1856 mfc_d2vmul[nn] += m_THM->m_mfmod->d2v(k, n) * mfc_dnjdmukp[in];
1857 }
1858 }
1859 }
1860
1861
1862 for (int k = 0; k < NN; ++k) {
1863 if (k == kp)
1864 ret[k * NN + kp] += 1.;
1865 ret[k * NN + kp] += -m_THM->m_exvolmod->df(kp, k) * ns[kp];
1866
1867 if (1 /*attrfl*/) {
1868 ret[k * NN + kp] += evc_d2fmul[evinds[k]];
1869
1870 ret[k * NN + kp] += mfc_d2vmul[mfinds[k]];
1871
1872 }
1873 }
1874 }
1875 }
1876
1877 return ret;
1878 }
1879
1880 bool ThermalModelRealGas::BroydenSolutionCriteriumRealGas::IsSolved(const std::vector<double>& x, const std::vector<double>& f, const std::vector<double>& xdelta) const
1881 {
1882 if (xdelta.size() == x.size()) {
1883 double maxdiff = 0.;
1884 for (size_t i = 0; i < xdelta.size(); ++i) {
1885 maxdiff = std::max(maxdiff, fabs(xdelta[i]));
1886 maxdiff = std::max(maxdiff, fabs(f[i]));
1887 }
1888 return (maxdiff < m_MaximumError);
1889 }
1890 else {
1891 return Broyden::BroydenSolutionCriterium::IsSolved(x, f, xdelta);
1892 }
1893 }
1894
1895 std::vector<double> ThermalModelRealGas::BroydenEquationsRealGasComponents::Equations(const std::vector<double> &x)
1896 {
1897 int NN = m_THM->Densities().size();
1898 vector<double> Ps(NN, 0.);
1899 for (int i = 0; i < NN; ++i) {
1900 Ps[i] = m_THM->TPS()->Particles()[i].Density(m_THM->Parameters(),
1901 IdealGasFunctions::Pressure,
1902 m_THM->UseWidth(),
1903 m_THM->ChemicalPotential(i) + x[m_THM->m_MapTodMuStar[i]]
1904 );
1905 }
1906
1907 vector<double> ns(NN, 0.);
1908 for (int i = 0; i < NN; ++i) {
1909 ns[i] = m_THM->TPS()->Particles()[i].Density(m_THM->Parameters(),
1910 IdealGasFunctions::ParticleDensity,
1911 m_THM->UseWidth(),
1912 m_THM->ChemicalPotential(i) + x[m_THM->m_MapTodMuStar[i]]
1913 );
1914 }
1915
1916 vector<double> np = m_THM->m_exvolmod->nsol(ns);
1917 m_THM->m_exvolmod->SetDensities(np);
1918 m_THM->m_mfmod->SetDensities(np);
1919
1920 vector<double> ret(m_N, 0.);
1921 for (size_t i = 0; i < ret.size(); ++i) {
1922 ret[i] = x[i];
1923 for (int j = 0; j < NN; ++j)
1924 ret[i] += -m_THM->m_exvolmod->df(j, m_THM->m_MapFromdMuStar[i]) * Ps[j];
1925 ret[i] += m_THM->m_mfmod->dv(m_THM->m_MapFromdMuStar[i]);
1926 }
1927 return ret;
1928 }
1929
1930 std::vector<double> ThermalModelRealGas::BroydenJacobianRealGasComponents::Jacobian(const std::vector<double> &x)
1931 {
1932 int NN = m_THM->m_densities.size();
1933 int NNdmu = m_THM->m_MapFromdMuStar.size();
1934
1935 MatrixXd densMatrix(NNdmu, NNdmu);
1936 VectorXd solVector(NNdmu), xVector(NNdmu);
1937
1938 std::vector<double> ret(NNdmu * NNdmu, 0.);
1939 {
1940 vector<double> Ps(NN, 0.);
1941 for (int i = 0; i < NN; ++i)
1942 Ps[i] = m_THM->TPS()->Particles()[i].Density(m_THM->Parameters(),
1943 IdealGasFunctions::Pressure,
1944 m_THM->UseWidth(),
1945 m_THM->ChemicalPotential(i) + x[m_THM->m_MapTodMuStar[i]]
1946 );
1947
1948 vector<double> ns(NN, 0.);
1949 for (int i = 0; i < NN; ++i)
1950 ns[i] = m_THM->TPS()->Particles()[i].Density(m_THM->Parameters(),
1951 IdealGasFunctions::ParticleDensity,
1952 m_THM->UseWidth(),
1953 m_THM->ChemicalPotential(i) + x[m_THM->m_MapTodMuStar[i]]
1954 );
1955
1956 vector<double> chi2s(NN, 0.);
1957 for (int i = 0; i < NN; ++i)
1958 chi2s[i] = m_THM->TPS()->Particles()[i].chiDimensionfull(2, m_THM->Parameters(),
1959 m_THM->UseWidth(),
1960 m_THM->ChemicalPotential(i) + x[m_THM->m_MapTodMuStar[i]]
1961 );
1962
1963 vector<double> evc_Ps(NNdmu, 0.);
1964 for (int i = 0; i < NN; ++i)
1965 evc_Ps[m_THM->m_MapTodMuStar[i]] += Ps[i];
1966 vector<double> ns_comps(NNdmu, 0.);
1967 for (int i = 0; i < NN; ++i)
1968 ns_comps[m_THM->m_MapTodMuStar[i]] += ns[i];
1969 vector<double> chi2_comps(NNdmu, 0.);
1970 for (int i = 0; i < NN; ++i)
1971 chi2_comps[m_THM->m_MapTodMuStar[i]] += chi2s[i];
1972
1973 vector<double> np = m_THM->m_exvolmod->nsol(ns);
1974 m_THM->m_exvolmod->SetDensities(np);
1975 m_THM->m_mfmod->SetDensities(np);
1976
1977 for (int i = 0; i < NNdmu; ++i) {
1978 for (int j = 0; j < NNdmu; ++j) {
1979 densMatrix(i, j) = 0.;
1980 if (i == j)
1981 densMatrix(i, j) += 1.;
1982
1983 densMatrix(i, j) += -m_THM->m_exvolmod->df(m_THM->m_MapFromdMuStar[i], m_THM->m_MapFromdMuStar[j])
1984 * ns_comps[i];
1985 }
1986 }
1987
1988
1989
1990 PartialPivLU<MatrixXd> decomp(densMatrix);
1991
1992 for (int kp = 0; kp < NNdmu; ++kp) {
1993
1994 for (int l = 0; l < NNdmu; ++l) {
1995 xVector[l] = 0.;
1996 if (l == kp) {
1997 xVector[l] = chi2_comps[l] * pow(xMath::GeVtoifm(), 3)
1998 * m_THM->m_exvolmod->f(m_THM->m_MapFromdMuStar[l]);
1999 }
2000 }
2001
2002 solVector = decomp.solve(xVector);
2003 for (int i = 0; i < NNdmu; ++i)
2004 if (solVector[i] > 1.) solVector[i] = 1.; // Stabilizer
2005
2006 vector<double> dnjdmukp(NNdmu, 0.);
2007 for (int j = 0; j < NNdmu; ++j) {
2008 dnjdmukp[j] = solVector[j];
2009 }
2010
2011 vector<double> evc_d2fmul(NNdmu, 0.);
2012 vector<double> mfc_d2vmul(NNdmu, 0.);
2013 for (int nn = 0; nn < NNdmu; ++nn) {
2014 for (int in = 0; in < NNdmu; ++in) {
2015 for (int inp = 0; inp < NNdmu; ++inp) {
2016 int n = m_THM->m_MapFromdMuStar[in];
2017 int np = m_THM->m_MapFromdMuStar[inp];
2018 int k = m_THM->m_MapFromdMuStar[nn];
2019 evc_d2fmul[nn] += -evc_Ps[in] * m_THM->m_exvolmod->d2f(n, k, np) * dnjdmukp[inp];
2020 }
2021 }
2022 }
2023
2024 for (int nn = 0; nn < NNdmu; ++nn) {
2025 for (int in = 0; in < NNdmu; ++in) {
2026 int n = m_THM->m_MapFromdMuStar[in];
2027 int k = m_THM->m_MapFromdMuStar[nn];
2028 mfc_d2vmul[nn] += m_THM->m_mfmod->d2v(k, n) * dnjdmukp[in];
2029 }
2030 }
2031
2032 for (int k = 0; k < NNdmu; ++k) {
2033 if (k == kp)
2034 ret[k * NNdmu + kp] += 1.;
2035 ret[k * NNdmu + kp] += -m_THM->m_exvolmod->df(m_THM->m_MapFromdMuStar[kp], m_THM->m_MapFromdMuStar[k])
2036 * ns_comps[kp];
2037
2038 ret[k * NNdmu + kp] += evc_d2fmul[k];
2039 ret[k * NNdmu + kp] += mfc_d2vmul[k];
2040 }
2041 }
2042 }
2043
2044 return ret;
2045 }
2046
2047}
ExcludedVolumeHelper.h
Contains some functions to deal with excluded volumes.
params
map< string, double > params
Definition NeutronStars-CSHRG.cpp:34
thermalfist::Broyden::BroydenSolutionCriterium::IsSolved
virtual bool IsSolved(const std::vector< double > &x, const std::vector< double > &f, const std::vector< double > &xdelta=std::vector< double >()) const
Definition Broyden.cpp:181
thermalfist::BroydenEquations::m_N
int m_N
The number of equations.
Definition Broyden.h:66
thermalfist::Broyden
Class implementing the Broyden method to solve a system of non-linear equations.
Definition Broyden.h:132
thermalfist::Broyden::MaxDifference
double MaxDifference() const
Definition Broyden.h:239
thermalfist::Broyden::Solve
virtual std::vector< double > Solve(const std::vector< double > &x0, BroydenSolutionCriterium *solcrit=NULL, int max_iterations=MAX_ITERS)
Definition Broyden.cpp:83
thermalfist::Broyden::Iterations
int Iterations() const
Definition Broyden.h:229
thermalfist::Broyden::MaxIterations
int MaxIterations() const
Definition Broyden.h:234
thermalfist::ExcludedVolumeModelMultiBase
Base class for multi-component excluded volume models.
Definition ExcludedVolumeModelsMulti.h:25
thermalfist::MeanFieldModelMultiBase
Base class for multi-component mean field models.
Definition MeanFieldModelsMulti.h:25
thermalfist::ThermalModelBase::CalculateSusceptibilityMatrix
virtual void CalculateSusceptibilityMatrix()
Calculates the conserved charges susceptibility matrix.
Definition ThermalModelBase.cpp:1292
thermalfist::ThermalModelBase::CalculateParticleChargeCorrelationMatrix
virtual void CalculateParticleChargeCorrelationMatrix()
Calculates the matrix of correlators between primordial (and also final) particle numbers and conserv...
Definition ThermalModelBase.cpp:1366
thermalfist::ThermalModelBase::FillChemicalPotentials
virtual void FillChemicalPotentials()
Sets the chemical potentials of all particles.
Definition ThermalModelBase.cpp:357
thermalfist::ThermalModelBase::SetChemicalPotentials
virtual void SetChemicalPotentials(const std::vector< double > &chem=std::vector< double >(0))
Sets the chemical potentials of all particles.
Definition ThermalModelBase.cpp:363
thermalfist::ThermalModelBase::ChemicalPotential
double ChemicalPotential(int i) const
Chemical potential of particle species i.
Definition ThermalModelBase.cpp:372
thermalfist::ThermalModelBase::ValidateCalculation
virtual void ValidateCalculation()
Checks whether issues have occured during the calculation of particle densities in the CalculateDensi...
Definition ThermalModelBase.cpp:633
thermalfist::ThermalModelBase::ComponentsNumber
int ComponentsNumber() const
Number of different particle species in the list.
Definition ThermalModelBase.h:66
thermalfist::ThermalModelBase::CalculateProxySusceptibilityMatrix
virtual void CalculateProxySusceptibilityMatrix()
Calculates the susceptibility matrix of conserved charges proxies.
Definition ThermalModelBase.cpp:1326
thermalfist::ThermalModelBase::ThermalModelBase
ThermalModelBase(ThermalParticleSystem *TPS, const ThermalModelParameters &params=ThermalModelParameters())
Construct a new ThermalModelBase object.
Definition ThermalModelBase.cpp:24
thermalfist::ThermalModelBase::UseWidth
bool UseWidth() const
Whether finite resonance widths are considered.
Definition ThermalModelBase.h:82
thermalfist::ThermalModelBase::CalculateDensities
virtual void CalculateDensities()
Calculates the primordial and total (after decays) densities of all species.
Definition ThermalModelBase.cpp:626
thermalfist::ThermalModelBase::Parameters
const ThermalModelParameters & Parameters() const
Definition ThermalModelBase.h:121
thermalfist::ThermalModelBase::GCE
@ GCE
Grand canonical ensemble.
Definition ThermalModelBase.h:36
thermalfist::ThermalModelBase::CalculateTwoParticleFluctuationsDecays
virtual void CalculateTwoParticleFluctuationsDecays()
Computes particle number correlations and fluctuations for all final-state particles which are marked...
Definition ThermalModelBase.cpp:969
thermalfist::ThermalModelRealGas::CalculateFluctuations
void CalculateFluctuations()
Calculate the fluctuations.
Definition ThermalModelRealGas.cpp:1596
thermalfist::ThermalModelRealGas::CalculateComponentsMap
void CalculateComponentsMap()
Partitions particles species into sets that have identical pair interactions.
Definition ThermalModelRealGas.cpp:82
thermalfist::ThermalModelRealGas::m_dMuStarIndices
std::vector< std::vector< int > > m_dMuStarIndices
Vector of component indices for each particle species.
Definition ThermalModelRealGas.h:364
thermalfist::ThermalModelRealGas::ThermalModelRealGas
ThermalModelRealGas(ThermalParticleSystem *TPS_, const ThermalModelParameters &params=ThermalModelParameters())
Construct a new ThermalModelRealGas object.
Definition ThermalModelRealGas.cpp:25
thermalfist::ThermalModelRealGas::SearchSingleSolutionUsingComponents
virtual std::vector< double > SearchSingleSolutionUsingComponents(const std::vector< double > &muStarInit)
Uses the Broyden method with a provided initial guess to determine the shifted chemical potentials by...
Definition ThermalModelRealGas.cpp:127
thermalfist::ThermalModelRealGas::m_ComponentMapCalculated
bool m_ComponentMapCalculated
Whether the mapping to components with the same VDW parameters has been calculated.
Definition ThermalModelRealGas.h:352
thermalfist::ThermalModelRealGas::m_MapTodMuStar
std::vector< int > m_MapTodMuStar
Mapping from particle species to dMuStar indices.
Definition ThermalModelRealGas.h:358
thermalfist::ThermalModelRealGas::SearchMultipleSolutions
std::vector< double > SearchMultipleSolutions(int iters=300)
Uses the Broyden method with different initial guesses to look for different possible solutions of th...
Definition ThermalModelRealGas.cpp:188
thermalfist::ThermalModelRealGas::SearchSingleSolution
virtual std::vector< double > SearchSingleSolution(const std::vector< double > &muStarInit)
Uses the Broyden method with a provided initial guess to determine the shifted chemical potentials by...
Definition ThermalModelRealGas.cpp:121
thermalfist::ThermalModelRealGas::SearchSingleSolutionDirect
virtual std::vector< double > SearchSingleSolutionDirect(const std::vector< double > &muStarInit)
Uses the Broyden method with a provided initial guess to determine the shifted chemical potentials by...
Definition ThermalModelRealGas.cpp:160
thermalfist::ThermalModelRealGas::CalculateEnergyDensity
virtual double CalculateEnergyDensity()
Calculate the energy density of the system.
Definition ThermalModelRealGas.cpp:1652
thermalfist::ThermalModelRealGas::m_chiarb
std::vector< double > m_chiarb
Vector of computed susceptibilities for a specified arbitraty charge.
Definition ThermalModelRealGas.h:267
thermalfist::ThermalModelRealGas::m_exvolmod
ExcludedVolumeModelMultiBase * m_exvolmod
Excluded volume model used in the real gas HRG model.
Definition ThermalModelRealGas.h:368
thermalfist::ThermalModelRealGas::CalculateChargeFluctuationsOld
virtual std::vector< double > CalculateChargeFluctuationsOld(const std::vector< double > &chgs, int order=4)
Calculate the charge fluctuations of the particle species (old method).
Definition ThermalModelRealGas.cpp:920
thermalfist::ThermalModelRealGas::m_MapFromdMuStar
std::vector< int > m_MapFromdMuStar
Mapping from dMuStar indices to particle species.
Definition ThermalModelRealGas.h:361
thermalfist::ThermalModelRealGas::CalculateTemperatureDerivatives
virtual void CalculateTemperatureDerivatives()
Calculate the temperature derivatives of the system.
Definition ThermalModelRealGas.cpp:1322
thermalfist::ThermalModelRealGas::CalculateTwoParticleCorrelations
void CalculateTwoParticleCorrelations()
Calculate the two-particle correlations of the particle species.
Definition ThermalModelRealGas.cpp:1230
thermalfist::ThermalModelRealGas::CalculatePrimordialDensities
virtual void CalculatePrimordialDensities()
Calculate the primordial densities of the particle species.
Definition ThermalModelRealGas.cpp:262
thermalfist::ThermalModelRealGas::MuShift
virtual double MuShift(int id) const
The shift in the chemical potential of particle species i due to the QvdW interactions.
Definition ThermalModelRealGas.cpp:1712
thermalfist::ThermalModelRealGas::CalculateEntropyDensity
virtual double CalculateEntropyDensity()
Calculate the entropy density of the system.
Definition ThermalModelRealGas.cpp:1670
thermalfist::ThermalModelRealGas::SetChemicalPotentials
virtual void SetChemicalPotentials(const std::vector< double > &chem=std::vector< double >(0))
Set the chemical potentials of the particle species.
Definition ThermalModelRealGas.cpp:75
thermalfist::ThermalModelRealGas::m_MuStar
std::vector< double > m_MuStar
Vector of the shifted chemical potentials.
Definition ThermalModelRealGas.h:355
thermalfist::ThermalModelRealGas::m_LastBroydenSuccessFlag
bool m_LastBroydenSuccessFlag
Whether Broyden's method was successfull.
Definition ThermalModelRealGas.h:349
thermalfist::ThermalModelRealGas::CalculateEnergyDensityDerivativeT
virtual double CalculateEnergyDensityDerivativeT()
Calculate the derivative of the energy density with respect to temperature.
Definition ThermalModelRealGas.cpp:1619
thermalfist::ThermalModelRealGas::m_DensitiesId
std::vector< double > m_DensitiesId
Vector of ideal gas densities with shifted chemical potentials.
Definition ThermalModelRealGas.h:340
thermalfist::ThermalModelRealGas::CalculatePressure
virtual double CalculatePressure()
Calculate the pressure of the system.
Definition ThermalModelRealGas.cpp:1687
thermalfist::ThermalModelRealGas::~ThermalModelRealGas
virtual ~ThermalModelRealGas(void)
Destroy the ThermalModelRealGas object.
Definition ThermalModelRealGas.cpp:45
thermalfist::ThermalModelRealGas::m_exvolmodideal
ExcludedVolumeModelMultiBase * m_exvolmodideal
Excluded volume model object in the ideal gas limit.
Definition ThermalModelRealGas.h:374
thermalfist::ThermalModelRealGas::FillChemicalPotentials
void FillChemicalPotentials()
Fill the chemical potentials of the particle species.
Definition ThermalModelRealGas.cpp:69
thermalfist::ThermalModelRealGas::m_scaldens
std::vector< double > m_scaldens
Vector of scalar densities. Not used.
Definition ThermalModelRealGas.h:343
thermalfist::ThermalModelRealGas::m_mfmod
MeanFieldModelMultiBase * m_mfmod
Mean field model used in the real gas HRG model.
Definition ThermalModelRealGas.h:371
thermalfist::ThermalModelRealGas::m_mfmodideal
MeanFieldModelMultiBase * m_mfmodideal
Mean field model object in the ideal gas limit.
Definition ThermalModelRealGas.h:377
thermalfist::ThermalModelRealGas::m_chi
std::vector< std::vector< double > > m_chi
Vector of computed susceptibilities values.
Definition ThermalModelRealGas.h:262
thermalfist::ThermalModelRealGas::CalculateChargeFluctuations
virtual std::vector< double > CalculateChargeFluctuations(const std::vector< double > &chgs, int order=4)
Calculate the charge fluctuations of the particle species.
Definition ThermalModelRealGas.cpp:284
thermalfist::ThermalModelRealGas::ParticleScalarDensity
virtual double ParticleScalarDensity(int part)
Calculate the scalar density of a particle species.
Definition ThermalModelRealGas.cpp:1706
thermalfist::ThermalModelRealGas::m_SearchMultipleSolutions
bool m_SearchMultipleSolutions
Whether multiple solutions are considered.
Definition ThermalModelRealGas.h:346
thermalfist::ThermalParticle
Class containing all information about a particle specie.
Definition ThermalParticle.h:64
thermalfist::ThermalParticle::Density
double Density(const ThermalModelParameters &params, IdealGasFunctions::Quantity type=IdealGasFunctions::ParticleDensity, bool useWidth=0, double mu=0.) const
Computes a specified ideal gas thermodynamic function.
Definition ThermalParticle.cpp:584
thermalfist::ThermalParticleSystem
Class containing the particle list.
Definition ThermalParticleSystem.h:30
thermalfist::xMath::GeVtoifm
constexpr double GeVtoifm()
A constant to transform GeV into fm .
Definition xMath.h:25
thermalfist::xMath::GeVtoifm3
constexpr double GeVtoifm3()
A constant to transform GeV into fm .
Definition xMath.h:31
thermalfist::xMath::mnucleon
constexpr double mnucleon()
Nucleon's mass. Value as in UrQMD.
Definition xMath.h:34
thermalfist
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition CosmicEoS.h:9
thermalfist::ThermalModelParameters
Structure containing all thermal parameters of the model.
Definition ThermalModelParameters.h:18
xMath.h
Contains some extra mathematical functions used in the code.
Generated by doxygen 1.13.2