Thermal-FIST 1.6
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
ThermalModelBase.cpp
Go to the documentation of this file.
1/*
2 * Thermal-FIST package
3 *
4 * Copyright (c) 2014-2025 Volodymyr Vovchenko
5 *
6 * GNU General Public License (GPLv3 or later)
7 */
8#include "HRGBase/ThermalModelBase.h"
9
10#include <cstdio>
11#include <algorithm>
12#include <cmath>
13#include <limits>
14#include <set>
15#include <sstream>
16
17#include <Eigen/Dense>
18
19#include "HRGBase/Utility.h"
20#include "HRGBase/ThermalParticleSystem.h"
21
22using namespace Eigen;
23
24using namespace std;
25
26namespace thermalfist {
27
28 namespace {
29 std::string ConservedChargeTag(const std::vector<int>& conserved)
30 {
31 const char labels[4] = {'B', 'Q', 'S', 'C'};
32 std::string tag;
33 for (int i = 0; i < 4; ++i) {
34 if (conserved[i] == 1)
35 tag.push_back(labels[i]);
36 }
37 if (tag.empty())
38 tag = "none";
39 return tag;
40 }
41
42 void WarnIllDefinedSoundSpeedOnce(const char* func,
43 const std::vector<int>& conserved,
44 double T,
45 const std::string& reason)
46 {
47 static std::set<std::string> warned;
48 std::ostringstream key;
49 key << func << "|" << ConservedChargeTag(conserved) << "|" << reason;
50 if (warned.insert(key.str()).second) {
51 std::cerr << "**WARNING** " << func << ": ill-defined for conserved set "
52 << ConservedChargeTag(conserved)
53 << " at T = " << T << " GeV. " << reason << std::endl;
54 }
55 }
56
57 bool IsMatrixInvertible(const MatrixXd& m, double threshold = 1e-14)
58 {
59 FullPivLU<MatrixXd> lu(m);
60 lu.setThreshold(threshold);
61 return lu.isInvertible();
62 }
63 } // namespace
64
65 ThermalModelBase::ThermalModelBase(ThermalParticleSystem *TPS_, const ThermalModelParameters& params) :
66 m_TPS(TPS_),
67 m_Parameters(params),
68 m_UseWidth(false),
69 m_PCE(false),
70 m_Calculated(false),
71 m_FeeddownCalculated(false),
72 m_TwoParticleCorrelationsCalculated(false),
73 m_FluctuationsCalculated(false),
74 m_TemperatureDerivativesCalculated(false),
75 m_GCECalculated(false),
76 m_NormBratio(false),
77 m_QuantumStats(true),
78 m_MaxDiff(0.),
79 m_useOpenMP(0),
80 m_IGFExtraConfig()
81 {
82 if (!Disclaimer::DisclaimerPrinted) {
83 Disclaimer::DisclaimerPrinted = Disclaimer::PrintDisclaimer();
84 }
85
86 m_QBgoal = 0.4;
87 m_SBgoal = 50.;
88 m_Chem.resize(m_TPS->Particles().size());
89 m_Volume = params.V;
90 m_densities.resize(m_TPS->Particles().size());
91 m_densitiestotal.resize(m_TPS->Particles().size());
92 m_densitiesbyfeeddown = std::vector< std::vector<double> >(ParticleDecayType::NumberOfDecayTypes, m_densitiestotal);
93
94 m_wprim.resize(m_TPS->Particles().size());
95 m_wtot.resize(m_TPS->Particles().size());
96 m_skewprim.resize(m_TPS->Particles().size());
97 m_skewtot.resize(m_TPS->Particles().size());
98 m_kurtprim.resize(m_TPS->Particles().size());
99 m_kurttot.resize(m_TPS->Particles().size());
100
101 m_ConstrainMuB = false;
102 m_ConstrainMuC = m_ConstrainMuQ = m_ConstrainMuS = true;
103
104 m_Susc.resize(4);
105 for (int i = 0; i < 4; ++i) m_Susc[i].resize(4);
106 m_dSuscdT = m_Susc;
107
108 m_NormBratio = false;
109
110 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
111 const ThermalParticle &tpart = m_TPS->Particles()[i];
112 for (size_t j = 0; j < tpart.Decays().size(); ++j) {
113 if (tpart.DecaysOriginal().size() == tpart.Decays().size() && tpart.Decays()[j].mBratio != tpart.DecaysOriginal()[j].mBratio)
114 m_NormBratio = true;
115 }
116 }
117
118 m_PrimCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(TPS()->ComponentsNumber(), 0.));
119 m_TotalCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(TPS()->ComponentsNumber(), 0.));
120 m_PrimChargesCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(4, 0.));
121 m_FinalChargesCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(4, 0.));
122
123 m_Ensemble = GCE;
124 m_InteractionModel = Ideal;
125
126 //SetStatistics(m_QuantumStats);
127 //SetCalculationType(IdealGasFunctions::Quadratures);
128 SetUseWidth(TPS()->ResonanceWidthIntegrationType());
129
130 ResetCalculatedFlags();
131
132 m_ValidityLog = "";
133 }
134
135
136 void ThermalModelBase::FillVirial(const std::vector<double>& /*ri*/)
137 {
138 }
139
140 void ThermalModelBase::SetUseWidth(bool useWidth)
141 {
142 if (!useWidth && m_TPS->ResonanceWidthIntegrationType() != ThermalParticle::ZeroWidth) {
143 m_TPS->SetResonanceWidthIntegrationType(ThermalParticle::ZeroWidth);
144 //m_TPS->ProcessDecays();
145 }
146 if (useWidth && m_TPS->ResonanceWidthIntegrationType() == ThermalParticle::ZeroWidth) {
147 m_TPS->SetResonanceWidthIntegrationType(ThermalParticle::BWTwoGamma);
148 }
149 m_UseWidth = useWidth;
150 }
151
152 void ThermalModelBase::SetUseWidth(ThermalParticle::ResonanceWidthIntegration type)
153 {
154 m_UseWidth = (type != ThermalParticle::ZeroWidth);
155 m_TPS->SetResonanceWidthIntegrationType(type);
156 }
157
158
159 void ThermalModelBase::SetNormBratio(bool normBratio) {
160 if (normBratio != m_NormBratio) {
161 m_NormBratio = normBratio;
162 if (m_NormBratio) {
163 m_TPS->NormalizeBranchingRatios();
164 }
165 else {
166 m_TPS->RestoreBranchingRatios();
167 }
168 }
169 }
170
171
172 void ThermalModelBase::ResetChemicalPotentials() {
173 m_Parameters.muS = m_Parameters.muB / 5.;
174 m_Parameters.muQ = -m_Parameters.muB / 50.;
175 m_Parameters.muC = 0.;
176 }
177
178
179 void ThermalModelBase::SetParameters(const ThermalModelParameters& params) {
180 m_Parameters = params;
181 m_Volume = m_Parameters.V;
182 FillChemicalPotentials();
183 ResetCalculatedFlags();
184 }
185
186 void ThermalModelBase::SetTemperature(double T)
187 {
188 m_Parameters.T = T;
189 ResetCalculatedFlags();
190 }
191
192 void ThermalModelBase::SetBaryonChemicalPotential(double muB)
193 {
194 m_Parameters.muB = muB;
195 FillChemicalPotentials();
196 ResetCalculatedFlags();
197 }
198
199 void ThermalModelBase::SetElectricChemicalPotential(double muQ)
200 {
201 m_Parameters.muQ = muQ;
202 FillChemicalPotentials();
203 ResetCalculatedFlags();
204 }
205
206 void ThermalModelBase::SetStrangenessChemicalPotential(double muS)
207 {
208 m_Parameters.muS = muS;
209 FillChemicalPotentials();
210 ResetCalculatedFlags();
211 }
212
213 void ThermalModelBase::SetCharmChemicalPotential(double muC)
214 {
215 m_Parameters.muC = muC;
216 FillChemicalPotentials();
217 ResetCalculatedFlags();
218 }
219
220 void ThermalModelBase::SetGammaS(double gammaS)
221 {
222 m_Parameters.gammaS = gammaS;
223 ResetCalculatedFlags();
224 }
225
226 void ThermalModelBase::SetGammaC(double gammaC)
227 {
228 m_Parameters.gammaC = gammaC;
229 ResetCalculatedFlags();
230 }
231
232 void ThermalModelBase::SetBaryonCharge(int B)
233 {
234 m_Parameters.B = B;
235 ResetCalculatedFlags();
236 }
237
238 void ThermalModelBase::SetElectricCharge(int Q)
239 {
240 m_Parameters.Q = Q;
241 ResetCalculatedFlags();
242 }
243
244 void ThermalModelBase::SetStrangeness(int S)
245 {
246 m_Parameters.S = S;
247 ResetCalculatedFlags();
248 }
249
250 void ThermalModelBase::SetCharm(int C)
251 {
252 m_Parameters.C = C;
253 ResetCalculatedFlags();
254 }
255
256 void ThermalModelBase::DisableMesonMesonVirial()
257 {
258 for (int i = 0; i < TPS()->ComponentsNumber(); ++i) {
259 for (int j = 0; j < TPS()->ComponentsNumber(); ++j) {
260 const ThermalParticle &part1 = TPS()->Particles()[i];
261 const ThermalParticle &part2 = TPS()->Particles()[j];
262 if (part1.BaryonCharge() == 0 && part2.BaryonCharge() == 0)
263 SetVirial(i, j, 0.);
264 }
265 }
266 }
267
268 void ThermalModelBase::DisableMesonMesonAttraction()
269 {
270 for (int i = 0; i < TPS()->ComponentsNumber(); ++i) {
271 for (int j = 0; j < TPS()->ComponentsNumber(); ++j) {
272 const ThermalParticle &part1 = TPS()->Particles()[i];
273 const ThermalParticle &part2 = TPS()->Particles()[j];
274 if (part1.BaryonCharge() == 0 && part2.BaryonCharge() == 0)
275 SetAttraction(i, j, 0.);
276 }
277 }
278 }
279
280 void ThermalModelBase::DisableMesonBaryonVirial()
281 {
282 for (int i = 0; i < TPS()->ComponentsNumber(); ++i) {
283 for (int j = 0; j < TPS()->ComponentsNumber(); ++j) {
284 const ThermalParticle &part1 = TPS()->Particles()[i];
285 const ThermalParticle &part2 = TPS()->Particles()[j];
286 if ((part1.BaryonCharge() == 0 && part2.BaryonCharge() != 0)
287 || (part1.BaryonCharge() != 0 && part2.BaryonCharge() == 0))
288 SetVirial(i, j, 0.);
289 }
290 }
291 }
292
293 void ThermalModelBase::DisableMesonBaryonAttraction()
294 {
295 for (int i = 0; i < TPS()->ComponentsNumber(); ++i) {
296 for (int j = 0; j < TPS()->ComponentsNumber(); ++j) {
297 const ThermalParticle &part1 = TPS()->Particles()[i];
298 const ThermalParticle &part2 = TPS()->Particles()[j];
299 if ((part1.BaryonCharge() == 0 && part2.BaryonCharge() != 0)
300 || (part1.BaryonCharge() != 0 && part2.BaryonCharge() == 0))
301 SetAttraction(i, j, 0.);
302 }
303 }
304 }
305
306 void ThermalModelBase::DisableBaryonBaryonVirial()
307 {
308 for (int i = 0; i < TPS()->ComponentsNumber(); ++i) {
309 for (int j = 0; j < TPS()->ComponentsNumber(); ++j) {
310 const ThermalParticle &part1 = TPS()->Particles()[i];
311 const ThermalParticle &part2 = TPS()->Particles()[j];
312 if ((part1.BaryonCharge() > 0 && part2.BaryonCharge() > 0)
313 || (part1.BaryonCharge() < 0 && part2.BaryonCharge() < 0))
314 SetVirial(i, j, 0.);
315 }
316 }
317 }
318
319 void ThermalModelBase::DisableBaryonBaryonAttraction()
320 {
321 for (int i = 0; i < TPS()->ComponentsNumber(); ++i) {
322 for (int j = 0; j < TPS()->ComponentsNumber(); ++j) {
323 const ThermalParticle &part1 = TPS()->Particles()[i];
324 const ThermalParticle &part2 = TPS()->Particles()[j];
325 if ((part1.BaryonCharge() > 0 && part2.BaryonCharge() > 0)
326 || (part1.BaryonCharge() < 0 && part2.BaryonCharge() < 0))
327 SetAttraction(i, j, 0.);
328 }
329 }
330 }
331
332 void ThermalModelBase::DisableBaryonAntiBaryonVirial()
333 {
334 for (int i = 0; i < TPS()->ComponentsNumber(); ++i) {
335 for (int j = 0; j < TPS()->ComponentsNumber(); ++j) {
336 const ThermalParticle &part1 = TPS()->Particles()[i];
337 const ThermalParticle &part2 = TPS()->Particles()[j];
338 if ((part1.BaryonCharge() > 0 && part2.BaryonCharge() < 0)
339 || (part1.BaryonCharge() < 0 && part2.BaryonCharge() > 0))
340 SetVirial(i, j, 0.);
341 }
342 }
343 }
344
345 void ThermalModelBase::DisableBaryonAntiBaryonAttraction()
346 {
347 for (int i = 0; i < TPS()->ComponentsNumber(); ++i) {
348 for (int j = 0; j < TPS()->ComponentsNumber(); ++j) {
349 const ThermalParticle &part1 = TPS()->Particles()[i];
350 const ThermalParticle &part2 = TPS()->Particles()[j];
351 if ((part1.BaryonCharge() > 0 && part2.BaryonCharge() < 0)
352 || (part1.BaryonCharge() < 0 && part2.BaryonCharge() > 0))
353 SetAttraction(i, j, 0.);
354 }
355 }
356 }
357
358 void ThermalModelBase::SetGammaq(double gammaq)
359 {
360 m_Parameters.gammaq = gammaq;
361 ResetCalculatedFlags();
362 }
363
364
365 void ThermalModelBase::ChangeTPS(ThermalParticleSystem *TPS_) {
366 m_TPS = TPS_;
367 m_Chem.resize(m_TPS->Particles().size());
368 m_densities.resize(m_TPS->Particles().size());
369 m_densitiestotal.resize(m_TPS->Particles().size());
370 m_densitiesbyfeeddown = std::vector< std::vector<double> >(ParticleDecayType::NumberOfDecayTypes, m_densitiestotal);
371 m_wprim.resize(m_TPS->Particles().size());
372 m_wtot.resize(m_TPS->Particles().size());
373 m_skewprim.resize(m_TPS->Particles().size());
374 m_skewtot.resize(m_TPS->Particles().size());
375 m_kurtprim.resize(m_TPS->Particles().size());
376 m_kurttot.resize(m_TPS->Particles().size());
377 m_PrimCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(TPS()->ComponentsNumber(), 0.));
378 m_TotalCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(TPS()->ComponentsNumber(), 0.));
379 m_PrimChargesCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(4, 0.));
380 m_FinalChargesCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(4, 0.));
381 ResetCalculatedFlags();
382 }
383
384 void ThermalModelBase::SetStatistics(bool stats) {
385 m_QuantumStats = stats;
386 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
387 m_TPS->Particle(i).UseStatistics(stats);
388 }
389
390 void ThermalModelBase::SetResonanceWidthIntegrationType(ThermalParticle::ResonanceWidthIntegration type)
391 {
392 if (!m_UseWidth) {
393 std::cerr << "**WARNING** ThermalModelBase::SetResonanceWidthIntegrationType: Using resonance widths is switched off!" << std::endl;
394 m_TPS->SetResonanceWidthIntegrationType(ThermalParticle::BWTwoGamma);
395 }
396 else
397 m_TPS->SetResonanceWidthIntegrationType(type);
398 }
399
400 void ThermalModelBase::FillChemicalPotentials() {
401 m_Chem.resize(m_TPS->Particles().size());
402 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
403 m_Chem[i] = m_TPS->Particles()[i].BaryonCharge() * m_Parameters.muB + m_TPS->Particles()[i].Strangeness() * m_Parameters.muS + m_TPS->Particles()[i].ElectricCharge() * m_Parameters.muQ + m_TPS->Particles()[i].Charm() * m_Parameters.muC;
404 }
405
406 void ThermalModelBase::SetChemicalPotentials(const std::vector<double>& chem)
407 {
408 if (chem.size() != m_TPS->Particles().size()) {
409 std::cerr << "**WARNING** " << m_TAG << "::SetChemicalPotentials(const std::vector<double> & chem): size of chem does not match number of hadrons in the list" << std::endl;
410 return;
411 }
412 m_Chem = chem;
413 }
414
415 double ThermalModelBase::ChemicalPotential(int i) const
416 {
417 if (i < 0 || i >= static_cast<int>(m_Chem.size())) {
418 throw std::out_of_range("ThermalModelBase::ChemicalPotential: i is out of bounds!");
419 }
420 return m_Chem[i];
421 }
422
423 void ThermalModelBase::SetChemicalPotential(int i, double chem)
424 {
425 if (i < 0 || i >= static_cast<int>(m_Chem.size())) {
426 throw std::out_of_range("ThermalModelBase::SetChemicalPotential: i is out of bounds!");
427 }
428 m_Chem[i] = chem;
429 }
430
431 double ThermalModelBase::FullIdealChemicalPotential(int i) const
432 {
433 if (i < 0 || i >= static_cast<int>(m_Chem.size())) {
434 throw std::out_of_range("ThermalModelBase::FullIdealChemicalPotential: i is out of bounds!");
435 }
436
437 double ret = ChemicalPotential(i);
438
439 ret += MuShift(i);
440
441 const ThermalParticle& part = m_TPS->Particles()[i];
442
443 if (!(m_Parameters.gammaq == 1.)) ret += log(m_Parameters.gammaq) * part.AbsoluteQuark() * m_Parameters.T;
444 if (!(m_Parameters.gammaS == 1. || part.AbsoluteStrangeness() == 0.)) ret += log(m_Parameters.gammaS) * part.AbsoluteStrangeness() * m_Parameters.T;
445 if (!(m_Parameters.gammaC == 1. || part.AbsoluteCharm() == 0.)) ret += log(m_Parameters.gammaC) * part.AbsoluteCharm() * m_Parameters.T;
446
447 return ret;
448 }
449
450 void ThermalModelBase::CalculateFeeddown() {
451 if (m_UseWidth && m_TPS->ResonanceWidthIntegrationType() == ThermalParticle::eBW) {
452 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
453 m_TPS->Particle(i).CalculateThermalBranchingRatios(m_Parameters, m_UseWidth, m_Chem[i] + MuShift(i));
454 }
455 m_TPS->ProcessDecays();
456 }
457
458 // Primordial
459 m_densitiesbyfeeddown[static_cast<int>(Feeddown::Primordial)] = m_densities;
460
461 // According to stability flags
462 int feed_index = static_cast<int>(Feeddown::StabilityFlag);
463 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
464 m_densitiestotal[i] = m_densities[i];
465 //const std::vector< std::pair<double, int> >& decayContributions = m_TPS->Particles()[i].DecayContributionsByFeeddown()[feed_index];
466 const ThermalParticleSystem::DecayContributionsToParticle& decayContributions = m_TPS->DecayContributionsByFeeddown()[feed_index][i];
467 for (size_t j = 0; j < decayContributions.size(); ++j)
468 if (i != decayContributions[j].second)
469 m_densitiestotal[i] += decayContributions[j].first * m_densities[decayContributions[j].second];
470 }
471
472 m_densitiesbyfeeddown[feed_index] = m_densitiestotal;
473
474 // Weak, EM, strong
475 for (feed_index = static_cast<int>(Feeddown::Weak); feed_index <= static_cast<int>(Feeddown::Strong); ++feed_index) {
476 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
477 m_densitiesbyfeeddown[feed_index][i] = m_densities[i];
478 //const std::vector< std::pair<double, int> >& decayContributions = m_TPS->Particles()[i].DecayContributionsByFeeddown()[feed_index];
479 const ThermalParticleSystem::DecayContributionsToParticle& decayContributions = m_TPS->DecayContributionsByFeeddown()[feed_index][i];
480 for (size_t j = 0; j < decayContributions.size(); ++j)
481 if (i != decayContributions[j].second)
482 m_densitiesbyfeeddown[feed_index][i] += decayContributions[j].first * m_densities[decayContributions[j].second];
483 }
484 }
485
486 m_FeeddownCalculated = true;
487 }
488
489
490 void ThermalModelBase::ConstrainChemicalPotentials(bool resetInitialValues)
491 {
492 if (resetInitialValues)
493 FixParameters();
494 else
495 FixParametersNoReset();
496 }
497
498 void ThermalModelBase::FixParameters() {
499 if (fabs(m_Parameters.muB) < 1e-6 && !m_ConstrainMuB) {
500 if (m_ConstrainMuS)
501 m_Parameters.muS = 0.;
502 if (m_ConstrainMuQ)
503 m_Parameters.muQ = 0.;
504 if (m_ConstrainMuC)
505 m_Parameters.muC = 0.;
506 FillChemicalPotentials();
507 CalculatePrimordialDensities();
508 return;
509 }
510 if (m_ConstrainMuB) {
511 m_Parameters.muB = xMath::mnucleon() / 2.;
512 }
513 double suppr = 10;
514 if (m_Parameters.muB > 0.150) suppr = 8.;
515 if (m_Parameters.muB > 0.300) suppr = 7.;
516 if (m_Parameters.muB > 0.450) suppr = 6.;
517 if (m_Parameters.muB > 0.600) suppr = 6.;
518 if (m_Parameters.muB > 0.750) suppr = 5.;
519 if (m_Parameters.muB > 0.900) suppr = 4.;
520 if (m_Parameters.muB > 1.000) suppr = 3.;
521 if (m_ConstrainMuS)
522 m_Parameters.muS = m_Parameters.muB / suppr;
523 if (m_ConstrainMuQ)
524 m_Parameters.muQ = -m_Parameters.muB / suppr / 10.;
525 if (m_ConstrainMuC)
526 m_Parameters.muC = -m_Parameters.muS;
527
528 FixParametersNoReset();
529 }
530
531 void ThermalModelBase::FixParametersNoReset() {
532 if (fabs(m_Parameters.muB) < 1e-6 && !m_ConstrainMuB) {
533 if (m_ConstrainMuS)
534 m_Parameters.muS = 0.;
535 if (m_ConstrainMuQ)
536 m_Parameters.muQ = 0.;
537 if (m_ConstrainMuC)
538 m_Parameters.muC = 0.;
539 //m_Parameters.muS = m_Parameters.muQ = m_Parameters.muC = 0.;
540 FillChemicalPotentials();
541 CalculatePrimordialDensities();
542 return;
543 }
544
545 m_ConstrainMuB &= m_TPS->hasBaryons();
546 m_ConstrainMuQ &= (m_TPS->hasCharged() && m_TPS->hasBaryons());
547 m_ConstrainMuS &= m_TPS->hasStrange();
548 m_ConstrainMuC &= m_TPS->hasCharmed();
549
550 vector<double> x22(4);
551 x22[0] = m_Parameters.muB;
552 x22[1] = m_Parameters.muQ;
553 x22[2] = m_Parameters.muS;
554 x22[3] = m_Parameters.muC;
555 vector<double> x2(4), xinit(4);
556 xinit[0] = x2[0] = m_Parameters.muB;
557 xinit[1] = x2[1] = m_Parameters.muQ;
558 xinit[2] = x2[2] = m_Parameters.muS;
559 xinit[3] = x2[3] = m_Parameters.muC;
560 int iter = 0, iterMAX = 2;
561 while (iter < iterMAX) {
562 BroydenEquationsChem eqs(this);
563 BroydenJacobianChem jaco(this);
564 BroydenChem broydn(this, &eqs, &jaco);
565 Broyden::BroydenSolutionCriterium crit(1.0E-8);
566 broydn.Solve(x22, &crit);
567 //std::cerr << "Broyden iters: " << broydn.Iterations() << std::endl;
568 break;
569 }
570 }
571
572 bool ThermalModelBase::FixChemicalPotentialsThroughDensities(double rhoB, double rhoQ, double rhoS, double rhoC,
573 double muBinit, double muQinit, double muSinit, double muCinit,
574 bool ConstrMuB, bool ConstrMuQ, bool ConstrMuS, bool ConstrMuC) {
575 double V = Parameters().V;
576 return SolveChemicalPotentials(
577 V * rhoB, V * rhoQ, V * rhoS, V * rhoC,
578 muBinit, muQinit, muSinit, muCinit,
579 ConstrMuB, ConstrMuQ, ConstrMuS, ConstrMuC);
580 }
581
582 bool ThermalModelBase::SolveChemicalPotentials(double totB, double totQ, double totS, double totC,
583 double muBinit, double muQinit, double muSinit, double muCinit,
584 bool ConstrMuB, bool ConstrMuQ, bool ConstrMuS, bool ConstrMuC) {
585 if (UsePartialChemicalEquilibrium()) {
586 std::cerr << "**WARNING** PCE enabled, cannot assume chemical equilibrium to do optimization..." << std::endl;
587 return false;
588 }
589
590 // TODO: Stability analysis for small densities/numbers
591 // if (abs(totB) < 1.e-6)
592 // totB = 0.;
593 // if (abs(totQ) < 1.e-6)
594 // totQ = 0.;
595 // if (abs(totS) < 1.e-6)
596 // totS = 0.;
597 // if (abs(totC) < 1.e-6)
598 // totC = 0.;
599
600 if (ConstrMuB)
601 m_Parameters.muB = muBinit;
602 if (ConstrMuS)
603 m_Parameters.muS = muSinit;
604 if (ConstrMuQ)
605 m_Parameters.muQ = muQinit;
606 if (ConstrMuC)
607 m_Parameters.muC = muCinit;
608 bool allzero = true;
609 allzero &= (totB == 0.0 && ConstrMuB) || (muBinit == 0 && !ConstrMuB);
610 allzero &= (totQ == 0.0 && ConstrMuQ) || (muQinit == 0 && !ConstrMuQ);
611 allzero &= (totS == 0.0 && ConstrMuS) || (muSinit == 0 && !ConstrMuS);
612 allzero &= (totC == 0.0 && ConstrMuC) || (muCinit == 0 && !ConstrMuC);
613
614
615 if (allzero) {
616 m_Parameters.muB = 0.;
617 m_Parameters.muS = 0.;
618 m_Parameters.muQ = 0.;
619 m_Parameters.muC = 0.;
620 FillChemicalPotentials();
621 CalculatePrimordialDensities();
622 return true;
623 }
624 vector<int> vConstr(4, 1);
625 vector<int> vType(4, 0);
626
627
628 vConstr[0] = m_TPS->hasBaryons() && ConstrMuB;
629 vConstr[1] = m_TPS->hasCharged() && ConstrMuQ;
630 vConstr[2] = m_TPS->hasStrange() && ConstrMuS;
631 vConstr[3] = m_TPS->hasCharmed() && ConstrMuC;
632
633 vType[0] = (int)(totB == 0.0);
634 vType[1] = (int)(totQ == 0.0);
635 vType[2] = (int)(totS == 0.0);
636 vType[3] = (int)(totC == 0.0);
637
638 vector<double> vTotals(4);
639 vTotals[0] = totB;
640 vTotals[1] = totQ;
641 vTotals[2] = totS;
642 vTotals[3] = totC;
643
644 vector<double> xin(4, 0.);
645 xin[0] = muBinit;
646 xin[1] = muQinit;
647 xin[2] = muSinit;
648 xin[3] = muCinit;
649
650 vector<double> xinactual;
651 for (int i = 0; i < 4; ++i) {
652 if (vConstr[i]) {
653 xinactual.push_back(xin[i]);
654 }
655 }
656
657 BroydenEquationsChemTotals eqs(vConstr, vType, vTotals, this);
658 BroydenJacobianChemTotals jaco(vConstr, vType, vTotals, this);
659 Broyden broydn(&eqs, &jaco);
660 Broyden::BroydenSolutionCriterium crit(1.0E-8);
661 broydn.Solve(xinactual, &crit);
662
663
664 //std::cerr << BaryonDensity() * Volume() << std::endl;
665
666 return (broydn.Iterations() < broydn.MaxIterations());
667 }
668
675
676 void ThermalModelBase::ValidateCalculation()
677 {
678 m_ValidityLog = "";
679
680 char cc[1000];
681
682 m_LastCalculationSuccessFlag = true;
683 for (size_t i = 0; i < m_densities.size(); ++i) {
684 if (m_densities[i] != m_densities[i]) {
685 m_LastCalculationSuccessFlag = false;
686
687 sprintf(cc, "Density for particle %d (%s) is NaN!\n", m_TPS->Particle(i).PdgId(), m_TPS->Particle(i).Name().c_str());
688 std::cerr << "**WARNING** " << cc << std::endl;
689
690
691 m_ValidityLog.append(cc);
692 }
693 //m_LastCalculationSuccessFlag &= (m_densities[i] == m_densities[i]);
694 }
695 }
696
697 std::vector<double> ThermalModelBase::CalculateChargeFluctuations(const std::vector<double>& /*chgs*/, int /*order*/, bool /*dimensionfull*/)
698 {
699 std::cerr << "**WARNING** " << m_TAG << "::CalculateChargeFluctuations(const std::vector<double>& chgs, int order) not implemented!" << std::endl;
700 return std::vector<double>();
701 }
702
703 std::vector<double> ThermalModelBase::CalculateGeneralizedSusceptibilities(const std::vector<std::vector<double>> &/*chgs*/)
704 {
705 throw std::runtime_error("ThermalModelBase::CalculateGeneralizedSusceptibilities: not implemented!");
706 }
707
708 double ThermalModelBase::CalculateHadronDensity() {
709 if (!m_Calculated) CalculateDensities();
710 double ret = 0.;
711 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
712 ret += m_densities[i];
713
714 return ret;
715 }
716
717 double ThermalModelBase::CalculateBaryonDensity() {
718 if (!m_Calculated) CalculateDensities();
719 double ret = 0.;
720 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
721 ret += m_TPS->Particles()[i].BaryonCharge() * m_densities[i];
722
723 return ret;
724 }
725
726 double ThermalModelBase::CalculateChargeDensity() {
727 if (!m_Calculated) CalculateDensities();
728 double ret = 0.;
729 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
730 ret += m_TPS->Particles()[i].ElectricCharge() * m_densities[i];
731
732 return ret;
733 }
734
735 double ThermalModelBase::CalculateStrangenessDensity() {
736 if (!m_Calculated) CalculateDensities();
737 double ret = 0.;
738 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
739 ret += m_TPS->Particles()[i].Strangeness() * m_densities[i];
740
741 return ret;
742 }
743
744 double ThermalModelBase::CalculateCharmDensity() {
745 if (!m_Calculated) CalculateDensities();
746 double ret = 0.;
747 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
748 ret += m_TPS->Particles()[i].Charm() * m_densities[i];
749 return ret;
750 }
751
752 double ThermalModelBase::CalculateAbsoluteBaryonDensity() {
753 if (!m_Calculated) CalculateDensities();
754 double ret = 0.;
755 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
756 ret += fabs((double)m_TPS->Particles()[i].BaryonCharge()) * m_densities[i];
757 return ret;
758 }
759
760 double ThermalModelBase::CalculateAbsoluteChargeDensity() {
761 if (!m_Calculated) CalculateDensities();
762 double ret = 0.;
763 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
764 ret += fabs((double)m_TPS->Particles()[i].ElectricCharge()) * m_densities[i];
765 return ret;
766 }
767
768 double ThermalModelBase::CalculateAbsoluteStrangenessDensity() {
769 if (!m_Calculated) CalculateDensities();
770 double ret = 0.;
771 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
772 ret += m_TPS->Particles()[i].AbsoluteStrangeness() * m_densities[i];
773 return ret;
774 }
775
776 double ThermalModelBase::CalculateAbsoluteStrangenessDensityModulo() {
777 if (!m_Calculated) CalculateDensities();
778 double ret = 0.;
779 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
780 ret += fabs((double)m_TPS->Particles()[i].Strangeness()) * m_densities[i];
781 return ret;
782 }
783
784 double ThermalModelBase::CalculateAbsoluteCharmDensity() {
785 if (!m_Calculated) CalculateDensities();
786 double ret = 0.;
787 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
788 ret += m_TPS->Particles()[i].AbsoluteCharm() * m_densities[i];
789 return ret;
790 }
791
792 double ThermalModelBase::CalculateAbsoluteCharmDensityModulo() {
793 if (!m_Calculated) CalculateDensities();
794 double ret = 0.;
795 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
796 ret += fabs((double)m_TPS->Particles()[i].Charm()) * m_densities[i];
797 return ret;
798 }
799
800 double ThermalModelBase::CalculateArbitraryChargeDensity() {
801 if (!m_Calculated) CalculateDensities();
802 double ret = 0.;
803 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i)
804 ret += m_TPS->Particles()[i].ArbitraryCharge() * m_densities[i];
805 return ret;
806 }
807
808
809 double ThermalModelBase::GetDensity(long long PDGID, const std::vector<double> *dens)
810 {
811 if (m_TPS->PdgToId(PDGID) != -1)
812 return dens->operator[](m_TPS->PdgToId(PDGID));
813
814 // 1 - Npart
815 if (PDGID == 1) return CalculateBaryonDensity();
816
817 // K0S or K0L
818 if (PDGID == 310 || PDGID == 130)
819 if (m_TPS->PdgToId(311) != -1 && m_TPS->PdgToId(-311) != -1)
820 return (dens->operator[](m_TPS->PdgToId(311)) + dens->operator[](m_TPS->PdgToId(-311))) / 2.;
821
822 // Id Pdg code has a trailing zero, try to construct a particle + anti-particle yield
823 if (PDGID % 10 == 0) {
824 long long tpdgid = PDGID / 10;
825 if (m_TPS->PdgToId(tpdgid) != -1 && m_TPS->PdgToId(-tpdgid) != -1)
826 return dens->operator[](m_TPS->PdgToId(tpdgid)) + dens->operator[](m_TPS->PdgToId(-tpdgid));
827 }
828
829 // 22122112 - nucleons
830 if (PDGID == 22122112 && m_TPS->PdgToId(2212) != -1 && m_TPS->PdgToId(2112) != -1)
831 return dens->operator[](m_TPS->PdgToId(2212)) + dens->operator[](m_TPS->PdgToId(2112));
832
833 std::cerr << "**WARNING** " << m_TAG << ": Density with PDG ID " << PDGID << " not found!" << std::endl;
834
835 return 0.;
836 }
837
838 double ThermalModelBase::GetDensity(long long PDGID, Feeddown::Type feeddown)
839 {
840 std::vector<double> *dens = NULL;
841 if (feeddown == Feeddown::Primordial)
842 dens = &m_densities;
843 else if (feeddown == Feeddown::StabilityFlag)
844 dens = &m_densitiestotal;
845 else if (static_cast<size_t>(feeddown) < m_densitiesbyfeeddown.size())
846 dens = &m_densitiesbyfeeddown[static_cast<int>(feeddown)];
847 else {
848 std::cerr << "**WARNING** " << m_TAG << ": GetDensity: Unknown feeddown: " << static_cast<int>(feeddown) << std::endl;
849 return 0.;
850 }
851
852 if (!m_Calculated)
853 CalculatePrimordialDensities();
854
855 if (feeddown != Feeddown::Primordial && !m_FeeddownCalculated)
856 CalculateFeeddown();
857
858 double ret = GetDensity(PDGID, dens);
859
860 // Weak decay contributions from K0S if this particle is not in the list
861 if (feeddown == Feeddown::Weak && m_TPS->PdgToId(310) == -1) {
862 // pi0
863 if (PDGID == 111) {
864 ret += 2. * 0.308 * GetDensity(310, dens);
865 }
866 // pi+,-
867 if (PDGID == 211 || PDGID == -211) {
868 ret += 0.692 * GetDensity(310, dens);
869 }
870 }
871
872 return ret;
873 }
874
875
876 std::vector<double> ThermalModelBase::GetIdealGasDensities() const {
877 std::vector<double> ret = m_densities;
878 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
879 ret[i] = m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::ParticleDensity, m_UseWidth, m_Chem[i]);
880 }
881 return ret;
882 }
883
884 void ThermalModelBase::ResetCalculatedFlags()
885 {
886 m_Calculated = false;
887 m_FeeddownCalculated = false;
888 m_TwoParticleCorrelationsCalculated = false;
889 m_SusceptibilitiesCalculated = false;
890 m_FluctuationsCalculated = false;
891 m_TemperatureDerivativesCalculated = false;
892 m_GCECalculated = false;
893 }
894
895 double ThermalModelBase::ConservedChargeDensity(ConservedCharge::Name chg)
896 {
897 if (chg == ConservedCharge::BaryonCharge)
898 return BaryonDensity();
899 if (chg == ConservedCharge::ElectricCharge)
900 return ElectricChargeDensity();
901 if (chg == ConservedCharge::StrangenessCharge)
902 return StrangenessDensity();
903 if (chg == ConservedCharge::CharmCharge)
904 return CharmDensity();
905 return 0.0;
906 }
907
908 double ThermalModelBase::AbsoluteConservedChargeDensity(ConservedCharge::Name chg)
909 {
910 if (chg == ConservedCharge::BaryonCharge)
911 return AbsoluteBaryonDensity();
912 if (chg == ConservedCharge::ElectricCharge)
913 return AbsoluteElectricChargeDensity();
914 if (chg == ConservedCharge::StrangenessCharge)
915 return AbsoluteStrangenessDensity();
916 if (chg == ConservedCharge::CharmCharge)
917 return AbsoluteCharmDensity();
918 return 0.0;
919 }
920
921 double ThermalModelBase::ConservedChargeDensitydT(ConservedCharge::Name chg)
922 {
923 if (!IsTemperatureDerivativesCalculated())
924 CalculateTemperatureDerivatives();
925
926 double ret = 0.0;
927 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
928 ret += m_TPS->Particles()[i].ConservedCharge(chg) * m_dndT[i];
929 }
930
931 return ret;
932 }
933
934 double ThermalModelBase::ChargedMultiplicity(int type)
935 {
936 if (!m_Calculated) CalculateDensities();
937 double ret = 0.0;
938 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
939 int tQ = m_TPS->Particles()[i].ElectricCharge();
940 bool fl = false;
941 if (type == 0 && tQ != 0)
942 fl = true;
943 if (type == 1 && tQ > 0)
944 fl = true;
945 if (type == -1 && tQ < 0)
946 fl = true;
947 if (fl)
948 ret += m_densities[i];
949 }
950 return ret * Volume();
951 }
952
953 double ThermalModelBase::ChargedScaledVariance(int type)
954 {
955 if (!m_FluctuationsCalculated) {
956 std::cerr << "**WARNING** " << m_TAG << ": ChargedScaledVariance(int): Fluctuations were not calculated\n";
957 return 1.;
958 }
959 double ret = 0.0;
960 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
961 int tQ = m_TPS->Particles()[i].ElectricCharge();
962 bool fl = false;
963 if (type == 0 && tQ != 0)
964 fl = true;
965 if (type == 1 && tQ > 0)
966 fl = true;
967 if (type == -1 && tQ < 0)
968 fl = true;
969 if (fl) {
970 for (int j = 0; j < m_TPS->ComponentsNumber(); ++j) {
971 int tQ2 = m_TPS->Particles()[j].ElectricCharge();
972 bool fl2 = false;
973 if (type == 0 && tQ2 != 0)
974 fl2 = true;
975 if (type == 1 && tQ2 > 0)
976 fl2 = true;
977 if (type == -1 && tQ2 < 0)
978 fl2 = true;
979
980 if (fl2) {
981 ret += m_PrimCorrel[i][j];
982 }
983 }
984 }
985 }
986 return ret * m_Parameters.T * Volume() / ChargedMultiplicity(type);
987 }
988
989 double ThermalModelBase::ChargedMultiplicityFinal(int type)
990 {
991 if (!m_Calculated) CalculateDensities();
992
993 int op = type;
994 if (type == -1)
995 op = 2;
996
997 double ret = 0.0;
998 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
999 ret += m_densities[i] * m_TPS->Particles()[i].Nch()[op];
1000 }
1001 return ret * Volume();
1002 }
1003
1004 double ThermalModelBase::ChargedScaledVarianceFinal(int type)
1005 {
1006 if (!m_FluctuationsCalculated) {
1007 std::cerr << "**WARNING** " << m_TAG << ": ChargedScaledVarianceFinal(int): Fluctuations were not calculated" << std::endl;
1008 return 1.;
1009 }
1010 int op = type;
1011 if (type == -1)
1012 op = 2;
1013 double ret = 0.0;
1014 for (int i = 0; i < m_TPS->ComponentsNumber(); ++i) {
1015 ret += m_densities[i] * Volume() * m_TPS->Particles()[i].DeltaNch()[op];
1016 for (int j = 0; j < m_TPS->ComponentsNumber(); ++j) {
1017 ret += m_PrimCorrel[i][j] * m_Parameters.T * Volume() * m_TPS->Particles()[i].Nch()[op] * m_TPS->Particles()[j].Nch()[op];
1018 }
1019 }
1020 return ret / ChargedMultiplicityFinal(type);
1021 }
1022
1023 void ThermalModelBase::CalculateTwoParticleCorrelations() {
1024 throw std::runtime_error("ThermalModelBase::CalculateTwoParticleCorrelations: Calculation of two-particle correlations and fluctuations is not implemented");
1025 }
1026
1027
1028 void ThermalModelBase::CalculateTwoParticleFluctuationsDecays()
1029 {
1030 // Decay contributions here are done according to Eq. (47) in nucl-th/0606036
1031
1032 int NN = m_densities.size();
1033
1034 // Fluctuations for all
1035 for (int i = 0; i < NN; ++i)
1036 //for(int j=0;j<NN;++j)
1037 {
1038 m_TotalCorrel[i][i] = m_PrimCorrel[i][i];
1039 //for (int r = 0; r < m_TPS->Particles()[i].DecayContributions().size(); ++r) {
1040 const ThermalParticleSystem::DecayContributionsToParticle& decayContributions = m_TPS->DecayContributionsByFeeddown()[Feeddown::StabilityFlag][i];
1041 for (size_t r = 0; r < decayContributions.size(); ++r) {
1042 int rr = decayContributions[r].second;
1043
1044 m_TotalCorrel[i][i] += 2. * m_PrimCorrel[i][rr] * decayContributions[r].first;
1045
1046 for (size_t r2 = 0; r2 < decayContributions.size(); ++r2) {
1047 int rr2 = decayContributions[r2].second;
1048 m_TotalCorrel[i][i] += m_PrimCorrel[rr][rr2] * decayContributions[r].first * decayContributions[r2].first;
1049 }
1050
1051 // Probabilistic decays
1052 m_TotalCorrel[i][i] += m_densities[rr] / m_Parameters.T * m_TPS->DecayCumulants()[i][r].first[1];
1053 //m_TotalCorrel[i][i] += m_densities[rr] / m_Parameters.T * m_TPS->Particles()[i].DecayCumulants()[r].first[1];
1054 //if (rr != i) { // && !m_TPS->Particles()[r].IsStable()) {
1055 // double nij = 0., ni = 0., nj = 0., dnij = 0.;
1056 // const auto& decayDistributions = TPS()->ResonanceFinalStatesDistributions()[rr];
1057 // for (size_t br = 0; br < decayDistributions.size(); ++br) {
1058 // nij += decayDistributions[br].first * decayDistributions[br].second[i] * decayDistributions[br].second[i];
1059 // ni += decayDistributions[br].first * decayDistributions[br].second[i];
1060 // nj += decayDistributions[br].first * decayDistributions[br].second[i];
1061 // }
1062 // dnij = nij - ni * nj;
1063 // m_TotalCorrel[i][i] += m_densities[rr] / Parameters().T * dnij;
1064 //}
1065 }
1066 }
1067
1068
1069 // Correlations only for stable
1070 for (int i = 0; i < NN; ++i) {
1071 if (m_TPS->Particles()[i].IsStable()) {
1072 for (int j = 0; j < NN; ++j) {
1073 if (j != i && m_TPS->Particles()[j].IsStable()) {
1074 m_TotalCorrel[i][j] = m_PrimCorrel[i][j];
1075
1076 const ThermalParticleSystem::DecayContributionsToParticle& decayContributionsI = m_TPS->DecayContributionsByFeeddown()[Feeddown::StabilityFlag][i];
1077 const ThermalParticleSystem::DecayContributionsToParticle& decayContributionsJ = m_TPS->DecayContributionsByFeeddown()[Feeddown::StabilityFlag][j];
1078
1079 for (size_t r = 0; r < decayContributionsJ.size(); ++r) {
1080 int rr = decayContributionsJ[r].second;
1081 m_TotalCorrel[i][j] += m_PrimCorrel[i][rr] * decayContributionsJ[r].first;
1082 }
1083
1084 for (size_t r = 0; r < decayContributionsI.size(); ++r) {
1085 int rr = decayContributionsI[r].second;
1086 m_TotalCorrel[i][j] += m_PrimCorrel[j][rr] * decayContributionsI[r].first;
1087 }
1088
1089 for (size_t r = 0; r < decayContributionsI.size(); ++r) {
1090 int rr = decayContributionsI[r].second;
1091
1092 for (size_t r2 = 0; r2 < decayContributionsJ.size(); ++r2) {
1093 int rr2 = decayContributionsJ[r2].second;
1094 m_TotalCorrel[i][j] += m_PrimCorrel[rr][rr2] * decayContributionsI[r].first * decayContributionsJ[r2].first;
1095 }
1096 }
1097
1098 // Contribution from probabilistic decays
1099 for (int r = 0; r < m_TPS->ComponentsNumber(); ++r) {
1100 if (r != i && r != j) { // && !m_TPS->Particles()[r].IsStable()) {
1101 double nij = 0., ni = 0., nj = 0., dnij = 0.;
1102 //const ThermalParticle &tpart = m_TPS->Particle(r);
1103 const ThermalParticleSystem::ResonanceFinalStatesDistribution &decayDistributions = m_TPS->ResonanceFinalStatesDistributions()[r];
1104 for (size_t br = 0; br < decayDistributions.size(); ++br) {
1105 nij += decayDistributions[br].first * decayDistributions[br].second[i] * decayDistributions[br].second[j];
1106 ni += decayDistributions[br].first * decayDistributions[br].second[i];
1107 nj += decayDistributions[br].first * decayDistributions[br].second[j];
1108 }
1109 dnij = nij - ni * nj;
1110 m_TotalCorrel[i][j] += m_densities[r] / m_Parameters.T * dnij;
1111 }
1112 }
1113 }
1114 }
1115 }
1116 }
1117
1118 for (int i = 0; i < NN; ++i) {
1119 m_wtot[i] = m_TotalCorrel[i][i];
1120 if (m_densitiestotal[i] > 0.) m_wtot[i] *= m_Parameters.T / m_densitiestotal[i];
1121 else m_wtot[i] = 1.;
1122 // Guard against NaN propagated from singular correlation matrices
1123 if (m_wtot[i] != m_wtot[i]) m_wtot[i] = 1.;
1124 }
1125 }
1126
1127 double ThermalModelBase::TwoParticleSusceptibilityPrimordial(int i, int j) const
1128 {
1129 if (!IsFluctuationsCalculated()) {
1130 throw std::runtime_error("ThermalModelBase::TwoParticleSusceptibilityPrimordial: fluctuations were not computed beforehand!");
1131 }
1132
1133 return m_PrimCorrel[i][j] / m_Parameters.T / m_Parameters.T / xMath::GeVtoifm() / xMath::GeVtoifm() / xMath::GeVtoifm();
1134 }
1135
1136 double ThermalModelBase::TwoParticleSusceptibilityPrimordialByPdg(long long id1, long long id2)
1137 {
1138 int i = TPS()->PdgToId(id1);
1139 int j = TPS()->PdgToId(id2);
1140
1141 if (i == -1) {
1142 std::cerr << "**WARNING** ThermalModelBase::TwoParticleSusceptibilityPrimordialByPdg: unknown pdg code " << id1 << std::endl;
1143 return 0.;
1144 }
1145 if (j == -1) {
1146 std::cerr << "**WARNING** ThermalModelBase::TwoParticleSusceptibilityPrimordialByPdg: unknown pdg code " << id2 << std::endl;
1147 return 0.;
1148 }
1149
1150 return TwoParticleSusceptibilityPrimordial(i, j);
1151 }
1152
1153 double ThermalModelBase::TwoParticleSusceptibilityTemperatureDerivativePrimordial(int i, int j) const {
1154 if (!IsFluctuationsCalculated() || !IsTemperatureDerivativesCalculated()) {
1155 throw std::runtime_error("ThermalModelBase::TwoParticleSusceptibilityPrimordial: temperature derivatives of fluctuations were not computed beforehand!");
1156 }
1157
1158 return m_PrimChi2sdT[i][j];
1159 }
1160
1161 double ThermalModelBase::TwoParticleSusceptibilityTemperatureDerivativePrimordialByPdg(long long id1, long long id2)
1162 {
1163 int i = TPS()->PdgToId(id1);
1164 int j = TPS()->PdgToId(id2);
1165
1166 if (i == -1) {
1167 std::cerr << "**WARNING** ThermalModelBase::TwoParticleSusceptibilityTemperatureDerivativePrimordialByPdg: unknown pdg code " << id1 << std::endl;
1168 return 0.;
1169 }
1170 if (j == -1) {
1171 std::cerr << "**WARNING** ThermalModelBase::TwoParticleSusceptibilityTemperatureDerivativePrimordialByPdg: unknown pdg code " << id2 << std::endl;
1172 return 0.;
1173 }
1174
1175 return TwoParticleSusceptibilityTemperatureDerivativePrimordial(i, j);
1176 }
1177
1178 double ThermalModelBase::NetParticleSusceptibilityPrimordialByPdg(long long id1, long long id2)
1179 {
1180 int i1 = TPS()->PdgToId(id1);
1181 int j1 = TPS()->PdgToId(id2);
1182
1183 if (i1 == -1) {
1184 std::cerr << "**WARNING** ThermalModelBase::NetParticleSusceptibilityPrimordialByPdg: unknown pdg code " << id1 << std::endl;
1185 return 0.;
1186 }
1187 if (j1 == -1) {
1188 std::cerr << "**WARNING** ThermalModelBase::NetParticleSusceptibilityPrimordialByPdg: unknown pdg code " << id2 << std::endl;
1189 return 0.;
1190 }
1191
1192 int i2 = TPS()->PdgToId(-id1);
1193 int j2 = TPS()->PdgToId(-id2);
1194
1195 // Both particles are neutral
1196 if (i2 == -1 && j2 == -1) {
1197 return TwoParticleSusceptibilityPrimordial(i1, j1);
1198 }
1199
1200 // First particle species is neutral
1201 if (i2 == -1) {
1202 return TwoParticleSusceptibilityPrimordial(i1, j1) - TwoParticleSusceptibilityPrimordial(i1, j2);
1203 }
1204
1205 // Second particle species is neutral
1206 if (j2 == -1) {
1207 return TwoParticleSusceptibilityPrimordial(i1, j1) - TwoParticleSusceptibilityPrimordial(i2, j1);
1208 }
1209
1210 // Both particles have antiparticles
1211 return TwoParticleSusceptibilityPrimordial(i1, j1) + TwoParticleSusceptibilityPrimordial(i2, j2) - TwoParticleSusceptibilityPrimordial(i1, j2) - TwoParticleSusceptibilityPrimordial(i2, j1);
1212 }
1213
1214 double ThermalModelBase::TwoParticleSusceptibilityFinal(int i, int j) const
1215 {
1216 if (!IsFluctuationsCalculated()) {
1217 throw std::runtime_error("ThermalModelBase::TwoParticleSusceptibilityFinal: fluctuations were not computed beforehand!");
1218 }
1219
1220 if (!m_TPS->Particle(i).IsStable() || !m_TPS->Particle(j).IsStable()) {
1221 int tid = i;
1222 if (!m_TPS->Particle(j).IsStable())
1223 tid = j;
1224 throw std::runtime_error("ThermalModelBase::TwoParticleSusceptibilityFinal: Particle " + std::to_string(tid) + " is not stable! Final correlations not computed for unstable particles.");
1225 }
1226
1227 return m_TotalCorrel[i][j] / m_Parameters.T / m_Parameters.T / xMath::GeVtoifm() / xMath::GeVtoifm() / xMath::GeVtoifm();
1228 }
1229
1230 double ThermalModelBase::TwoParticleSusceptibilityFinalByPdg(long long id1, long long id2)
1231 {
1232 int i = TPS()->PdgToId(id1);
1233 int j = TPS()->PdgToId(id2);
1234
1235 if (i == -1) {
1236 std::cerr << "**WARNING** ThermalModelBase::TwoParticleSusceptibilityFinalByPdg: unknown pdg code " << id1 << std::endl;
1237 return 0.;
1238 }
1239 if (j == -1) {
1240 std::cerr << "**WARNING** ThermalModelBase::TwoParticleSusceptibilityFinalByPdg: unknown pdg code " << id2 << std::endl;
1241 return 0.;
1242 }
1243
1244 return TwoParticleSusceptibilityFinal(i, j);
1245 }
1246
1247 double ThermalModelBase::NetParticleSusceptibilityFinalByPdg(long long id1, long long id2)
1248 {
1249 int i1 = TPS()->PdgToId(id1);
1250 int j1 = TPS()->PdgToId(id2);
1251
1252 if (i1 == -1) {
1253 std::cerr << "**WARNING** ThermalModelBase::NetParticleSusceptibilityFinalByPdg: unknown pdg code " << id1 << std::endl;
1254 return 0.;
1255 }
1256 if (j1 == -1) {
1257 std::cerr << "**WARNING** ThermalModelBase::NetParticleSusceptibilityFinalByPdg: unknown pdg code " << id2 << std::endl;
1258 return 0.;
1259 }
1260
1261 int i2 = TPS()->PdgToId(-id1);
1262 int j2 = TPS()->PdgToId(-id2);
1263
1264 // Both particles are neutral
1265 if (i2 == -1 && j2 == -1) {
1266 return TwoParticleSusceptibilityFinal(i1, j1);
1267 }
1268
1269 // First particle species is neutral
1270 if (i2 == -1) {
1271 return TwoParticleSusceptibilityFinal(i1, j1) - TwoParticleSusceptibilityFinal(i1, j2);
1272 }
1273
1274 // Second particle species is neutral
1275 if (j2 == -1) {
1276 return TwoParticleSusceptibilityFinal(i1, j1) - TwoParticleSusceptibilityFinal(i2, j1);
1277 }
1278
1279 // Both particles have antiparticles
1280 return TwoParticleSusceptibilityFinal(i1, j1) + TwoParticleSusceptibilityFinal(i2, j2) - TwoParticleSusceptibilityFinal(i1, j2) - TwoParticleSusceptibilityFinal(i2, j1);
1281 }
1282
1283 double ThermalModelBase::PrimordialParticleChargeSusceptibility(int i, ConservedCharge::Name chg) const
1284 {
1285 if (!IsFluctuationsCalculated()) {
1286 throw std::runtime_error("ThermalModelBase::PrimordialParticleChargeSusceptibility: fluctuations were not computed beforehand!");
1287 }
1288
1289 return m_PrimChargesCorrel[i][static_cast<int>(chg)] / m_Parameters.T / m_Parameters.T / xMath::GeVtoifm() / xMath::GeVtoifm() / xMath::GeVtoifm();
1290 }
1291
1292 double ThermalModelBase::PrimordialParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg)
1293 {
1294 int i = TPS()->PdgToId(id1);
1295 if (i == -1) {
1296 std::cerr << "**WARNING** ThermalModelBase::PrimordialParticleChargeSusceptibilityByPdg: unknown pdg code " << id1 << std::endl;
1297 return 0.;
1298 }
1299
1300 return PrimordialParticleChargeSusceptibility(i, chg);
1301 }
1302
1303 double ThermalModelBase::PrimordialNetParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg)
1304 {
1305 int i1 = TPS()->PdgToId(id1);
1306 if (i1 == -1) {
1307 std::cerr << "**WARNING** ThermalModelBase::PrimordialNetParticleChargeSusceptibilityByPdg: unknown pdg code " << id1 << std::endl;
1308 return 0.;
1309 }
1310
1311 int i2 = TPS()->PdgToId(-id1);
1312 if (i2 == -1)
1313 return PrimordialParticleChargeSusceptibility(i1, chg);
1314
1315 return PrimordialParticleChargeSusceptibility(i1, chg) - PrimordialParticleChargeSusceptibility(i2, chg);
1316 }
1317
1318 double ThermalModelBase::FinalParticleChargeSusceptibility(int i, ConservedCharge::Name chg) const
1319 {
1320 if (!IsFluctuationsCalculated()) {
1321 throw std::runtime_error("ThermalModelBase::FinalParticleChargeSusceptibility: fluctuations were not computed beforehand!");
1322 }
1323
1324 return m_FinalChargesCorrel[i][static_cast<int>(chg)] / m_Parameters.T / m_Parameters.T / xMath::GeVtoifm() / xMath::GeVtoifm() / xMath::GeVtoifm();
1325 }
1326
1327 double ThermalModelBase::FinalParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg)
1328 {
1329 int i = TPS()->PdgToId(id1);
1330 if (i == -1) {
1331 std::cerr << "**WARNING** ThermalModelBase::FinalParticleChargeSusceptibilityByPdg: unknown pdg code " << id1 << std::endl;
1332 return 0.;
1333 }
1334
1335 return FinalParticleChargeSusceptibility(i, chg);
1336 }
1337
1338 double ThermalModelBase::FinalNetParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg)
1339 {
1340 int i1 = TPS()->PdgToId(id1);
1341 if (i1 == -1) {
1342 std::cerr << "**WARNING** ThermalModelBase::FinalNetParticleChargeSusceptibilityByPdg: unknown pdg code " << id1 << std::endl;
1343 return 0.;
1344 }
1345
1346 int i2 = TPS()->PdgToId(-id1);
1347 if (i2 == -1)
1348 return FinalParticleChargeSusceptibility(i1, chg);
1349
1350 return FinalParticleChargeSusceptibility(i1, chg) - FinalParticleChargeSusceptibility(i2, chg);
1351 }
1352
1353 void ThermalModelBase::CalculateSusceptibilityMatrix()
1354 {
1355 // Ensure two-particle correlations are calculated first
1356 if (!m_TwoParticleCorrelationsCalculated)
1357 CalculateTwoParticleCorrelations();
1358
1359 m_Susc.resize(4);
1360 for (int i = 0; i < 4; ++i) m_Susc[i].resize(4);
1361 m_dSuscdT = m_Susc;
1362
1363 for (int i = 0; i < 4; ++i) {
1364 for (int j = 0; j < 4; ++j) {
1365 m_Susc[i][j] = 0.;
1366 m_dSuscdT[i][j] = 0.;
1367 for (size_t k = 0; k < m_PrimCorrel.size(); ++k) {
1368 int c1 = 0;
1369 if (i == 0) c1 = m_TPS->Particles()[k].BaryonCharge();
1370 if (i == 1) c1 = m_TPS->Particles()[k].ElectricCharge();
1371 if (i == 2) c1 = m_TPS->Particles()[k].Strangeness();
1372 if (i == 3) c1 = m_TPS->Particles()[k].Charm();
1373 for (size_t kp = 0; kp < m_PrimCorrel.size(); ++kp) {
1374 int c2 = 0;
1375 if (j == 0) c2 = m_TPS->Particles()[kp].BaryonCharge();
1376 if (j == 1) c2 = m_TPS->Particles()[kp].ElectricCharge();
1377 if (j == 2) c2 = m_TPS->Particles()[kp].Strangeness();
1378 if (j == 3) c2 = m_TPS->Particles()[kp].Charm();
1379 m_Susc[i][j] += c1 * c2 * m_PrimCorrel[k][kp];
1380
1381 if (IsTemperatureDerivativesCalculated())
1382 m_dSuscdT[i][j] += c1 * c2 * m_PrimChi2sdT[k][kp];
1383 }
1384 }
1385 m_Susc[i][j] = m_Susc[i][j] / m_Parameters.T / m_Parameters.T / xMath::GeVtoifm() / xMath::GeVtoifm() / xMath::GeVtoifm();
1386 }
1387 }
1388 m_SusceptibilitiesCalculated = true;
1389 }
1390
1391 void ThermalModelBase::CalculateProxySusceptibilityMatrix()
1392 {
1393 m_ProxySusc.resize(4);
1394 for (int i = 0; i < 4; ++i)
1395 m_ProxySusc[i].resize(4);
1396
1397 // Up to 3, no charm here yet
1398 for (int i = 0; i < 3; ++i) {
1399 for (int j = 0; j < 3; ++j) {
1400 m_ProxySusc[i][j] = 0.;
1401 for (size_t k = 0; k < m_TotalCorrel.size(); ++k) {
1402 if (m_TPS->Particles()[k].IsStable()) {
1403 int c1 = 0;
1404 //if (i == 0) c1 = m_TPS->Particles()[k].BaryonCharge();
1405 if (i == 0) c1 = 1 * (m_TPS->Particles()[k].PdgId() == 2212) - 1 * (m_TPS->Particles()[k].PdgId() == -2212);
1406 if (i == 1) c1 = m_TPS->Particles()[k].ElectricCharge();
1407 //if (i == 1) c1 = 1 * (m_TPS->Particles()[k].PdgId() == 211) - 1 * (m_TPS->Particles()[k].PdgId() == -211);
1408 if (i == 2) c1 = 1 * (m_TPS->Particles()[k].PdgId() == 321) - 1 * (m_TPS->Particles()[k].PdgId() == -321);
1409 if (i == 3) c1 = m_TPS->Particles()[k].Charm();
1410 for (size_t kp = 0; kp < m_TotalCorrel.size(); ++kp) {
1411 if (m_TPS->Particles()[kp].IsStable()) {
1412 int c2 = 0;
1413 //if (j == 0) c2 = m_TPS->Particles()[kp].BaryonCharge();
1414 if (j == 0) c2 = 1 * (m_TPS->Particles()[kp].PdgId() == 2212) - 1 * (m_TPS->Particles()[kp].PdgId() == -2212);
1415 if (j == 1) c2 = m_TPS->Particles()[kp].ElectricCharge();
1416 //if (j == 1) c2 = 1 * (m_TPS->Particles()[kp].PdgId() == 211) - 1 * (m_TPS->Particles()[kp].PdgId() == -211);
1417 if (j == 2) c2 = 1 * (m_TPS->Particles()[kp].PdgId() == 321) - 1 * (m_TPS->Particles()[kp].PdgId() == -321);
1418 if (j == 3) c2 = m_TPS->Particles()[kp].Charm();
1419 m_ProxySusc[i][j] += c1 * c2 * m_TotalCorrel[k][kp];
1420 }
1421 }
1422 }
1423 }
1424 m_ProxySusc[i][j] = m_ProxySusc[i][j] / m_Parameters.T / m_Parameters.T / xMath::GeVtoifm() / xMath::GeVtoifm() / xMath::GeVtoifm();
1425 }
1426 }
1427 //std::cerr << "chi2netp/chi2skellam = " << m_ProxySusc[0][0] / (m_densitiestotal[m_TPS->PdgToId(2212)] + m_densitiestotal[m_TPS->PdgToId(-2212)]) * pow(m_Parameters.T * xMath::GeVtoifm(), 3) << std::endl;
1428 //std::cerr << "chi2netpi/chi2skellam = " << m_ProxySusc[1][1] / (m_densitiestotal[m_TPS->PdgToId(211)] + m_densitiestotal[m_TPS->PdgToId(-211)]) * pow(m_Parameters.T * xMath::GeVtoifm(), 3) << std::endl;
1429 }
1430
1431 void ThermalModelBase::CalculateParticleChargeCorrelationMatrix()
1432 {
1433 m_PrimChargesCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(4, 0.));
1434 m_FinalChargesCorrel = std::vector< std::vector<double> >(TPS()->ComponentsNumber(), std::vector<double>(4, 0.));
1435
1436 for (size_t i = 0; i < TPS()->ComponentsNumber(); ++i) {
1437 int c1 = 1;
1438 for (int chg = 0; chg < 4; ++chg) {
1439 for (size_t j = 0; j < TPS()->ComponentsNumber(); ++j) {
1440 int c2 = TPS()->Particle(j).ConservedCharge((ConservedCharge::Name)chg);
1441 m_PrimChargesCorrel[i][chg] += c1 * c2 * m_PrimCorrel[i][j];
1442 }
1443 }
1444 }
1445
1446 for (size_t i = 0; i < TPS()->ComponentsNumber(); ++i) {
1447 if (m_TPS->Particles()[i].IsStable()) {
1448 int c1 = 1;
1449 for (int chg = 0; chg < 4; ++chg) {
1450 for (size_t j = 0; j < TPS()->ComponentsNumber(); ++j) {
1451 if (m_TPS->Particles()[j].IsStable()) {
1452 int c2 = TPS()->Particle(j).ConservedCharge((ConservedCharge::Name)chg);
1453 m_FinalChargesCorrel[i][chg] += c1 * c2 * m_TotalCorrel[i][j];
1454 }
1455 }
1456 }
1457 }
1458 }
1459 }
1460
1461 std::vector<double> ThermalModelBase::PartialPressures() {
1462 if (!IsCalculated())
1463 CalculatePrimordialDensities();
1464
1465 std::vector<double> ret(5, 0.);
1466 for (size_t i = 0; i < TPS()->ComponentsNumber(); ++i) {
1467 const ThermalParticle& tpart = TPS()->Particle(i);
1468 int ind = 0;
1469 if (tpart.BaryonCharge() == 1)
1470 ind = 1;
1471 if (tpart.BaryonCharge() == -1)
1472 ind = 2;
1473 if (tpart.BaryonCharge() > 1)
1474 ind = 3;
1475 if (tpart.BaryonCharge() < -1)
1476 ind = 4;
1477 ret[ind] += m_TPS->Particles()[i].Density(m_Parameters, IdealGasFunctions::Pressure, m_UseWidth, FullIdealChemicalPotential(i));
1478 }
1479
1480 return ret;
1481 }
1482
1483 void ThermalModelBase::CalculateFluctuations() {
1484 std::cerr << "**WARNING** " << m_TAG << ": Calculation of fluctuations is not implemented" << std::endl;
1485 }
1486
1487 void ThermalModelBase::CalculateTemperatureDerivatives() {
1488 std::cerr << "**WARNING** " << m_TAG << ": Calculation of temperature derivatives is not implemented" << std::endl;
1489 }
1490
1491 std::vector<double> ThermalModelBase::BroydenEquationsChem::Equations(const std::vector<double>& x)
1492 {
1493 std::vector<double> ret(m_N, 0.);
1494
1495 int i1 = 0;
1496 if (m_THM->ConstrainMuB()) { m_THM->SetBaryonChemicalPotential(x[i1]); i1++; }
1497 if (m_THM->ConstrainMuQ()) { m_THM->SetElectricChemicalPotential(x[i1]); i1++; }
1498 if (m_THM->ConstrainMuS()) { m_THM->SetStrangenessChemicalPotential(x[i1]); i1++; }
1499 if (m_THM->ConstrainMuC()) { m_THM->SetCharmChemicalPotential(x[i1]); i1++; }
1500 m_THM->FillChemicalPotentials();
1501 m_THM->CalculatePrimordialDensities();
1502
1503 i1 = 0;
1504
1505 // Baryon charge
1506 if (m_THM->ConstrainMuB()) {
1507 double fBd = m_THM->CalculateBaryonDensity();
1508 double fSd = m_THM->CalculateEntropyDensity();
1509
1510 ret[i1] = (fBd / fSd - 1. / m_THM->SoverB()) * m_THM->SoverB();
1511
1512 i1++;
1513 }
1514
1515 // Electric charge
1516 if (m_THM->ConstrainMuQ()) {
1517 double fBd = m_THM->CalculateBaryonDensity();
1518 double fQd = m_THM->CalculateChargeDensity();
1519
1520 // Update: remove division by Q/B to allow for charge neutrality
1521 ret[i1] = (fQd / fBd - m_THM->QoverB());
1522 //ret[i1] = (fQd / fBd - m_THM->QoverB()) / m_THM->QoverB();
1523
1524 i1++;
1525 }
1526
1527
1528 // Strangeness
1529 if (m_THM->ConstrainMuS()) {
1530 double fSd = m_THM->CalculateStrangenessDensity();
1531 double fASd = m_THM->CalculateAbsoluteStrangenessDensity();
1532 if (fASd < 1.e-25)
1533 fASd = m_THM->CalculateAbsoluteStrangenessDensityModulo();
1534
1535 ret[i1] = fSd / fASd;
1536
1537 i1++;
1538 }
1539
1540
1541 // Charm
1542 if (m_THM->ConstrainMuC()) {
1543 double fCd = m_THM->CalculateCharmDensity();
1544 double fACd = m_THM->CalculateAbsoluteCharmDensity();
1545 if (fACd < 1.e-25)
1546 fACd = m_THM->CalculateAbsoluteCharmDensityModulo();
1547
1548 ret[i1] = fCd / fACd;
1549
1550 i1++;
1551 }
1552
1553 return ret;
1554 }
1555
1556 std::vector<double> ThermalModelBase::BroydenJacobianChem::Jacobian(const std::vector<double>& x)
1557 {
1558 int i1 = 0;
1559 // Analytic calculations of Jacobian not yet supported if entropy per baryon is involved
1560 if (m_THM->ConstrainMuB()) {
1561 throw std::runtime_error("ThermalModelBase::ConstrainChemicalPotentials: analytic calculation of the Jacobian not supported if muB is constrained");
1562 }
1563
1564 if (m_THM->ConstrainMuQ()) { m_THM->SetElectricChemicalPotential(x[i1]); i1++; }
1565 if (m_THM->ConstrainMuS()) { m_THM->SetStrangenessChemicalPotential(x[i1]); i1++; }
1566 if (m_THM->ConstrainMuC()) { m_THM->SetCharmChemicalPotential(x[i1]); i1++; }
1567 m_THM->FillChemicalPotentials();
1568 m_THM->CalculatePrimordialDensities();
1569
1570 double fBd = m_THM->CalculateBaryonDensity();
1571 double fQd = m_THM->CalculateChargeDensity();
1572 double fSd = m_THM->CalculateStrangenessDensity();
1573 double fASd = m_THM->CalculateAbsoluteStrangenessDensity();
1574 if (fASd < 1.e-25) {
1575 fASd = m_THM->CalculateAbsoluteStrangenessDensityModulo();
1576 }
1577 double fCd = m_THM->CalculateCharmDensity();
1578 double fACd = m_THM->CalculateAbsoluteCharmDensity();
1579 if (fACd < 1.e-25) {
1580 fACd = m_THM->CalculateAbsoluteCharmDensityModulo();
1581 }
1582
1583 vector<double> chi2s;
1584 chi2s.resize(m_THM->Densities().size());
1585 for (size_t i = 0; i < chi2s.size(); ++i) {
1586 chi2s[i] = m_THM->TPS()->Particle(i).chiDimensionfull(2, m_THM->Parameters(), m_THM->UseWidth(), m_THM->ChemicalPotential(i) + m_THM->MuShift(i))
1587 * xMath::GeVtoifm3();
1588 }
1589
1590 int NNN = 0;
1591 if (m_THM->ConstrainMuQ()) NNN++;
1592 if (m_THM->ConstrainMuS()) NNN++;
1593 if (m_THM->ConstrainMuC()) NNN++;
1594 MatrixXd ret(NNN, NNN);
1595
1596 i1 = 0;
1597 // Electric charge derivatives
1598 if (m_THM->ConstrainMuQ()) {
1599 int i2 = 0;
1600
1601 double d1 = 0., d2 = 0.;
1602
1603 if (m_THM->ConstrainMuQ()) {
1604 d1 = 0.;
1605 for (size_t i = 0; i < chi2s.size(); ++i)
1606 d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i];
1607
1608 d2 = 0.;
1609 for (size_t i = 0; i < chi2s.size(); ++i)
1610 d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i];
1611
1612 // Update: remove division by Q/B to allow for charge neutrality
1613 ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2);
1614 //ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2) / m_THM->QoverB();
1615
1616 i2++;
1617 }
1618
1619
1620 if (m_THM->ConstrainMuS()) {
1621 d1 = 0.;
1622 for (size_t i = 0; i < chi2s.size(); ++i)
1623 d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i];
1624
1625 d2 = 0.;
1626 for (size_t i = 0; i < chi2s.size(); ++i)
1627 d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i];
1628
1629 // Update: remove division by Q/B to allow for charge neutrality
1630 ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2);
1631 //ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2) / m_THM->QoverB();
1632
1633 i2++;
1634 }
1635
1636
1637 if (m_THM->ConstrainMuC()) {
1638 d1 = 0.;
1639 for (size_t i = 0; i < chi2s.size(); ++i)
1640 d1 += m_THM->TPS()->Particle(i).ElectricCharge() * m_THM->TPS()->Particle(i).Charm() * chi2s[i];
1641
1642 d2 = 0.;
1643 for (size_t i = 0; i < chi2s.size(); ++i)
1644 d2 += m_THM->TPS()->Particle(i).BaryonCharge() * m_THM->TPS()->Particle(i).Charm() * chi2s[i];
1645
1646 // Update: remove division by Q/B to allow for charge neutrality
1647 ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2);
1648 //ret(i1, i2) = (d1 / fBd - fQd / fBd / fBd * d2) / m_THM->QoverB();
1649
1650 i2++;
1651 }
1652
1653 i1++;
1654 }
1655
1656
1657 // Strangeness derivatives
1658 if (m_THM->ConstrainMuS()) {
1659 int i2 = 0;
1660
1661 double d1 = 0., d2 = 0.;
1662
1663 if (m_THM->ConstrainMuQ()) {
1664 d1 = 0.;
1665 for (size_t i = 0; i < chi2s.size(); ++i)
1666 d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i];
1667
1668 d2 = 0.;
1669 for (size_t i = 0; i < chi2s.size(); ++i)
1670 d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i];
1671
1672 ret(i1, i2) = d1 / fASd - fSd / fASd / fASd * d2;
1673
1674 i2++;
1675 }
1676
1677
1678 if (m_THM->ConstrainMuS()) {
1679 d1 = 0.;
1680 for (size_t i = 0; i < chi2s.size(); ++i)
1681 d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i];
1682
1683 d2 = 0.;
1684 for (size_t i = 0; i < chi2s.size(); ++i)
1685 d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i];
1686
1687 ret(i1, i2) = d1 / fASd - fSd / fASd / fASd * d2;
1688
1689 i2++;
1690 }
1691
1692
1693 if (m_THM->ConstrainMuC()) {
1694 d1 = 0.;
1695 for (size_t i = 0; i < chi2s.size(); ++i)
1696 d1 += m_THM->TPS()->Particle(i).Strangeness() * m_THM->TPS()->Particle(i).Charm() * chi2s[i];
1697
1698 d2 = 0.;
1699 for (size_t i = 0; i < chi2s.size(); ++i)
1700 d2 += m_THM->TPS()->Particle(i).AbsoluteStrangeness() * m_THM->TPS()->Particle(i).Charm() * chi2s[i];
1701
1702 ret(i1, i2) = d1 / fASd - fSd / fASd / fASd * d2;
1703
1704 i2++;
1705 }
1706
1707 i1++;
1708 }
1709
1710
1711 // Charm derivatives
1712 if (m_THM->ConstrainMuC()) {
1713 int i2 = 0;
1714
1715 double d1 = 0., d2 = 0.;
1716
1717 if (m_THM->ConstrainMuQ()) {
1718 d1 = 0.;
1719 for (size_t i = 0; i < chi2s.size(); ++i)
1720 d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).ElectricCharge() * chi2s[i];
1721
1722 d2 = 0.;
1723 for (size_t i = 0; i < chi2s.size(); ++i)
1724 d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).ElectricCharge() *chi2s[i];
1725
1726 ret(i1, i2) = d1 / fACd - fCd / fACd / fACd * d2;
1727
1728 i2++;
1729 }
1730
1731
1732 if (m_THM->ConstrainMuS()) {
1733 d1 = 0.;
1734 for (size_t i = 0; i < chi2s.size(); ++i)
1735 d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i];
1736
1737 d2 = 0.;
1738 for (size_t i = 0; i < chi2s.size(); ++i)
1739 d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).Strangeness() * chi2s[i];
1740
1741 ret(i1, i2) = d1 / fACd - fCd / fACd / fACd * d2;
1742
1743 i2++;
1744 }
1745
1746
1747 if (m_THM->ConstrainMuC()) {
1748 d1 = 0.;
1749 for (size_t i = 0; i < chi2s.size(); ++i)
1750 d1 += m_THM->TPS()->Particle(i).Charm() * m_THM->TPS()->Particle(i).Charm() * chi2s[i];
1751
1752 d2 = 0.;
1753 for (size_t i = 0; i < chi2s.size(); ++i)
1754 d2 += m_THM->TPS()->Particle(i).AbsoluteCharm() * m_THM->TPS()->Particle(i).Charm() * chi2s[i];
1755
1756 ret(i1, i2) = d1 / fACd - fCd / fACd / fACd * d2;
1757
1758 i2++;
1759 }
1760
1761 i1++;
1762 }
1763
1764 std::vector<double> retVec(NNN*NNN, 0);
1765 for (int i = 0; i < NNN; ++i)
1766 for (int j = 0; j < NNN; ++j)
1767 retVec[i*NNN + j] = ret(i, j);
1768
1769
1770 return retVec;
1771 }
1772
1773 std::vector<double> ThermalModelBase::BroydenChem::Solve(const std::vector<double>& x0, BroydenSolutionCriterium * solcrit, int max_iterations)
1774 {
1775 if (m_Equations == NULL) {
1776 throw std::runtime_error("Broyden::Solve: Equations to solve not specified!");
1777 }
1778
1779 int NNN = 0;
1780 std::vector<double> xcur;
1781 if (m_THM->ConstrainMuB()) { xcur.push_back(x0[0]); NNN++; }
1782 if (m_THM->ConstrainMuQ()) { xcur.push_back(x0[1]); NNN++; }
1783 if (m_THM->ConstrainMuS()) { xcur.push_back(x0[2]); NNN++; }
1784 if (m_THM->ConstrainMuC()) { xcur.push_back(x0[3]); NNN++; }
1785 if (NNN == 0) {
1786 m_THM->FillChemicalPotentials();
1787 return xcur;
1788 }
1789
1790 m_Equations->SetDimension(NNN);
1791
1792 m_MaxIterations = max_iterations;
1793
1794 BroydenSolutionCriterium *SolutionCriterium = solcrit;
1795 bool UseDefaultSolutionCriterium = false;
1796 if (SolutionCriterium == NULL) {
1797 SolutionCriterium = new BroydenSolutionCriterium(TOL);
1798 UseDefaultSolutionCriterium = true;
1799 }
1800
1801 BroydenJacobian *JacobianInUse = m_Jacobian;
1802 bool UseDefaultJacobian = false;
1803 if (JacobianInUse == NULL || m_THM->ConstrainMuB()) {
1804 JacobianInUse = new BroydenJacobian(m_Equations);
1805 UseDefaultJacobian = true;
1806 }
1807 m_Iterations = 0;
1808 double &maxdiff = m_MaxDifference;
1809 int N = m_Equations->Dimension();
1810
1811
1812
1813 std::vector<double> tmpvec, xdeltavec = xcur;
1814 VectorXd xold(N), xnew(N), xdelta(N);
1815 VectorXd fold(N), fnew(N), fdelta(N);
1816
1817 xold = VectorXd::Map(&xcur[0], xcur.size());
1818
1819 MatrixXd Jac = Eigen::Map< Matrix<double, Dynamic, Dynamic, RowMajor> >(&JacobianInUse->Jacobian(xcur)[0], N, N);
1820
1821 bool constrmuB = m_THM->ConstrainMuB();
1822 bool constrmuQ = m_THM->ConstrainMuQ();
1823 bool constrmuS = m_THM->ConstrainMuS();
1824 bool constrmuC = m_THM->ConstrainMuC();
1825 bool repeat = false;
1826 NNN = 0;
1827 if (m_THM->ConstrainMuB()) {
1828 for (int j = 0; j < Jac.rows(); ++j)
1829 if (Jac(NNN, j) > 1.e8) { repeat = true; m_THM->ConstrainMuB(false); }
1830 double S = m_THM->CalculateEntropyDensity();
1831 double nB = m_THM->CalculateBaryonDensity();
1832 if (abs(S) < 1.e-25 || abs(nB) < 1.e-25) { repeat = true; m_THM->ConstrainMuB(false); }
1833 NNN++;
1834 }
1835 if (m_THM->ConstrainMuQ()) {
1836 for (int j = 0; j < Jac.rows(); ++j)
1837 if (Jac(NNN, j) > 1.e8) { repeat = true; m_THM->ConstrainMuQ(false); }
1838 double nQ = m_THM->CalculateChargeDensity();
1839 double nB = m_THM->CalculateBaryonDensity();
1840 if (abs(nQ) < 1.e-25 || abs(nB) < 1.e-25) { repeat = true; m_THM->ConstrainMuQ(false); }
1841 NNN++;
1842 }
1843 if (m_THM->ConstrainMuS()) {
1844 for (int j = 0; j < Jac.rows(); ++j)
1845 if (Jac(NNN, j) > 1.e8) { repeat = true; m_THM->ConstrainMuS(false); }
1846
1847 double nS = m_THM->CalculateAbsoluteStrangenessDensity();
1848 if (abs(nS) < 1.e-25) { nS = m_THM->CalculateAbsoluteStrangenessDensityModulo(); }
1849 if (abs(nS) < 1.e-25) { repeat = true; m_THM->ConstrainMuS(false); }
1850 NNN++;
1851 }
1852 if (m_THM->ConstrainMuC()) {
1853 for (int j = 0; j < Jac.rows(); ++j)
1854 if (Jac(NNN, j) > 1.e8) { repeat = true; m_THM->ConstrainMuC(false); }
1855 double nC = m_THM->CalculateAbsoluteCharmDensity();
1856 if (abs(nC) < 1.e-25) { nC = m_THM->CalculateAbsoluteCharmDensityModulo(); }
1857 if (abs(nC) < 1.e-25) { repeat = true; m_THM->ConstrainMuC(false); }
1858 NNN++;
1859 }
1860 if (repeat) {
1861 std::vector<double> ret = Solve(x0, solcrit, max_iterations);
1862 m_THM->ConstrainMuQ(constrmuB);
1863 m_THM->ConstrainMuQ(constrmuQ);
1864 m_THM->ConstrainMuS(constrmuS);
1865 m_THM->ConstrainMuC(constrmuC);
1866 return ret;
1867 }
1868
1869 if (Jac.determinant() == 0.0)
1870 {
1871 std::cerr << "**WARNING** Singular Jacobian in Broyden::Solve" << std::endl;
1872 return xcur;
1873 }
1874
1875 MatrixXd Jinv = Jac.inverse();
1876 tmpvec = m_Equations->Equations(xcur);
1877 fold = VectorXd::Map(&tmpvec[0], tmpvec.size());
1878
1879 for (m_Iterations = 1; m_Iterations < max_iterations; ++m_Iterations) {
1880 xnew = xold - Jinv * fold;
1881
1882 VectorXd::Map(&xcur[0], xcur.size()) = xnew;
1883
1884 tmpvec = m_Equations->Equations(xcur);
1885 fnew = VectorXd::Map(&tmpvec[0], tmpvec.size());
1886
1887
1888 maxdiff = 0.;
1889 for (size_t i = 0; i < xcur.size(); ++i) {
1890 maxdiff = std::max(maxdiff, fabs(fnew[i]));
1891 }
1892
1893 xdelta = xnew - xold;
1894 fdelta = fnew - fold;
1895
1896 VectorXd::Map(&xdeltavec[0], xdeltavec.size()) = xdelta;
1897
1898 if (SolutionCriterium->IsSolved(xcur, tmpvec, xdeltavec))
1899 break;
1900
1901 if (!m_UseNewton) // Use Broyden's method
1902 {
1903
1904 double norm = 0.;
1905 for (int i = 0; i < N; ++i)
1906 for (int j = 0; j < N; ++j)
1907 norm += xdelta[i] * Jinv(i, j) * fdelta[j];
1908 VectorXd p1(N);
1909 p1 = (xdelta - Jinv * fdelta);
1910 for (int i = 0; i < N; ++i) p1[i] *= 1. / norm;
1911 Jinv = Jinv + (p1 * xdelta.transpose()) * Jinv;
1912 }
1913 else // Use Newton's method
1914 {
1915 Jac = Eigen::Map< Matrix<double, Dynamic, Dynamic, RowMajor> >(&JacobianInUse->Jacobian(xcur)[0], N, N);
1916 Jinv = Jac.inverse();
1917 }
1918
1919 xold = xnew;
1920 fold = fnew;
1921 }
1922
1923 if (m_Iterations == max_iterations) {
1924 std::cerr << "**WARNING** Reached maximum number of iterations in Broyden procedure" << std::endl;
1925 }
1926
1927 if (UseDefaultSolutionCriterium) {
1928 delete SolutionCriterium;
1929 SolutionCriterium = NULL;
1930 }
1931 if (UseDefaultJacobian) {
1932 delete JacobianInUse;
1933 JacobianInUse = NULL;
1934 }
1935 return xcur;
1936 }
1937
1938
1939 ThermalModelBase::BroydenEquationsChemTotals::BroydenEquationsChemTotals(const std::vector<int>& vConstr, const std::vector<int>& vType, const std::vector<double>& vTotals, ThermalModelBase * model) :
1940 BroydenEquations(), m_Constr(vConstr), m_Type(vType), m_Totals(vTotals), m_THM(model)
1941 {
1942 m_N = 0;
1943 for (size_t i = 0; i < m_Constr.size(); ++i)
1944 m_N += m_Constr[i];
1945 }
1946
1947 std::vector<double> ThermalModelBase::BroydenEquationsChemTotals::Equations(const std::vector<double>& x)
1948 {
1949 std::vector<double> ret(m_N, 0.);
1950
1951 int i1 = 0;
1952 for (int i = 0; i < 4; ++i) {
1953 if (m_Constr[i]) {
1954 if (i == 0) m_THM->SetBaryonChemicalPotential(x[i1]);
1955 if (i == 1) m_THM->SetElectricChemicalPotential(x[i1]);
1956 if (i == 2) m_THM->SetStrangenessChemicalPotential(x[i1]);
1957 if (i == 3) m_THM->SetCharmChemicalPotential(x[i1]);
1958 i1++;
1959 }
1960 }
1961 m_THM->FillChemicalPotentials();
1962 m_THM->CalculatePrimordialDensities();
1963
1964 vector<double> dens(4, 0.), absdens(4, 0.);
1965 if (m_Constr[0]) {
1966 dens[0] = m_THM->CalculateBaryonDensity();
1967 absdens[0] = m_THM->CalculateAbsoluteBaryonDensity();
1968 }
1969 if (m_Constr[1]) {
1970 dens[1] = m_THM->CalculateChargeDensity();
1971 absdens[1] = m_THM->CalculateAbsoluteChargeDensity();
1972 }
1973 if (m_Constr[2]) {
1974 dens[2] = m_THM->CalculateStrangenessDensity();
1975 absdens[2] = m_THM->CalculateAbsoluteStrangenessDensity();
1976 if (absdens[2] < 1.e-25) {
1977 absdens[2] = m_THM->CalculateAbsoluteStrangenessDensityModulo();
1978 }
1979 }
1980 if (m_Constr[3]) {
1981 dens[3] = m_THM->CalculateCharmDensity();
1982 absdens[3] = m_THM->CalculateAbsoluteCharmDensity();
1983 if (absdens[3] < 1.e-25) {
1984 absdens[3] = m_THM->CalculateAbsoluteCharmDensityModulo();
1985 }
1986 }
1987
1988 i1 = 0;
1989
1990 for (int i = 0; i < 4; ++i) {
1991 if (m_Constr[i]) {
1992 if (m_Type[i] == 0)
1993 ret[i1] = (dens[i] * m_THM->Parameters().V - m_Totals[i]) / m_Totals[i];
1994 else
1995 ret[i1] = dens[i] / absdens[i];
1996 i1++;
1997 }
1998 }
1999
2000 return ret;
2001 }
2002
2003 std::vector<double> ThermalModelBase::BroydenJacobianChemTotals::Jacobian(const std::vector<double>& x)
2004 {
2005 int NNN = 0;
2006 for (int i = 0; i < 4; ++i) NNN += m_Constr[i];
2007
2008 int i1 = 0;
2009
2010 for (int i = 0; i < 4; ++i) {
2011 if (m_Constr[i]) {
2012 if (i == 0) m_THM->SetBaryonChemicalPotential(x[i1]);
2013 if (i == 1) m_THM->SetElectricChemicalPotential(x[i1]);
2014 if (i == 2) m_THM->SetStrangenessChemicalPotential(x[i1]);
2015 if (i == 3) m_THM->SetCharmChemicalPotential(x[i1]);
2016 i1++;
2017 }
2018 }
2019
2020 m_THM->FillChemicalPotentials();
2021 m_THM->CalculatePrimordialDensities();
2022
2023 vector<double> dens(4, 0.), absdens(4, 0.);
2024 if (m_Constr[0]) {
2025 dens[0] = m_THM->CalculateBaryonDensity();
2026 absdens[0] = m_THM->CalculateAbsoluteBaryonDensity();
2027 }
2028 if (m_Constr[1]) {
2029 dens[1] = m_THM->CalculateChargeDensity();
2030 absdens[1] = m_THM->CalculateAbsoluteChargeDensity();
2031 }
2032 if (m_Constr[2]) {
2033 dens[2] = m_THM->CalculateStrangenessDensity();
2034 absdens[2] = m_THM->CalculateAbsoluteStrangenessDensity();
2035 if (absdens[2] < 1.e-25) {
2036 absdens[2] = m_THM->CalculateAbsoluteStrangenessDensityModulo();
2037 }
2038 }
2039 if (m_Constr[3]) {
2040 dens[3] = m_THM->CalculateCharmDensity();
2041 absdens[3] = m_THM->CalculateAbsoluteCharmDensity();
2042 if (absdens[3] < 1.e-25) {
2043 absdens[3] = m_THM->CalculateAbsoluteCharmDensityModulo();
2044 }
2045 }
2046
2047 vector<double> chi2s;
2048 chi2s.resize(m_THM->Densities().size());
2049 for (size_t i = 0; i < chi2s.size(); ++i) {
2050 chi2s[i] = m_THM->TPS()->Particle(i).chiDimensionfull(2, m_THM->Parameters(), m_THM->UseWidth(), m_THM->ChemicalPotential(i) + m_THM->MuShift(i))
2051 * xMath::GeVtoifm3();
2052 }
2053
2054 vector< vector<double> > deriv(4, vector<double>(4)), derivabs(4, vector<double>(4));
2055 for (int i = 0; i < 4; ++i)
2056 for (int j = 0; j < 4; ++j) {
2057 deriv[i][j] = 0.;
2058 for (size_t part = 0; part < chi2s.size(); ++part)
2059 deriv[i][j] += m_THM->TPS()->Particles()[part].GetCharge(i) * m_THM->TPS()->Particles()[part].GetCharge(j) * chi2s[part];
2060
2061 derivabs[i][j] = 0.;
2062 for (size_t part = 0; part < chi2s.size(); ++part)
2063 derivabs[i][j] += m_THM->TPS()->Particles()[part].GetAbsCharge(i) * m_THM->TPS()->Particles()[part].GetCharge(j) * chi2s[part];
2064 }
2065
2066
2067 MatrixXd ret(NNN, NNN);
2068
2069 i1 = 0;
2070
2071 for (int i = 0; i < 4; ++i) {
2072 if (m_Constr[i]) {
2073 int i2 = 0;
2074 for (int j = 0; j < 4; ++j)
2075 if (m_Constr[j]) {
2076 ret(i1, i2) = 0.;
2077 if (m_Type[i] == 0)
2078 ret(i1, i2) = deriv[i][j] * m_THM->Parameters().V / m_Totals[i];
2079 else
2080 ret(i1, i2) = deriv[i][j] / absdens[i] - dens[i] / absdens[i] / absdens[i] * derivabs[i][j];
2081 i2++;
2082 }
2083 i1++;
2084 }
2085 }
2086
2087 std::vector<double> retVec(NNN*NNN, 0);
2088 for (int i = 0; i < NNN; ++i)
2089 for (int j = 0; j < NNN; ++j)
2090 retVec[i*NNN + j] = ret(i, j);
2091
2092
2093 return retVec;
2094 }
2095
2096 void ThermalModelBase::SetDensityModelForParticleSpecies(int i, GeneralizedDensity *density_model)
2097 {
2098 if (i >= 0 && i < ComponentsNumber()) {
2099 TPS()->Particle(i).SetGeneralizedDensity(density_model);
2100 } else {
2101 std::cerr << "**WARNING** ThermalModelBase::SetDensityModelForParticleSpecies(): Particle id " << i << " is outside the range!" << std::endl;
2102 }
2103 }
2104
2105 void ThermalModelBase::SetDensityModelForParticleSpeciesByPdg(long long PDGID, GeneralizedDensity *density_model)
2106 {
2107 int id = TPS()->PdgToId(PDGID);
2108 if (id != -1) {
2109 TPS()->Particle(id).SetGeneralizedDensity(density_model);
2110 } else {
2111 std::cerr << "**WARNING** ThermalModelBase::SetDensityModelForParticleSpeciesByPdg(): Pdg code " << PDGID << " is outside the range!" << std::endl;
2112 }
2113 }
2114
2115 void ThermalModelBase::ClearDensityModels() {
2116 for (auto& particle : TPS()->Particles())
2117 particle.ClearGeneralizedDensity();
2118 }
2119
2120 void ThermalModelBase::SetMagneticField(double B, int lmax) {
2121 m_IGFExtraConfig.MagneticField.B = B;
2122 if (lmax >= 0)
2123 m_IGFExtraConfig.MagneticField.lmax = lmax;
2124 for (auto& particle : TPS()->Particles())
2125 particle.SetMagneticField(B, m_IGFExtraConfig.MagneticField.lmax);
2126 RecomputeThresholdsDueToMagneticField();
2127 ResetCalculatedFlags();
2128 }
2129
2130 void ThermalModelBase::RecomputeThresholdsDueToMagneticField() {
2131 const double &B = m_IGFExtraConfig.MagneticField.B;
2132 for (int ipart = 0; ipart < TPS()->ComponentsNumber(); ++ipart) {
2133 ThermalParticle &part = TPS()->Particle(ipart);
2134 if (!part.ZeroWidthEnforced() || part.ElectricCharge()==0) {
2135 if (B == 0.) {
2136 part.CalculateAndSetDynamicalThreshold();
2137 part.FillCoefficientsDynamical();
2138 } else if (!part.ZeroWidthEnforced()) {
2139 double szmax = (part.Degeneracy() - 1.) / 2.;
2140 if (szmax > 0.5) {
2141 double mmin = sqrt(2. * abs(part.ElectricCharge()) * B * (szmax - 0.5));
2142 part.CalculateAndSetDynamicalThreshold();
2143 if (mmin > part.DecayThresholdMassDynamical()) {
2144 part.SetDecayThresholdMassDynamical(mmin);
2145 }
2146 part.FillCoefficientsDynamical();
2147 }
2148 }
2149 }
2150 }
2151 }
2152
2153 void ThermalModelBase::ClearMagneticField() {
2154 m_IGFExtraConfig.MagneticField.B = 0.;
2155 for (auto& particle : TPS()->Particles())
2156 particle.ClearMagneticField();
2157 ResetCalculatedFlags();
2158 }
2159
2160 double ThermalModelBase::Susc(ConservedCharge::Name i, ConservedCharge::Name j) const {
2161 assert(IsSusceptibilitiesCalculated());
2162 return m_Susc[i][j];
2163 }
2164
2165 double ThermalModelBase::dSuscdT(ConservedCharge::Name i, ConservedCharge::Name j) const {
2166 assert(IsSusceptibilitiesCalculated() && IsTemperatureDerivativesCalculated());
2167 return m_dSuscdT[i][j];
2168 }
2169
2170 double ThermalModelBase::ProxySusc(ConservedCharge::Name i, ConservedCharge::Name j) const {
2171 assert(IsFluctuationsCalculated());
2172 return m_ProxySusc[i][j];
2173 }
2174
2175 double ThermalModelBase::GetdndT(int i) const {
2176 if (IsTemperatureDerivativesCalculated() && i >= 0 && i < ComponentsNumber())
2177 return m_dndT[i];
2178 else {
2179 std::cerr << "**WARNING** ThermalModelBase::GetdndT(): Temperature derivatives are not calculated or particle id " << i << " is outside the range!" << std::endl;
2180 return 0.;
2181 }
2182 return 0.;
2183 }
2184
2185 double ThermalModelBase::SusceptibilityDimensionfull(ConservedCharge::Name i, ConservedCharge::Name j) const {
2186 double ret = 0.;
2187 for (size_t k = 0; k < m_PrimCorrel.size(); ++k) {
2188 int c1 = m_TPS->Particles()[k].ConservedCharge(i);
2189 for (size_t kp = 0; kp < m_PrimCorrel.size(); ++kp) {
2190 int c2 = m_TPS->Particles()[kp].ConservedCharge(j);
2191 ret += c1 * c2 * m_PrimCorrel[k][kp];
2192 }
2193 }
2194 return ret / xMath::GeVtoifm() / xMath::GeVtoifm() / xMath::GeVtoifm();
2195 }
2196
2197 double ThermalModelBase::CalculateHeatCapacityMu()
2198 {
2199 double dedT = CalculateEnergyDensityDerivativeT(); // fm^-3
2200
2201 // T * ds/dT = (de/dT - mu_i d n_i / d T)
2202 double ret = dedT;
2203 for(int i = 0; i < TPS()->ComponentsNumber(); ++i) {
2204 ret -= m_Chem[i] * m_dndT[i];
2205 }
2206 return ret;
2207 }
2208
2209 // Auxiliary function to compute the susceptibility matrix
2210 MatrixXd GetSusceptibilityMatrix(ThermalModelBase* model, const vector<int>& ConservedDensities) {
2211 int N = 0;
2212 for (int i = 0; i < 4; ++i)
2213 if (ConservedDensities[i] == 1)
2214 N++;
2215
2216 assert(N > 0);
2217
2218 MatrixXd ret(N, N);
2219
2220 // Pi(i,j) = d P / dmui dmuj
2221 {
2222 int i1 = 0;
2223 for(int i = 0; i < 4; ++i) {
2224 if (ConservedDensities[i] == 0)
2225 continue;
2226 int i2 = 0;
2227 for(int j = 0; j < 4; ++j) {
2228 if (ConservedDensities[j] == 0)
2229 continue;
2230 ret(i1, i2) = model->SusceptibilityDimensionfull((ConservedCharge::Name)i, (ConservedCharge::Name)j) * xMath::GeVtoifm3();
2231 i2++;
2232 }
2233 i1++;
2234 }
2235 }
2236
2237 return ret;
2238 }
2239
2240 // Auxiliary function to compute the Hessian matrix of pressure
2241 MatrixXd GetPressureHessianMatrix(ThermalModelBase* model, const vector<int>& ConservedDensities) {
2242 int N = 1;
2243 for (int i = 0; i < 4; ++i)
2244 if (ConservedDensities[i] == 1)
2245 N++;
2246
2247 MatrixXd ret(N, N);
2248
2249 // Compute the Hessian matrix of pressure
2250 // Pi(0,0) = d^2 P / d T^2 = ds/dT
2251 ret(0, 0) = model->CalculateEntropyDensityDerivativeT(); // fm^-3 GeV^-1
2252
2253 // Pi(0,i) = Pi(i,0) = d^2 P / dT dmui = d ni / d T
2254 {
2255 int i1 = 0;
2256 for(int i = 0; i < 4; ++i) {
2257 if (ConservedDensities[i] == 0)
2258 continue;
2259 ret(0, i1 + 1) = ret(i1 + 1, 0) = model->ConservedChargeDensitydT((ConservedCharge::Name)i);
2260 i1++;
2261 }
2262 }
2263
2264 // Pi(i,j) = d^2 P / dmui dmuj
2265 {
2266 int i1 = 0;
2267 for(int i = 0; i < 4; ++i) {
2268 if (ConservedDensities[i] == 0)
2269 continue;
2270 int i2 = 0;
2271 for(int j = 0; j < 4; ++j) {
2272 if (ConservedDensities[j] == 0)
2273 continue;
2274 ret(i1 + 1, i2 + 1) = model->SusceptibilityDimensionfull((ConservedCharge::Name)i, (ConservedCharge::Name)j) * xMath::GeVtoifm3();
2275 i2++;
2276 }
2277 i1++;
2278 }
2279 }
2280
2281 return ret;
2282 }
2283
2284 vector<vector<double>> ThermalModelBase::PressureHessian()
2285 {
2286 // Ensure susceptibilities and temperature derivatives are calculated
2287 if (!m_FluctuationsCalculated)
2288 CalculateFluctuations();
2289 if (!m_TemperatureDerivativesCalculated)
2290 CalculateTemperatureDerivatives();
2291
2292 // Use the existing helper with all conserved charges enabled
2293 vector<int> allCharges = {1, 1, 1, 1}; // B, Q, S, C
2294 MatrixXd H_eigen = GetPressureHessianMatrix(this, allCharges);
2295
2296 // Convert Eigen matrix to vector<vector<double>>
2297 int N = H_eigen.rows();
2298 vector<vector<double>> H(N, vector<double>(N, 0.0));
2299 for (int i = 0; i < N; ++i) {
2300 for (int j = 0; j < N; ++j) {
2301 H[i][j] = H_eigen(i, j);
2302 }
2303 }
2304
2305 return H;
2306 }
2307
2308 bool ThermalModelBase::IsThermodynamicallyStable()
2309 {
2310 // Ensure susceptibilities and temperature derivatives are calculated
2311 if (!m_FluctuationsCalculated)
2312 CalculateFluctuations();
2313 if (!m_TemperatureDerivativesCalculated)
2314 CalculateTemperatureDerivatives();
2315
2316 // Get the Hessian matrix using the existing helper
2317 vector<int> allCharges = {1, 1, 1, 1}; // B, Q, S, C
2318 MatrixXd H = GetPressureHessianMatrix(this, allCharges);
2319
2320 // Use SelfAdjointEigenSolver since Hessian is symmetric
2321 Eigen::SelfAdjointEigenSolver<MatrixXd> solver(H);
2322 if (solver.info() != Eigen::Success) {
2323 return false; // Eigenvalue computation failed
2324 }
2325
2326 // Check if all eigenvalues are non-negative (positive semi-definite)
2327 const auto& eigenvalues = solver.eigenvalues();
2328 for (int i = 0; i < eigenvalues.size(); ++i) {
2329 if (eigenvalues(i) < 0.0) {
2330 return false;
2331 }
2332 }
2333
2334 return true;
2335 }
2336
2337 double ThermalModelBase::CalculateAdiabaticSpeedOfSoundSquared(bool rhoBconst, bool rhoQconst, bool rhoSconst, bool rhoCconst) {
2338 double ret = 0.;
2339 const double nanv = std::numeric_limits<double>::quiet_NaN();
2340 const double eps = 1e-14;
2341
2342 double T = Parameters().T;
2343
2344 std::vector<double> mus = {Parameters().muB, Parameters().muQ, Parameters().muS, Parameters().muC};
2345 std::vector<int> ConservedCharges = {static_cast<int>(rhoBconst), static_cast<int>(rhoQconst), static_cast<int>(rhoSconst), static_cast<int>(rhoCconst)};
2346
2347 if (!TPS()->hasBaryons())
2348 ConservedCharges[0] = 0;
2349 if (!TPS()->hasCharged())
2350 ConservedCharges[1] = 0;
2351 if (!TPS()->hasStrange())
2352 ConservedCharges[2] = 0;
2353 if (!TPS()->hasCharmed())
2354 ConservedCharges[3] = 0;
2355
2356 // Drop conserved charge constraints that are vacuous:
2357 // if the absolute density is zero, the constraint is trivially satisfied
2358 // and including it would make the susceptibility/Hessian matrix singular
2359 // (e.g. strangeness at T=0 when mu_S < m_K)
2360 for (int i = 0; i < 4; ++i) {
2361 if (ConservedCharges[i] == 1
2362 && AbsoluteConservedChargeDensity((ConservedCharge::Name)i) < eps)
2363 ConservedCharges[i] = 0;
2364 }
2365
2366 // Zero chemical potentials
2367 {
2368 bool AllMuZero = true;
2369 for(int i = 0; i < ConservedCharges.size(); ++i)
2370 AllMuZero &= (ConservedCharges[i] == 0 || mus[i] == 0.0);
2371 if (AllMuZero) {
2372 if (T == 0.) {
2373 WarnIllDefinedSoundSpeedOnce(
2374 "CalculateAdiabaticSpeedOfSoundSquared", ConservedCharges, T,
2375 "all relevant chemical potentials are zero (0/0 limit at T=0).");
2376 return nanv;
2377 }
2378 // At mu = 0 cs2 = s(T) / e'(T)
2379 // Calculate temperature derivatives if not already
2380 if (!IsTemperatureDerivativesCalculated())
2381 CalculateTemperatureDerivatives();
2382 return EntropyDensity() / CalculateEnergyDensityDerivativeT();
2383 //return EntropyDensity() / HeatCapacityMu();
2384 }
2385 }
2386
2387 int Ndens = 0;
2388 for(int i = 0; i < 4; ++i)
2389 if (ConservedCharges[i] == 1)
2390 Ndens++;
2391 if (Ndens == 0) {
2392 WarnIllDefinedSoundSpeedOnce(
2393 "CalculateAdiabaticSpeedOfSoundSquared", ConservedCharges, T,
2394 "no conserved densities selected.");
2395 return nanv;
2396 }
2397
2398 // Compute fluctuations if not already
2399 if (!IsFluctuationsCalculated())
2400 CalculateFluctuations();
2401
2402 // Zero temperature
2403 if (T == 0.) {
2404
2405 MatrixXd chi2matr = GetSusceptibilityMatrix(this, ConservedCharges);
2406 if (!IsMatrixInvertible(chi2matr)) {
2407 WarnIllDefinedSoundSpeedOnce(
2408 "CalculateAdiabaticSpeedOfSoundSquared", ConservedCharges, T,
2409 "susceptibility matrix is singular at T=0.");
2410 return nanv;
2411 }
2412 MatrixXd chi2inv = chi2matr.inverse();
2413 ret = 0.;
2414 int i1 = 0, i2 = 0;
2415 for(int i = 0; i < 4; ++i) {
2416 if (ConservedCharges[i] == 0)
2417 continue;
2418 i2 = 0;
2419 for(int j = 0; j < 4; ++j) {
2420 if (ConservedCharges[j] == 0)
2421 continue;
2422
2423 auto chg1 = (ConservedCharge::Name)i;
2424 auto chg2 = (ConservedCharge::Name)j;
2425 ret += ConservedChargeDensity(chg1) * ConservedChargeDensity(chg2) * chi2inv(i1, i2);
2426 i2++;
2427 }
2428 i1++;
2429 }
2430
2431 const double ep = EnergyDensity() + Pressure();
2432 if (std::abs(ep) < eps) {
2433 WarnIllDefinedSoundSpeedOnce(
2434 "CalculateAdiabaticSpeedOfSoundSquared", ConservedCharges, T,
2435 "EnergyDensity + Pressure is zero.");
2436 return nanv;
2437 }
2438 ret /= ep;
2439 return ret;
2440 }
2441
2442 // General case
2443
2444 // Calculate temperature derivatives if not already
2445 if (!IsTemperatureDerivativesCalculated())
2446 CalculateTemperatureDerivatives();
2447
2448 // double dedT = HeatCapacityMu(); // fm^-3
2449 // Compute the Hessian matrix of pressure
2450 MatrixXd PiMatr = GetPressureHessianMatrix(this, ConservedCharges);
2451
2452 // Compute the inverse of the Hessian matrix
2453 if (!IsMatrixInvertible(PiMatr)) {
2454 WarnIllDefinedSoundSpeedOnce(
2455 "CalculateAdiabaticSpeedOfSoundSquared", ConservedCharges, T,
2456 "pressure Hessian matrix is singular.");
2457 return nanv;
2458 }
2459 MatrixXd PiMatrInv = PiMatr.inverse();
2460 // Prepare the vector x
2461 double s = EntropyDensity();
2462 vector<double> x = {s};
2463 for(int i = 0; i < ConservedCharges.size(); ++i)
2464 if (ConservedCharges[i] == 1)
2465 x.push_back(ConservedChargeDensity((ConservedCharge::Name)i));
2466
2467 ret = 0.;
2468 for(int i = 0; i < x.size(); ++i)
2469 for(int j = 0; j < x.size(); ++j)
2470 ret += x[i] * x[j] * PiMatrInv(i,j);
2471 const double ep = EnergyDensity() + Pressure();
2472 if (std::abs(ep) < eps) {
2473 WarnIllDefinedSoundSpeedOnce(
2474 "CalculateAdiabaticSpeedOfSoundSquared", ConservedCharges, T,
2475 "EnergyDensity + Pressure is zero.");
2476 return nanv;
2477 }
2478 ret /= ep;
2479
2480 return ret;
2481 }
2482
2483 double ThermalModelBase::CalculateIsothermalSpeedOfSoundSquared(bool rhoBconst, bool rhoQconst, bool rhoSconst, bool rhoCconst) {
2484 double ret = 0.;
2485 const double nanv = std::numeric_limits<double>::quiet_NaN();
2486 const double eps = 1e-14;
2487
2488 double T = Parameters().T;
2489
2490 std::vector<double> mus = {Parameters().muB, Parameters().muQ, Parameters().muS, Parameters().muC};
2491 std::vector<int> ConservedCharges = {static_cast<int>(rhoBconst), static_cast<int>(rhoQconst), static_cast<int>(rhoSconst), static_cast<int>(rhoCconst)};
2492
2493 if (!TPS()->hasBaryons())
2494 ConservedCharges[0] = 0;
2495 if (!TPS()->hasCharged())
2496 ConservedCharges[1] = 0;
2497 if (!TPS()->hasStrange())
2498 ConservedCharges[2] = 0;
2499 if (!TPS()->hasCharmed())
2500 ConservedCharges[3] = 0;
2501
2502 // Drop conserved charge constraints that are vacuous:
2503 // if the absolute density is zero, the constraint is trivially satisfied
2504 // and including it would make the susceptibility matrix singular
2505 // (e.g. strangeness at T=0 when mu_S < m_K)
2506 for (int i = 0; i < 4; ++i) {
2507 if (ConservedCharges[i] == 1
2508 && AbsoluteConservedChargeDensity((ConservedCharge::Name)i) < eps)
2509 ConservedCharges[i] = 0;
2510 }
2511
2512 int Ndens = 0;
2513 for (int i = 0; i < 4; ++i)
2514 if (ConservedCharges[i] == 1)
2515 Ndens++;
2516 if (Ndens == 0) {
2517 WarnIllDefinedSoundSpeedOnce(
2518 "CalculateIsothermalSpeedOfSoundSquared", ConservedCharges, T,
2519 "no conserved densities selected.");
2520 return nanv;
2521 }
2522
2523 // If all chemical potentials are zero, c_T is undefined (return 0.)
2524 {
2525 bool AllMuZero = true;
2526 for(int i = 0; i < ConservedCharges.size(); ++i)
2527 AllMuZero &= (ConservedCharges[i] == 0 || mus[i] == 0.0);
2528 if (AllMuZero) {
2529 return 0.;
2530 }
2531 }
2532
2533 auto chi2Matr = GetSusceptibilityMatrix(this, ConservedCharges);
2534 if (!IsMatrixInvertible(chi2Matr)) {
2535 WarnIllDefinedSoundSpeedOnce(
2536 "CalculateIsothermalSpeedOfSoundSquared", ConservedCharges, T,
2537 "susceptibility matrix is singular.");
2538 return nanv;
2539 }
2540 // Compute the inverse of the susceptibility matrix
2541 MatrixXd chi2inv= chi2Matr.inverse();
2542
2543 double numerator = 0.;
2544 double denominator = 0.;
2545 int i1 = 0, i2 = 0;
2546 for(int i = 0; i < 4; ++i) {
2547 if (ConservedCharges[i] == 0)
2548 continue;
2549 denominator += mus[i] * ConservedChargeDensity((ConservedCharge::Name)i);
2550 i2 = 0;
2551 for(int j = 0; j < 4; ++j) {
2552 if (ConservedCharges[j] == 0)
2553 continue;
2554
2555 auto chg1 = (ConservedCharge::Name)i;
2556 auto chg2 = (ConservedCharge::Name)j;
2557
2558 numerator += ConservedChargeDensity(chg1) * ConservedChargeDensity(chg2) * chi2inv(i1, i2);
2559 denominator += T * ConservedChargeDensitydT(chg1) * chi2inv(i1, i2) * ConservedChargeDensity(chg2);
2560
2561 ret += ConservedChargeDensity(chg1) * ConservedChargeDensity(chg2) * chi2inv(i1, i2);
2562 i2++;
2563 }
2564 i1++;
2565 }
2566
2567 if (std::abs(denominator) < eps) {
2568 WarnIllDefinedSoundSpeedOnce(
2569 "CalculateIsothermalSpeedOfSoundSquared", ConservedCharges, T,
2570 "denominator is zero.");
2571 return nanv;
2572 }
2573 return numerator / denominator;
2574 }
2575
2576 double ThermalModelBase::CalculateHeatCapacity(bool rhoBconst, bool rhoQconst, bool rhoSconst, bool rhoCconst) {
2577 double T = Parameters().T;
2578
2579 if (T == 0.) {
2580 return 0.;
2581 }
2582
2583 std::vector<double> mus = {Parameters().muB, Parameters().muQ, Parameters().muS, Parameters().muC};
2584 std::vector<int> ConservedCharges = {static_cast<int>(rhoBconst), static_cast<int>(rhoQconst), static_cast<int>(rhoSconst), static_cast<int>(rhoCconst)};
2585
2586 if (!TPS()->hasBaryons())
2587 ConservedCharges[0] = 0;
2588 if (!TPS()->hasCharged())
2589 ConservedCharges[1] = 0;
2590 if (!TPS()->hasStrange())
2591 ConservedCharges[2] = 0;
2592 if (!TPS()->hasCharmed())
2593 ConservedCharges[3] = 0;
2594
2595 // Drop conserved charge constraints that are vacuous:
2596 // if the absolute density is zero, the constraint is trivially satisfied
2597 // and including it would make the Hessian matrix singular
2598 // (e.g. strangeness at T=0 when mu_S < m_K)
2599 {
2600 const double eps = 1e-14;
2601 for (int i = 0; i < 4; ++i) {
2602 if (ConservedCharges[i] == 1
2603 && AbsoluteConservedChargeDensity((ConservedCharge::Name)i) < eps)
2604 ConservedCharges[i] = 0;
2605 }
2606 }
2607
2608 // Calculate the numerator
2609 double TdsdT = HeatCapacityMu();
2610
2611 // If all chemical potentials are zero, return TdsT
2612 {
2613 bool AllMuZero = true;
2614 for(int i = 0; i < ConservedCharges.size(); ++i)
2615 AllMuZero &= (ConservedCharges[i] == 0 || mus[i] == 0.0);
2616 if (AllMuZero) {
2617 return TdsdT;
2618 }
2619 }
2620
2621 if (!IsFluctuationsCalculated())
2622 CalculateFluctuations();
2623
2624 // Calculate the denominator
2625 double denominator = 1.;
2626 auto PiMatr = GetPressureHessianMatrix(this, ConservedCharges);
2627 auto PiMatrInv = PiMatr.inverse();
2628 int N = PiMatr.rows();
2629 for(int i = 1; i < N; ++i) {
2630 denominator -= PiMatr(0, i) * PiMatrInv(i, 0);
2631 }
2632
2633 double heatCapacity = TdsdT / denominator;
2634 return heatCapacity;
2635 }
2636
2637} // namespace thermalfist
params
map< string, double > params
Definition NeutronStars-CSHRG.cpp:34
Utility.h
Contains some helper functions.
thermalfist::Broyden::BroydenSolutionCriterium
Sub-class where it is determined whether the required accuracy is achieved in the Broyden's method.
Definition Broyden.h:144
thermalfist::BroydenEquations
Abstract class which defines the system of non-linear equations to be solved by the Broyden's method.
Definition Broyden.h:32
thermalfist::Broyden
Class implementing the Broyden method to solve a system of non-linear equations.
Definition Broyden.h:132
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::BroydenJacobian
Class which implements calculation of the Jacobian needed for the Broyden's method.
Definition Broyden.h:77
thermalfist::Disclaimer::PrintDisclaimer
static bool PrintDisclaimer()
Definition Utility.cpp:49
thermalfist::Disclaimer::DisclaimerPrinted
static bool DisclaimerPrinted
Definition Utility.h:28
thermalfist::GeneralizedDensity
Implements the possibility of a generalized calculation of the densities. For example,...
Definition IdealGasFunctions.h:982
thermalfist::ThermalModelBase
Abstract base class for an HRG model implementation.
Definition ThermalModelBase.h:29
thermalfist::ThermalModelBase::CalculateSusceptibilityMatrix
virtual void CalculateSusceptibilityMatrix()
Calculates the conserved charges susceptibility matrix.
Definition ThermalModelBase.cpp:1353
thermalfist::ThermalModelBase::FullIdealChemicalPotential
virtual double FullIdealChemicalPotential(int i) const
Chemical potential entering the ideal gas expressions of particle species i.
Definition ThermalModelBase.cpp:431
thermalfist::ThermalModelBase::FixParametersNoReset
virtual void FixParametersNoReset()
Method which actually implements ConstrainChemicalPotentialsNoReset() (for backward compatibility).
Definition ThermalModelBase.cpp:531
thermalfist::ThermalModelBase::ChangeTPS
virtual void ChangeTPS(ThermalParticleSystem *TPS)
Change the particle list.
Definition ThermalModelBase.cpp:365
thermalfist::ThermalModelBase::SetGammaC
virtual void SetGammaC(double gammaC)
Set the charm quark fugacity factor.
Definition ThermalModelBase.cpp:226
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::SetTemperature
virtual void SetTemperature(double T)
Set the temperature.
Definition ThermalModelBase.cpp:186
thermalfist::ThermalModelBase::CalculateFluctuations
virtual void CalculateFluctuations()
Computes the fluctuation observables.
Definition ThermalModelBase.cpp:1483
thermalfist::ThermalModelBase::UsePartialChemicalEquilibrium
bool UsePartialChemicalEquilibrium()
Whether partial chemical equilibrium with additional chemical potentials is used.
Definition ThermalModelBase.h:484
thermalfist::ThermalModelBase::SolveChemicalPotentials
virtual bool SolveChemicalPotentials(double totB=0., double totQ=0., double totS=0., double totC=0., double muBinit=0., double muQinit=0., double muSinit=0., double muCinit=0., bool ConstrMuB=true, bool ConstrMuQ=true, bool ConstrMuS=true, bool ConstrMuC=true)
The procedure which calculates the chemical potentials which reproduce the specified total baryon,...
Definition ThermalModelBase.cpp:582
thermalfist::ThermalModelBase::TwoParticleSusceptibilityPrimordialByPdg
virtual double TwoParticleSusceptibilityPrimordialByPdg(long long id1, long long id2)
Returns the computed primordial particle number (cross-)susceptibility for particles with pdg codes ...
Definition ThermalModelBase.cpp:1136
thermalfist::ThermalModelBase::CalculateChargeFluctuations
virtual std::vector< double > CalculateChargeFluctuations(const std::vector< double > &chgs, int order=4, bool dimensionfull=false)
Calculates fluctuations (diagonal susceptibilities) of an arbitrary "conserved" charge.
Definition ThermalModelBase.cpp:697
thermalfist::ThermalModelBase::SetVirial
virtual void SetVirial(int, int, double)
Set the excluded volume coefficient .
Definition ThermalModelBase.h:239
thermalfist::ThermalModelBase::FillChemicalPotentials
virtual void FillChemicalPotentials()
Sets the chemical potentials of all particles.
Definition ThermalModelBase.cpp:400
thermalfist::ThermalModelBase::CalculateTemperatureDerivatives
virtual void CalculateTemperatureDerivatives()
Computes the temperature derivatives of densities, shifted chemical potentials, and primordial hadron...
Definition ThermalModelBase.cpp:1487
thermalfist::ThermalModelBase::ConstrainChemicalPotentials
void ConstrainChemicalPotentials(bool resetInitialValues=true)
Constrains the chemical potentials by the conservation laws imposed.
Definition ThermalModelBase.cpp:490
thermalfist::ThermalModelBase::FixParameters
virtual void FixParameters()
Method which actually implements ConstrainChemicalPotentials() (for backward compatibility).
Definition ThermalModelBase.cpp:498
thermalfist::ThermalModelBase::TwoParticleSusceptibilityFinalByPdg
virtual double TwoParticleSusceptibilityFinalByPdg(long long id1, long long id2)
Returns the computed final particle number (cross-)susceptibility for particles with pdg codes id1 a...
Definition ThermalModelBase.cpp:1230
thermalfist::ThermalModelBase::TwoParticleSusceptibilityFinal
virtual double TwoParticleSusceptibilityFinal(int i, int j) const
Returns the computed final particle number (cross-)susceptibility for particles with ids i and j....
Definition ThermalModelBase.cpp:1214
thermalfist::ThermalModelBase::SetGammaS
virtual void SetGammaS(double gammaS)
Set the strange quark fugacity factor.
Definition ThermalModelBase.cpp:220
thermalfist::ThermalModelBase::SetCharm
virtual void SetCharm(int C)
Set the total charm (for canonical ensemble only)
Definition ThermalModelBase.cpp:250
thermalfist::ThermalModelBase::FinalNetParticleChargeSusceptibilityByPdg
virtual double FinalNetParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg)
Returns the computed final (after decays) net particle vs conserved charge (cross-)susceptibility fo...
Definition ThermalModelBase.cpp:1338
thermalfist::ThermalModelBase::SetResonanceWidthIntegrationType
void SetResonanceWidthIntegrationType(ThermalParticle::ResonanceWidthIntegration type)
Set the ThermalParticle::ResonanceWidthIntegration scheme for all particles. Calls the corresponding ...
Definition ThermalModelBase.cpp:390
thermalfist::ThermalModelBase::CalculateFeeddown
virtual void CalculateFeeddown()
Calculates the total densities which include feeddown contributions.
Definition ThermalModelBase.cpp:450
thermalfist::ThermalModelBase::TwoParticleSusceptibilityTemperatureDerivativePrimordial
virtual double TwoParticleSusceptibilityTemperatureDerivativePrimordial(int i, int j) const
Returns the computed temperature derivative of the primordial particle number (cross-)susceptibility ...
Definition ThermalModelBase.cpp:1153
thermalfist::ThermalModelBase::SetUseWidth
void SetUseWidth(bool useWidth)
Sets whether finite resonance widths are used. Deprecated.
Definition ThermalModelBase.cpp:140
thermalfist::ThermalModelBase::FinalParticleChargeSusceptibilityByPdg
virtual double FinalParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg)
Returns the computed final (after decays) particle vs conserved charge (cross-)susceptibility for pa...
Definition ThermalModelBase.cpp:1327
thermalfist::ThermalModelBase::DisableBaryonAntiBaryonAttraction
virtual void DisableBaryonAntiBaryonAttraction()
Definition ThermalModelBase.cpp:345
thermalfist::ThermalModelBase::FinalParticleChargeSusceptibility
virtual double FinalParticleChargeSusceptibility(int i, ConservedCharge::Name chg) const
Returns the computed final (after decays) particle vs conserved charge (cross-)susceptibility for pa...
Definition ThermalModelBase.cpp:1318
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::SetElectricCharge
virtual void SetElectricCharge(int Q)
Set the total electric charge (for canonical ensemble only)
Definition ThermalModelBase.cpp:238
thermalfist::ThermalModelBase::SetBaryonChemicalPotential
virtual void SetBaryonChemicalPotential(double muB)
Set the baryon chemical potential.
Definition ThermalModelBase.cpp:192
thermalfist::ThermalModelBase::CalculateTwoParticleCorrelations
virtual void CalculateTwoParticleCorrelations()
Computes the fluctuations and correlations of the primordial particle numbers.
Definition ThermalModelBase.cpp:1023
thermalfist::ThermalModelBase::SusceptibilityDimensionfull
virtual double SusceptibilityDimensionfull(ConservedCharge::Name i, ConservedCharge::Name j) const
A 2nd order susceptibility of conserved charges.
Definition ThermalModelBase.cpp:2185
thermalfist::ThermalModelBase::SetNormBratio
void SetNormBratio(bool normBratio)
Whether branching ratios are renormalized to 100%.
Definition ThermalModelBase.cpp:159
thermalfist::ThermalModelBase::SetParameters
virtual void SetParameters(const ThermalModelParameters &params)
The thermal parameters.
Definition ThermalModelBase.cpp:179
thermalfist::ThermalModelBase::SetElectricChemicalPotential
virtual void SetElectricChemicalPotential(double muQ)
Set the electric chemical potential.
Definition ThermalModelBase.cpp:199
thermalfist::ThermalModelBase::PrimordialNetParticleChargeSusceptibilityByPdg
virtual double PrimordialNetParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg)
Returns the computed primordial net particle vs conserved charge (cross-)susceptibility for particle...
Definition ThermalModelBase.cpp:1303
thermalfist::ThermalModelBase::TwoParticleSusceptibilityTemperatureDerivativePrimordialByPdg
virtual double TwoParticleSusceptibilityTemperatureDerivativePrimordialByPdg(long long id1, long long id2)
Returns the computed temperature derivative of the primordial particle number (cross-)susceptibility ...
Definition ThermalModelBase.cpp:1161
thermalfist::ThermalModelBase::NetParticleSusceptibilityFinalByPdg
virtual double NetParticleSusceptibilityFinalByPdg(long long id1, long long id2)
Returns the computed (cross-)susceptibility between final net-particle numbers for pdg codes id1 and...
Definition ThermalModelBase.cpp:1247
thermalfist::ThermalModelBase::SetChemicalPotential
virtual void SetChemicalPotential(int i, double chem)
Sets the chemical potential of particle species i.
Definition ThermalModelBase.cpp:423
thermalfist::ThermalModelBase::FixChemicalPotentialsThroughDensities
virtual bool FixChemicalPotentialsThroughDensities(double rhoB=0., double rhoQ=0., double rhoS=0., double rhoC=0., double muBinit=0., double muQinit=0., double muSinit=0., double muCinit=0., bool ConstrMuB=true, bool ConstrMuQ=true, bool ConstrMuS=true, bool ConstrMuC=true)
The procedure which calculates the chemical potentials which reproduce the specified baryon,...
Definition ThermalModelBase.cpp:572
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::SetBaryonCharge
virtual void SetBaryonCharge(int B)
Set the total baryon number (for canonical ensemble only)
Definition ThermalModelBase.cpp:232
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::PrimordialParticleChargeSusceptibility
virtual double PrimordialParticleChargeSusceptibility(int i, ConservedCharge::Name chg) const
Returns the computed primordial particle vs conserved charge (cross-)susceptibility for particle wit...
Definition ThermalModelBase.cpp:1283
thermalfist::ThermalModelBase::CalculatePrimordialDensities
virtual void CalculatePrimordialDensities()=0
Calculates the primordial densities of all species.
thermalfist::ThermalModelBase::SetStrangenessChemicalPotential
virtual void SetStrangenessChemicalPotential(double muS)
Set the strangeness chemical potential.
Definition ThermalModelBase.cpp:206
thermalfist::ThermalModelBase::TwoParticleSusceptibilityPrimordial
virtual double TwoParticleSusceptibilityPrimordial(int i, int j) const
Returns the computed primordial particle number (cross-)susceptibility for particles with ids i and ...
Definition ThermalModelBase.cpp:1127
thermalfist::ThermalModelBase::CalculateDensities
virtual void CalculateDensities()
Calculates the primordial and total (after decays) densities of all species.
Definition ThermalModelBase.cpp:669
thermalfist::ThermalModelBase::FillVirial
virtual void FillVirial(const std::vector< double > &ri=std::vector< double >(0))
Fills the excluded volume coefficients based on the provided radii parameters for all species.
Definition ThermalModelBase.cpp:136
thermalfist::ThermalModelBase::SetGammaq
virtual void SetGammaq(double gammaq)
Set the light quark fugacity factor.
Definition ThermalModelBase.cpp:358
thermalfist::ThermalModelBase::SetStrangeness
virtual void SetStrangeness(int S)
Set the total strangeness (for canonical ensemble only)
Definition ThermalModelBase.cpp:244
thermalfist::ThermalModelBase::SetAttraction
virtual void SetAttraction(int, int, double)
Set the vdW mean field attraction coefficient .
Definition ThermalModelBase.h:251
thermalfist::ThermalModelBase::PrimordialParticleChargeSusceptibilityByPdg
virtual double PrimordialParticleChargeSusceptibilityByPdg(long long id1, ConservedCharge::Name chg)
Returns the computed primordial particle vs conserved charge (cross-)susceptibility for particle wit...
Definition ThermalModelBase.cpp:1292
thermalfist::ThermalModelBase::Parameters
const ThermalModelParameters & Parameters() const
Definition ThermalModelBase.h:121
thermalfist::ThermalModelBase::NetParticleSusceptibilityPrimordialByPdg
virtual double NetParticleSusceptibilityPrimordialByPdg(long long id1, long long id2)
Returns the computed (cross-)susceptibility between primordial net-particle numbers for pdg codes id...
Definition ThermalModelBase.cpp:1178
thermalfist::ThermalModelBase::SetCharmChemicalPotential
virtual void SetCharmChemicalPotential(double muC)
Set the charm chemical potential.
Definition ThermalModelBase.cpp:213
thermalfist::ThermalModelBase::SetStatistics
virtual void SetStatistics(bool stats)
Definition ThermalModelBase.cpp:384
thermalfist::ThermalModelBase::Volume
double Volume() const
System volume (fm )
Definition ThermalModelBase.h:517
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::ThermalParticle
Class containing all information about a particle specie.
Definition ThermalParticle.h:64
thermalfist::ThermalParticle::AbsoluteCharm
double AbsoluteCharm() const
Absolute charm quark content |s|.
Definition ThermalParticle.h:448
thermalfist::ThermalParticle::BaryonCharge
int BaryonCharge() const
Particle's baryon number.
Definition ThermalParticle.h:396
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::ThermalParticle::Decays
const ParticleDecaysVector & Decays() const
A vector of particle's decays.
Definition ThermalParticle.h:579
thermalfist::ThermalParticle::ResonanceWidthIntegration
ResonanceWidthIntegration
Treatment of finite resonance widths.
Definition ThermalParticle.h:82
thermalfist::ThermalParticle::eBW
@ eBW
Energy-dependent Breit-Wigner scheme (eBW)
Definition ThermalParticle.h:87
thermalfist::ThermalParticle::BWTwoGamma
@ BWTwoGamma
Energy-independent Breit-Wigner in +-2\Gamma interval.
Definition ThermalParticle.h:84
thermalfist::ThermalParticle::ZeroWidth
@ ZeroWidth
Zero-width approximation.
Definition ThermalParticle.h:83
thermalfist::ThermalParticle::AbsoluteStrangeness
double AbsoluteStrangeness() const
Absolute strange quark content |s|.
Definition ThermalParticle.h:442
thermalfist::ThermalParticle::AbsoluteQuark
double AbsoluteQuark() const
Absolute light quark content |u,d|.
Definition ThermalParticle.h:436
thermalfist::ThermalParticle::DecaysOriginal
const ParticleDecaysVector & DecaysOriginal() const
A backup copy of particle's decays.
Definition ThermalParticle.h:598
thermalfist::ThermalParticleSystem
Class containing the particle list.
Definition ThermalParticleSystem.h:30
thermalfist::ThermalParticleSystem::ResonanceFinalStatesDistribution
std::vector< std::pair< double, std::vector< int > > > ResonanceFinalStatesDistribution
Definition ThermalParticleSystem.h:183
thermalfist::ThermalParticleSystem::DecayContributionsToParticle
std::vector< SingleDecayContribution > DecayContributionsToParticle
A vector of SingleDecayContribution where each element corresponds to a certain resonance species.
Definition ThermalParticleSystem.h:146
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::GetSusceptibilityMatrix
MatrixXd GetSusceptibilityMatrix(ThermalModelBase *model, const vector< int > &ConservedDensities)
Definition ThermalModelBase.cpp:2210
thermalfist::GetPressureHessianMatrix
MatrixXd GetPressureHessianMatrix(ThermalModelBase *model, const vector< int > &ConservedDensities)
Definition ThermalModelBase.cpp:2241
thermalfist::ConservedCharge::Name
Name
Set of all conserved charges considered.
Definition ThermalParticle.h:46
thermalfist::ConservedCharge::ElectricCharge
@ ElectricCharge
Electric charge.
Definition ThermalParticle.h:48
thermalfist::Feeddown::Strong
@ Strong
Feeddown from strong decays.
Definition ParticleDecay.h:44
thermalfist::Feeddown::Weak
@ Weak
Feeddown from strong, electromagnetic, and weak decays.
Definition ParticleDecay.h:42
thermalfist::Feeddown::StabilityFlag
@ StabilityFlag
Feeddown from all particles marked as unstable.
Definition ParticleDecay.h:41
thermalfist::Feeddown::Primordial
@ Primordial
No feeddown.
Definition ParticleDecay.h:40
thermalfist::ParticleDecayType::NumberOfDecayTypes
static const int NumberOfDecayTypes
Definition ParticleDecay.h:71
thermalfist::ThermalModelParameters
Structure containing all thermal parameters of the model.
Definition ThermalModelParameters.h:18
Generated by doxygen 1.13.2