external_packages/hydro_from_external_file/FluidcellStatistic.cpp

0001 // Copyright @ Chun Shen 2014
0002 #include <iostream>
0003 #include <sstream>
0004 #include <fstream>
0005 #include <cmath>
0006 #include <iomanip>
0007 #include <string>
0008 #include <cstdlib>
0009
0010 #include "./FluidcellStatistic.h"
0011 #include "./SurfaceFinder.h"
0012
0013 using namespace std;
0014
0015 FluidcellStatistic::FluidcellStatistic(void* hydroinfo_ptr_in,
0016                                        ParameterReader* paraRdr_in) {
0017     paraRdr = paraRdr_in;
0018     hydro_type = paraRdr->getVal("hydro_type");
0019
0020     grid_tau0 = 0.0;
0021     grid_tauf = 0.0;
0022     grid_x0 = 0.0;
0023     grid_y0 = 0.0;
0024     if (hydro_type == 0) {
0025 #ifdef USE_HDF5
0026         hydroinfo_ptr = (HydroinfoH5*) hydroinfo_ptr_in;
0027         grid_tau0 = hydroinfo_ptr->getHydrogridTau0();
0028         grid_tauf = hydroinfo_ptr->getHydrogridTaumax();
0029         grid_x0 = - hydroinfo_ptr->getHydrogridXmax();
0030         grid_y0 = - hydroinfo_ptr->getHydrogridYmax();
0031 #endif
0032     } else {
0033         hydroinfo_MUSIC_ptr = (Hydroinfo_MUSIC*) hydroinfo_ptr_in;
0034         grid_tau0 = hydroinfo_MUSIC_ptr->get_hydro_tau0();
0035         grid_tauf = hydroinfo_MUSIC_ptr->get_hydro_tau_max();
0036         grid_x0 = (- hydroinfo_MUSIC_ptr->get_hydro_x_max()
0037                    + hydroinfo_MUSIC_ptr->get_hydro_dx());
0038         grid_y0 = grid_x0;
0039     }
0040     T_dec = paraRdr->getVal("T_cut");
0041     grid_dt = paraRdr->getVal("grid_dt");
0042     grid_dx = paraRdr->getVal("grid_dx");
0043     grid_dy = paraRdr->getVal("grid_dy");
0044     hbarC = 0.19733;
0045 }
0046
0047 FluidcellStatistic::~FluidcellStatistic() {}
0048
0049 void FluidcellStatistic::checkFreezeoutSurface(double Tdec) {
0050     cout << "check the position of the freeze-out surface "
0051          << "at (T = " << Tdec << " GeV) ... " << endl;
0052
0053     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0054     int nx = static_cast<int>(abs(2*grid_x0)/grid_dx) + 1;
0055     int ny = static_cast<int>(abs(2*grid_y0)/grid_dy) + 1;
0056
0057     double tau_local;
0058     hydrofluidCell* fluidCellptr = new hydrofluidCell();
0059     ofstream output;
0060     output.open("results/checkFreezeoutSurface.dat", std::ofstream::app);
0061
0062     for (int itime = 0; itime < ntime; itime++) {
0063         // loop over time evolution
0064         tau_local = grid_tau0 + itime*grid_dt;
0065         for (int i = 0; i < nx; i++) {
0066             // loops over the transverse plane
0067             double x_local = grid_x0 + i*grid_dx;
0068             for (int j = 0; j < ny; j++) {
0069                 double y_local = grid_y0 + j*grid_dy;
0070                 //hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0071                 //                            fluidCellptr);
0072                 double temp_local = fluidCellptr->temperature;
0073                 if (fabs(temp_local - Tdec) < 0.001) {
0074                     output << tau_local << "   " << x_local << "   "
0075                            << y_local << "   "
0076                            << sqrt(x_local*x_local + y_local*y_local) << endl;
0077                 }
0078             }
0079         }
0080     }
0081     output.close();
0082     delete fluidCellptr;
0083     return;
0084 }
0085
0086 void FluidcellStatistic::output_temperature_vs_tau() {
0087     cout << "output temperature vs tau ..." << endl;
0088
0089     int nt = static_cast<int>((grid_tauf - grid_tau0)/grid_dt + 1);
0090     int nx = static_cast<int>(abs(2*grid_x0)/grid_dx) + 1;
0091     int ny = static_cast<int>(abs(2*grid_y0)/grid_dy) + 1;
0092
0093     ofstream output;
0094     output.open("results/tau_vs_T.dat");
0095
0096     hydrofluidCell* fluidCellptr = new hydrofluidCell;
0097     for (int it = 0; it < nt; it++) {
0098         double tau_local = grid_tau0 + it*grid_dt;
0099         double avgTemp, stdTemp;
0100         int numCell = 0;
0101         double tempSum = 0.0;
0102         double tempSumsq = 0.0;
0103         double tempMax = 0.0;
0104         for (int i = 0; i < nx; i++) {
0105             // loops over the transverse plane
0106             double x_local = grid_x0 + i*grid_dx;
0107             for (int j = 0; j < ny; j++) {
0108                 double y_local = grid_y0 + j*grid_dy;
0109                 // get hydro information
0110                 if (hydro_type == 0) {
0111 #ifdef USE_HDF5
0112                     hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0113                                                 fluidCellptr);
0114 #endif
0115                 } else {
0116                     hydroinfo_MUSIC_ptr->getHydroValues(
0117                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0118                 }
0119                 double temp_local = fluidCellptr->temperature;
0120                 if (temp_local > T_dec) {
0121                     numCell += 1;
0122                     tempSum += temp_local;
0123                     tempSumsq += temp_local*temp_local;
0124                     if (tempMax < temp_local)
0125                         tempMax = temp_local;
0126                 }
0127             }
0128         }
0129         if (numCell == 0) {
0130             avgTemp = 0.0;
0131             stdTemp = 0.0;
0132         } else {
0133             avgTemp = tempSum/numCell;
0134             stdTemp = sqrt((tempSumsq - tempSum*tempSum/numCell)/(numCell - 1));
0135         }
0136         output << tau_local << "   " << avgTemp << "   " << stdTemp
0137                << "   " << tempMax << endl;
0138     }
0139     delete fluidCellptr;
0140     return;
0141 }
0142
0143 void FluidcellStatistic::output_flowvelocity_vs_tau() {
0144     cout << "output u^tau vs tau ... " << endl;
0145
0146     int nt = static_cast<int>((grid_tauf - grid_tau0)/grid_dt + 1);
0147     int nx = static_cast<int>(abs(2.*grid_x0)/grid_dx) + 1;
0148     int ny = static_cast<int>(abs(2.*grid_y0)/grid_dy) + 1;
0149
0150     hydrofluidCell* fluidCellptr = new hydrofluidCell;
0151
0152     ofstream output;
0153     output.open("results/tau_vs_v.dat", std::ofstream::app);
0154     output << "# tau (fm)  V4 (fm^4)  <u^tau>  theta" << endl;
0155
0156     for (int it = 1; it < nt; it++) {
0157         double tau_local = grid_tau0 + it*grid_dt;
0158         double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
0159         double V4 = 0.0;
0160         double avg_utau = 0.0;
0161         double avg_theta = 0.0;
0162         int count = 0;
0163         for (int i = 0; i < nx; i++) {
0164             // loops over the transverse plane
0165             double x_local = grid_x0 + i*grid_dx;
0166             for (int j = 0; j < ny; j++) {
0167                 double y_local = grid_y0 + j*grid_dy;
0168                 if (hydro_type == 0) {
0169 #ifdef USE_HDF5
0170                     hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0171                                                 fluidCellptr);
0172 #endif
0173                 } else {
0174                     hydroinfo_MUSIC_ptr->getHydroValues(
0175                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0176                 }
0177                 double temp_local = fluidCellptr->temperature;
0178                 if (temp_local > T_dec) {
0179                     // inside freeze-out surface
0180                     double vx_local = fluidCellptr->vx;
0181                     double vy_local = fluidCellptr->vy;
0182                     double v_perp = sqrt(vx_local*vx_local + vy_local*vy_local);
0183                     double gamma = 1./sqrt(1. - v_perp*v_perp);
0184
0185                     double theta = compute_local_expansion_rate(
0186                                                 tau_local, x_local, y_local);
0187                     avg_utau += gamma*volume_element;
0188                     avg_theta += theta*volume_element;
0189                     V4 += volume_element;
0190                     count++;
0191                 }
0192             }
0193         }
0194         if (count < 1) {
0195             avg_utau = 1.0;
0196             avg_theta = 0.0;
0197             V4 = 0.0;
0198         } else {
0199             avg_utau = avg_utau/V4;
0200             avg_theta = avg_theta/V4;
0201         }
0202         output << tau_local << "   " << V4 << "  " << avg_utau << "  "
0203                << avg_theta << endl;
0204     }
0205     return;
0206 }
0207
0208 void FluidcellStatistic::output_temperature_vs_avg_utau() {
0209     cout << "output utau vs T ... " << endl;
0210
0211     int nt = static_cast<int>((grid_tauf - grid_tau0)/grid_dt + 1);
0212     int nx = static_cast<int>(abs(2*grid_x0)/grid_dx) + 1;
0213     int ny = static_cast<int>(abs(2*grid_y0)/grid_dy) + 1;
0214
0215     ofstream output("results/T_vs_v.dat");
0216     output << "# T (GeV)  V (fm^4)  <u^tau>  theta" << endl;
0217
0218     int n_bin = 41;
0219     double T_max = 0.5;
0220     double dT = (T_max - T_dec)/(n_bin - 1);
0221     double *T_bin_avg = new double[n_bin];
0222     double *V4 = new double[n_bin];
0223     double *avg_utau = new double[n_bin];
0224     double *avg_theta = new double[n_bin];
0225     for (int i = 0; i < n_bin; i++) {
0226         T_bin_avg[i] = 0.0;
0227         V4[i] = 0.0;
0228         avg_utau[i] = 0.0;
0229         avg_theta[i] = 0.0;
0230     }
0231
0232     hydrofluidCell* fluidCellptr = new hydrofluidCell;
0233     for (int it = 1; it < nt; it++) {
0234         double tau_local = grid_tau0 + it*grid_dt;
0235         double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
0236         for (int i = 0; i < nx; i++) {
0237             // loops over the transverse plane
0238             double x_local = grid_x0 + i*grid_dx;
0239             for (int j = 0; j < ny; j++) {
0240                 double y_local = grid_y0 + j*grid_dy;
0241
0242                 // get hydro information
0243                 if (hydro_type == 0) {
0244 #ifdef USE_HDF5
0245                     hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0246                                                 fluidCellptr);
0247 #endif
0248                 } else {
0249                     hydroinfo_MUSIC_ptr->getHydroValues(
0250                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0251                 }
0252                 double temp_local = fluidCellptr->temperature;
0253                 if (temp_local > T_dec) {
0254                     double vx_local = fluidCellptr->vx;
0255                     double vy_local = fluidCellptr->vy;
0256                     double utau_local = 1./sqrt(1. - vx_local*vx_local
0257                                                 - vy_local*vy_local);
0258                     double theta_local = compute_local_expansion_rate(
0259                                                 tau_local, x_local, y_local);
0260                     int T_idx = static_cast<int>((temp_local - T_dec)/dT);
0261                     if (T_idx >= 0 && T_idx < n_bin-1) {
0262                         T_bin_avg[T_idx] += temp_local*volume_element;
0263                         avg_utau[T_idx] += utau_local*volume_element;
0264                         avg_theta[T_idx] += theta_local*volume_element;
0265                         V4[T_idx] += volume_element;
0266                     }
0267                 }
0268             }
0269         }
0270     }
0271     for (int i = 0; i < n_bin; i++) {
0272         double T_bin, utau_bin, theta_bin;
0273         if (fabs(T_bin_avg[i]) < 1e-10) {
0274             T_bin = T_dec + i*dT;
0275             utau_bin = 1.0;
0276             theta_bin = 0.0;
0277         } else {
0278             T_bin = T_bin_avg[i]/(V4[i] + 1e-15);
0279             utau_bin = avg_utau[i]/(V4[i] + 1e-15);
0280             theta_bin = avg_theta[i]/(V4[i] + 1e-15);
0281         }
0282         output << scientific << setw(18) << setprecision(8)
0283                << T_bin << "   " << V4[i] << "  " << utau_bin << "  "
0284                << theta_bin << endl;
0285     }
0286     output.close();
0287
0288     delete[] T_bin_avg;
0289     delete[] V4;
0290     delete[] avg_utau;
0291     delete[] avg_theta;
0292     delete fluidCellptr;
0293
0294     return;
0295 }
0296
0297 void FluidcellStatistic::output_momentum_anisotropy_vs_tau() {
0298     cout << "output momentum anisotropy vs tau ... " << endl;
0299
0300     int nt = static_cast<int>((grid_tauf - grid_tau0)/grid_dt + 1);
0301     int nx = static_cast<int>(abs(2*grid_x0)/grid_dx) + 1;
0302     int ny = static_cast<int>(abs(2*grid_y0)/grid_dy) + 1;
0303
0304     hydrofluidCell* fluidCellptr = new hydrofluidCell;
0305     ofstream output;
0306     output.open("results/tau_vs_epsP.dat");
0307
0308     for (int it = 0; it < nt; it++) {
0309         double tau_local = grid_tau0 + it*grid_dt;
0310         double TxxSum = 0.0;
0311         double TyySum = 0.0;
0312         double epsP = 0.0;
0313         for (int i = 0; i < nx; i++) {
0314             // loops over the transverse plane
0315             double x_local = grid_x0 + i*grid_dx;
0316             for(int j = 0; j < ny; j++) {
0317                 double y_local = grid_y0 + j*grid_dy;
0318                 // get hydro information
0319                 if (hydro_type == 0) {
0320 #ifdef USE_HDF5
0321                     hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0322                                                 fluidCellptr);
0323 #endif
0324                 } else {
0325                     hydroinfo_MUSIC_ptr->getHydroValues(
0326                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0327                 }
0328                 double temp_local = fluidCellptr->temperature;
0329                 if (temp_local > T_dec) {
0330                     double e_local = fluidCellptr->ed;
0331                     double p_local = fluidCellptr->pressure;
0332                     double vx_local = fluidCellptr->vx;
0333                     double vy_local = fluidCellptr->vy;
0334                     double v_perp = sqrt(vx_local*vx_local
0335                                          + vy_local*vy_local);
0336                     double gamma = 1./sqrt(1. - v_perp*v_perp);
0337                     double ux_local = gamma*vx_local;
0338                     double uy_local = gamma*vy_local;
0339                     double pi_xx = fluidCellptr->pi[1][1];
0340                     double pi_yy = fluidCellptr->pi[2][2];
0341                     TxxSum += (
0342                         (e_local+p_local)*ux_local*ux_local + p_local + pi_xx);
0343                     TyySum += (
0344                         (e_local+p_local)*uy_local*uy_local + p_local + pi_yy);
0345                 }
0346             }
0347         }
0348         if (fabs(TxxSum)+ fabs(TyySum) < 1e-10)
0349             epsP = 0.0;
0350         else
0351             epsP = (TxxSum - TyySum)/(TxxSum + TyySum);
0352
0353         output << tau_local-0.6 << "   " << epsP << endl;
0354     }
0355     return;
0356 }
0357
0358 void FluidcellStatistic::outputTempasTauvsX() {
0359     cout << "output 2D contour plot for temperature as function of "
0360          << "tau and x ... " << endl;
0361
0362     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0363     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
0364
0365     double tau_local;
0366     hydrofluidCell* fluidCellptr = new hydrofluidCell();
0367     ofstream output;
0368     output.open("results/TempasTauvsX.dat");
0369
0370     for (int itime = 0; itime < ntime; itime++) {
0371         // loop over time evolution
0372         tau_local = grid_tau0 + itime*grid_dt;
0373         for (int i = 0; i < nx; i++) {
0374             // loops over the transverse plane
0375             double x_local = grid_x0 + i*grid_dx;
0376             double y_local = 0.0;
0377             // get hydro information
0378             if (hydro_type == 0) {
0379 #ifdef USE_HDF5
0380                 hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0381                                             fluidCellptr);
0382 #endif
0383             } else {
0384                 hydroinfo_MUSIC_ptr->getHydroValues(
0385                         x_local, y_local, 0.0, tau_local, fluidCellptr);
0386             }
0387             double temp_local = fluidCellptr->temperature;
0388             if (temp_local > 0.05)
0389                 output << temp_local << "   " ;
0390             else
0391                 output << 0.0 << "   " ;
0392         }
0393         output << endl;
0394     }
0395     output.close();
0396     delete fluidCellptr;
0397     return;
0398 }
0399
0400
0401 void FluidcellStatistic::outputKnudersonNumberasTauvsX() {
0402     cout << "output Knudersen Number as a function of tau and x ..." << endl;
0403
0404     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0405     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx);
0406
0407     double eps = 1e-15;
0408
0409     hydrofluidCell* fluidCellptr = new hydrofluidCell;
0410     ofstream output;
0411     output.open("results/KnudsenNumberasTauvsX.dat");
0412
0413     for (int itime = 1; itime < ntime; itime++) {
0414         // loop over time evolution
0415         double tau_local = grid_tau0 + itime*grid_dt;
0416         for (int i = 0; i < nx; i++) {
0417             // loops over the transverse plane
0418             double x_local = grid_x0 + i*grid_dx;
0419             double y_local = 0.0;
0420             double theta = compute_local_expansion_rate(tau_local,
0421                                                         x_local, y_local);
0422
0423             if (hydro_type == 0) {
0424 #ifdef USE_HDF5
0425                 hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0426                                             fluidCellptr);
0427 #endif
0428             } else {
0429                 hydroinfo_MUSIC_ptr->getHydroValues(
0430                         x_local, y_local, 0.0, tau_local, fluidCellptr);
0431             }
0432
0433             double eta_s = 0.08;
0434             double L_micro = (5.*eta_s
0435                               /(fabs(fluidCellptr->temperature) + eps))*hbarC;
0436             double L_macro = 1/(fabs(theta) + eps);
0437             double Knudsen = L_micro/L_macro;
0438
0439             output << Knudsen << "    ";
0440         }
0441         output << endl;
0442     }
0443     output.close();
0444     delete fluidCellptr;
0445     return;
0446 }
0447
0448 double FluidcellStatistic::compute_local_expansion_rate(
0449                         double tau_local, double x_local, double y_local) {
0450     // this function computes the local expansion rate at the given
0451     // space-time point (tau_loca, x_local, y_local)
0452     // theta = \patial_mu u^\mu + u^tau/tau
0453
0454     hydrofluidCell* fluidCellptr = new hydrofluidCell;
0455     hydrofluidCell* fluidCellptrt1 = new hydrofluidCell;
0456     hydrofluidCell* fluidCellptrt2 = new hydrofluidCell;
0457     hydrofluidCell* fluidCellptrx1 = new hydrofluidCell;
0458     hydrofluidCell* fluidCellptrx2 = new hydrofluidCell;
0459     hydrofluidCell* fluidCellptry1 = new hydrofluidCell;
0460     hydrofluidCell* fluidCellptry2 = new hydrofluidCell;
0461
0462     if (hydro_type == 0) {
0463 #ifdef USE_HDF5
0464         hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0465                                     fluidCellptr);
0466         hydroinfo_ptr->getHydroinfo(tau_local-grid_dt, x_local, y_local,
0467                                     fluidCellptrt1);
0468         hydroinfo_ptr->getHydroinfo(tau_local+grid_dt, x_local, y_local,
0469                                     fluidCellptrt2);
0470         hydroinfo_ptr->getHydroinfo(tau_local, x_local-grid_dx, y_local,
0471                                     fluidCellptrx1);
0472         hydroinfo_ptr->getHydroinfo(tau_local, x_local+grid_dx, y_local,
0473                                     fluidCellptrx2);
0474         hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local-grid_dy,
0475                                     fluidCellptry1);
0476         hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local+grid_dy,
0477                                     fluidCellptry2);
0478 #endif
0479     } else {
0480         hydroinfo_MUSIC_ptr->getHydroValues(
0481                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0482         hydroinfo_MUSIC_ptr->getHydroValues(
0483                 x_local, y_local, 0.0, tau_local - grid_dt, fluidCellptrt1);
0484         hydroinfo_MUSIC_ptr->getHydroValues(
0485                 x_local, y_local, 0.0, tau_local + grid_dt, fluidCellptrt2);
0486         hydroinfo_MUSIC_ptr->getHydroValues(
0487                 x_local - grid_dx, y_local, 0.0, tau_local, fluidCellptrx1);
0488         hydroinfo_MUSIC_ptr->getHydroValues(
0489                 x_local + grid_dx, y_local, 0.0, tau_local, fluidCellptrx2);
0490         hydroinfo_MUSIC_ptr->getHydroValues(
0491                 x_local, y_local - grid_dy, 0.0, tau_local, fluidCellptry1);
0492         hydroinfo_MUSIC_ptr->getHydroValues(
0493                 x_local, y_local + grid_dy, 0.0, tau_local, fluidCellptry2);
0494     }
0495
0496     double u0 = 1./sqrt(1. - fluidCellptr->vx*fluidCellptr->vx
0497                         + fluidCellptr->vy*fluidCellptr->vy);
0498     double u0t1 = 1./sqrt(1. - fluidCellptrt1->vx*fluidCellptrt1->vx
0499                           + fluidCellptrt1->vy*fluidCellptrt1->vy);
0500     double u0t2 = 1./sqrt(1. - fluidCellptrt2->vx*fluidCellptrt2->vx
0501                           + fluidCellptrt2->vy*fluidCellptrt2->vy);
0502     double u1x1 = (fluidCellptrx1->vx
0503                    /sqrt(1. - fluidCellptrx1->vx*fluidCellptrx1->vx
0504                          + fluidCellptrx1->vy*fluidCellptrx1->vy));
0505     double u1x2 = (fluidCellptrx2->vx
0506                    /sqrt(1. - fluidCellptrx2->vx*fluidCellptrx2->vx
0507                          + fluidCellptrx2->vy*fluidCellptrx2->vy));
0508     double u2y1 = (fluidCellptry1->vy
0509                    /sqrt(1. - fluidCellptry1->vx*fluidCellptry1->vx
0510                          + fluidCellptry1->vy*fluidCellptry1->vy));
0511     double u2y2 = (fluidCellptry2->vy
0512                    /sqrt(1. - fluidCellptry2->vx*fluidCellptry2->vx
0513                          + fluidCellptry2->vy*fluidCellptry2->vy));
0514
0515     double d0u0 = (u0t2 - u0t1)/2./grid_dt;
0516     double d1u1 = (u1x2 - u1x1)/2./grid_dx;
0517     double d2u2 = (u2y2 - u2y1)/2./grid_dy;
0518     double theta = (d0u0 + d1u1 + d2u2 + u0/tau_local);
0519
0520     delete fluidCellptr;
0521     delete fluidCellptrt1;
0522     delete fluidCellptrt2;
0523     delete fluidCellptrx1;
0524     delete fluidCellptrx2;
0525     delete fluidCellptry1;
0526     delete fluidCellptry2;
0527     return(theta);
0528 }
0529
0530 void FluidcellStatistic::outputinverseReynoldsNumberasTauvsX() {
0531     cout << "output inverse Reynolds number as a function of "
0532          << "tau and x ... " << endl;
0533
0534     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0535     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
0536
0537     double MAX = 1000.;
0538
0539     hydrofluidCell* fluidCellptr = new hydrofluidCell();
0540     ofstream output;
0541     output.open("results/inverseReynoldsNumberasTauvsX.dat");
0542
0543     for (int itime = 0; itime < ntime; itime++) {
0544         // loop over time evolution
0545         double tau_local = grid_tau0 + itime*grid_dt;
0546         for(int i = 0; i < nx; i++) {
0547             // loops over the transverse plane
0548             double x_local = grid_x0 + i*grid_dx;
0549             double y_local = 0.0;
0550
0551             if (hydro_type == 0) {
0552 #ifdef USE_HDF5
0553                 hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0554                                             fluidCellptr);
0555 #endif
0556             } else {
0557                 hydroinfo_MUSIC_ptr->getHydroValues(
0558                         x_local, y_local, 0.0, tau_local, fluidCellptr);
0559             }
0560             double pi2 = (
0561                     fluidCellptr->pi[0][0]*fluidCellptr->pi[0][0]
0562                     + fluidCellptr->pi[1][1]*fluidCellptr->pi[1][1]
0563                     + fluidCellptr->pi[2][2]*fluidCellptr->pi[2][2]
0564                     + fluidCellptr->pi[3][3]*fluidCellptr->pi[3][3]
0565                     - 2.*(fluidCellptr->pi[0][1]*fluidCellptr->pi[0][1]
0566                           + fluidCellptr->pi[0][2]*fluidCellptr->pi[0][2]
0567                           + fluidCellptr->pi[0][3]*fluidCellptr->pi[0][3])
0568                     + 2.*(fluidCellptr->pi[1][2]*fluidCellptr->pi[1][2]
0569                           + fluidCellptr->pi[1][3]*fluidCellptr->pi[1][3]
0570                           + fluidCellptr->pi[2][3]*fluidCellptr->pi[2][3]));
0571
0572             double inverseReynold;
0573
0574             if (pi2 >= 0)
0575                 inverseReynold = sqrt(pi2)/(fluidCellptr->ed
0576                                             + fluidCellptr->pressure + 1e-15);
0577             else
0578                 inverseReynold = MAX;
0579
0580             output << inverseReynold << "    " ;
0581         }
0582         output << endl;
0583     }
0584     output.close();
0585     delete fluidCellptr;
0586     return;
0587 }
0588
0589 void FluidcellStatistic::outputBulkinverseReynoldsNumberasTauvsX() {
0590     cout << "output bulk inverse Reynolds number as a function of "
0591          << "tau and x ... " << endl;
0592
0593     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0594     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
0595
0596     hydrofluidCell* fluidCellptr = new hydrofluidCell();
0597     ofstream output;
0598     output.open("results/inverseReynoldsNumberasTauvsX.dat");
0599
0600     for (int itime = 0 ; itime < ntime; itime++) {
0601         // loop over time evolution
0602         double tau_local = grid_tau0 + itime*grid_dt;
0603         for (int i = 0; i < nx; i++) {
0604             // loops over the transverse plane
0605             double x_local = grid_x0 + i*grid_dx;
0606             double y_local = 0.0;
0607
0608             if (hydro_type == 0) {
0609 #ifdef USE_HDF5
0610                 hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0611                                             fluidCellptr);
0612 #endif
0613             } else {
0614                 hydroinfo_MUSIC_ptr->getHydroValues(
0615                         x_local, y_local, 0.0, tau_local, fluidCellptr);
0616             }
0617
0618             double inverseReynold;
0619             inverseReynold = fabs(fluidCellptr->bulkPi)/fluidCellptr->pressure;
0620             output << inverseReynold << "    " ;
0621         }
0622         output << endl;
0623     }
0624     output.close();
0625     delete fluidCellptr;
0626     return;
0627 }
0628
0629 void FluidcellStatistic::analysis_hydro_volume_for_photon(double T_cut) {
0630     double V_3 = calculate_hypersurface_3volume(T_cut);
0631     double V_4 = calculate_spacetime_4volume(T_cut);
0632     double average_tau = calculate_average_tau(T_cut);
0633     double average_T4 = calculate_average_temperature4(T_cut);
0634     double average_GammaT =
0635         calculate_average_integrated_photonRate_parameterization(T_cut);
0636     stringstream output;
0637     output << "results/volume_info_for_photon_Tcut_" << T_cut << ".dat";
0638     ofstream of(output.str().c_str());
0639     of << "# V_4  <tau>  V_3  <T_4>  <GammaT> " << endl;
0640     of << scientific << setw(18) << setprecision(8)
0641        << V_4 << "  " << average_tau << "  "
0642        << V_3 << "  " << average_T4 << "  " << average_GammaT << endl;
0643     of.close();
0644 }
0645
0646 double FluidcellStatistic::calculate_spacetime_4volume(double T_cut) {
0647     // this function calculates the space-time 4 volume of the medium
0648     // inside a give temperature T_cut [GeV]
0649     // the output volume V_4 is in [fm^4]
0650     // deta = 1
0651
0652     cout << "compute 4-volume inside T = " << T_cut << " GeV ..." << endl;
0653
0654     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0655     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
0656     int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy) + 1;
0657
0658     hydrofluidCell* fluidCellptr = new hydrofluidCell();
0659
0660     double volume = 0.0;
0661     for (int itime = 0; itime < ntime; itime++) {
0662         // loop over time evolution
0663         double tau_local = grid_tau0 + itime*grid_dt;
0664         double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
0665         for (int i = 0; i < nx; i++) {
0666             // loops over the transverse plane
0667             double x_local = grid_x0 + i*grid_dx;
0668             for (int j = 0; j < ny; j++) {
0669                 double y_local = grid_y0 + j*grid_dy;
0670                 // get hydro information
0671                 if (hydro_type == 0) {
0672 #ifdef USE_HDF5
0673                     hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0674                                                 fluidCellptr);
0675 #endif
0676                 } else {
0677                     hydroinfo_MUSIC_ptr->getHydroValues(
0678                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0679                 }
0680                 double T_local = fluidCellptr->temperature;  // GeV
0681                 if (T_local > T_cut) {
0682                     volume += volume_element;
0683                 }
0684             }
0685         }
0686     }
0687     delete fluidCellptr;
0688     return(volume);
0689 }
0690
0691 double FluidcellStatistic::calculate_average_tau(double T_cut) {
0692     // this function calculates the average <\tau> of the medium
0693     // inside a give temperature T_cut [GeV]
0694     // the output <tau> is in [fm]
0695
0696     cout << "compute <tau> ... " << endl;
0697
0698     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0699     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
0700     int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy) + 1;
0701
0702     hydrofluidCell* fluidCellptr = new hydrofluidCell();
0703
0704     double average_tau = 0.0;
0705     double volume = 0.0;
0706     for (int itime = 0; itime < ntime; itime++) {
0707         // loop over time evolution
0708         double tau_local = grid_tau0 + itime*grid_dt;
0709         double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
0710         for (int i = 0; i < nx; i++) {
0711             // loops over the transverse plane
0712             double x_local = grid_x0 + i*grid_dx;
0713             for (int j = 0; j < ny; j++) {
0714                 double y_local = grid_y0 + j*grid_dy;
0715                 // get hydro information
0716                 if (hydro_type == 0) {
0717 #ifdef USE_HDF5
0718                     hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0719                                                 fluidCellptr);
0720 #endif
0721                 } else {
0722                     hydroinfo_MUSIC_ptr->getHydroValues(
0723                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0724                 }
0725                 double T_local = fluidCellptr->temperature;  // GeV
0726                 if (T_local > T_cut) {
0727                     volume += volume_element;
0728                     average_tau += tau_local*volume_element;
0729                 }
0730             }
0731         }
0732     }
0733     average_tau /= volume;
0734     delete fluidCellptr;
0735     return(average_tau);
0736 }
0737
0738 double FluidcellStatistic::calculate_average_temperature4(double T_cut) {
0739     // this function calculates the average <T^4> of the medium
0740     // inside a give temperature T_cut [GeV]
0741     // the output <T^4> is in [GeV^4]
0742
0743     cout << "compute <T^4> ..." << endl;
0744     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0745     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
0746     int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy) + 1;
0747
0748     hydrofluidCell* fluidCellptr = new hydrofluidCell();
0749
0750     double average_T4 = 0.0;
0751     double volume = 0.0;
0752     for (int itime = 0; itime < ntime; itime++) {
0753         // loop over time evolution
0754         double tau_local = grid_tau0 + itime*grid_dt;
0755         double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
0756         for (int i = 0; i < nx; i++) {
0757             // loops over the transverse plane
0758             double x_local = grid_x0 + i*grid_dx;
0759             for (int j = 0; j < ny; j++) {
0760                 double y_local = grid_y0 + j*grid_dy;
0761                 // get hydro information
0762                 if (hydro_type == 0) {
0763 #ifdef USE_HDF5
0764                     hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0765                                                 fluidCellptr);
0766 #endif
0767                 } else {
0768                     hydroinfo_MUSIC_ptr->getHydroValues(
0769                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0770                 }
0771                 double T_local = fluidCellptr->temperature;  // GeV
0772                 if (T_local > T_cut) {
0773                     volume += volume_element;
0774                     average_T4 += (
0775                             T_local*T_local*T_local*T_local*volume_element);
0776                 }
0777             }
0778         }
0779     }
0780     average_T4 /= volume;
0781     delete fluidCellptr;
0782     return(average_T4);
0783 }
0784
0785 double FluidcellStatistic::
0786     calculate_average_integrated_photonRate_parameterization(double T_cut) {
0787     // this function calculates the average Gamma(T) of the medium
0788     // Gamma(T) = 1.746 T^4.095       for T > Tsw = 0.18
0789     //          = 11663.923 T^9.024   for T < Tsw
0790     // inside a give temperature T_cut [GeV]
0791     // the output <Gamma(T)> is in [1/fm^4]
0792
0793     cout << "compute <Gamma(T)> ..." << endl;
0794     int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
0795     int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
0796     int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy) + 1;
0797
0798     hydrofluidCell* fluidCellptr = new hydrofluidCell();
0799
0800     double average_GammaT = 0.0;
0801     double volume = 0.0;
0802     double Tsw = 0.18;   // switching temperature for QGP rates to HG rates
0803     for (int itime = 0; itime < ntime; itime++) {
0804         // loop over time evolution
0805         double tau_local = grid_tau0 + itime*grid_dt;
0806         double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
0807         for (int i = 0; i < nx; i++) {
0808             // loops over the transverse plane
0809             double x_local = grid_x0 + i*grid_dx;
0810             for (int j = 0; j < ny; j++) {
0811                 double y_local = grid_y0 + j*grid_dy;
0812                 // get hydro information
0813                 if (hydro_type == 0) {
0814 #ifdef USE_HDF5
0815                     hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
0816                                                 fluidCellptr);
0817 #endif
0818                 } else {
0819                     hydroinfo_MUSIC_ptr->getHydroValues(
0820                             x_local, y_local, 0.0, tau_local, fluidCellptr);
0821                 }
0822                 double T_local = fluidCellptr->temperature;  // GeV
0823                 if (T_local > T_cut) {
0824                     volume += volume_element;
0825                     double GammaT = 0.0;
0826                     if (T_local > Tsw) {
0827                         GammaT = 1.746*pow(T_local, 4.095);
0828                     } else {
0829                         GammaT = 11663.923*pow(T_local, 9.024);
0830                     }
0831                     average_GammaT += GammaT*volume_element;
0832                 }
0833             }
0834         }
0835     }
0836     average_GammaT /= volume;
0837     delete fluidCellptr;
0838     return(average_GammaT);
0839 }
0840
0841 double FluidcellStatistic::calculate_hypersurface_3volume(double T_cut) {
0842     // this function calculates the surface area of an isothermal hyper-surface
0843     // at a give temperature T_cut [GeV]
0844     // the output volume V_3 = int u^\mu d^3\sigma_\mu is in [fm^3]
0845     // deta = 1
0846
0847     cout << "compute 3-volume at T = " << T_cut << " GeV ... " << endl;
0848     // first find the hyper-surface
0849     void* hydro_ptr = NULL;
0850     if (hydro_type == 1) {
0851 #ifdef USE_HDF5
0852         hydro_ptr = hydroinfo_ptr;
0853 #endif
0854     } else {
0855         hydro_ptr = hydroinfo_MUSIC_ptr;
0856     }
0857     SurfaceFinder surfFinder(hydro_ptr, paraRdr, T_cut);
0858     surfFinder.Find_full_hypersurface();
0859
0860     // load the hyper-surface file
0861     ifstream surf("hyper_surface_2+1d.dat");
0862     if (!surf.is_open()) {
0863         cout << "FluidcellStatistic::calculate_hypersurface_3volume:"
0864              << "can not open the file hyper_surface_2+1d.dat!" << endl;
0865         exit(1);
0866     }
0867
0868     double volume = 0.0;
0869     double tau, x, y, da0, da1, da2, vx, vy, T;
0870     surf >> tau >> x >> y >> da0 >> da1 >> da2 >> T >> vx >> vy;
0871     while (!surf.eof()) {
0872         double u0 = 1./sqrt(1. - vx*vx - vy*vy);
0873         double ux = u0*vx;
0874         double uy = u0*vy;
0875         volume += tau*(u0*da0 + ux*da1 + uy*da2);
0876         surf >> tau >> x >> y >> da0 >> da1 >> da2 >> T >> vx >> vy;
0877     }
0878     return(volume);
0879 }
0880