cpc2-chi2-vs-T.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 the fit to ALICE 2.76 TeV data, 0-5% centrality, as in 1512.08046
// Four variants of HRG model: 
// 1. Ideal HRG: <config> = 0
// 2. Diagonal EV-HRG with bag model parametrization r = r_p * (m/m_p)^1/3, where r_p = 0.5 is proton radius parameter (as in 1512.08046): <config> = 1
// 3. Diagonal EV-HRG with constant radius parameter r = 0.3 fm for all baryons and r = 0 for all mesons (as in 1201.0693): <config> = 2
// 4. QvdW-HRG with a and b for baryons only, fixed to nuclear ground state (as in 1609.03975): <config> = 3
// 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("#Fitting 2.76 TeV ALICE data at \\mu = 0 in Id-HRG model\n");
 
    modeltype = "Id-HRG";
  }
  else if (config == 1) // EV-HRG, r = 0.3 fm for baryons, r = 0 for mesons, no light nuclei, as in 1512.08046
  {
    model = new ThermalModelEVDiagonal(&TPS);
    double rad = 0.3;
    for (int i = 0; i < model->TPS()->ComponentsNumber(); ++i) {
      if (model->TPS()->Particle(i).BaryonCharge() == 0) model->SetRadius(i, 0.);
      else model->SetRadius(i, rad);
    }
 
    printf("#Fitting 2.76 TeV ALICE data at \\mu = 0 in EV-HRG model with r = %lf fm for baryons, and r = 0 for mesons\n", rad);
 
    modeltype = "EV-HRG-TwoComponent";
  }
  else if (config == 2) // EV-HRG, Bag Model with r_p = 0.5 fm, as in 1512.08046
  {
    model = new ThermalModelEVDiagonal(&TPS);
    double radProton = 0.5;
    for (int i = 0; i < model->TPS()->ComponentsNumber(); ++i) {
      model->SetRadius(i, radProton * pow(model->TPS()->Particle(i).Mass() / 0.938, 1/3.));
    }
 
    printf("#Fitting 2.76 TeV ALICE data at \\mu = 0 in Bag Model EV-HRG model with proton r = %lf fm\n", radProton);
 
    modeltype = "EV-HRG-BagModel";
  }
  else if (config == 3) // 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 repulsion
          model->SetVirial(i, j, b);     // Set QvdW attraction
        }
        else // no vdW interactions for meson-meson, meson-baryon or baryon-anitbaryon pairs
        {
          model->SetAttraction(i, j, 0.);
          model->SetVirial(i, j, 0.);
        }
      }
    }
 
    printf("#Fitting 2.76 TeV ALICE data 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 fitter
  ThermalModelFit fitter(model);
 
  // By default T, muB, and R parameters are fitted, while others (gammaS etc.) are fixed
  // Initial parameters values are taken from those currently set in the ThermalModel object
  // Here we do not fit muB, which is set to zero above
  fitter.SetParameterFitFlag("muB", false); // Do not fit muB
 
  
 
  // R is fitted by default
  // We can still specify the initial value, the initial delta used by minuit,
  // and the lower and upper limits using the SetParameter function
  double Rinit  = 10.0;
  double Rdelta = 1.0;
  double Rmin   = 0.0;
  double Rmax   = 30.0;
  fitter.SetParameter("R", Rinit, Rdelta, Rmin, Rmax);
 
  // Load the experimental data
  vector<FittedQuantity> quantities = ThermalModelFit::loadExpDataFromFile(string(ThermalFIST_INPUT_FOLDER) + "/data/ALICE-PbPb2.76TeV-0-5-1512.08046.dat");
  fitter.SetQuantities(quantities);
 
 
  printf("%15s%15s%15s%15s\n",
    "T[MeV]",   // Temperature in MeV
    "R[fm]",    // System radius in fm
    "chi2",     // chi_2
    "chi2_dof"  // Reduced chi2
  );
 
  
  // Prepare for output to file
  char tmpc[1000];
  sprintf(tmpc, "cpc2.%s.ALICE2_76.chi2.TDep.out", modeltype.c_str());
  FILE *fout = fopen(tmpc, "w");
 
  fprintf(fout, "%15s%15s%15s%15s\n",
    "T[MeV]",   // Temperature in MeV
    "R[fm]",    // System radius in fm
    "chi2",     // chi_2
    "chi2_dof"  // Reduced chi2
  );
 
  double wt1 = get_wall_time(); // Timing
 
  int iters = 0; // Number of data points
 
  // Temperature interval, in GeV
  double Tmin = 0.100;
  double Tmax = 0.2501;
  double dT   = 0.002;
 
  if (config == 0)
    dT = 0.001;
 
  if (config == 2) {
    dT = 0.005;
    Tmax = 0.4001;
  }
 
  for (double T = Tmin; T <= Tmax; T += dT) {
    // We also do not fit T, but fix it at each iteration to a given value
    fitter.SetParameterFitFlag("T", false);
    fitter.SetParameterValue("T", T); // Set the temperature
 
    ThermalModelFitParameters result = fitter.PerformFit(false);  // We still have to fit the radius, the argument suppresses the output during minimization  
 
    double Rfit = result.R.value;
    double chi2 = result.chi2;
    double chi2dof = result.chi2ndf;
 
    printf("%15lf%15lf%15lf%15lf\n", T * 1000., Rfit, chi2, chi2 / (result.ndf - 1.));
 
    fprintf(fout, "%15lf%15lf%15lf%15lf\n", T * 1000., Rfit, chi2, chi2 / (result.ndf - 1.));
 
    iters++;
 
    fflush(stdout);
 
  }
 
  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;
}