PCE-Saha-LHC.cpp

Go to the documentation of this file.

/*
 * Thermal-FIST package
 * 
 * Copyright (c) 2014-2020 Volodymyr Vovchenko
 *
 * GNU General Public License (GPLv3 or later)
 */
#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 "HRGPCE.h"

#include "ThermalFISTConfig.h"

using namespace std;

#ifdef ThermalFIST_USENAMESPACE
using namespace thermalfist;
#endif

// This is an example of doing PCE-HRG model calculations at the LHC energies using Thermal-FIST
// Usage: PCE-Saha-LHC
int main(int argc, char *argv[])
{
    // The default particle list. As of version 1.3 this is PDG2020 list including light nuclei
    ThermalParticleSystem parts(ThermalFIST_DEFAULT_LIST_FILE);

    // To reproduce arXiv:1903.10024 use the PDG2014 list
    //ThermalParticleSystem TPS(string(ThermalFIST_INPUT_FOLDER) + "/list/PDG2014/list-withnuclei.dat"); 

    // Use ideal HRG model
    ThermalModelIdeal model(&parts);
    
    // PCE-HRG model
    ThermalModelPCE modelpce(&model);
    modelpce.UseSahaForNuclei(true); // Light nuclei evaluated using the Saha equation (arXiv:1903.10024)
    modelpce.FreezeLonglivedResonances(false); // All strongly decaying resonance are in partial equilibrium

    // Chemical freeze-out conditions: 2.76 TeV 0-10% Pb-Pb collisions
    ThermalModelParameters params_chemical_freezeout;
    params_chemical_freezeout.T = 0.155; // Temperature in GeV
    params_chemical_freezeout.muB = 0.;
    params_chemical_freezeout.V = 4700.; // Volume in fm^3

    model.SetParameters(params_chemical_freezeout);
    model.FillChemicalPotentials(); // Fills chemical potentials for all species at Tch

    // Set the chemical freeze-out as an "initial" condition for PCE
    modelpce.SetChemicalFreezeout(params_chemical_freezeout);

    // The list of chemical potentials for output, coded by the pdg code
    vector<long long> pdgcodes_stable;
    pdgcodes_stable.push_back(211);  // pions (pi+)
    pdgcodes_stable.push_back(321);  // kaons (K+)
    pdgcodes_stable.push_back(2212); // protons (p+)
    pdgcodes_stable.push_back(3122); // Lambdas
    pdgcodes_stable.push_back(3222); // Sigma+
    pdgcodes_stable.push_back(3312); // Xi-
    pdgcodes_stable.push_back(3334); // Omega

    // The list of yield ratios to output
    vector< pair<long long, long long> > ratios;
    // First the nuclei
    ratios.push_back(make_pair(1000010020, 2212));  // d/p
    ratios.push_back(make_pair(1000020030, 2212)); // He3/p
    ratios.push_back(make_pair(1000010030, 2212)); // H3/p
    ratios.push_back(make_pair(1000020040, 2212)); // He4/p
    ratios.push_back(make_pair(1010010030, 2212)); // Hypertriton/p
    ratios.push_back(make_pair(1010010040, 2212)); // HyperHydrogen4/p
    // Now the resonances
    ratios.push_back(make_pair(313, -321));    // K^*0 / K^-
    ratios.push_back(make_pair(113, 211));     // rho^0/ pi^+
    ratios.push_back(make_pair(3124, 3122));   // \Lambda(1520)/\Lambda
    ratios.push_back(make_pair(9010221, 211)); // f0(980) / pi^+
    ratios.push_back(make_pair(2224, 2212));   // \Delta(1232)++/p

    // Preparing the output files
    // The file to output the parameters (volume, entropy, chemical potentials)
    FILE* fout_params = fopen("PCE.LHC.Parameters.dat", "w");
    fprintf(fout_params, "%15s %15s %15s ", "T[MeV]", "V/Vch", "S/Sch");
    for (int i = 0; i < pdgcodes_stable.size(); ++i) {
        fprintf(fout_params, "%15s ", ("mu_" + string(parts.ParticleByPDG(pdgcodes_stable[i]).Name())).c_str());
    }
    fprintf(fout_params, "\n");

    // The file to output the yield ratios
    FILE* fout_ratios = fopen("PCE.LHC.Ratios.dat", "w");
    fprintf(fout_ratios, "%15s ", "T[MeV]");
    for (int i = 0; i < ratios.size(); ++i) {
        fprintf(fout_ratios, "%15s ", (parts.ParticleByPDG(ratios[i].first).Name() + "/" + parts.ParticleByPDG(ratios[i].second).Name()).c_str());
    }
    fprintf(fout_ratios, "\n");

    // The temperature scan
    double T0 = params_chemical_freezeout.T;
    double dT = 0.001;   // steps of 1 MeV
    double Tmin = 0.070; // Down to 70 MeV

    // Store the value of the total entropy at the chemical freeze-out
    double entropy_chemical_freezeout = modelpce.ThermalModel()->EntropyDensity() * params_chemical_freezeout.V;

    // Loop over temperatures
    for (double T = T0; T >= Tmin - 1.e-9; T -= dT) {
        printf("T = %lf MeV\n", T * 1.e3);

        // Compute the PCE chemical potentials at a given temperature
        modelpce.CalculatePCE(T);

        // Output the parameters at the current temperature
        fprintf(fout_params, "%15lf %15lf %15lf ",
            T * 1.e3,
            modelpce.ThermalModel()->Volume() / params_chemical_freezeout.V,
            modelpce.ThermalModel()->EntropyDensity() * modelpce.ThermalModel()->Volume() / entropy_chemical_freezeout
            );

        for (int i = 0; i < pdgcodes_stable.size(); ++i) {
            fprintf(fout_params, "%15lf ",
                modelpce.ChemicalPotentials()[ parts.PdgToId(pdgcodes_stable[i]) ]
            );
        }
        fprintf(fout_params, "\n");

        // Output the yield ratios at the current temperature
        fprintf(fout_ratios, "%15lf ", T * 1.e3);
        for (int i = 0; i < ratios.size(); ++i) {
            fprintf(fout_ratios, "%15E ",
                modelpce.ThermalModel()->GetYield(ratios[i].first, Feeddown::Electromagnetic) /
                modelpce.ThermalModel()->GetYield(ratios[i].second, Feeddown::Electromagnetic));
        }
        fprintf(fout_ratios, "\n");
        
    }

    fclose(fout_params);
    fclose(fout_ratios);

    return 0;
}