Thermal-FIST  1.3
Package for hadron resonance gas model applications
CalculationTmu.cpp
Go to the documentation of this file.
1 /*
2  * Thermal-FIST package
3  *
4  * Copyright (c) 2014-2018 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 
15 #include "HRGBase.h"
16 #include "HRGEV.h"
17 #include "HRGFit.h"
18 #include "HRGVDW.h"
19 
20 #include "ThermalFISTConfig.h"
21 
22 using namespace std;
23 
24 #ifdef ThermalFIST_USENAMESPACE
25 using namespace thermalfist;
26 #endif
27 
28 // Calculates equation of state, particle yields and fluctuations at a given set of T-mu values
29 // Usage: CalculationTmu
30 int main(int argc, char *argv[])
31 {
32  int ModelType = 1; // 0 - Ideal HRG, 1 - QvdW HRG (from 1609.03975)
33  if (argc > 1)
34  ModelType = atoi(argv[1]);
35 
36  std::string prefix = "QvdW-HRG";
37  if (ModelType != 1)
38  prefix = "IdealHRG";
39 
40  // Fill the T-mu values where calculations should be performed
41  vector<double> Tvalues, muvalues;
42 
43  // Here done by hand
44  // Alternatively one can read those from external file, or populate in a loop, etc.
45  // Note that all energy units are in GeV!
46  // 1
47  Tvalues.push_back(0.100); muvalues.push_back(0.600);
48  // 2
49  Tvalues.push_back(0.130); muvalues.push_back(0.500);
50  // 3
51  Tvalues.push_back(0.160); muvalues.push_back(0.000);
52 
53 
54  // Create the hadron list instance and read the list from file
55 
56 
57  //ThermalParticleSystem TPS(string(ThermalFIST_INPUT_FOLDER) + "/list/thermus23mod/list.dat"); // <-- modified THERMUS-2.3 list
58  ThermalParticleSystem TPS(string(ThermalFIST_INPUT_FOLDER) + "/list/PDG2014/list.dat"); // <-- Default list, no light nuclei
59  //ThermalParticleSystem TPS(string(ThermalFIST_INPUT_FOLDER) + "/list/PDG2014/list-withnuclei.dat"); // <-- Default list, with light nuclei
60 
61  // Create the ThermalModel instance
62  // Choose the class which fits the required variant of HRG model
63  //ThermalModelIdeal model(&TPS);
64  //ThermalModelCanonical model(&TPS);
65  ThermalModelBase *model;
66  if (ModelType == 1) // QvdW-HRG
67  {
68  model = new ThermalModelVDWFull(&TPS);
69 
70  // Set the QvdW interaction parameters
71  // As in 1609.03975, here we consider (anti)baryon-(anti)baryon interactions only
72  double a = 0.329;
73  double b = 3.42;
74 
75  // Loop over all hadron-hadron pairs to set a and b for each of these pairs
76  for (int i1 = 0; i1 < model->TPS()->Particles().size(); ++i1) {
77  for (int i2 = 0; i2 < model->TPS()->Particles().size(); ++i2) {
78  const ThermalParticle &part1 = model->TPS()->Particles()[i1];
79  const ThermalParticle &part2 = model->TPS()->Particles()[i2];
80 
81  int B1 = part1.BaryonCharge();
82  int B2 = part2.BaryonCharge();
83 
84  // Or use pdgid's to identify the two particles
85  //int pdgid1 = part1.PdgId();
86  //int pdgid2 = part2.PdgId();
87 
88  // Meson-meson, meson-baryon, baryon-antibaryon non-interacting
89  if (!(B1 * B2 > 0)) {
90 
91  model->SetVirial(i1, i2, 0.); // No repulsion
92 
93  model->SetAttraction(i1, i2, 0.); // No attraction
94  continue;
95  }
96  else {
97  // BB excluded volume
98  model->SetVirial(i1, i2, b);
99 
100  // BB attraction
101  model->SetAttraction(i1, i2, a);
102  }
103  }
104  }
105  }
106  else {
107  model = new ThermalModelIdeal(&TPS);
108  }
109 
110  // Use (or not) finite resonance width
111  model->SetUseWidth(true);
112 
113  // Include (or not) quantum statistics
114  model->SetStatistics(true);
115 
116  // Output, here on screen, to write into file use, e.g., fprintf
117  printf("%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s\n",
118  "T[GeV]", "muB[GeV]",
119  "P[GeV/fm3]", "e[GeV/fm3]", "s[fm-3]",
120  "<K+>", "<pi+>", "<K+>/<pi+>",
121  "w[K+]", "w[pi+]",
122  "<N->", "w[N-]",
123  "chi3B/chi2B", "chi4B/chi2B",
124  "chi3Q/chi2Q", "chi4Q/chi2Q",
125  "chi3S/chi2S", "chi4S/chi2S");
126 
127  // The same output to file
128  std::string filename = prefix + ".CalculationTmu.dat";
129  FILE *f = fopen(filename.c_str(), "w");
130  fprintf(f, "%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s%15s\n",
131  "T[GeV]", "muB[GeV]",
132  "P[GeV/fm3]", "e[GeV/fm3]", "s[fm-3]",
133  "<K+>", "<pi+>", "<K+>/<pi+>",
134  "w[K+]", "w[pi+]",
135  "<N->", "w[N-]",
136  "chi3B/chi2B", "chi4B/chi2B",
137  "chi3Q/chi2Q", "chi4Q/chi2Q",
138  "chi3S/chi2S", "chi4S/chi2S");
139 
140  // Iterate over all the T-muB pair values
141  for (int i = 0; i < Tvalues.size(); ++i) {
142  double T = Tvalues[i];
143  double muB = muvalues[i];
144 
145 
146  // Set temperature and baryon chemical potential
147  model->SetTemperature(T);
148  model->SetBaryonChemicalPotential(muB);
149 
150  // Constrain muB from strangeness neutrality condition
151  model->ConstrainMuS(true);
152  // Alternatively set the muS value manually
153  // model->ConstrainMuS(false);
154  // model->SetStrangenessChemicalPotential(0.);
155 
156  // Constrain muq from Q/B = 0.4 condition
157  model->ConstrainMuQ(true);
158  model->SetQoverB(0.4);
159  // Alternatively set the muQ value manually
160  //model->ConstrainMuQ(false);
161  //model->SetElectricChemicalPotential(0.);
162 
163  // Chemical non-equilbrium parameters
164  model->SetGammaq(1.);
165  model->SetGammaS(1.);
166 
167  // Set volume
168  model->SetVolumeRadius(3.); // System radius R in fm, volume is V = (4/3) * \pi * R^3
169  //model->SetVolume(5000.); //<-- Alternative, volume V in fm^3
170 
171 
172  // Determine muS and/or muQ from constraints, if there are any
173  model->FixParameters();
174 
175  // Calculate all hadron densities, both primordial and final
176  model->CalculateDensities();
177 
178  // Calculate fluctuations
179  model->CalculateFluctuations();
180 
181  // Equation of state parameters
182  double p = model->CalculatePressure(); // Pressure in GeV/fm3
183  double e = model->CalculateEnergyDensity(); // Energy density in GeV/fm3
184  double s = model->CalculateEntropyDensity(); // Entropy density in fm-3
185 
186  // Calculate final yields,
187  // Usage: model->GetDensity(pdgid, feeddown)
188  // pdgid -- PDG code for the desired particle species
189  // feeddown: 0 - primordial, 1 - final, 2 - final with additional feeddown from weak decays
190  // yield = density * volume
191  double yieldKplus = model->GetDensity(321, Feeddown::StabilityFlag) * model->Volume();
192  double yieldpiplus = model->GetDensity(211, Feeddown::StabilityFlag) * model->Volume();
193 
194  // Scaled variance of final state particle number fluctuations
195  double wKplus = model->ScaledVarianceTotal( model->TPS()->PdgToId(321) );
196  double wpiplus = model->ScaledVarianceTotal( model->TPS()->PdgToId(211) );
197 
198  // Charged particle mean multplicities, after decays
199  // Argument: 0 - all charged, 1 - positively charged, 2 - negatively charged
200  double Nch = model->ChargedMultiplicityFinal(0);
201  double Nplus = model->ChargedMultiplicityFinal(1);
202  double Nminus = model->ChargedMultiplicityFinal(2);
203 
204  // Scaled variance for charged particle multplicity distribution, after decays
205  double wNch = model->ChargedScaledVarianceFinal(0);
206  double wNplus = model->ChargedScaledVarianceFinal(1);
207  double wNminus = model->ChargedScaledVarianceFinal(2);
208 
209  // Higher-order fluctuations of conserved charges B, Q, S
210 
211  // Array of charges of all particles in the list
212  // E.g., if the ith particle has baryon charge 1, then chargesB[i] = 1 etc.
213  vector<double> chargesB(model->Densities().size()), chargesQ(model->Densities().size()), chargesS(model->Densities().size());
214 
215  // Array with the values of the calculated susceptibilities
216  vector<double> chchis;
217 
218  // Baryon number
219  for (int i = 0; i < model->TPS()->Particles().size(); ++i) {
220  chargesB[i] = model->TPS()->Particles()[i].BaryonCharge();
221  }
222 
223  // Calculation of susceptibilities chi1-chi4
224  chchis = model->CalculateChargeFluctuations(chargesB, 4);
225  double chi2B = chchis[1];
226  double chi3B = chchis[2];
227  double chi4B = chchis[3];
228 
229  // Electric charge, same procedure
230  for (int i = 0; i < model->TPS()->Particles().size(); ++i) {
231  chargesQ[i] = model->TPS()->Particles()[i].ElectricCharge();
232  }
233 
234  chchis = model->CalculateChargeFluctuations(chargesQ, 4);
235  double chi2Q = chchis[1];
236  double chi3Q = chchis[2];
237  double chi4Q = chchis[3];
238 
239  // Strangeness
240  for (int i = 0; i < model->TPS()->Particles().size(); ++i) {
241  chargesS[i] = model->TPS()->Particles()[i].Strangeness();
242  }
243 
244  chchis = model->CalculateChargeFluctuations(chargesS, 4);
245  double chi2S = chchis[1];
246  double chi3S = chchis[2];
247  double chi4S = chchis[3];
248 
249 
250  printf("%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf\n",
251  T,
252  muB,
253  p,
254  e,
255  s,
256  yieldKplus,
257  yieldpiplus,
258  yieldKplus/yieldpiplus,
259  wKplus,
260  wpiplus,
261  Nminus,
262  wNminus,
263  chi3B / chi2B,
264  chi4B / chi2B,
265  chi3Q / chi2Q,
266  chi4Q / chi2Q,
267  chi3S / chi2S,
268  chi4S / chi2S);
269 
270  fprintf(f, "%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf\n",
271  T,
272  muB,
273  p,
274  e,
275  s,
276  yieldKplus,
277  yieldpiplus,
278  yieldKplus / yieldpiplus,
279  wKplus,
280  wpiplus,
281  Nminus,
282  wNminus,
283  chi3B / chi2B,
284  chi4B / chi2B,
285  chi3Q / chi2Q,
286  chi4Q / chi2Q,
287  chi3S / chi2S,
288  chi4S / chi2S);
289  }
290 
291  fclose(f);
292 
293  delete model;
294 
295  return 0;
296 }
297 
Abstract base class for an HRG model implementation.
const std::vector< double > & Densities() const
virtual void SetGammaS(double gammaS)
Set the strange quark fugacity factor.
virtual void SetBaryonChemicalPotential(double muB)
Set the baryon chemical potential.
virtual void SetStatistics(bool stats)
virtual void FixParameters()
Method which actually implements ConstrainChemicalPotentials() (for backward compatibility).
void SetUseWidth(bool useWidth)
Sets whether finite resonance widths are used. Deprecated.
int PdgToId(long long pdgid)
Transforms PDG ID to a 0-based particle id number.
Class containing the particle list.
virtual void SetTemperature(double T)
Set the temperature.
virtual void SetGammaq(double gammaq)
Set the light quark fugacity factor.
void SetQoverB(double QB)
The electric-to-baryon charge ratio to be used to constrain the electric chemical potential...
virtual void SetVirial(int, int, double)
Set the excluded volume coefficient .
double ChargedMultiplicityFinal(int type=0)
Multiplicity of charged particles including the feeddown contributions in accordance with the stabili...
int BaryonCharge() const
Particle&#39;s baryon number.
Class containing all information about a particle specie.
double ChargedScaledVarianceFinal(int type=0)
Scaled variance of charged particles including the feeddown contributions in accordance with the stab...
double GetDensity(long long PDGID, int feeddown)
Same as GetDensity(int,Feeddown::Type)
const std::vector< ThermalParticle > & Particles() const
Returns the vector of all particle species.
virtual double CalculateEnergyDensity()=0
virtual double CalculateEntropyDensity()=0
virtual void CalculateDensities()
Calculates the primordial and total (after decays) densities of all species.
virtual double CalculatePressure()=0
Implementation of the equation of state functions.
virtual void CalculateFluctuations()
Computes the fluctuation observables.
virtual std::vector< double > CalculateChargeFluctuations(const std::vector< double > &chgs, int order=4)
Calculates fluctuations (diagonal susceptibilities) of an arbitrary "conserved" charge.
virtual void SetAttraction(int, int, double)
Set the vdW mean field attraction coefficient .
Class implementing the Ideal HRG model.
double ScaledVarianceTotal(int id) const
Scaled variance of final particle number fluctuations for particle species id.
ThermalModelVDW ThermalModelVDWFull
For backward compatibility.
int main(int argc, char *argv[])
The main namespace where all classes and functions of the Thermal-FIST library reside.
double Volume() const
System volume (fm )
ThermalParticleSystem * TPS()
void SetVolumeRadius(double R)
Sets the system radius.