Thermal-FIST 1.5
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
CosmicTrajectory.cpp
Go to the documentation of this file.
1/*
2 * Thermal-FIST package
3 *
4 * Copyright (c) 2019-2025 Volodymyr Vovchenko
5 *
6 * GNU General Public License (GPLv3 or later)
7 */
8#include <string.h>
9#include <fstream>
10#include <iostream>
11#include <iomanip>
12#include <ctime>
13#include <cstdio>
14#include <vector>
15#include <string>
16
17#include "HRGBase.h"
18#include "HRGEV.h"
19#include "HRGVDW.h"
20#include "HRGFit.h"
21
22#include "ThermalFISTConfig.h"
23
24#include "CosmicEos/CosmicEoS.h"
25
26
27using namespace std;
28using namespace thermalfist;
29
30
31// Compute chemical potentials as a function of temperature in the Early Universe
32// Reference: V. Vovchenko et al., Phys. Rev. Lett. 126 (2021) 012701, https://inspirehep.net/literature/1815329
33int main(int argc, char *argv[])
34{
35 // Baryon asymmetry
36 double b = 8.6e-11;
37 // Electric charge asymmetry
38 double q = 0.0;
39
40 // Electron flavor asymmetry (read from command line)
41 double le = 0.0;
42 if (argc > 1)
43 le = atof(argv[1]);
44
45 // Muon flavor asymmetry (read from command line)
46 double lmu = 0.0;
47 if (argc > 2)
48 lmu = atof(argv[2]);
49
50 // Tau flavor asymmetry (read from command line)
51 double ltau = 0.0;
52 if (argc > 3)
53 ltau = atof(argv[3]);
54
55 // If very small, set to zero
56 if (fabs(le) < 1.e-13)
57 le = 0.0;
58 if (fabs(lmu) < 1.e-13)
59 lmu = 0.0;
60 if (fabs(ltau) < 1.e-13)
61 ltau = 0.0;
62
63
64 // Particle list for the HRG model
65 ThermalParticleSystem parts(string(ThermalFIST_INPUT_FOLDER) + "/list/PDG2020/list.dat");
66
67 // HRG model
68 ThermalModelIdeal modelHRG(&parts);
69
70 // For interactin HRG use e.g.
71 // ThermalModelVDW model(&parts);
72 // double a = 0.329, b = 3.42;
73 // SetVDWHRGInteractionParameters(&model, a, b);
74
75 // Resonance widths
76 bool useWidth = false;
77 modelHRG.SetUseWidth(useWidth);
78 if (useWidth)
79 modelHRG.SetUseWidth(ThermalParticle::eBWconstBR);
80
81 // Quantum statistics
82 bool useQStats = true;
83 modelHRG.SetStatistics(useQStats);
84
85 // Effective mass model for pions (needed for BEC)
86 bool interactingpions = true;
87
88 // Cosmic EoS object
89 CosmicEoS cosmos(&modelHRG, interactingpions);
90
91
92
93
94 cout << "le + lmu = " << setw(15) << le + lmu << endl;
95 cout << "le - lmu = " << setw(15) << le - lmu << endl;
96 cout << "le = " << setw(15) << le << endl;
97 cout << "lmu = " << setw(15) << lmu << endl;
98 cout << "ltau = " << setw(15) << ltau << endl;
99
100 cosmos.SetAsymmetries(vector<double>({ b, 0., le, lmu, ltau }));
101
102 // Range of temperatures
103 double Tmin = 0.010; // GeV
104 double Tmax = 0.180; // GeV
105 double dT = 0.001; // GeV
106 vector<double> Temps;
107 for (double tT = Tmin; tT <= Tmax + 0.1 * dT; tT += dT) {
108 Temps.push_back(tT);
109 }
110
111 // Output file (write le + lmu since it's the main ingredient for pion condersation)
112 string filename = "CosmicTrajectory";
113 {
114 char cc[500];
115 sprintf(cc, "le+lmu.%lf", le + lmu);
116 filename += "." + string(cc);
117 }
118 filename += ".dat";
119
120 ofstream fout(filename);
121
122 // We will output the trajectories and the various properties of the EoS
123 fout << setw(15) << "T[MeV]" << " ";
124 fout << setw(15) << "muB[MeV]" << " ";
125 fout << setw(15) << "muQ[MeV]" << " ";
126 fout << setw(15) << "mue[MeV]" << " ";
127 fout << setw(15) << "mum[MeV]" << " ";
128 fout << setw(15) << "mut[MeV]" << " ";
129 fout << setw(15) << "pion_bec" << " ";
130 fout << setw(15) << "nI/T3" << " ";
131 fout << setw(15) << "pT4" << " ";
132 fout << setw(15) << "eT4" << " ";
133 fout << setw(15) << "IT4" << " ";
134 fout << setw(15) << "pT4_QCD" << " ";
135 fout << setw(15) << "eT4_QCD" << " ";
136 fout << setw(15) << "IT4_QCD" << " ";
137 fout << setw(15) << "sT3" << " ";
138 fout << setw(15) << "sT3_QCD" << " ";
139 fout << setw(15) << "IT4_e" << " ";
140 fout << setw(15) << "IT4_mu" << " ";
141 fout << setw(15) << "IT4_tau" << " ";
142 fout << setw(15) << "nQ_QCD/T3" << " ";
143 fout << setw(15) << "npi_QCD/T3" << " ";
144 fout << setw(15) << "ne/T3" << " ";
145 fout << setw(15) << "nmu/T3" << " ";
146 fout << setw(15) << "ntau/T3" << " ";
147 fout << endl;
148
149 // On-screen
150 cout << setw(15) << "T[MeV]" << " ";
151 cout << setw(15) << "muB[MeV]" << " ";
152 cout << setw(15) << "muQ[MeV]" << " ";
153 cout << setw(15) << "mue[MeV]" << " ";
154 cout << setw(15) << "mum[MeV]" << " ";
155 cout << setw(15) << "mut[MeV]" << " ";
156 cout << setw(15) << "pion_bec" << endl;
157
158 // Timer
159 double wt1 = get_wall_time();
160 int iters = 0; // number of iterations
161
162
163 // Initial guess for the chemical potentials (muB, muQ, mue, mumu, mutau)
164 vector<double> prev = vector<double>({ 0.700, -1.e-7, -1.e-7, -1.e-7, -1.e-7 });
165
166 // Loop over the temperatures
167 for (auto&& T : Temps) {
168 // Solve the equations for the chemical potentials
169 vector<double> chems = cosmos.SolveChemicalPotentials(T, prev);
170 iters++;
171 // Will be the initial guess for the next temperature
172 prev = chems;
173
174 // Check for pion condensation
175 if (!interactingpions && abs(chems[1]) > 0.139)
176 break;
177
178 cout << setw(15) << T * 1.e3 << " ";
179 cout << setw(15) << chems[0] * 1.e3 << " ";
180 cout << setw(15) << chems[1] * 1.e3 << " ";
181 cout << setw(15) << chems[2] * 1.e3 << " ";
182 cout << setw(15) << chems[3] * 1.e3 << " ";
183 cout << setw(15) << chems[4] * 1.e3 << " ";
184 cout << setw(15) << cosmos.InPionCondensedPhase() << endl;
185
186
187 // Write to file
188 fout << setw(15) << T * 1.e3 << " ";
189 fout << setw(15) << chems[0] * 1.e3 << " ";
190 fout << setw(15) << chems[1] * 1.e3 << " ";
191 fout << setw(15) << chems[2] * 1.e3 << " ";
192 fout << setw(15) << chems[3] * 1.e3 << " ";
193 fout << setw(15) << chems[4] * 1.e3 << " ";
194 fout << setw(15) << cosmos.InPionCondensedPhase() << " ";
195
196 fout << setw(15) << cosmos.IsospinChargeDensity() / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
197
198 fout << setw(15) << cosmos.Pressure() / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
199 fout << setw(15) << cosmos.EnergyDensity() / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
200 fout << setw(15) << (cosmos.EnergyDensity() - 3. * cosmos.Pressure()) / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
201 fout << setw(15) << cosmos.PressureHRG() / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
202 fout << setw(15) << cosmos.EnergyDensityHRG() / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
203 fout << setw(15) << (cosmos.EnergyDensityHRG() - 3. * cosmos.PressureHRG()) / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
204 fout << setw(15) << cosmos.EntropyDensity() / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
205 fout << setw(15) << cosmos.EntropyDensityHRG() / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
206
207 //fout << setw(15) << (cosmos.EnergyDensity() - 3. * cosmos.Pressure()) / pow(T, 4) / pow(xMath::GeVtoifm(), 3);
208 fout << setw(15) << (cosmos.EnergyDensityChargedLepton(0) - 3. * cosmos.PressureChargedLepton(0)) / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
209 fout << setw(15) << (cosmos.EnergyDensityChargedLepton(1) - 3. * cosmos.PressureChargedLepton(1)) / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
210 fout << setw(15) << (cosmos.EnergyDensityChargedLepton(2) - 3. * cosmos.PressureChargedLepton(2)) / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
211
212 fout << setw(15) << cosmos.ElectricChargeDensityHRG() / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
213 fout << setw(15) << (cosmos.ElectricChargeDensityHRG() - cosmos.HRGModel()->ElectricChargeDensity()) / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
214 fout << setw(15) << cosmos.NetDensityChargedLepton(0) / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
215 fout << setw(15) << cosmos.NetDensityChargedLepton(1) / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
216 fout << setw(15) << cosmos.NetDensityChargedLepton(2) / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
217 fout << endl;
218 }
219
220 // Optional cross-check (run backwards in temperature and make sure the trajectory is the same)
221 if (0) {
222 for(int iT = Temps.size() - 1; 0 && iT >= 0; iT--) {
223 double T = Temps[iT];
224 vector<double> chems = cosmos.SolveChemicalPotentials(T, prev);
225 iters++;
226 prev = chems;
227
228 if (!interactingpions && abs(chems[1]) > 0.139)
229 break;
230
231 cout << setw(15) << T * 1.e3 << " ";
232 cout << setw(15) << chems[0] * 1.e3 << " ";
233 cout << setw(15) << chems[1] * 1.e3 << " ";
234 cout << setw(15) << chems[2] * 1.e3 << " ";
235 cout << setw(15) << chems[3] * 1.e3 << " ";
236 cout << setw(15) << chems[4] * 1.e3 << " ";
237 cout << setw(15) << cosmos.InPionCondensedPhase() << endl;
238
239 fout << setw(15) << T * 1.e3 << " ";
240 fout << setw(15) << chems[0] * 1.e3 << " ";
241 fout << setw(15) << chems[1] * 1.e3 << " ";
242 fout << setw(15) << chems[2] * 1.e3 << " ";
243 fout << setw(15) << chems[3] * 1.e3 << " ";
244 fout << setw(15) << chems[4] * 1.e3 << " ";
245 fout << setw(15) << cosmos.InPionCondensedPhase() << " ";
246 fout << setw(15) << cosmos.Pressure() / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
247 fout << setw(15) << cosmos.EnergyDensity() / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
248 fout << setw(15) << (cosmos.EnergyDensity() - 3. * cosmos.Pressure()) / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
249 fout << setw(15) << cosmos.PressureHRG() / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
250 fout << setw(15) << cosmos.EnergyDensityHRG() / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
251 fout << setw(15) << (cosmos.EnergyDensityHRG() - 3. * cosmos.PressureHRG()) / pow(T, 4) / pow(xMath::GeVtoifm(), 3) << " ";
252 fout << setw(15) << cosmos.EntropyDensity() / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
253 fout << setw(15) << cosmos.EntropyDensityHRG() / pow(T, 3) / pow(xMath::GeVtoifm(), 3) << " ";
254 fout << endl;
255 }
256 }
257
258 double wt2 = get_wall_time();
259
260 printf("%30s %lf s\n", "Running time:", (wt2 - wt1));
261 printf("%30s %lf s\n", "Time per single calculation:", (wt2 - wt1) / iters);
262
263 return 0;
264}
265
main
int main(int argc, char *argv[])
Definition CosmicTrajectory.cpp:33
thermalfist::CosmicEoS
Class implementing cosmological equation of state as a mixture of HRG model, ideal gases of leptons a...
Definition CosmicEoS.h:57
thermalfist::CosmicEoS::PressureHRG
double PressureHRG()
Calculates the partial pressure of the HRG part.
Definition CosmicEoS.cpp:160
thermalfist::CosmicEoS::SolveChemicalPotentials
std::vector< double > SolveChemicalPotentials(double T, const std::vector< double > &muInit=std::vector< double >())
Calculates the values of the chemical potential (B,Q,{L}) that satisfy the given asymmetry constraint...
Definition CosmicEoS.cpp:287
thermalfist::CosmicEoS::EnergyDensity
double EnergyDensity()
Calculates the total energy density.
Definition CosmicEoS.cpp:167
thermalfist::CosmicEoS::EntropyDensityHRG
double EntropyDensityHRG()
Calculates the entropy density of the HRG part.
Definition CosmicEoS.cpp:134
thermalfist::CosmicEoS::HRGModel
ThermalModelBase * HRGModel() const
Gets the pointer to the HRG model object.
Definition CosmicEoS.h:78
thermalfist::CosmicEoS::PressureChargedLepton
double PressureChargedLepton(int iL)
Calculates the partial pressure of charged lepton flavor.
Definition CosmicEoS.cpp:203
thermalfist::CosmicEoS::EnergyDensityHRG
double EnergyDensityHRG()
Calculates the energy density of the HRG part.
Definition CosmicEoS.cpp:185
thermalfist::CosmicEoS::InPionCondensedPhase
bool InPionCondensedPhase() const
Checks if the system has non-zero BEC of pions.
Definition CosmicEoS.cpp:338
thermalfist::CosmicEoS::ElectricChargeDensityHRG
double ElectricChargeDensityHRG(bool absolute=false)
Calculates the electric charge density of the HRG part.
Definition CosmicEoS.cpp:245
thermalfist::CosmicEoS::EnergyDensityChargedLepton
double EnergyDensityChargedLepton(int iL)
Calculates the energy density of charged lepton flavor.
Definition CosmicEoS.cpp:212
thermalfist::CosmicEoS::NetDensityChargedLepton
double NetDensityChargedLepton(int iL)
Calculates the net density of charged lepton flavor.
Definition CosmicEoS.cpp:192
thermalfist::CosmicEoS::SetAsymmetries
void SetAsymmetries(const std::vector< double > &asymmetries)
Sets all asymmetries at once.
Definition CosmicEoS.h:194
thermalfist::CosmicEoS::Pressure
double Pressure()
Calculates the total pressure.
Definition CosmicEoS.cpp:141
thermalfist::CosmicEoS::EntropyDensity
double EntropyDensity()
Calculates the total entropy density.
Definition CosmicEoS.cpp:116
thermalfist::CosmicEoS::IsospinChargeDensity
double IsospinChargeDensity(bool absolute=false)
Calculates the isospin charge density.
Definition CosmicEoS.cpp:254
thermalfist::ThermalModelBase::SetUseWidth
void SetUseWidth(bool useWidth)
Sets whether finite resonance widths are used. Deprecated.
Definition ThermalModelBase.cpp:98
thermalfist::ThermalModelBase::SetStatistics
virtual void SetStatistics(bool stats)
Definition ThermalModelBase.cpp:341
thermalfist::ThermalModelIdeal
Class implementing the Ideal HRG model.
Definition ThermalModelIdeal.h:20
thermalfist::ThermalParticle::eBWconstBR
@ eBWconstBR
Energy-dependent Breit-Wigner scheme (eBW) with constant branching ratios when evaluating feeddown.
Definition ThermalParticle.h:88
thermalfist::ThermalParticleSystem
Class containing the particle list.
Definition ThermalParticleSystem.h:30
thermalfist::xMath::GeVtoifm
constexpr double GeVtoifm()
A constant to transform GeV into fm .
Definition xMath.h:25
thermalfist
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition CosmicEoS.h:9
thermalfist::get_wall_time
double get_wall_time()
Definition Utility.cpp:191
Generated by doxygen 1.13.2