1808.05245.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Canonical statistical model for light nuclei at the LHC\n",
    "\n",
    "This notebook implements the canonical statistical model calculations of various light (anti-)nuclei-to-proton ratios as a function charged pion multiplicity at the LHC, as reported in the paper:\n",
    "\n",
    "- V. Vovchenko, B. Dönigus, H. Stoecker, *Multiplicity dependence of light nuclei production at LHC energies in the canonical statistical model*, [Phys. Lett. B **785**, 171 (2018)](https://doi.org/10.1016/j.physletb.2018.08.041), [arXiv:1808.05245 [hep-ph]](https://arxiv.org/abs/1808.05245)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "###############################################################################\n",
      "#                                                                             #\n",
      "# This is Thermal-FIST version 1.2.2                                          #\n",
      "#                                                                             #\n",
      "# Copyright (c) 2019 Volodymyr Vovchenko <vovchenko@fias.uni-frankfurt.de>    #\n",
      "#                                                                             #\n",
      "# Distributed under the GNU General Public License 3.0 (GPLv3 or later)       #\n",
      "#                                                                             #\n",
      "# Please cite when using this code:                                           #\n",
      "# V. Vovchenko, H. Stoecker, arXiv:1901.05249 [nucl-th]                       #\n",
      "#                                                                             #\n",
      "# The latest version is available at https://github.com/vlvovch/Thermal-FIST  #\n",
      "#                                                                             #\n",
      "###############################################################################\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#include \"inc/loadFISTlibs.h\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using namespace std;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "using namespace thermalfist;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bool useWidth  = true;\n",
    "bool useQStats = false; // Maxwel-Boltzmann statistics used in 1808.05245"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "// Particle list\n",
    "ThermalParticleSystem parts(string(INPUT_FOLDER) + \"/list/PDG2014/list-withnuclei.dat\");\n",
    "//ThermalParticleSystem parts(string(INPUT_FOLDER) + \"/list/PDG2014/list.dat\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "// Fill the list of the considered yield ratios\n",
    "\n",
    "vector<int> pdgs1, pdgs2;\n",
    "vector<string> names1, names2;\n",
    "\n",
    "names1.push_back(\"d\");\n",
    "names2.push_back(\"p\");\n",
    "pdgs1.push_back(1000010020);\n",
    "pdgs2.push_back(2212);\n",
    "\n",
    "names1.push_back(\"He3\");\n",
    "names2.push_back(\"p\");\n",
    "pdgs1.push_back(1000020030);\n",
    "pdgs2.push_back(2212);\n",
    "\n",
    "names1.push_back(\"H3La\");\n",
    "names2.push_back(\"p\");\n",
    "pdgs1.push_back(1010010030);\n",
    "pdgs2.push_back(2212);\n",
    "\n",
    "names1.push_back(\"He4\");\n",
    "names2.push_back(\"p\");\n",
    "pdgs1.push_back(1000020040);\n",
    "pdgs2.push_back(2212);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "// Computes limiting GCE values\n",
    "void ComputeGCERatios(vector<double>& ratiosGCE, double Tch = 0.155) \n",
    "{\n",
    "    ratiosGCE.resize(0);\n",
    "    \n",
    "    ThermalModelIdeal modelgce(&parts);\n",
    "\n",
    "    if (useWidth)\n",
    "        modelgce.SetUseWidth(ThermalParticle::eBW);\n",
    "    else\n",
    "        modelgce.SetUseWidth(ThermalParticle::ZeroWidth);\n",
    "    modelgce.SetStatistics(useQStats);\n",
    "    // Include quantum statistics for mesons only\n",
    "    if (useQStats){\n",
    "        for (int i = 0; i < modelgce.TPS()->Particles().size(); ++i) {\n",
    "            ThermalParticle &part = modelgce.TPS()->Particle(i);\n",
    "            if (part.BaryonCharge() != 0)\n",
    "                part.UseStatistics(false);\n",
    "        }\n",
    "    }\n",
    "\n",
    "    modelgce.SetTemperature(Tch);\n",
    "    modelgce.SetBaryonChemicalPotential(0.);\n",
    "    modelgce.SetElectricChemicalPotential(0.);\n",
    "    modelgce.SetStrangenessChemicalPotential(0.);\n",
    "\n",
    "    modelgce.FillChemicalPotentials();\n",
    "\n",
    "    modelgce.CalculateDensities();\n",
    "\n",
    "    for (int i = 0; i < names1.size(); ++i) {\n",
    "        ratiosGCE.push_back(modelgce.GetDensity(pdgs1[i], 1) / modelgce.GetDensity(pdgs2[i], 1));\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "// Prepare the output (on the screen and to a file)\n",
    "\n",
    "vector<std::ostream*> outs;\n",
    "outs.push_back(&std::cout);\n",
    "\n",
    "outs.push_back(NULL);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "// Routine performing the (logarithmic scale) scan over Vc and dNpi/dy.\n",
    "// Parameter k determines Vc = k dV/dy\n",
    "\n",
    "void PerformScan(ThermalModelBase *model, \n",
    "                   double k,\n",
    "                   double Tch,\n",
    "                   const string& filename)\n",
    "{\n",
    "    ofstream fout(filename);\n",
    "    outs[1] = &fout;\n",
    "    \n",
    "    for(std::ostream* out : outs) {\n",
    "        *out << setw(15) << \"dNpipm/dy\";\n",
    "        *out << setw(15) << \"dV/dy[fm3]\";\n",
    "        *out << setw(15) << \"Vc[fm3]\";\n",
    "\n",
    "        for (int i = 0; i < names1.size(); ++i) {\n",
    "            *out << setw(15) << (names1[i] + \"/\" + names2[i]);\n",
    "        }\n",
    "\n",
    "        *out << endl;\n",
    "    }\n",
    "    \n",
    "    vector<double> ratiosGCE;\n",
    "    ComputeGCERatios(ratiosGCE, Tch);\n",
    "    \n",
    "    model->SetTemperature(Tch);\n",
    "    \n",
    "    double Vmin = 1.;\n",
    "    double Vmax = 15000.;\n",
    "    int iters = 30;\n",
    "    double logVmin = log10(Vmin), logVmax = log10(Vmax);\n",
    "    double dlogV = (logVmax - logVmin) / (iters - 1);\n",
    "    \n",
    "    for (double logV = logVmin; logV <= logVmax + 1.e-3; logV += dlogV) {\n",
    "        double V = pow(10., logV);\n",
    "\n",
    "        model->SetVolume(V);\n",
    "        model->SetCanonicalVolume(k * V);\n",
    "        model->CalculateDensities();\n",
    "\n",
    "        for(std::ostream* out : outs) {\n",
    "            *out << setw(15) << 2. * model->GetDensity(211, 1) * V; // Charged pions dN/dy\n",
    "            *out << setw(15) << V; // dVdy\n",
    "            *out << setw(15) << k * V; // Vc\n",
    "\n",
    "            for (int i = 0; i < names1.size(); ++i) {\n",
    "                *out << setw(15) << model->GetDensity(pdgs1[i], 1) / model->GetDensity(pdgs2[i], 1);\n",
    "            }\n",
    "\n",
    "            *out << endl;\n",
    "        }\n",
    "    }\n",
    "    \n",
    "    for(std::ostream* out : outs) {\n",
    "        *out << setw(15) << \"HML\";\n",
    "        *out << setw(15) << \"HML\";\n",
    "        *out << setw(15) << \"HML\";\n",
    "\n",
    "        for (int i = 0; i < ratiosGCE.size(); ++i) {\n",
    "            *out << setw(15) << ratiosGCE[i];\n",
    "        }\n",
    "\n",
    "        *out << endl;\n",
    "    }\n",
    "    \n",
    "    fout.close();\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "// Prepare the thermal model for scans\n",
    "ThermalModelCanonical modelce(&parts);\n",
    "\n",
    "if (useWidth)\n",
    "    modelce.SetUseWidth(ThermalParticle::eBW);\n",
    "else\n",
    "    modelce.SetUseWidth(ThermalParticle::ZeroWidth);\n",
    "\n",
    "modelce.SetStatistics(useQStats);\n",
    "// Include quantum statistics for mesons only\n",
    "if (useQStats){\n",
    "    for (int i = 0; i < modelce.TPS()->Particles().size(); ++i) {\n",
    "        ThermalParticle &part = modelce.TPS()->Particle(i);\n",
    "        if (part.BaryonCharge() != 0)\n",
    "            part.UseStatistics(false);\n",
    "    }\n",
    "}\n",
    "\n",
    "modelce.SetBaryonChemicalPotential(0.);\n",
    "modelce.SetElectricChemicalPotential(0.);\n",
    "modelce.SetStrangenessChemicalPotential(0.);\n",
    "\n",
    "modelce.SetBaryonCharge(0);\n",
    "modelce.SetElectricCharge(0);\n",
    "modelce.SetStrangeness(0);\n",
    "\n",
    "modelce.FillChemicalPotentials();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      dNpipm/dy     dV/dy[fm3]        Vc[fm3]            d/p          He3/p         H3La/p          He4/p\n",
      "      0.0766118              1              1    4.39316e-05    7.18969e-10    2.74001e-10    1.04936e-14\n",
      "       0.111174        1.39317        1.39317     5.3271e-05    1.07856e-09    4.52602e-10    1.83853e-14\n",
      "       0.163537        1.94091        1.94091    6.65589e-05    1.72585e-09     7.7999e-10    3.51208e-14\n",
      "        0.24466        2.70402        2.70402    8.55916e-05    2.93138e-09    1.39649e-09    7.33834e-14\n",
      "       0.373351        3.76715        3.76715    0.000113073    5.25294e-09    2.58693e-09    1.67317e-13\n",
      "       0.581907        5.24826        5.24826    0.000153105    9.87053e-09    4.93896e-09    4.13383e-13\n",
      "       0.924051        7.31171        7.31171    0.000211909    1.93339e-08    9.67768e-09    1.09583e-12\n",
      "        1.48162        10.1864        10.1864    0.000298724    3.92104e-08    1.93529e-08    3.07929e-12\n",
      "        2.36001        14.1914        14.1914    0.000426539    8.15386e-08    3.91534e-08    9.02968e-12\n",
      "        3.67016         19.771         19.771    0.000611603    1.71281e-07    7.90859e-08    2.70485e-11\n",
      "         5.5338        27.5443        27.5443    0.000869954    3.55908e-07     1.5672e-07    8.04872e-11\n",
      "        8.13402        38.3738        38.3738     0.00120938    7.13181e-07    2.98575e-07    2.30143e-10\n",
      "        11.7495        53.4611        53.4611     0.00161717    1.34026e-06    5.34787e-07    6.09339e-10\n",
      "        16.7693        74.4802        74.4802     0.00205072    2.29784e-06    8.79415e-07    1.43661e-09\n",
      "        23.7381        103.763        103.763      0.0024505    3.52682e-06    1.30482e-06    2.92047e-09\n",
      "        33.4264         144.56         144.56     0.00277631     4.8482e-06    1.74736e-06    5.06374e-09\n",
      "        46.9118        201.396        201.396     0.00302335    6.08032e-06     2.1491e-06    7.59963e-09\n",
      "        65.6923        280.578        280.578     0.00320516    7.12554e-06    2.48297e-06    1.01659e-08\n",
      "        91.8522        390.892        390.892     0.00333743    7.96229e-06    2.74623e-06    1.24926e-08\n",
      "        128.294        544.578        544.578     0.00343314    8.60849e-06    2.94734e-06    1.44523e-08\n",
      "        179.063        758.688        758.688     0.00350218    9.09614e-06    3.09792e-06    1.60233e-08\n",
      "        249.789        1056.98        1056.98     0.00355191    9.45855e-06     3.2092e-06    1.72415e-08\n",
      "        348.322        1472.55        1472.55     0.00358769    9.72511e-06    3.29073e-06    1.81646e-08\n",
      "        485.599         2051.5         2051.5     0.00361337    9.91966e-06    3.35005e-06    1.88529e-08\n",
      "         676.85        2858.09        2858.09     0.00363182     1.0061e-05    3.39308e-06    1.93606e-08\n",
      "        943.273        3981.79        3981.79     0.00364517    1.01636e-05    3.42423e-06    1.97331e-08\n",
      "        1314.43         5547.3         5547.3     0.00365478    1.02376e-05    3.44667e-06    2.00035e-08\n",
      "        1831.53        7728.32        7728.32     0.00366163     1.0291e-05    3.46297e-06    2.01999e-08\n",
      "           2552        10766.8        10766.8     0.00366656    1.03301e-05    3.47493e-06    2.03443e-08\n",
      "        3555.81          15000          15000     0.00366989    1.03569e-05      3.483e-06    2.04435e-08\n",
      "            HML            HML            HML     0.00367894    1.04273e-05    3.50421e-06    2.07054e-08\n"
     ]
    }
   ],
   "source": [
    "// Tch = 155 MeV, Vc = dV/dy\n",
    "PerformScan(&modelce, 1., 0.155, \"out/1808.05245.Vc1.Tch155.dat\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      dNpipm/dy     dV/dy[fm3]        Vc[fm3]            d/p          He3/p         H3La/p          He4/p\n",
      "      0.0928902              1              3     9.3132e-05    3.49993e-09    1.68812e-09    9.41377e-14\n",
      "        0.14269        1.39317         4.1795     0.00012402    6.36975e-09    3.15818e-09    2.20307e-13\n",
      "       0.223818        1.94091        5.82274    0.000169141    1.21337e-08    6.08185e-09     5.5709e-13\n",
      "       0.356972        2.70402        8.11205    0.000235558    2.40468e-08     1.2002e-08    1.50615e-12\n",
      "       0.572334        3.76715        11.3014    0.000333638    4.92144e-08    2.41157e-08    4.29739e-12\n",
      "       0.906276        5.24826        15.7448    0.000477543    1.02861e-07    4.88451e-08    1.27191e-11\n",
      "        1.39585        7.31171        21.9351    0.000684023     2.1588e-07    9.82739e-08    3.81465e-11\n",
      "          2.086        10.1864        30.5593    0.000967632    4.44788e-07    1.92794e-07    1.12526e-10\n",
      "        3.04719        14.1914        42.5742     0.00133117    8.76227e-07    3.61178e-07    3.15399e-10\n",
      "        4.38301         19.771         59.313     0.00175298    1.60448e-06    6.31401e-07    8.08614e-10\n",
      "        6.23721        27.5443        82.6328      0.0021822    2.66048e-06    1.00659e-06    1.82491e-09\n",
      "        8.81228        38.3738        115.121     0.00256137    3.94058e-06    1.44485e-06    3.52905e-09\n",
      "        12.3939        53.4611        160.383     0.00286169    5.25027e-06    1.87952e-06    5.83331e-09\n",
      "        17.3803        74.4802        223.441     0.00308655    6.43004e-06    2.26147e-06    8.41585e-09\n",
      "        24.3252        103.763         311.29     0.00325125    7.40974e-06    2.57276e-06    1.09292e-08\n",
      "        33.9994         144.56        433.679     0.00337083    8.18383e-06    2.81539e-06    1.31486e-08\n",
      "        47.4763        201.396        604.187     0.00345725     8.7767e-06    2.99939e-06    1.49853e-08\n",
      "        66.2512        280.578        841.734     0.00351955    9.22166e-06    3.13652e-06    1.64403e-08\n",
      "        92.4072        390.892        1172.68     0.00356441    9.55113e-06    3.23755e-06    1.75595e-08\n",
      "        128.847        544.578        1633.73     0.00359667    9.79282e-06    3.31139e-06    1.84028e-08\n",
      "        179.615        758.688        2276.06     0.00361981     9.9689e-06    3.36505e-06     1.9029e-08\n",
      "        250.342        1056.98        3170.93     0.00363647    1.00968e-05    3.40394e-06    1.94901e-08\n",
      "        348.869        1472.55        4417.64     0.00364854    1.01895e-05    3.43208e-06    1.98275e-08\n",
      "        486.129         2051.5        6154.51     0.00365717    1.02561e-05    3.45233e-06    2.00717e-08\n",
      "        677.369        2858.09        8574.26     0.00366336    1.03046e-05    3.46714e-06    2.02501e-08\n",
      "        943.822        3981.79        11945.4     0.00366777    1.03398e-05    3.47788e-06    2.03801e-08\n",
      "        1315.05         5547.3        16641.9     0.00367065    1.03631e-05    3.48479e-06    2.04662e-08\n",
      "        1832.21        7728.32          23185     0.00367256    1.03783e-05    3.48915e-06    2.05229e-08\n",
      "        2552.65        10766.8        32300.5     0.00367432    1.03916e-05    3.49308e-06    2.05728e-08\n",
      "        3556.33          15000          45000     0.00367604    1.04045e-05    3.49707e-06     2.0621e-08\n",
      "            HML            HML            HML     0.00367894    1.04273e-05    3.50421e-06    2.07054e-08\n"
     ]
    }
   ],
   "source": [
    "// Tch = 155 MeV, Vc = 3dV/dy\n",
    "PerformScan(&modelce, 3., 0.155, \"out/1808.05245.Vc3.Tch155.dat\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      dNpipm/dy     dV/dy[fm3]        Vc[fm3]            d/p          He3/p         H3La/p          He4/p\n",
      "       0.155843              1              1    0.000121572    6.93548e-09    3.58469e-09    3.20659e-13\n",
      "       0.228961        1.39317        1.39317    0.000154338    1.14019e-08    6.27757e-09    6.35122e-13\n",
      "       0.342221        1.94091        1.94091    0.000201594    1.99102e-08    1.14085e-08    1.38302e-12\n",
      "       0.522145        2.70402        2.70402    0.000270353     3.6667e-08    2.14172e-08    3.29281e-12\n",
      "       0.814642        3.76715        3.76715    0.000371248    7.07102e-08    4.13367e-08    8.48472e-12\n",
      "        1.29622        5.24826        5.24826    0.000520053    1.41691e-07    8.15378e-08    2.33455e-11\n",
      "        2.08113        7.31171        7.31171    0.000738742    2.91922e-07    1.62884e-07    6.74105e-11\n",
      "        3.30782        10.1864        10.1864     0.00105366    6.08188e-07    3.24761e-07    1.99431e-10\n",
      "        5.11027        14.1914        14.1914     0.00148636    1.24955e-06    6.32542e-07    5.84656e-10\n",
      "        7.63879         19.771         19.771      0.0020341    2.45001e-06    1.17113e-06    1.62753e-09\n",
      "        11.1291        27.5443        27.5443     0.00264894    4.42527e-06    2.00202e-06    4.10234e-09\n",
      "        15.9458        38.3738        38.3738     0.00324412    7.16104e-06    3.08779e-06    8.97286e-09\n",
      "        22.6192        53.4611        53.4611     0.00374446    1.02991e-05    4.27187e-06     1.6655e-08\n",
      "        31.8946        74.4802        74.4802      0.0041281    1.33584e-05    5.37905e-06    2.64166e-08\n",
      "        44.8047        103.763        103.763      0.0044108    1.60238e-05    6.31307e-06    3.68112e-08\n",
      "        62.7835         144.56         144.56     0.00461624    1.81911e-05    7.05484e-06    4.65604e-08\n",
      "        87.8263        201.396        201.396     0.00476468    1.98808e-05    7.62346e-06    5.49551e-08\n",
      "        122.712        280.578        280.578     0.00487164    2.11635e-05    8.05001e-06    6.17843e-08\n",
      "        171.311        390.892        390.892      0.0049486    2.21206e-05    8.36555e-06    6.71318e-08\n",
      "        239.016        544.578        544.578     0.00500393    2.28264e-05    8.59684e-06    7.12116e-08\n",
      "        333.339        758.688        758.688      0.0050437    2.33427e-05    8.76533e-06     7.4269e-08\n",
      "        464.748        1056.98        1056.98     0.00507223    2.37183e-05    8.88749e-06    7.65307e-08\n",
      "        647.827        1472.55        1472.55     0.00509267    2.39903e-05    8.97575e-06    7.81891e-08\n",
      "        902.876         2051.5         2051.5     0.00510743    2.41869e-05    9.03939e-06     7.9399e-08\n",
      "        1258.18        2858.09        2858.09     0.00511817    2.43294e-05    9.08549e-06    8.02809e-08\n",
      "        1753.17        3981.79        3981.79     0.00512598    2.44337e-05     9.1195e-06    8.09289e-08\n",
      "        2442.83         5547.3         5547.3     0.00513136    2.45075e-05    9.14368e-06    8.13877e-08\n",
      "        3403.72        7728.32        7728.32     0.00513464    2.45545e-05    9.15871e-06    8.16805e-08\n",
      "        4742.41        10766.8        10766.8     0.00513692    2.45869e-05    9.16853e-06    8.18844e-08\n",
      "        6607.36          15000          15000     0.00513932    2.46185e-05    9.17828e-06    8.20843e-08\n",
      "            HML            HML            HML     0.00514494    2.46926e-05    9.20297e-06    8.25483e-08\n"
     ]
    }
   ],
   "source": [
    "// Tch = 170 MeV, Vc = dV/dy\n",
    "PerformScan(&modelce, 1., 0.170, \"out/1808.05245.Vc1.Tch170.dat\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "C++14",
   "language": "C++14",
   "name": "xcpp14"
  },
  "language_info": {
   "codemirror_mode": "text/x-c++src",
   "file_extension": ".cpp",
   "mimetype": "text/x-c++src",
   "name": "c++",
   "version": "14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}