vhlle/src/eos.cpp

0001 #include <string>
0002 #include <iostream>
0003 #include <fstream>
0004 #include <cstdlib>
0005 #include "eos.h"
0006 #include "inc.h"
0007 #include <TGraph.h>
0008
0009 using namespace std;
0010
0011 const double bagp = pow(247.19 / 197.32, 4) / gevtofm;
0012 const double bagt = pow(247.19 / 197.32, 4) / gevtofm;
0013 const double deg = 16.0 + 3.0 * 12.0 * (7.0 / 8.0);
0014
0015 // EoS choise --> MOVED to Makefile
0016 //#define TABLE // Laine, etc
0017 //#define SIMPLE  // p=e/3
0018
0019 double EoS::s(double e, double nb, double nq, double ns) {
0020   double T, mub, muq, mus, p;
0021   eos(e, nb, nq, ns, T, mub, muq, mus, p);
0022   if (T > 0.0)
0023     return (e + p - mub * nb - muq * nq - mus * ns) / T;
0024   else
0025     return 0.;
0026 }
0027
0028 EoSs::EoSs(string fname, int ncols) {
0029 #if defined TABLE || defined LAINE_CFO
0030
0031   int edat = 10000;
0032   double* e = new double[edat];
0033   double* pGrid = new double[edat];
0034   double* tpGrid = new double[edat];
0035   double* muGrid = new double[edat];
0036
0037   ifstream finput(fname.c_str(), ios::in);
0038   // ofstream fout("debug.txt");
0039   if (!finput) {
0040     cerr << "can't open input file \"" << fname.c_str() << "\"" << endl;
0041     exit(1);
0042   }
0043   edat = 0;
0044   while (!finput.eof()) {
0045     if (ncols == 3) {
0046       finput >> e[edat] >> pGrid[edat] >> tpGrid[edat];
0047       muGrid[edat] = 0.;
0048     } else {
0049       finput >> e[edat] >> pGrid[edat] >> tpGrid[edat] >> muGrid[edat];
0050     }
0051     if (pGrid[edat] < 0.) pGrid[edat] = 0.;
0052     edat++;
0053   }
0054   finput.close();
0055
0056   gp = new TGraph(edat, e, pGrid);
0057   gT = new TGraph(edat, e, tpGrid);
0058   gmu = new TGraph(edat, e, muGrid);
0059
0060 #elif defined SIMPLE
0061 // nothing
0062 #endif
0063 }
0064
0065 EoSs::~EoSs(void) {}
0066
0067 double EoSs::p(double e) {
0068 #if defined TABLE
0069   return gp->Eval(e);
0070 #elif defined SIMPLE
0071   return e / 3.;
0072 #endif
0073 }
0074
0075 double EoSs::dpe(double e) {
0076 #if defined TABLE
0077   return (gp->Eval(e * 1.1) - gp->Eval(e)) / (0.1 * e);
0078 #elif defined SIMPLE
0079   return 1. / 3.;
0080 #endif
0081 }
0082
0083 double EoSs::t(double e) {
0084 #if defined TABLE
0085   return gT->Eval(e);
0086 #elif defined SIMPLE
0087   const double cnst =
0088       (16 + 0.5 * 21.0 * 2.5) * pow(C_PI, 2) / 30.0 / pow(0.197326968, 3);
0089   return e > 0. ? 1.0 * pow(e / cnst, 0.25) : 0.;
0090 #endif
0091 }
0092
0093 double EoSs::mu(double e) {
0094 #if defined TABLE
0095   return gmu->Eval(e);
0096 #elif defined SIMPLE
0097   return 0.;
0098 #endif
0099 }