vhlle/src/icGlauber.cpp

0001 #include <fstream>
0002 #include <iomanip>
0003 #include <iomanip>
0004 #include <TF1.h>
0005 #include <TF2.h>
0006 #include <TGraph.h>
0007
0008 #include "fld.h"
0009 #include "eos.h"
0010 #include "icGlauber.h"
0011 #include "inc.h"
0012
0013 using namespace std;
0014
0015 // --------------------------------------------
0016 //   Initial state from optical Glauber model
0017 // --------------------------------------------
0018
0019 // Au nucleus parameters for optical Glauber
0020 const double A = 197.0;    // mass number
0021 const double Ra = 6.37;    // radius
0022 const double dlt = 0.54;   // diffuseness
0023 const double sigma = 4.0;  // NN cross section in fm^2
0024
0025 const int nphi = 301;
0026
0027 ICGlauber::ICGlauber(double e, double impactPar, double _tau0) {
0028   epsilon = e;
0029   b = impactPar;
0030   tau0 = _tau0;
0031 }
0032
0033 ICGlauber::~ICGlauber(void) {}
0034
0035 double ICGlauber::eProfile(double x, double y) {
0036  prms[0] = sqrt((x + b / 2.0) * (x + b / 2.0) + y * y);
0037  iff->SetParameters(prms);
0038  const double tpp = iff->Integral(-3.0 * Ra, 3.0 * Ra, 1.0e-9);
0039  prms[0] = sqrt((x - b / 2.0) * (x - b / 2.0) + y * y);
0040  iff->SetParameters(prms);
0041  const double tmm = iff->Integral(-3.0 * Ra, 3.0 * Ra, 1.0e-9);
0042  return epsilon *
0043         pow(1. / rho0 * (tpp * (1.0 - pow((1.0 - sigma * tmm / A), A)) +
0044                          tmm * (1.0 - pow((1.0 - sigma * tpp / A), A))),
0045             1.0);
0046 }
0047
0048 void ICGlauber::findRPhi(void) {
0049   _rphi = new double[nphi];
0050   for (int iphi = 0; iphi < nphi; iphi++) {
0051     double phi = iphi * C_PI * 2. / (nphi - 1);
0052     double r = 0., r1 = 0., r2 = 2. * Ra;
0053     while (fabs((r2 - r1) / r2) > 0.001 && r2 > 0.001) {
0054       r = 0.5 * (r1 + r2);
0055       if (eProfile(r * cos(phi), r * sin(phi)) > 0.5)
0056         r1 = r;
0057       else
0058         r2 = r;
0059     }
0060     _rphi[iphi] = r;
0061   }
0062 }
0063
0064
0065 double ICGlauber::rPhi(double phi) {
0066   const double cpi = C_PI;
0067   phi = phi - 2. * cpi * floor(phi / 2. / cpi);
0068   int iphi = (int)(phi / (2. * cpi) * (nphi - 1));
0069   int iphi1 = iphi + 1;
0070   if (iphi1 == nphi) iphi = nphi - 2;
0071   return _rphi[iphi] * (1. - (phi / (2. * cpi) * (nphi - 1) - iphi)) +
0072          _rphi[iphi1] * (phi / (2. * cpi) * (nphi - 1) - iphi);
0073 }
0074
0075
0076 void ICGlauber::setIC(Fluid *f, EoS *eos) {
0077   double e, nb, nq, vx = 0., vy = 0., vz = 0.;
0078   Cell *c;
0079   ofstream fvel("velocity_debug.txt");
0080
0081   TF2 *ff = 0;
0082   double prms2[2], intgr2;
0083   cout << "finding normalization constant\n";
0084   ff = new TF2("ThicknessF", this, &ICGlauber::Thickness, -3.0 * Ra, 3.0 * Ra,
0085                -3.0 * Ra, 3.0 * Ra, 2, "IC", "Thickness");
0086   prms2[0] = Ra;
0087   prms2[1] = dlt;
0088   ff->SetParameters(prms2);
0089   intgr2 = ff->Integral(-3.0 * Ra, 3.0 * Ra, -3.0 * Ra, 3.0 * Ra, 1.0e-9);
0090   if (intgr2 == 0.0) {
0091     cerr << "IC::setICGlauber Error! ff->Integral == 0; Return -1\n";
0092     delete ff;
0093     exit(1);
0094   }
0095   delete ff;
0096   cout << "a = " << A / intgr2 << endl;
0097   prms[1] = A / intgr2;
0098   prms[2] = Ra;
0099   prms[3] = dlt;
0100   iff = new TF1("WoodSaxonDF", this, &ICGlauber::WoodSaxon, -3.0 * Ra, 3.0 * Ra, 4,
0101                 "IC", "WoodSaxon");
0102   prms[0] = 0.0;
0103   iff->SetParameters(prms);
0104   const double tpp = iff->Integral(-3.0 * Ra, 3.0 * Ra, 1.0e-9);
0105   rho0 = 2.0 * tpp * (1.0 - pow((1.0 - sigma * tpp / A), A));
0106
0107   findRPhi();  // fill in R(phi) table
0108   cout << "R(phi) =  ";
0109   for (int jj = 0; jj < 5; jj++) cout << rPhi(jj * C_PI / 2.) << "  ";  // test
0110   cout << endl;
0111   //--------------
0112   double avv_num = 0., avv_den = 0.;
0113   double Etotal = 0.0;
0114
0115   for (int ix = 0; ix < f->getNX(); ix++)
0116     for (int iy = 0; iy < f->getNY(); iy++)
0117       for (int iz = 0; iz < f->getNZ(); iz++) {
0118         c = f->getCell(ix, iy, iz);
0119         double x = f->getX(ix);
0120         double y = f->getY(iy);
0121         double eta = f->getZ(iz);
0122         double etaFactor;
0123         double eta1 = fabs(eta) < 1.3 ? 0.0 : fabs(eta) - 1.3;
0124         etaFactor = exp(-eta1 * eta1 / 2.1 / 2.1) * (fabs(eta) < 5.3 ? 1.0 : 0.0);
0125         e = eProfile(x, y) * etaFactor;
0126         if (e < 0.5) e = 0.0;
0127         vx = vy = 0.0;
0128         nb = nq = 0.0;
0129         vz = 0.0;
0130
0131       avv_num += sqrt(vx * vx + vy * vy) * e;
0132       avv_den += e;
0133
0134         c->setPrimVar(eos, tau0, e, nb, nq, 0., vx, vy, vz);
0135         double _p = eos->p(e, nb, nq, 0.);
0136         const double gamma2 = 1.0 / (1.0 - vx * vx - vy * vy - vz * vz);
0137         Etotal +=
0138             ((e + _p) * gamma2 * (cosh(eta) + vz * sinh(eta)) - _p * cosh(eta));
0139         c->saveQprev();
0140
0141         if (e > 0.) c->setAllM(1.);
0142       }
0143   fvel.close();
0144   cout << "average initial flow = " << avv_num / avv_den << endl;
0145   cout << "total energy = " << Etotal *f->getDx() * f->getDy() * f->getDz() *
0146                                    tau0 << endl;
0147 }
0148
0149 double ICGlauber::Thickness(double *x, double *p) {
0150   // p[0]: Ra radius; p[1]: delta = 0.54fm
0151   double intgrl, prms[4];
0152   TF1 *iff = 0;
0153   iff = new TF1("WoodSaxonDF", this, &ICGlauber::WoodSaxon, -3.0 * p[0], 3.0 * p[0], 4,
0154                 "IC", "WoodSaxon");
0155   prms[0] = sqrt(x[0] * x[0] + x[1] * x[1]);  //
0156   prms[1] = 1.0;  // normalization parameter which must be found.
0157   prms[2] = p[0];
0158   prms[3] = p[1];
0159   iff->SetParameters(prms);
0160   intgrl = iff->Integral(-3.0 * p[0], 3.0 * p[0], 1.0e-9);
0161   if (iff) delete iff;
0162   return intgrl;
0163 }
0164
0165 double ICGlauber::WoodSaxon(double *x, double *p) {
0166   return p[1] / (exp((sqrt(x[0] * x[0] + p[0] * p[0]) - p[2]) / p[3]) + 1.0);
0167 }