Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Hydroinfo_MUSIC.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Hydroinfo_MUSIC.cpp
1 // Hydroinfo_MUSIC.cpp is a part of the MARTINI event generator.
2 // Copyright (C) 2009-2010 Bjoern Schenke.
3 // MARTINI is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // This file contains routines to read in hydro data from files and functions
7 // that return interpolated data at a given space-time point
8 
9 #include <iostream>
10 #include <fstream>
11 #include <iomanip>
12 #include <cstdlib>
13 #include <cmath>
14 #include <sstream>
15 #include <vector>
16 #include <string>
17 
18 #include "./Hydroinfo_MUSIC.h"
19 
20 using namespace std;
21 
22 Hydroinfo_MUSIC::Hydroinfo_MUSIC() {
23  verbose_ = 9;
24  hbarC = 0.19733;
25  boost_invariant = false;
26 }
27 
28 Hydroinfo_MUSIC::~Hydroinfo_MUSIC() {
29  clean_hydro_event();
30 }
31 
32 void Hydroinfo_MUSIC::clean_hydro_event() {
33  if (boost_invariant) {
34  lattice_2D.clear();
35  } else {
36  lattice_3D.clear();
37  lattice_3D_ideal.clear();
38  }
39 }
40 
41 void Hydroinfo_MUSIC::readHydroData(int whichHydro, int nskip_tau_in,
42  string input_filename_in,
43  string hydro_ideal_filename_in,
44  string hydro_shear_filename_in,
45  string hydro_bulk_filename_in) {
46  // all hydro data is stored in tau steps (not t)
47  // evolution is converted to tau when accessing the hydro data
48  lattice_2D.clear();
49  lattice_3D.clear();
50  lattice_3D_ideal.clear();
51 
52  input_filename = input_filename_in;
53  hydro_ideal_filename = hydro_ideal_filename_in;
54  hydro_shear_filename = hydro_shear_filename_in;
55  hydro_bulk_filename = hydro_bulk_filename_in;
56 
57  // read in setups of the hydro simulation
58  hydroWhichHydro = whichHydro;
59  if (hydroWhichHydro < 10) {
60  ostringstream config_file;
61  config_file << input_filename;
62  ifstream configuration;
63  configuration.open(config_file.str().c_str(), ios::in);
64  if (!configuration) {
65  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
66  << "Unable to open file: " << config_file.str() << endl;
67  exit(1);
68  }
69  string temp1;
70  string temp_name;
71  while (!configuration.eof()) {
72  getline(configuration, temp1);
73  stringstream ss(temp1);
74  ss >> temp_name;
75 
76  // read in grid information
77  if (temp_name == "Initial_time_tau_0") {
78  ss >> hydroTau0;
79  } else if (temp_name == "Delta_Tau") {
80  ss >> hydroDtau;
81  } else if (temp_name == "X_grid_size_in_fm") {
82  double temp;
83  ss >> temp;
84  hydroXmax = temp/2.;
85  } else if (temp_name == "Grid_size_in_x") {
86  ss >> ixmax;
87  } else if (temp_name == "Eta_grid_size") {
88  double temp;
89  ss >> temp;
90  hydro_eta_max = temp/2.;
91  } else if (temp_name == "Grid_size_in_eta") {
92  ss >> ietamax;
93  } else if (temp_name == "output_evolution_every_N_timesteps") {
94  ss >> nskip_tau;
95  } else if (temp_name == "output_evolution_every_N_x") {
96  ss >> nskip_x;
97  } else if (temp_name == "output_evolution_every_N_eta") {
98  ss >> nskip_eta;
99  }
100  // read in additioinal information
101  if (temp_name == "Include_Rhob_Yes_1_No_0") {
102  ss >> turn_on_rhob;
103  } else if (temp_name == "Include_Shear_Visc_Yes_1_No_0") {
104  ss >> turn_on_shear;
105  } else if (temp_name == "Include_Bulk_Visc_Yes_1_No_0") {
106  ss >> turn_on_bulk;
107  } else if (temp_name == "turn_on_baryon_diffusion") {
108  ss >> turn_on_diff;
109  }
110  }
111  configuration.close();
112  hydroDx = 2.*hydroXmax/(ixmax - 1.);
113  hydroDeta = 2.*hydro_eta_max/(static_cast<double>(ietamax));
114  }
115 
116  use_tau_eta_coordinate = 1;
117  if (use_tau_eta_coordinate == 0) {
118  cout << "Hydroinfo_MUSIC:: Warning hydro grid is set to "
119  << "cartesian coordinates, please make sure this is correct!"
120  << endl;
121  }
122 
123  if (whichHydro == 6) {
124  // 3+1D MUSIC hydro (Schenke, Jeon, Gale)
125  cout << "Using 3+1D Jeon Schenke hydro reading data ..." << endl;
126  boost_invariant = false;
127 
128  ixmax = static_cast<int>(2.*hydroXmax/hydroDx + 0.001);
129  ietamax = static_cast<int>(2.*hydro_eta_max/hydroDeta + 0.001);
130 
131  // read in temperature, QGP fraction , flow velocity
132  // The name of the evolution file: evolution_name
133  string evolution_name = hydro_ideal_filename;
134  cout << "Evolution file name = " << evolution_name << endl;
135  ifstream fin;
136  fin.open(evolution_name.c_str(), ios::in);
137  if (!fin) {
138  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
139  << "Unable to open file: " << evolution_name << endl;
140  exit(1);
141  }
142 
143  double T, vx, vy, vz, QGPfrac;
144  fluidCell_3D newCell;
145  int ik = 0;
146  while (!fin.eof()) {
147  ik++;
148  fin >> T;
149  fin >> QGPfrac;
150  fin >> vx;
151  fin >> vy;
152  fin >> vz;
153  newCell.temperature = T;
154  newCell.vx = vx;
155  newCell.vy = vy;
156  newCell.vz = vz;
157 
158  lattice_3D.push_back(newCell);
159  if (ik%50000 == 0)
160  cout << "o" << flush;
161  }
162  cout << ik << endl;
163  fin.close();
164  } else if (whichHydro == 8) {
165  // event-by-event (2+1)-d MUSIC hydro from JF
166  // there are two slices in medium in eta_s
167  // one at eta_s = -15. and the other at eta_s = 0.0
168  // only the medium at middle rapidity will be kept in the memory
169  boost_invariant = true;
170  cout << "Reading event-by-event hydro evolution data "
171  << "from standard MUSIC ..." << endl;
172 
173  ixmax = static_cast<int>(2.*hydroXmax/hydroDx + 0.001);
174  ietamax = 1;
175  nskip_tau = nskip_tau_in;
176 
177  hydroDx *= nskip_x;
178  hydroDtau *= nskip_tau;
179  hydroDeta *= nskip_eta;
180 
181  int n_eta = 2; // there are two slices in eta_s
182  // number of fluid cell in the transverse plane
183  int num_fluid_cell_trans = ixmax*ixmax;
184 
185  // read in hydro evolution
186  string evolution_name = hydro_ideal_filename;
187  string evolution_name_Wmunu = hydro_shear_filename;
188  string evolution_name_Pi = hydro_bulk_filename;
189 
190  std::FILE *fin;
191  string evolution_file_name = evolution_name;
192  cout << "Evolution file name = " << evolution_file_name << endl;
193  fin = std::fopen(evolution_file_name.c_str(), "rb");
194 
195  if (fin == NULL) {
196  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
197  << "Unable to open file: " << evolution_file_name << endl;
198  exit(1);
199  }
200 
201  std::FILE *fin1 = NULL;
202  if (turn_on_shear == 1) {
203  fin1 = std::fopen(evolution_name_Wmunu.c_str(), "rb");
204  if (fin1 == NULL) {
205  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
206  << "Unable to open file: " << evolution_name_Wmunu << endl;
207  exit(1);
208  }
209  }
210 
211  std::FILE *fin2 = NULL;
212  if (turn_on_bulk == 1) {
213  fin2 = std::fopen(evolution_name_Pi.c_str(), "rb");
214  if (fin2 == NULL) {
215  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
216  << "Unable to open file: " << evolution_name_Pi << endl;
217  exit(1);
218  }
219  }
220 
221  int ik = 0;
222  fluidCell_2D newCell;
223  float T, QGPfrac, vx, vy, vz;
224  double ux, uy, ueta;
225  float pi00 = 0.0;
226  float pi01 = 0.0;
227  float pi02 = 0.0;
228  float pi03 = 0.0;
229  float pi11 = 0.0;
230  float pi12 = 0.0;
231  float pi13 = 0.0;
232  float pi22 = 0.0;
233  float pi23 = 0.0;
234  float pi33 = 0.0;;
235  float bulkPi = 0.0;
236  float e_plus_P = 1e-15;
237  float cs2 = 0.0;
238  int size = sizeof(float);
239  while (true) {
240  int status = 0;
241  status = std::fread(&T, size, 1, fin);
242  status += std::fread(&QGPfrac, size, 1, fin);
243  status += std::fread(&vx, size, 1, fin);
244  status += std::fread(&vy, size, 1, fin);
245  status += std::fread(&vz, size, 1, fin);
246 
247  if (status != 5) { // this is the end of file
248  break;
249  }
250 
251  int status_pi = 0;
252  if (turn_on_shear == 1) {
253  status_pi = std::fread(&pi00, size, 1, fin1);
254  status_pi += std::fread(&pi01, size, 1, fin1);
255  status_pi += std::fread(&pi02, size, 1, fin1);
256  status_pi += std::fread(&pi03, size, 1, fin1);
257  status_pi += std::fread(&pi11, size, 1, fin1);
258  status_pi += std::fread(&pi12, size, 1, fin1);
259  status_pi += std::fread(&pi13, size, 1, fin1);
260  status_pi += std::fread(&pi22, size, 1, fin1);
261  status_pi += std::fread(&pi23, size, 1, fin1);
262  status_pi += std::fread(&pi33, size, 1, fin1);
263 
264  if (status_pi != 10) {
265  cout << "Error:Hydroinfo_MUSIC::readHydroData: "
266  << "Wmunu file does not have the same number of "
267  << "fluid cells as the ideal file!" << endl;
268  exit(1);
269  }
270  }
271 
272  int status_bulkPi = 0;
273  if (turn_on_bulk == 1) {
274  status_bulkPi = std::fread(&bulkPi, size, 1, fin2);
275  status_bulkPi += std::fread(&e_plus_P, size, 1, fin2);
276  status_bulkPi += std::fread(&cs2, size, 1, fin2);
277 
278  if (status_bulkPi != 3) {
279  cout << "Error:Hydroinfo_MUSIC::readHydroData: "
280  << "bulkPi file does not have the same number of "
281  << "fluid cells as the ideal file!" << endl;
282  exit(1);
283  }
284  }
285 
286  int ieta_idx = static_cast<int>(ik/num_fluid_cell_trans) % n_eta;
287  int itau_idx = static_cast<int>(ik/(num_fluid_cell_trans*n_eta));
288  ik++;
289  if (itau_idx%nskip_tau != 0) // skip in tau
290  continue;
291 
292  // print out tau information
293  double tau_local = hydroTau0 + itau_idx*hydroDtau/nskip_tau;
294  if ((ik-1)%(num_fluid_cell_trans*n_eta) == 0) {
295  cout << "read in tau frame: " << itau_idx
296  << " tau_local = " << setprecision(3) << tau_local
297  << " fm ..."<< endl;
298  }
299 
300  if (ieta_idx == (n_eta-1)) {
301  // store the hydro medium at eta_s = 0.0
302  double v2 = vx*vx + vy*vy + vz*vz;
303  if (v2 > 1.0) {
304  cerr << "[Hydroinfo_MUSIC::readHydroData:] Error: "
305  << "v > 1! vx = " << vx << ", vy = " << vy
306  << ", vz = " << vz << ", T = " << T << endl;
307  if (T > 0.01) {
308  exit(1);
309  } else {
310  v2 = 0.0;
311  }
312  }
313  double gamma = 1./sqrt(1. - v2);
314  ux = gamma*vx;
315  uy = gamma*vy;
316  ueta = gamma*vz; // assuming eta = 0
317 
318  newCell.temperature = T;
319  // convert vx and vy to longitudinal co-moving frame
320  newCell.ux = ux;
321  newCell.uy = uy;
322  newCell.ueta = ueta;
323 
324  // pi^\mu\nu tensor
325  newCell.pi00 = pi00;
326  newCell.pi01 = pi01;
327  newCell.pi02 = pi02;
328  newCell.pi11 = pi11;
329  newCell.pi12 = pi12;
330  newCell.pi22 = pi22;
331  newCell.pi33 = pi33;
332 
333  // bulk pressure
334  if (T > 0.18) {
335  // QGP phase prefactor is divided out here
336  newCell.bulkPi = bulkPi/(15.*(1./3. - cs2)*e_plus_P);
337  } else {
338  newCell.bulkPi = bulkPi; // [1/fm^4]
339  }
340  lattice_2D.push_back(newCell);
341  }
342  }
343  std::fclose(fin);
344  if (turn_on_shear == 1) {
345  std::fclose(fin1);
346  }
347  if (turn_on_bulk == 1) {
348  std::fclose(fin2);
349  }
350  cout << endl;
351  cout << "number of fluid cells: " << lattice_2D.size() << endl;
352  } else if (whichHydro == 9) {
353  // event-by-event (2+1)-d MUSIC hydro
354  // the output medium is at middle rapidity
355  boost_invariant = true;
356  cout << "Reading event-by-event hydro evolution data "
357  << "from (2+1)D MUSIC ..." << endl;
358 
359  ietamax = 1;
360 
361  hydroDx *= nskip_x;
362  hydroDtau *= nskip_tau;
363  hydroDeta *= nskip_eta;
364 
365  nskip_tau = nskip_tau_in;
366  ixmax = static_cast<int>(2.*hydroXmax/hydroDx + 0.001);
367 
368  int n_eta = 1;
369  // number of fluid cell in the transverse plane
370  int num_fluid_cell_trans = ixmax*ixmax;
371 
372  // read in hydro evolution
373  string evolution_name = hydro_ideal_filename;
374  string evolution_name_Wmunu = hydro_shear_filename;
375  string evolution_name_Pi = hydro_bulk_filename;
376 
377  std::FILE *fin;
378  string evolution_file_name = evolution_name;
379  cout << "Evolution file name = " << evolution_file_name << endl;
380  fin = std::fopen(evolution_file_name.c_str(), "rb");
381 
382  if (fin == NULL) {
383  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
384  << "Unable to open file: " << evolution_file_name << endl;
385  exit(1);
386  }
387 
388  std::FILE *fin1 = NULL;
389  if (turn_on_shear == 1) {
390  fin1 = std::fopen(evolution_name_Wmunu.c_str(), "rb");
391  if (fin1 == NULL) {
392  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
393  << "Unable to open file: " << evolution_name_Wmunu << endl;
394  exit(1);
395  }
396  }
397 
398  std::FILE *fin2 = NULL;
399  if (turn_on_bulk == 1) {
400  fin2 = std::fopen(evolution_name_Pi.c_str(), "rb");
401  if (fin2 == NULL) {
402  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
403  << "Unable to open file: " << evolution_name_Pi << endl;
404  exit(1);
405  }
406  }
407 
408  int ik = 0;
409  fluidCell_2D newCell;
410  double T, QGPfrac, ux, uy, ueta;
411  double vx, vy, vz;
412  float pi00 = 0.0;
413  float pi01 = 0.0;
414  float pi02 = 0.0;
415  float pi03 = 0.0;
416  float pi11 = 0.0;
417  float pi12 = 0.0;
418  float pi13 = 0.0;
419  float pi22 = 0.0;
420  float pi23 = 0.0;
421  float pi33 = 0.0;;
422  float bulkPi = 0.0;
423  float e_plus_P = 1e-15;
424  float cs2 = 0.0;
425  int size = sizeof(double);
426  while (true) {
427  int status = 0;
428  status = std::fread(&T, size, 1, fin);
429  status += std::fread(&QGPfrac, size, 1, fin);
430  status += std::fread(&vx, size, 1, fin);
431  status += std::fread(&vy, size, 1, fin);
432  status += std::fread(&vz, size, 1, fin);
433  if (status != 5) { // this is the end of file
434  break;
435  }
436 
437  double v2 = vx*vx + vy*vy + vz*vz;
438  if (v2 > 1.) {
439  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
440  << "v > 1! vx = " << vx << ", vy = " << vy
441  << ", vz = " << vz << endl;
442  exit(1);
443  }
444  double gamma = 1./sqrt(1. - v2);
445  ux = vx*gamma;
446  uy = vy*gamma;
447  ueta = vz*gamma; // assuming at the eta = 0
448 
449  int status_pi = 0;
450  if (turn_on_shear == 1) {
451  status_pi = std::fread(&pi00, size, 1, fin1);
452  status_pi += std::fread(&pi01, size, 1, fin1);
453  status_pi += std::fread(&pi02, size, 1, fin1);
454  status_pi += std::fread(&pi03, size, 1, fin1);
455  status_pi += std::fread(&pi11, size, 1, fin1);
456  status_pi += std::fread(&pi12, size, 1, fin1);
457  status_pi += std::fread(&pi13, size, 1, fin1);
458  status_pi += std::fread(&pi22, size, 1, fin1);
459  status_pi += std::fread(&pi23, size, 1, fin1);
460  status_pi += std::fread(&pi33, size, 1, fin1);
461 
462  if (status_pi != 10) {
463  cout << "Error:Hydroinfo_MUSIC::readHydroData: "
464  << "Wmunu file does not have the same number of "
465  << "fluid cells as the ideal file!" << endl;
466  exit(1);
467  }
468  }
469 
470  int status_bulkPi = 0;
471  if (turn_on_bulk == 1) {
472  status_bulkPi = std::fread(&bulkPi, size, 1, fin2);
473  status_bulkPi += std::fread(&e_plus_P, size, 1, fin2);
474  status_bulkPi += std::fread(&cs2, size, 1, fin2);
475 
476  if (status_bulkPi != 3) {
477  cout << "Error:Hydroinfo_MUSIC::readHydroData: "
478  << "bulkPi file does not have the same number of "
479  << "fluid cells as the ideal file!" << endl;
480  exit(1);
481  }
482  }
483 
484  int ieta_idx = static_cast<int>(ik/num_fluid_cell_trans) % n_eta;
485  int itau_idx = static_cast<int>(ik/(num_fluid_cell_trans*n_eta));
486  ik++;
487  if (itau_idx%nskip_tau != 0) // skip in tau
488  continue;
489 
490  // print out tau information
491  double tau_local = hydroTau0 + itau_idx*hydroDtau/nskip_tau;
492  if ((ik-1)%(num_fluid_cell_trans*n_eta) == 0) {
493  cout << "read in tau frame: " << itau_idx
494  << " tau_local = " << setprecision(3) << tau_local
495  << " fm ..."<< endl;
496  }
497 
498  if (ieta_idx == (n_eta-1)) {
499  newCell.temperature = T;
500  newCell.ux = ux;
501  newCell.uy = uy;
502  newCell.ueta = ueta;
503 
504  // pi^\mu\nu tensor
505  newCell.pi00 = pi00;
506  newCell.pi01 = pi01;
507  newCell.pi02 = pi02;
508  newCell.pi11 = pi11;
509  newCell.pi12 = pi12;
510  newCell.pi22 = pi22;
511  newCell.pi33 = pi33;
512 
513  // bulk pressure
514  if (T > 0.18) {
515  // QGP phase prefactor is divided out here
516  newCell.bulkPi = bulkPi/(15.*(1./3. - cs2)*e_plus_P);
517  } else {
518  newCell.bulkPi = bulkPi; // [1/fm^4]
519  }
520  lattice_2D.push_back(newCell);
521  }
522  }
523  std::fclose(fin);
524  if (turn_on_shear == 1) {
525  std::fclose(fin1);
526  }
527  if (turn_on_bulk == 1) {
528  std::fclose(fin2);
529  }
530  cout << endl;
531  cout << "number of fluid cells: " << lattice_2D.size() << endl;
532  } else if (whichHydro == 10 || whichHydro == 11) {
533  // new MUSIC hydro (no regular grid)
534  cout << "Using 3+1D new MUSIC hydro reading data ..." << endl;
535  if (whichHydro == 10) {
536  boost_invariant = false;
537  } else {
538  boost_invariant = true;
539  }
540 
541  // read in temperature and flow velocity
542  // The name of the evolution file: evolution_name
543  string evolution_name = hydro_ideal_filename;
544  cout << "Evolution file name = " << evolution_name << endl;
545  std::FILE *fin;
546  fin = std::fopen(evolution_name.c_str(), "rb");
547  if (fin == NULL) {
548  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
549  << "Unable to open file: " << evolution_name << endl;
550  exit(1);
551  }
552 
553  float header[16];
554  int status = std::fread(&header, sizeof(float), 16, fin);
555  if (status == 0) {
556  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
557  << "Can not read the evolution file header" << endl;
558  exit(1);
559  }
560 
561  hydroTau0 = header[0];
562  hydroDtau = header[1];
563  ixmax = static_cast<int>(header[2]);
564  hydroDx = header[3];
565  hydroXmax = std::abs(header[4]);
566  ietamax = static_cast<int>(header[8]);
567  hydroDeta = header[9];
568  hydro_eta_max = std::abs(header[10]);
569  turn_on_rhob = static_cast<int>(header[11]);
570  turn_on_shear = static_cast<int>(header[12]);
571  turn_on_bulk = static_cast<int>(header[13]);
572  turn_on_diff = static_cast<int>(header[14]);
573  const int nVar_per_cell = static_cast<int>(header[15]);
574 
575  float cell_info[nVar_per_cell];
576 
577  int itau_max = 0;
578  fluidCell_3D_ideal zeroCell;
579  zeroCell.itau = 0;
580  zeroCell.ix = 0;
581  zeroCell.iy = 0;
582  zeroCell.ieta = 0;
583  zeroCell.temperature = 0.;
584  zeroCell.ed = 0.;
585  zeroCell.pressure = 0.;
586  zeroCell.ux = 0.;
587  zeroCell.uy = 0.;
588  zeroCell.uz = 0.;
589  lattice_3D_ideal.push_back(zeroCell);
590  int ik = 0;
591  while (true) {
592  status = 0;
593  status = std::fread(&cell_info, sizeof(float), nVar_per_cell, fin);
594  if (status == 0) break;
595  if (status != nVar_per_cell) {
596  cerr << "[Hydroinfo_MUSIC::readHydroData]: ERROR: "
597  << "the evolution file format is not correct" << endl;
598  exit(1);
599  }
600 
601  if (itau_max < static_cast<int>(cell_info[0]))
602  itau_max = static_cast<int>(cell_info[0]);
603  fluidCell_3D_ideal newCell;
604  newCell.itau = static_cast<int>(cell_info[0]);
605  newCell.ix = static_cast<int>(cell_info[1]);
606  newCell.iy = static_cast<int>(cell_info[2]);
607  newCell.ieta = static_cast<int>(cell_info[3]);
608  newCell.temperature = cell_info[6];
609  newCell.ed = cell_info[4];
610  newCell.pressure = cell_info[5];
611  newCell.ux = cell_info[8];
612  newCell.uy = cell_info[9];
613  newCell.uz = cell_info[10];
614  lattice_3D_ideal.push_back(newCell);
615  ik++;
616  if (verbose_ > 7) {
617  if (ik%50000 == 0)
618  cout << "o" << flush;
619  }
620  }
621  cout << endl;
622  std::fclose(fin);
623  itaumax = itau_max;
624  // create the index map
625  long long ncells = (itaumax + 1)*ixmax*ixmax*ietamax;
626  idx_map_.resize(ncells, 0);
627  for (int i = 0; i < lattice_3D_ideal.size(); i++) {
628  const auto cell_i = lattice_3D_ideal[i];
629  long long cell_idx = (
630  ( (cell_i.itau*ietamax + cell_i.ieta)*ixmax
631  + cell_i.iy)*ixmax + cell_i.ix);
632  idx_map_[cell_idx] = i;
633  }
634  hydroTauMax = hydroTau0 + hydroDtau*itaumax;
635  } else {
636  cout << "Hydroinfo_MUSIC:: This option is obsolete! whichHydro = "
637  << whichHydro << endl;
638  exit(1);
639  }
640 
641  // One final step for easy automation of MARTINI:
642  // hydroTauMax is reset for the case where writing to evolution.dat
643  // ended early (due to all cells freezing out):
644  if (whichHydro == 6) {
645  hydroTauMax = (
646  hydroTau0 + hydroDtau*static_cast<int>(
647  static_cast<double>(lattice_3D.size())
648  /((2.*hydroXmax/hydroDx+1.)*(2.*hydroXmax/hydroDx+1.)
649  *2.*(hydro_eta_max/hydroDeta))));
650  itaumax = static_cast<int>((hydroTauMax-hydroTau0)/hydroDtau+0.001);
651  }
652  if (whichHydro == 8 || whichHydro == 9) {
653  hydroTauMax = (
654  hydroTau0 + hydroDtau*static_cast<int>(
655  static_cast<double>(lattice_2D.size())
656  /((2.*hydroXmax/hydroDx)*(2.*hydroXmax/hydroDx)) - 1));
657  itaumax = static_cast<int>((hydroTauMax - hydroTau0)/hydroDtau);
658  }
659 
660  cout << "hydro_tau0 = " << hydroTau0 << " fm"<< endl;
661  cout << "hydro_tau_max = " << hydroTauMax << " fm" << endl;
662  cout << "hydry_dtau = " << hydroDtau << " fm" << endl;
663  cout << "hydro_Xmax = " << hydroXmax << " fm" << endl;
664  cout << "hydro_dx = " << hydroDx << " fm" << endl;
665  cout << "hydro_eta_max = " << hydro_eta_max << endl;
666  cout << "hydro_deta = " << hydroDeta << endl;
667 }
668 
669 
670 void Hydroinfo_MUSIC::getHydroValues(double x, double y,
671  double z, double t, hydrofluidCell* info) {
672 // For interpolation of evolution files in tau-eta coordinates. Only the
673 // reading of MUSIC's evolution_xyeta.dat file is implemented here.
674 // For simplicity, hydro_eta_max refers to MUSIC's eta_size, and similarly for
675 // hydroDeta; however, x, y, z, and t are as usual to stay compatible with
676 // MARTINI.
677  double tau, eta;
678  if (use_tau_eta_coordinate == 1) {
679  if (t*t > z*z) {
680  tau = sqrt(t*t-z*z);
681  eta = 0.5*log((t+z)/(t-z));
682  } else {
683  tau = 0.;
684  eta = 0.;
685  }
686  } else {
687  // if the medium is given in cartesian coordinates
688  // set tau and eta to t and z
689  tau = t;
690  eta = z;
691  }
692 
693  int ieta = floor((hydro_eta_max+eta)/hydroDeta + 0.0001);
694  int itau = floor((tau-hydroTau0)/hydroDtau + 0.0001);
695  int ix = floor((hydroXmax+x)/hydroDx + 0.0001);
696  int iy = floor((hydroXmax+y)/hydroDx + 0.0001);
697 
698  double xfrac = (x - (static_cast<double>(ix)*hydroDx - hydroXmax))/hydroDx;
699  double yfrac = (y - (static_cast<double>(iy)*hydroDx - hydroXmax))/hydroDx;
700  double etafrac = (eta/hydroDeta - static_cast<double>(ieta)
701  + 0.5*static_cast<double>(ietamax));
702  double taufrac = (tau - hydroTau0)/hydroDtau - static_cast<double>(itau);
703 
704  if (boost_invariant) {
705  ieta = 0;
706  etafrac = 0.;
707  }
708 
709  if (ix < 0 || ix >= ixmax) {
710  if (verbose_ > 7) {
711  cout << "[Hydroinfo_MUSIC::getHydroValues]: "
712  << "WARNING - x out of range x=" << x
713  << ", ix=" << ix << ", ixmax=" << ixmax << endl;
714  cout << "x=" << x << " y=" << y << " eta=" << eta
715  << " ix=" << ix << " iy=" << iy << " ieta=" << ieta << endl;
716  cout << "t=" << t << " tau=" << tau
717  << " itau=" << itau << " itaumax=" << itaumax << endl;
718  }
719 
720  info->temperature = 0.0;
721  info->vx = 0.0;
722  info->vy = 0.0;
723  info->vz = 0.0;
724  return;
725  }
726  if (iy < 0 || iy >= ixmax) {
727  if (verbose_ > 7) {
728  cout << "[Hydroinfo_MUSIC::getHydroValues]: "
729  << "WARNING - y out of range, y=" << y << ", iy=" << iy
730  << ", iymax=" << ixmax << endl;
731  cout << "x=" << x << " y=" << y << " eta=" << eta
732  << " ix=" << ix << " iy=" << iy << " ieta=" << ieta << endl;
733  cout << "t=" << t << " tau=" << tau
734  << " itau=" << itau << " itaumax=" << itaumax << endl;
735  }
736 
737  info->temperature = 0.0;
738  info->vx = 0.0;
739  info->vy = 0.0;
740  info->vz = 0.0;
741  return;
742  }
743  if (itau < 0 || itau > itaumax) {
744  if (verbose_ > 7) {
745  cout << "[Hydroinfo_MUSIC::getHydroValues]: WARNING - "
746  << "tau out of range, itau=" << itau
747  << ", itaumax=" << itaumax << endl;
748  cout << "[MARTINI:Hydroinfo_MUSIC::getHydroValues]: tau= " << tau
749  << ", hydroTauMax = " << hydroTauMax
750  << ", hydroDtau = " << hydroDtau << endl;
751  }
752 
753  info->temperature = 0.0;
754  info->vx = 0.0;
755  info->vy = 0.0;
756  info->vz = 0.0;
757  return;
758  }
759  if (ieta < 0 || ieta >= ietamax) {
760  if (verbose_ > 7) {
761  cout << "[Hydroinfo_MUSIC::getHydroValues]: WARNING - "
762  << "eta out of range, ieta=" << ieta
763  << ", ietamax=" << ietamax
764  << endl;
765  }
766  info->temperature = 0.0;
767  info->vx = 0.0;
768  info->vy = 0.0;
769  info->vz = 0.0;
770  return;
771  }
772 
773  // The array of positions on the 4-dimensional rectangle:
774  int position[2][2][2][2];
775  for (int ipx = 0; ipx < 2; ipx++) {
776  int px;
777  if (ipx == 0 || ix == ixmax-1) {
778  px = ix;
779  } else {
780  px = ix + 1;
781  }
782  for (int ipy = 0; ipy < 2; ipy++) {
783  int py;
784  if (ipy == 0 || iy == ixmax-1) {
785  py = iy;
786  } else {
787  py = iy + 1;
788  }
789  for (int ipeta = 0; ipeta < 2; ipeta++) {
790  int peta;
791  if (ipeta == 0 || ieta == ietamax-1) {
792  peta = ieta;
793  } else {
794  peta = ieta + 1;
795  }
796  for (int iptau = 0; iptau < 2; iptau++) {
797  int ptau;
798  if (iptau == 0 || itau == itaumax) {
799  ptau = itau;
800  } else {
801  ptau = itau + 1;
802  }
803  position[ipx][ipy][ipeta][iptau] = (
804  px + ixmax*(py + ixmax*(peta + ietamax*ptau)));
805  }
806  }
807  }
808  }
809 
810  // And now, the interpolation:
811  double T = 0.0;
812  double ed = 0.;
813  double p = 0.;
814  double vx = 0.0;
815  double vy = 0.0;
816  double vz = 0.0;
817  double ux = 0.0;
818  double uy = 0.0;
819  double uz = 0.0;
820  double ueta = 0.0;
821  double pi00 = 0.0;
822  double pi01 = 0.0;
823  double pi02 = 0.0;
824  double pi03 = 0.0;
825  double pi11 = 0.0;
826  double pi12 = 0.0;
827  double pi13 = 0.0;
828  double pi22 = 0.0;
829  double pi23 = 0.0;
830  double pi33 = 0.0;
831  double bulkPi = 0.0;
832 
833  fluidCell_2D *HydroCell_2D_ptr1, *HydroCell_2D_ptr2;
834  fluidCell_3D *HydroCell_3D_ptr1, *HydroCell_3D_ptr2;
835  fluidCell_3D_ideal *HydroCell_3D_ideal_ptr1, *HydroCell_3D_ideal_ptr2;
836  for (int iptau = 0; iptau < 2; iptau++) {
837  double taufactor;
838  if (iptau == 0)
839  taufactor = 1. - taufrac;
840  else
841  taufactor = taufrac;
842  for (int ipeta = 0; ipeta < 2; ipeta++) {
843  double etafactor;
844  if (ipeta == 0)
845  etafactor = 1. - etafrac;
846  else
847  etafactor = etafrac;
848  for (int ipy = 0; ipy < 2; ipy++) {
849  double yfactor;
850  if (ipy == 0)
851  yfactor = 1. - yfrac;
852  else
853  yfactor = yfrac;
854 
855  double prefrac = yfactor*etafactor*taufactor;
856 
857  if (hydroWhichHydro == 8 || hydroWhichHydro == 9) {
858  HydroCell_2D_ptr1 = (
859  &lattice_2D[position[0][ipy][ipeta][iptau]]);
860  HydroCell_2D_ptr2 = (
861  &lattice_2D[position[1][ipy][ipeta][iptau]]);
862  T += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->temperature
863  + xfrac*HydroCell_2D_ptr2->temperature);
864  ux += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->ux
865  + xfrac*HydroCell_2D_ptr2->ux);
866  uy += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->uy
867  + xfrac*HydroCell_2D_ptr2->uy);
868  ueta += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->ueta
869  + xfrac*HydroCell_2D_ptr2->ueta);
870  pi00 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi00
871  + xfrac*HydroCell_2D_ptr2->pi00);
872  pi01 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi01
873  + xfrac*HydroCell_2D_ptr2->pi01);
874  pi02 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi02
875  + xfrac*HydroCell_2D_ptr2->pi02);
876  pi11 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi11
877  + xfrac*HydroCell_2D_ptr2->pi11);
878  pi12 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi12
879  + xfrac*HydroCell_2D_ptr2->pi12);
880  pi22 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi22
881  + xfrac*HydroCell_2D_ptr2->pi22);
882  pi33 += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->pi33
883  + xfrac*HydroCell_2D_ptr2->pi33);
884  bulkPi += prefrac*((1. - xfrac)*HydroCell_2D_ptr1->bulkPi
885  + xfrac*HydroCell_2D_ptr2->bulkPi);
886  } else if (hydroWhichHydro == 6) {
887  HydroCell_3D_ptr1 = (
888  &lattice_3D[position[0][ipy][ipeta][iptau]]);
889  HydroCell_3D_ptr2 = (
890  &lattice_3D[position[1][ipy][ipeta][iptau]]);
891  T += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->temperature
892  + xfrac*HydroCell_3D_ptr2->temperature);
893  vx += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->vx
894  + xfrac*HydroCell_3D_ptr2->vx);
895  vy += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->vy
896  + xfrac*HydroCell_3D_ptr2->vy);
897  vz += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->vz
898  + xfrac*HydroCell_3D_ptr2->vz);
899  pi00 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi00
900  + xfrac*HydroCell_3D_ptr2->pi00);
901  pi01 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi01
902  + xfrac*HydroCell_3D_ptr2->pi01);
903  pi02 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi02
904  + xfrac*HydroCell_3D_ptr2->pi02);
905  pi03 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi03
906  + xfrac*HydroCell_3D_ptr2->pi03);
907  pi11 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi11
908  + xfrac*HydroCell_3D_ptr2->pi11);
909  pi12 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi12
910  + xfrac*HydroCell_3D_ptr2->pi12);
911  pi13 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi13
912  + xfrac*HydroCell_3D_ptr2->pi13);
913  pi22 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi22
914  + xfrac*HydroCell_3D_ptr2->pi22);
915  pi23 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi23
916  + xfrac*HydroCell_3D_ptr2->pi23);
917  pi33 += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->pi33
918  + xfrac*HydroCell_3D_ptr2->pi33);
919  bulkPi += prefrac*((1. - xfrac)*HydroCell_3D_ptr1->bulkPi
920  + xfrac*HydroCell_3D_ptr2->bulkPi);
921  } else if (hydroWhichHydro == 10 || hydroWhichHydro == 11) {
922  HydroCell_3D_ideal_ptr1 = (
923  &lattice_3D_ideal[idx_map_[position[0][ipy][ipeta][iptau]]]);
924  HydroCell_3D_ideal_ptr2 = (
925  &lattice_3D_ideal[idx_map_[position[1][ipy][ipeta][iptau]]]);
926  T += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->temperature
927  + xfrac*HydroCell_3D_ideal_ptr2->temperature);
928  ed += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->ed
929  + xfrac*HydroCell_3D_ideal_ptr2->ed);
930  p += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->pressure
931  + xfrac*HydroCell_3D_ideal_ptr2->pressure);
932  ux += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->ux
933  + xfrac*HydroCell_3D_ideal_ptr2->ux);
934  uy += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->uy
935  + xfrac*HydroCell_3D_ideal_ptr2->uy);
936  uz += prefrac*((1. - xfrac)*HydroCell_3D_ideal_ptr1->uz
937  + xfrac*HydroCell_3D_ideal_ptr2->uz);
938  }
939  }
940  }
941  }
942 
943  if (hydroWhichHydro == 10) {
944  double ut = sqrt(1. + ux*ux + uy*uy + uz*uz);
945  vx = ux/ut;
946  vy = uy/ut;
947  vz = uz/ut;
948  } else if (hydroWhichHydro == 11) {
949  double eta_local = 0.5*log((t + z)/(t - z));
950  double sinh_eta = sinh(eta_local);
951  double cosh_eta = cosh(eta_local);
952  ueta = uz;
953  double utau = sqrt(1. + ux*ux + uy*uy + ueta*ueta);
954  uz = utau*sinh_eta + ueta*cosh_eta;
955  double ut = utau*cosh_eta + ueta*sinh_eta;
956  vx = ux/ut;
957  vy = uy/ut;
958  vz = uz/ut;
959  } else if (hydroWhichHydro == 8 || hydroWhichHydro == 9) {
960  double eta_local = 0.5*log((t + z)/(t - z));
961  double sinh_eta = sinh(eta_local);
962  double cosh_eta = cosh(eta_local);
963  double utau = sqrt(1. + ux*ux + uy*uy + ueta*ueta);
964  uz = utau*sinh_eta + ueta*cosh_eta;
965  double ut = utau*cosh_eta + ueta*sinh_eta;
966  vx = ux/ut;
967  vy = uy/ut;
968  vz = uz/ut;
969  }
970 
971  info->temperature = T;
972  info->vx = vx;
973  info->vy = vy;
974  info->vz = vz;
975 
976  if (hydroWhichHydro < 10) {
977  info->ed = 1.0; // pi's are already divided by e+P
978  info->sd = 0.0;
979  info->pressure = 0.0;
980  } else {
981  info->ed = ed;
982  info->sd = (ed + p)/(T + 1e-16);
983  info->pressure = p;
984  }
985 
986  info->pi[0][0] = pi00;
987  info->pi[0][1] = pi01;
988  info->pi[0][2] = pi02;
989  info->pi[0][3] = pi03;
990  info->pi[1][0] = pi01;
991  info->pi[1][1] = pi11;
992  info->pi[1][2] = pi12;
993  info->pi[1][3] = pi13;
994  info->pi[2][0] = pi02;
995  info->pi[2][1] = pi12;
996  info->pi[2][2] = pi22;
997  info->pi[2][3] = pi23;
998  info->pi[3][0] = pi03;
999  info->pi[3][1] = pi13;
1000  info->pi[3][2] = pi23;
1001  info->pi[3][3] = pi33;
1002 
1003  info->bulkPi = bulkPi;
1004  return;
1005 }
1006 
1007 
1008 void Hydroinfo_MUSIC::output_temperature_evolution(string filename_base) {
1009  hydrofluidCell *hydroInfo = new hydrofluidCell;
1010  for (int i = 0; i < itaumax; i++) {
1011  double tau = hydroTau0 + i*hydroDtau;
1012  ostringstream filename;
1013  filename << filename_base << "_tau_" << tau << ".dat";
1014  ofstream temp_evo(filename.str().c_str());
1015  for (int ix = 0; ix < ixmax; ix++) {
1016  double x_local = -hydroXmax + ix*hydroDx;
1017  for (int iy = 0; iy < ixmax; iy++) {
1018  double y_local = -hydroXmax + iy*hydroDx;
1019  getHydroValues(x_local, y_local, 0.0, tau, hydroInfo);
1020  double temp_local = hydroInfo->temperature;
1021  temp_evo << scientific << setw(16) << setprecision(8)
1022  << temp_local << " ";
1023  }
1024  temp_evo << endl;
1025  }
1026  temp_evo.close();
1027  }
1028  delete hydroInfo;
1029 }
1030 
1031 
1032 void Hydroinfo_MUSIC::update_grid_info(
1033  double tau0, double tau_max, double dtau,
1034  double x_max, double dx, double eta_max, double deta) {
1035  hydroTau0 = tau0;
1036  hydroTauMax = tau_max;
1037  hydroDtau = dtau;
1038  hydroXmax = x_max;
1039  hydroDx = dx;
1040  hydro_eta_max = eta_max;
1041  hydroDeta = deta;
1042  if (hydroWhichHydro == 8) {
1043  itaumax = static_cast<int>((tau_max-tau0)/dtau+0.001);
1044  ixmax = static_cast<int>(2*x_max/dx+0.001);
1045  ietamax = static_cast<int>(2*eta_max/deta+0.001);
1046  }
1047  if (hydroWhichHydro == 6) {
1048  itaumax = static_cast<int>((tau_max-tau0)/dtau+0.001);
1049  ixmax = static_cast<int>(2*x_max/dx+0.001);
1050  ietamax = static_cast<int>(2*eta_max/deta+0.001);
1051  }
1052 }