Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FluidcellStatistic.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FluidcellStatistic.cpp
1 // Copyright @ Chun Shen 2014
2 #include <iostream>
3 #include <sstream>
4 #include <fstream>
5 #include <cmath>
6 #include <iomanip>
7 #include <string>
8 #include <cstdlib>
9 
10 #include "./FluidcellStatistic.h"
11 #include "./SurfaceFinder.h"
12 
13 using namespace std;
14 
15 FluidcellStatistic::FluidcellStatistic(void* hydroinfo_ptr_in,
16  ParameterReader* paraRdr_in) {
17  paraRdr = paraRdr_in;
18  hydro_type = paraRdr->getVal("hydro_type");
19 
20  grid_tau0 = 0.0;
21  grid_tauf = 0.0;
22  grid_x0 = 0.0;
23  grid_y0 = 0.0;
24  if (hydro_type == 0) {
25 #ifdef USE_HDF5
26  hydroinfo_ptr = (HydroinfoH5*) hydroinfo_ptr_in;
27  grid_tau0 = hydroinfo_ptr->getHydrogridTau0();
28  grid_tauf = hydroinfo_ptr->getHydrogridTaumax();
29  grid_x0 = - hydroinfo_ptr->getHydrogridXmax();
30  grid_y0 = - hydroinfo_ptr->getHydrogridYmax();
31 #endif
32  } else {
33  hydroinfo_MUSIC_ptr = (Hydroinfo_MUSIC*) hydroinfo_ptr_in;
34  grid_tau0 = hydroinfo_MUSIC_ptr->get_hydro_tau0();
35  grid_tauf = hydroinfo_MUSIC_ptr->get_hydro_tau_max();
36  grid_x0 = (- hydroinfo_MUSIC_ptr->get_hydro_x_max()
37  + hydroinfo_MUSIC_ptr->get_hydro_dx());
38  grid_y0 = grid_x0;
39  }
40  T_dec = paraRdr->getVal("T_cut");
41  grid_dt = paraRdr->getVal("grid_dt");
42  grid_dx = paraRdr->getVal("grid_dx");
43  grid_dy = paraRdr->getVal("grid_dy");
44  hbarC = 0.19733;
45 }
46 
47 FluidcellStatistic::~FluidcellStatistic() {}
48 
49 void FluidcellStatistic::checkFreezeoutSurface(double Tdec) {
50  cout << "check the position of the freeze-out surface "
51  << "at (T = " << Tdec << " GeV) ... " << endl;
52 
53  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
54  int nx = static_cast<int>(abs(2*grid_x0)/grid_dx) + 1;
55  int ny = static_cast<int>(abs(2*grid_y0)/grid_dy) + 1;
56 
57  double tau_local;
58  hydrofluidCell* fluidCellptr = new hydrofluidCell();
59  ofstream output;
60  output.open("results/checkFreezeoutSurface.dat", std::ofstream::app);
61 
62  for (int itime = 0; itime < ntime; itime++) {
63  // loop over time evolution
64  tau_local = grid_tau0 + itime*grid_dt;
65  for (int i = 0; i < nx; i++) {
66  // loops over the transverse plane
67  double x_local = grid_x0 + i*grid_dx;
68  for (int j = 0; j < ny; j++) {
69  double y_local = grid_y0 + j*grid_dy;
70  //hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
71  // fluidCellptr);
72  double temp_local = fluidCellptr->temperature;
73  if (fabs(temp_local - Tdec) < 0.001) {
74  output << tau_local << " " << x_local << " "
75  << y_local << " "
76  << sqrt(x_local*x_local + y_local*y_local) << endl;
77  }
78  }
79  }
80  }
81  output.close();
82  delete fluidCellptr;
83  return;
84 }
85 
86 void FluidcellStatistic::output_temperature_vs_tau() {
87  cout << "output temperature vs tau ..." << endl;
88 
89  int nt = static_cast<int>((grid_tauf - grid_tau0)/grid_dt + 1);
90  int nx = static_cast<int>(abs(2*grid_x0)/grid_dx) + 1;
91  int ny = static_cast<int>(abs(2*grid_y0)/grid_dy) + 1;
92 
93  ofstream output;
94  output.open("results/tau_vs_T.dat");
95 
96  hydrofluidCell* fluidCellptr = new hydrofluidCell;
97  for (int it = 0; it < nt; it++) {
98  double tau_local = grid_tau0 + it*grid_dt;
99  double avgTemp, stdTemp;
100  int numCell = 0;
101  double tempSum = 0.0;
102  double tempSumsq = 0.0;
103  double tempMax = 0.0;
104  for (int i = 0; i < nx; i++) {
105  // loops over the transverse plane
106  double x_local = grid_x0 + i*grid_dx;
107  for (int j = 0; j < ny; j++) {
108  double y_local = grid_y0 + j*grid_dy;
109  // get hydro information
110  if (hydro_type == 0) {
111 #ifdef USE_HDF5
112  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
113  fluidCellptr);
114 #endif
115  } else {
116  hydroinfo_MUSIC_ptr->getHydroValues(
117  x_local, y_local, 0.0, tau_local, fluidCellptr);
118  }
119  double temp_local = fluidCellptr->temperature;
120  if (temp_local > T_dec) {
121  numCell += 1;
122  tempSum += temp_local;
123  tempSumsq += temp_local*temp_local;
124  if (tempMax < temp_local)
125  tempMax = temp_local;
126  }
127  }
128  }
129  if (numCell == 0) {
130  avgTemp = 0.0;
131  stdTemp = 0.0;
132  } else {
133  avgTemp = tempSum/numCell;
134  stdTemp = sqrt((tempSumsq - tempSum*tempSum/numCell)/(numCell - 1));
135  }
136  output << tau_local << " " << avgTemp << " " << stdTemp
137  << " " << tempMax << endl;
138  }
139  delete fluidCellptr;
140  return;
141 }
142 
143 void FluidcellStatistic::output_flowvelocity_vs_tau() {
144  cout << "output u^tau vs tau ... " << endl;
145 
146  int nt = static_cast<int>((grid_tauf - grid_tau0)/grid_dt + 1);
147  int nx = static_cast<int>(abs(2.*grid_x0)/grid_dx) + 1;
148  int ny = static_cast<int>(abs(2.*grid_y0)/grid_dy) + 1;
149 
150  hydrofluidCell* fluidCellptr = new hydrofluidCell;
151 
152  ofstream output;
153  output.open("results/tau_vs_v.dat", std::ofstream::app);
154  output << "# tau (fm) V4 (fm^4) <u^tau> theta" << endl;
155 
156  for (int it = 1; it < nt; it++) {
157  double tau_local = grid_tau0 + it*grid_dt;
158  double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
159  double V4 = 0.0;
160  double avg_utau = 0.0;
161  double avg_theta = 0.0;
162  int count = 0;
163  for (int i = 0; i < nx; i++) {
164  // loops over the transverse plane
165  double x_local = grid_x0 + i*grid_dx;
166  for (int j = 0; j < ny; j++) {
167  double y_local = grid_y0 + j*grid_dy;
168  if (hydro_type == 0) {
169 #ifdef USE_HDF5
170  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
171  fluidCellptr);
172 #endif
173  } else {
174  hydroinfo_MUSIC_ptr->getHydroValues(
175  x_local, y_local, 0.0, tau_local, fluidCellptr);
176  }
177  double temp_local = fluidCellptr->temperature;
178  if (temp_local > T_dec) {
179  // inside freeze-out surface
180  double vx_local = fluidCellptr->vx;
181  double vy_local = fluidCellptr->vy;
182  double v_perp = sqrt(vx_local*vx_local + vy_local*vy_local);
183  double gamma = 1./sqrt(1. - v_perp*v_perp);
184 
185  double theta = compute_local_expansion_rate(
186  tau_local, x_local, y_local);
187  avg_utau += gamma*volume_element;
188  avg_theta += theta*volume_element;
189  V4 += volume_element;
190  count++;
191  }
192  }
193  }
194  if (count < 1) {
195  avg_utau = 1.0;
196  avg_theta = 0.0;
197  V4 = 0.0;
198  } else {
199  avg_utau = avg_utau/V4;
200  avg_theta = avg_theta/V4;
201  }
202  output << tau_local << " " << V4 << " " << avg_utau << " "
203  << avg_theta << endl;
204  }
205  return;
206 }
207 
208 void FluidcellStatistic::output_temperature_vs_avg_utau() {
209  cout << "output utau vs T ... " << endl;
210 
211  int nt = static_cast<int>((grid_tauf - grid_tau0)/grid_dt + 1);
212  int nx = static_cast<int>(abs(2*grid_x0)/grid_dx) + 1;
213  int ny = static_cast<int>(abs(2*grid_y0)/grid_dy) + 1;
214 
215  ofstream output("results/T_vs_v.dat");
216  output << "# T (GeV) V (fm^4) <u^tau> theta" << endl;
217 
218  int n_bin = 41;
219  double T_max = 0.5;
220  double dT = (T_max - T_dec)/(n_bin - 1);
221  double *T_bin_avg = new double[n_bin];
222  double *V4 = new double[n_bin];
223  double *avg_utau = new double[n_bin];
224  double *avg_theta = new double[n_bin];
225  for (int i = 0; i < n_bin; i++) {
226  T_bin_avg[i] = 0.0;
227  V4[i] = 0.0;
228  avg_utau[i] = 0.0;
229  avg_theta[i] = 0.0;
230  }
231 
232  hydrofluidCell* fluidCellptr = new hydrofluidCell;
233  for (int it = 1; it < nt; it++) {
234  double tau_local = grid_tau0 + it*grid_dt;
235  double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
236  for (int i = 0; i < nx; i++) {
237  // loops over the transverse plane
238  double x_local = grid_x0 + i*grid_dx;
239  for (int j = 0; j < ny; j++) {
240  double y_local = grid_y0 + j*grid_dy;
241 
242  // get hydro information
243  if (hydro_type == 0) {
244 #ifdef USE_HDF5
245  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
246  fluidCellptr);
247 #endif
248  } else {
249  hydroinfo_MUSIC_ptr->getHydroValues(
250  x_local, y_local, 0.0, tau_local, fluidCellptr);
251  }
252  double temp_local = fluidCellptr->temperature;
253  if (temp_local > T_dec) {
254  double vx_local = fluidCellptr->vx;
255  double vy_local = fluidCellptr->vy;
256  double utau_local = 1./sqrt(1. - vx_local*vx_local
257  - vy_local*vy_local);
258  double theta_local = compute_local_expansion_rate(
259  tau_local, x_local, y_local);
260  int T_idx = static_cast<int>((temp_local - T_dec)/dT);
261  if (T_idx >= 0 && T_idx < n_bin-1) {
262  T_bin_avg[T_idx] += temp_local*volume_element;
263  avg_utau[T_idx] += utau_local*volume_element;
264  avg_theta[T_idx] += theta_local*volume_element;
265  V4[T_idx] += volume_element;
266  }
267  }
268  }
269  }
270  }
271  for (int i = 0; i < n_bin; i++) {
272  double T_bin, utau_bin, theta_bin;
273  if (fabs(T_bin_avg[i]) < 1e-10) {
274  T_bin = T_dec + i*dT;
275  utau_bin = 1.0;
276  theta_bin = 0.0;
277  } else {
278  T_bin = T_bin_avg[i]/(V4[i] + 1e-15);
279  utau_bin = avg_utau[i]/(V4[i] + 1e-15);
280  theta_bin = avg_theta[i]/(V4[i] + 1e-15);
281  }
282  output << scientific << setw(18) << setprecision(8)
283  << T_bin << " " << V4[i] << " " << utau_bin << " "
284  << theta_bin << endl;
285  }
286  output.close();
287 
288  delete[] T_bin_avg;
289  delete[] V4;
290  delete[] avg_utau;
291  delete[] avg_theta;
292  delete fluidCellptr;
293 
294  return;
295 }
296 
297 void FluidcellStatistic::output_momentum_anisotropy_vs_tau() {
298  cout << "output momentum anisotropy vs tau ... " << endl;
299 
300  int nt = static_cast<int>((grid_tauf - grid_tau0)/grid_dt + 1);
301  int nx = static_cast<int>(abs(2*grid_x0)/grid_dx) + 1;
302  int ny = static_cast<int>(abs(2*grid_y0)/grid_dy) + 1;
303 
304  hydrofluidCell* fluidCellptr = new hydrofluidCell;
305  ofstream output;
306  output.open("results/tau_vs_epsP.dat");
307 
308  for (int it = 0; it < nt; it++) {
309  double tau_local = grid_tau0 + it*grid_dt;
310  double TxxSum = 0.0;
311  double TyySum = 0.0;
312  double epsP = 0.0;
313  for (int i = 0; i < nx; i++) {
314  // loops over the transverse plane
315  double x_local = grid_x0 + i*grid_dx;
316  for(int j = 0; j < ny; j++) {
317  double y_local = grid_y0 + j*grid_dy;
318  // get hydro information
319  if (hydro_type == 0) {
320 #ifdef USE_HDF5
321  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
322  fluidCellptr);
323 #endif
324  } else {
325  hydroinfo_MUSIC_ptr->getHydroValues(
326  x_local, y_local, 0.0, tau_local, fluidCellptr);
327  }
328  double temp_local = fluidCellptr->temperature;
329  if (temp_local > T_dec) {
330  double e_local = fluidCellptr->ed;
331  double p_local = fluidCellptr->pressure;
332  double vx_local = fluidCellptr->vx;
333  double vy_local = fluidCellptr->vy;
334  double v_perp = sqrt(vx_local*vx_local
335  + vy_local*vy_local);
336  double gamma = 1./sqrt(1. - v_perp*v_perp);
337  double ux_local = gamma*vx_local;
338  double uy_local = gamma*vy_local;
339  double pi_xx = fluidCellptr->pi[1][1];
340  double pi_yy = fluidCellptr->pi[2][2];
341  TxxSum += (
342  (e_local+p_local)*ux_local*ux_local + p_local + pi_xx);
343  TyySum += (
344  (e_local+p_local)*uy_local*uy_local + p_local + pi_yy);
345  }
346  }
347  }
348  if (fabs(TxxSum)+ fabs(TyySum) < 1e-10)
349  epsP = 0.0;
350  else
351  epsP = (TxxSum - TyySum)/(TxxSum + TyySum);
352 
353  output << tau_local-0.6 << " " << epsP << endl;
354  }
355  return;
356 }
357 
358 void FluidcellStatistic::outputTempasTauvsX() {
359  cout << "output 2D contour plot for temperature as function of "
360  << "tau and x ... " << endl;
361 
362  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
363  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
364 
365  double tau_local;
366  hydrofluidCell* fluidCellptr = new hydrofluidCell();
367  ofstream output;
368  output.open("results/TempasTauvsX.dat");
369 
370  for (int itime = 0; itime < ntime; itime++) {
371  // loop over time evolution
372  tau_local = grid_tau0 + itime*grid_dt;
373  for (int i = 0; i < nx; i++) {
374  // loops over the transverse plane
375  double x_local = grid_x0 + i*grid_dx;
376  double y_local = 0.0;
377  // get hydro information
378  if (hydro_type == 0) {
379 #ifdef USE_HDF5
380  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
381  fluidCellptr);
382 #endif
383  } else {
384  hydroinfo_MUSIC_ptr->getHydroValues(
385  x_local, y_local, 0.0, tau_local, fluidCellptr);
386  }
387  double temp_local = fluidCellptr->temperature;
388  if (temp_local > 0.05)
389  output << temp_local << " " ;
390  else
391  output << 0.0 << " " ;
392  }
393  output << endl;
394  }
395  output.close();
396  delete fluidCellptr;
397  return;
398 }
399 
400 
401 void FluidcellStatistic::outputKnudersonNumberasTauvsX() {
402  cout << "output Knudersen Number as a function of tau and x ..." << endl;
403 
404  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
405  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx);
406 
407  double eps = 1e-15;
408 
409  hydrofluidCell* fluidCellptr = new hydrofluidCell;
410  ofstream output;
411  output.open("results/KnudsenNumberasTauvsX.dat");
412 
413  for (int itime = 1; itime < ntime; itime++) {
414  // loop over time evolution
415  double tau_local = grid_tau0 + itime*grid_dt;
416  for (int i = 0; i < nx; i++) {
417  // loops over the transverse plane
418  double x_local = grid_x0 + i*grid_dx;
419  double y_local = 0.0;
420  double theta = compute_local_expansion_rate(tau_local,
421  x_local, y_local);
422 
423  if (hydro_type == 0) {
424 #ifdef USE_HDF5
425  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
426  fluidCellptr);
427 #endif
428  } else {
429  hydroinfo_MUSIC_ptr->getHydroValues(
430  x_local, y_local, 0.0, tau_local, fluidCellptr);
431  }
432 
433  double eta_s = 0.08;
434  double L_micro = (5.*eta_s
435  /(fabs(fluidCellptr->temperature) + eps))*hbarC;
436  double L_macro = 1/(fabs(theta) + eps);
437  double Knudsen = L_micro/L_macro;
438 
439  output << Knudsen << " ";
440  }
441  output << endl;
442  }
443  output.close();
444  delete fluidCellptr;
445  return;
446 }
447 
448 double FluidcellStatistic::compute_local_expansion_rate(
449  double tau_local, double x_local, double y_local) {
450  // this function computes the local expansion rate at the given
451  // space-time point (tau_loca, x_local, y_local)
452  // theta = \patial_mu u^\mu + u^tau/tau
453 
454  hydrofluidCell* fluidCellptr = new hydrofluidCell;
455  hydrofluidCell* fluidCellptrt1 = new hydrofluidCell;
456  hydrofluidCell* fluidCellptrt2 = new hydrofluidCell;
457  hydrofluidCell* fluidCellptrx1 = new hydrofluidCell;
458  hydrofluidCell* fluidCellptrx2 = new hydrofluidCell;
459  hydrofluidCell* fluidCellptry1 = new hydrofluidCell;
460  hydrofluidCell* fluidCellptry2 = new hydrofluidCell;
461 
462  if (hydro_type == 0) {
463 #ifdef USE_HDF5
464  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
465  fluidCellptr);
466  hydroinfo_ptr->getHydroinfo(tau_local-grid_dt, x_local, y_local,
467  fluidCellptrt1);
468  hydroinfo_ptr->getHydroinfo(tau_local+grid_dt, x_local, y_local,
469  fluidCellptrt2);
470  hydroinfo_ptr->getHydroinfo(tau_local, x_local-grid_dx, y_local,
471  fluidCellptrx1);
472  hydroinfo_ptr->getHydroinfo(tau_local, x_local+grid_dx, y_local,
473  fluidCellptrx2);
474  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local-grid_dy,
475  fluidCellptry1);
476  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local+grid_dy,
477  fluidCellptry2);
478 #endif
479  } else {
480  hydroinfo_MUSIC_ptr->getHydroValues(
481  x_local, y_local, 0.0, tau_local, fluidCellptr);
482  hydroinfo_MUSIC_ptr->getHydroValues(
483  x_local, y_local, 0.0, tau_local - grid_dt, fluidCellptrt1);
484  hydroinfo_MUSIC_ptr->getHydroValues(
485  x_local, y_local, 0.0, tau_local + grid_dt, fluidCellptrt2);
486  hydroinfo_MUSIC_ptr->getHydroValues(
487  x_local - grid_dx, y_local, 0.0, tau_local, fluidCellptrx1);
488  hydroinfo_MUSIC_ptr->getHydroValues(
489  x_local + grid_dx, y_local, 0.0, tau_local, fluidCellptrx2);
490  hydroinfo_MUSIC_ptr->getHydroValues(
491  x_local, y_local - grid_dy, 0.0, tau_local, fluidCellptry1);
492  hydroinfo_MUSIC_ptr->getHydroValues(
493  x_local, y_local + grid_dy, 0.0, tau_local, fluidCellptry2);
494  }
495 
496  double u0 = 1./sqrt(1. - fluidCellptr->vx*fluidCellptr->vx
497  + fluidCellptr->vy*fluidCellptr->vy);
498  double u0t1 = 1./sqrt(1. - fluidCellptrt1->vx*fluidCellptrt1->vx
499  + fluidCellptrt1->vy*fluidCellptrt1->vy);
500  double u0t2 = 1./sqrt(1. - fluidCellptrt2->vx*fluidCellptrt2->vx
501  + fluidCellptrt2->vy*fluidCellptrt2->vy);
502  double u1x1 = (fluidCellptrx1->vx
503  /sqrt(1. - fluidCellptrx1->vx*fluidCellptrx1->vx
504  + fluidCellptrx1->vy*fluidCellptrx1->vy));
505  double u1x2 = (fluidCellptrx2->vx
506  /sqrt(1. - fluidCellptrx2->vx*fluidCellptrx2->vx
507  + fluidCellptrx2->vy*fluidCellptrx2->vy));
508  double u2y1 = (fluidCellptry1->vy
509  /sqrt(1. - fluidCellptry1->vx*fluidCellptry1->vx
510  + fluidCellptry1->vy*fluidCellptry1->vy));
511  double u2y2 = (fluidCellptry2->vy
512  /sqrt(1. - fluidCellptry2->vx*fluidCellptry2->vx
513  + fluidCellptry2->vy*fluidCellptry2->vy));
514 
515  double d0u0 = (u0t2 - u0t1)/2./grid_dt;
516  double d1u1 = (u1x2 - u1x1)/2./grid_dx;
517  double d2u2 = (u2y2 - u2y1)/2./grid_dy;
518  double theta = (d0u0 + d1u1 + d2u2 + u0/tau_local);
519 
520  delete fluidCellptr;
521  delete fluidCellptrt1;
522  delete fluidCellptrt2;
523  delete fluidCellptrx1;
524  delete fluidCellptrx2;
525  delete fluidCellptry1;
526  delete fluidCellptry2;
527  return(theta);
528 }
529 
530 void FluidcellStatistic::outputinverseReynoldsNumberasTauvsX() {
531  cout << "output inverse Reynolds number as a function of "
532  << "tau and x ... " << endl;
533 
534  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
535  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
536 
537  double MAX = 1000.;
538 
539  hydrofluidCell* fluidCellptr = new hydrofluidCell();
540  ofstream output;
541  output.open("results/inverseReynoldsNumberasTauvsX.dat");
542 
543  for (int itime = 0; itime < ntime; itime++) {
544  // loop over time evolution
545  double tau_local = grid_tau0 + itime*grid_dt;
546  for(int i = 0; i < nx; i++) {
547  // loops over the transverse plane
548  double x_local = grid_x0 + i*grid_dx;
549  double y_local = 0.0;
550 
551  if (hydro_type == 0) {
552 #ifdef USE_HDF5
553  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
554  fluidCellptr);
555 #endif
556  } else {
557  hydroinfo_MUSIC_ptr->getHydroValues(
558  x_local, y_local, 0.0, tau_local, fluidCellptr);
559  }
560  double pi2 = (
561  fluidCellptr->pi[0][0]*fluidCellptr->pi[0][0]
562  + fluidCellptr->pi[1][1]*fluidCellptr->pi[1][1]
563  + fluidCellptr->pi[2][2]*fluidCellptr->pi[2][2]
564  + fluidCellptr->pi[3][3]*fluidCellptr->pi[3][3]
565  - 2.*(fluidCellptr->pi[0][1]*fluidCellptr->pi[0][1]
566  + fluidCellptr->pi[0][2]*fluidCellptr->pi[0][2]
567  + fluidCellptr->pi[0][3]*fluidCellptr->pi[0][3])
568  + 2.*(fluidCellptr->pi[1][2]*fluidCellptr->pi[1][2]
569  + fluidCellptr->pi[1][3]*fluidCellptr->pi[1][3]
570  + fluidCellptr->pi[2][3]*fluidCellptr->pi[2][3]));
571 
572  double inverseReynold;
573 
574  if (pi2 >= 0)
575  inverseReynold = sqrt(pi2)/(fluidCellptr->ed
576  + fluidCellptr->pressure + 1e-15);
577  else
578  inverseReynold = MAX;
579 
580  output << inverseReynold << " " ;
581  }
582  output << endl;
583  }
584  output.close();
585  delete fluidCellptr;
586  return;
587 }
588 
589 void FluidcellStatistic::outputBulkinverseReynoldsNumberasTauvsX() {
590  cout << "output bulk inverse Reynolds number as a function of "
591  << "tau and x ... " << endl;
592 
593  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
594  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
595 
596  hydrofluidCell* fluidCellptr = new hydrofluidCell();
597  ofstream output;
598  output.open("results/inverseReynoldsNumberasTauvsX.dat");
599 
600  for (int itime = 0 ; itime < ntime; itime++) {
601  // loop over time evolution
602  double tau_local = grid_tau0 + itime*grid_dt;
603  for (int i = 0; i < nx; i++) {
604  // loops over the transverse plane
605  double x_local = grid_x0 + i*grid_dx;
606  double y_local = 0.0;
607 
608  if (hydro_type == 0) {
609 #ifdef USE_HDF5
610  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
611  fluidCellptr);
612 #endif
613  } else {
614  hydroinfo_MUSIC_ptr->getHydroValues(
615  x_local, y_local, 0.0, tau_local, fluidCellptr);
616  }
617 
618  double inverseReynold;
619  inverseReynold = fabs(fluidCellptr->bulkPi)/fluidCellptr->pressure;
620  output << inverseReynold << " " ;
621  }
622  output << endl;
623  }
624  output.close();
625  delete fluidCellptr;
626  return;
627 }
628 
629 void FluidcellStatistic::analysis_hydro_volume_for_photon(double T_cut) {
630  double V_3 = calculate_hypersurface_3volume(T_cut);
631  double V_4 = calculate_spacetime_4volume(T_cut);
632  double average_tau = calculate_average_tau(T_cut);
633  double average_T4 = calculate_average_temperature4(T_cut);
634  double average_GammaT =
635  calculate_average_integrated_photonRate_parameterization(T_cut);
636  stringstream output;
637  output << "results/volume_info_for_photon_Tcut_" << T_cut << ".dat";
638  ofstream of(output.str().c_str());
639  of << "# V_4 <tau> V_3 <T_4> <GammaT> " << endl;
640  of << scientific << setw(18) << setprecision(8)
641  << V_4 << " " << average_tau << " "
642  << V_3 << " " << average_T4 << " " << average_GammaT << endl;
643  of.close();
644 }
645 
646 double FluidcellStatistic::calculate_spacetime_4volume(double T_cut) {
647  // this function calculates the space-time 4 volume of the medium
648  // inside a give temperature T_cut [GeV]
649  // the output volume V_4 is in [fm^4]
650  // deta = 1
651 
652  cout << "compute 4-volume inside T = " << T_cut << " GeV ..." << endl;
653 
654  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
655  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
656  int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy) + 1;
657 
658  hydrofluidCell* fluidCellptr = new hydrofluidCell();
659 
660  double volume = 0.0;
661  for (int itime = 0; itime < ntime; itime++) {
662  // loop over time evolution
663  double tau_local = grid_tau0 + itime*grid_dt;
664  double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
665  for (int i = 0; i < nx; i++) {
666  // loops over the transverse plane
667  double x_local = grid_x0 + i*grid_dx;
668  for (int j = 0; j < ny; j++) {
669  double y_local = grid_y0 + j*grid_dy;
670  // get hydro information
671  if (hydro_type == 0) {
672 #ifdef USE_HDF5
673  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
674  fluidCellptr);
675 #endif
676  } else {
677  hydroinfo_MUSIC_ptr->getHydroValues(
678  x_local, y_local, 0.0, tau_local, fluidCellptr);
679  }
680  double T_local = fluidCellptr->temperature; // GeV
681  if (T_local > T_cut) {
682  volume += volume_element;
683  }
684  }
685  }
686  }
687  delete fluidCellptr;
688  return(volume);
689 }
690 
691 double FluidcellStatistic::calculate_average_tau(double T_cut) {
692  // this function calculates the average <\tau> of the medium
693  // inside a give temperature T_cut [GeV]
694  // the output <tau> is in [fm]
695 
696  cout << "compute <tau> ... " << endl;
697 
698  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
699  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
700  int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy) + 1;
701 
702  hydrofluidCell* fluidCellptr = new hydrofluidCell();
703 
704  double average_tau = 0.0;
705  double volume = 0.0;
706  for (int itime = 0; itime < ntime; itime++) {
707  // loop over time evolution
708  double tau_local = grid_tau0 + itime*grid_dt;
709  double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
710  for (int i = 0; i < nx; i++) {
711  // loops over the transverse plane
712  double x_local = grid_x0 + i*grid_dx;
713  for (int j = 0; j < ny; j++) {
714  double y_local = grid_y0 + j*grid_dy;
715  // get hydro information
716  if (hydro_type == 0) {
717 #ifdef USE_HDF5
718  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
719  fluidCellptr);
720 #endif
721  } else {
722  hydroinfo_MUSIC_ptr->getHydroValues(
723  x_local, y_local, 0.0, tau_local, fluidCellptr);
724  }
725  double T_local = fluidCellptr->temperature; // GeV
726  if (T_local > T_cut) {
727  volume += volume_element;
728  average_tau += tau_local*volume_element;
729  }
730  }
731  }
732  }
733  average_tau /= volume;
734  delete fluidCellptr;
735  return(average_tau);
736 }
737 
738 double FluidcellStatistic::calculate_average_temperature4(double T_cut) {
739  // this function calculates the average <T^4> of the medium
740  // inside a give temperature T_cut [GeV]
741  // the output <T^4> is in [GeV^4]
742 
743  cout << "compute <T^4> ..." << endl;
744  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
745  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
746  int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy) + 1;
747 
748  hydrofluidCell* fluidCellptr = new hydrofluidCell();
749 
750  double average_T4 = 0.0;
751  double volume = 0.0;
752  for (int itime = 0; itime < ntime; itime++) {
753  // loop over time evolution
754  double tau_local = grid_tau0 + itime*grid_dt;
755  double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
756  for (int i = 0; i < nx; i++) {
757  // loops over the transverse plane
758  double x_local = grid_x0 + i*grid_dx;
759  for (int j = 0; j < ny; j++) {
760  double y_local = grid_y0 + j*grid_dy;
761  // get hydro information
762  if (hydro_type == 0) {
763 #ifdef USE_HDF5
764  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
765  fluidCellptr);
766 #endif
767  } else {
768  hydroinfo_MUSIC_ptr->getHydroValues(
769  x_local, y_local, 0.0, tau_local, fluidCellptr);
770  }
771  double T_local = fluidCellptr->temperature; // GeV
772  if (T_local > T_cut) {
773  volume += volume_element;
774  average_T4 += (
775  T_local*T_local*T_local*T_local*volume_element);
776  }
777  }
778  }
779  }
780  average_T4 /= volume;
781  delete fluidCellptr;
782  return(average_T4);
783 }
784 
785 double FluidcellStatistic::
786  calculate_average_integrated_photonRate_parameterization(double T_cut) {
787  // this function calculates the average Gamma(T) of the medium
788  // Gamma(T) = 1.746 T^4.095 for T > Tsw = 0.18
789  // = 11663.923 T^9.024 for T < Tsw
790  // inside a give temperature T_cut [GeV]
791  // the output <Gamma(T)> is in [1/fm^4]
792 
793  cout << "compute <Gamma(T)> ..." << endl;
794  int ntime = static_cast<int>((grid_tauf - grid_tau0)/grid_dt) + 1;
795  int nx = static_cast<int>(fabs(2.*grid_x0)/grid_dx) + 1;
796  int ny = static_cast<int>(fabs(2.*grid_y0)/grid_dy) + 1;
797 
798  hydrofluidCell* fluidCellptr = new hydrofluidCell();
799 
800  double average_GammaT = 0.0;
801  double volume = 0.0;
802  double Tsw = 0.18; // switching temperature for QGP rates to HG rates
803  for (int itime = 0; itime < ntime; itime++) {
804  // loop over time evolution
805  double tau_local = grid_tau0 + itime*grid_dt;
806  double volume_element = tau_local*grid_dt*grid_dx*grid_dy;
807  for (int i = 0; i < nx; i++) {
808  // loops over the transverse plane
809  double x_local = grid_x0 + i*grid_dx;
810  for (int j = 0; j < ny; j++) {
811  double y_local = grid_y0 + j*grid_dy;
812  // get hydro information
813  if (hydro_type == 0) {
814 #ifdef USE_HDF5
815  hydroinfo_ptr->getHydroinfo(tau_local, x_local, y_local,
816  fluidCellptr);
817 #endif
818  } else {
819  hydroinfo_MUSIC_ptr->getHydroValues(
820  x_local, y_local, 0.0, tau_local, fluidCellptr);
821  }
822  double T_local = fluidCellptr->temperature; // GeV
823  if (T_local > T_cut) {
824  volume += volume_element;
825  double GammaT = 0.0;
826  if (T_local > Tsw) {
827  GammaT = 1.746*pow(T_local, 4.095);
828  } else {
829  GammaT = 11663.923*pow(T_local, 9.024);
830  }
831  average_GammaT += GammaT*volume_element;
832  }
833  }
834  }
835  }
836  average_GammaT /= volume;
837  delete fluidCellptr;
838  return(average_GammaT);
839 }
840 
841 double FluidcellStatistic::calculate_hypersurface_3volume(double T_cut) {
842  // this function calculates the surface area of an isothermal hyper-surface
843  // at a give temperature T_cut [GeV]
844  // the output volume V_3 = int u^\mu d^3\sigma_\mu is in [fm^3]
845  // deta = 1
846 
847  cout << "compute 3-volume at T = " << T_cut << " GeV ... " << endl;
848  // first find the hyper-surface
849  void* hydro_ptr = NULL;
850  if (hydro_type == 1) {
851 #ifdef USE_HDF5
852  hydro_ptr = hydroinfo_ptr;
853 #endif
854  } else {
855  hydro_ptr = hydroinfo_MUSIC_ptr;
856  }
857  SurfaceFinder surfFinder(hydro_ptr, paraRdr, T_cut);
858  surfFinder.Find_full_hypersurface();
859 
860  // load the hyper-surface file
861  ifstream surf("hyper_surface_2+1d.dat");
862  if (!surf.is_open()) {
863  cout << "FluidcellStatistic::calculate_hypersurface_3volume:"
864  << "can not open the file hyper_surface_2+1d.dat!" << endl;
865  exit(1);
866  }
867 
868  double volume = 0.0;
869  double tau, x, y, da0, da1, da2, vx, vy, T;
870  surf >> tau >> x >> y >> da0 >> da1 >> da2 >> T >> vx >> vy;
871  while (!surf.eof()) {
872  double u0 = 1./sqrt(1. - vx*vx - vy*vy);
873  double ux = u0*vx;
874  double uy = u0*vy;
875  volume += tau*(u0*da0 + ux*da1 + uy*da2);
876  surf >> tau >> x >> y >> da0 >> da1 >> da2 >> T >> vx >> vy;
877  }
878  return(volume);
879 }
880