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