cpc1-HRG-TDep.cpp

Go to the documentation of this file.

#include <string.h>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <ctime>
#include <cstdio>
 
#include "HRGBase.h"
#include "HRGEV.h"
#include "HRGFit.h"
#include "HRGVDW.h"
 
#include "ThermalFISTConfig.h"
 
using namespace std;
 
#ifdef ThermalFIST_USENAMESPACE
using namespace thermalfist;
#endif
 
 
// Temperature dependence of HRG thermodynamics at \mu = 0
// Three variants of the HRG model: 
// 1. Ideal HRG: <config> = 0
// 2. EV-HRG with constant radius parameter r = 0.3 fm for all hadrons (as in 1412.5478): <config> = 1
// 3. QvdW-HRG with a and b for baryons only, fixed to nuclear ground state (as in 1609.03975): <config> = 2
// Usage: cpc1HRGTDep <config>
int main(int argc, char *argv[])
{
  // Particle list file
  // Here we will use the list from THERMUS-2.3, for comparing the results with THERMUS-2.3
  string listname = string(ThermalFIST_INPUT_FOLDER) + "/list/thermus23/list.dat";
 
  // Alternative: use the default PDG2014 list
  //string listname = string(ThermalFIST_INPUT_FOLDER) + "/list/PDG2014/list.dat";
  
  // Create the hadron list instance and read the list from file
  ThermalParticleSystem TPS(listname);
 
  // Which variant of the HRG model to use
  int config = 0;
 
  // Read config from command line
  if (argc > 1)
    config = atoi(argv[1]);
 
  
  string modeltype; // For output
  
 
  // Pointer to the thermal model instance used in calculations
  ThermalModelBase *model;
 
  
 
  if (config == 0) // Ideal HRG
  {
    model = new ThermalModelIdeal(&TPS);
    
    printf("#Calculating thermodynamics at \\mu = 0 in Id-HRG model\n");
 
    modeltype = "Id-HRG";
  }
  else if (config == 1) // EV-HRG, r = 0.3 fm, to reproduce 1412.5478
  {
    model = new ThermalModelEVDiagonal(&TPS);
    // Set r = 0.3 fm for each hadron in the list
    double rad = 0.3;
    for (int i = 0; i < model->TPS()->ComponentsNumber(); ++i)
      model->SetRadius(i, rad);
 
    printf("#Calculating thermodynamics at \\mu = 0 in EV-HRG model with r = %lf fm\n", rad);
 
    modeltype = "EV-HRG";
  }
  else if (config == 2) // QvdW-HRG, to reproduce 1609.03975
  {
    model = new ThermalModelVDWFull(&TPS);
 
    // vdW parameters, for baryon-baryon, antibaryon-antibaryon ONLY, otherwise zero
    double a = 0.329; // In GeV*fm3
    double b = 3.42;  // In fm3
 
    // Iterate over all pairs of hadron species and set a and b
    for (int i = 0; i < model->TPS()->ComponentsNumber(); ++i) {
      for (int j = 0; j < model->TPS()->ComponentsNumber(); ++j) {
        int B1 = model->TPS()->Particle(i).BaryonCharge(); // Baryon number of 1st species
        int B2 = model->TPS()->Particle(j).BaryonCharge(); // Baryon number of 2nd species
 
        if ((B1 > 0 && B2 > 0) || (B1 < 0 && B2 < 0))  // baryon-baryon or antibaryon-antibaryon
        {
          model->SetAttraction(i, j, a); // Set QvdW attraction
          model->SetVirial(i, j, b);     // Set QvdW repulsion
        }
        else // no vdW interactions for meson-meson, meson-baryon or baryon-anitbaryon pairs
        {
          model->SetAttraction(i, j, 0.);
          model->SetVirial(i, j, 0.);
        }
      }
    }
 
    printf("#Calculating thermodynamics at \\mu = 0 in QvdW-HRG model\n");
 
    modeltype = "QvdW-HRG";
  }
  else // Ideal HRG by default
  {
    model = new ThermalModelIdeal(&TPS);
 
    modeltype = "Id-HRG";
  }
 
  // Use quantum statistics
  model->SetStatistics(true);
  //model->SetStatistics(false);
 
  // Use mass integration over Breit-Wigner shapes in +-2Gamma interval, as in THERMUS
  model->SetUseWidth(ThermalParticle::BWTwoGamma);
  //model->SetUseWidth(ThermalParticle::ZeroWidth);
 
  // Set chemical potentials to zero
  model->SetBaryonChemicalPotential(0.0);
  model->SetElectricChemicalPotential(0.0);
  model->SetStrangenessChemicalPotential(0.0);
  model->SetCharmChemicalPotential(0.0);
  model->FillChemicalPotentials();
 
  
  // Prepare output
  char tmpc[1000];
  sprintf(tmpc, "cpc1.%s.TDep.out", modeltype.c_str());
  FILE *fout = fopen(tmpc, "w");
 
  // Output to screen
  printf("%15s%15s%15s%15s%15s%15s%15s\n",
    "T[MeV]", // Temperature
    "p/T^4",  // Scaled pressure
    "e/T^4",  // Scaled energy density
    "s/T^3",   // Scaled entropy density
    "chi2B",   // Baryon number susceptibility
    "chi4B",   // 4th order baryon cumulant
    "chi2B-chi4B"   // Difference of 2nd and 4th order baryon susceptibilities
  );
 
  // Output to file
  fprintf(fout, "%15s%15s%15s%15s%15s%15s%15s\n",
    "T[MeV]", // Temperature
    "p/T^4",  // Scaled pressure
    "e/T^4",  // Scaled energy density
    "s/T^3",   // Scaled entropy density
    "chi2B",   // Baryon number susceptibility
    "chi4B",   // 4th order baryon cumulant
    "chi2B-chi4B"   // Difference of 2nd and 4th order baryon susceptibilities
  );
 
 
  double wt1 = get_wall_time(); // Timing
 
  int iters = 0; // Number of data points
 
  // Temperature interval, in GeV
  double Tmin = 0.020;
  double Tmax = 0.2001;
  double dT   = 0.001;
 
  for (double T = Tmin; T <= Tmax; T += dT) {
    model->SetTemperature(T);
 
    // Calculates densities, solves all necessary transcendental equations, if necessary
    model->CalculateDensities();
 
    // Output temperature, scale to get the unit of MeV
    printf("%15lf", T * 1000.);
    fprintf(fout, "%15lf", T * 1000.);
 
    // Pressure, in [GeV/fm3]
    double p = model->CalculatePressure();
    // Scaled pressure, pow(xMath::GeVtoifm(), 3) needed to make it dimensionless
    double pT4 = p / pow(T, 4) / pow(xMath::GeVtoifm(), 3);
    printf("%15lf", pT4);
    fprintf(fout, "%15lf", pT4);
 
    // Energy density
    double e = model->CalculateEnergyDensity();
    double eT4 = e / pow(T, 4) / pow(xMath::GeVtoifm(), 3);
    printf("%15lf", eT4);
    fprintf(fout, "%15lf", eT4);
 
    // Entropy density
    double s = model->CalculateEntropyDensity();
    double sT3 = s / pow(T, 3) / pow(xMath::GeVtoifm(), 3);
    printf("%15lf", sT3);
    fprintf(fout, "%15lf", sT3);
 
 
    // Baryon number fluctuations
 
    // Vector containing baryon charge of each species
    vector<double> baryon_charges;
    for (int i = 0; i < model->TPS()->ComponentsNumber(); ++i) {
      baryon_charges.push_back(static_cast<double>(model->TPS()->Particle(i).BaryonCharge()));
    }
 
    // Calculate baryon number cumulants up to 4th order
    vector<double> chiB = model->CalculateChargeFluctuations(baryon_charges, 4);
 
    if (chiB.size() >= 4) {
      // chi2B
      printf("%15E", chiB[1]);
      fprintf(fout, "%15E", chiB[1]);
 
      // chi4B
      printf("%15E", chiB[3]);
      fprintf(fout, "%15E", chiB[3]);
 
      // chi2B - chi4B
      printf("%15E", chiB[1] - chiB[3]);
      fprintf(fout, "%15E", chiB[1] - chiB[3]);
    }
 
    printf("\n");
    fprintf(fout, "\n");
 
    iters++;
 
  }
 
  fclose(fout);
 
  delete model;
  
 
  double wt2 = get_wall_time();
 
  printf("%30s %lf s\n", "Running time:", (wt2 - wt1));
  printf("%30s %lf s\n", "Time per single calculation:", (wt2 - wt1) / iters);
 
  return 0;
}