trento/scripts/global-time-converter.py

0001 #!/usr/bin/env python3
0002 import numpy as np, h5py, math
0003 import sys, os, warnings
0004 import matplotlib.pyplot as plt
0005
0006 def help():
0007     """
0008   This script transforms the initial condition from
0009   (tau, x, y, eta) frame to (t, x, y, z) frame. This
0010   conversion assumes longitudinal Bjorken expansion
0011   and transverse still.
0012
0013   The entropy density in longitudinal-local frame at
0014   constant proper time tau0 is:
0015     s(tau0, x, y, eta) = s0(x,y,eta)/tau0
0016                        = ns0(x,y,eta) <-- trento3d
0017
0018   Assuming Bjorken expansion for a short proper time,
0019     s(tau, x, y, eta) = 1/tau*s0(x,y,eta)
0020                       = tau0/tau*ns0(x,y,eta) -- (1)
0021
0022   Expression (1) is transformed into (t, x, y, z) via,
0023     s'(t0, x, y, z) = 1/tau*ns0(x, y, eta)
0024                 tau = sqrt(t0^2 - z^2)
0025                 eta = 1/2*ln[(t0+z)/(t0-z)]
0026   where tau0 is absorbed into normalization of ns0
0027   (see --normalization option of the trento3d model).
0028
0029   Usage:
0030     {:s} trento3d-hdf5-output t0 [list-of-event-id-to-convert]
0031
0032   For example, to convert all events to t0=1fm/c:
0033     {:s} ic.hdf5 1.0
0034   To convert only events #2 and #3 to t0=1fm/c:
0035     {:s} ic.hdf5 1.0 2 3
0036 """
0037     print(help.__doc__.format(__file__, __file__, __file__))
0038
0039 tiny = 1e-9
0040
0041 def convert(dataset, t0):
0042     field = dataset.value
0043     Nx, Ny, Neta = field.shape
0044     deta = dataset.attrs['deta']
0045     dxy = dataset.attrs['dxy']
0046     Lx, Ly, Leta = Nx*dxy/2., Ny*dxy/2., (Neta-1)*deta/2.
0047     x = np.linspace(-Lx, Lx, Nx)
0048     y = np.linspace(-Ly, Ly, Ny)
0049     eta = np.linspace(-Leta, Leta, Neta)
0050     # set maximum z-range, about +/- 6-unit of rapidity
0051     z_max = t0*math.tanh(6.)
0052     zarray = np.linspace(-z_max, z_max, Neta)
0053
0054     new_field = np.zeros_like(field)
0055
0056     for iz, z in enumerate(zarray):
0057         eta = 0.5*np.log((t0+z)/(t0-z))
0058         if np.abs(eta) > Leta:
0059             continue
0060         residue, index = math.modf((eta+Leta)/deta)
0061         index = int(index)
0062         new_field[:,:,iz] = (  field[:,:,index]  *(1.-residue) \
0063                              + field[:,:,index+1]*residue)  \
0064                              / np.sqrt(t0**2-z**2)
0065     return x, y, zarray, new_field
0066
0067 def plot(x, y, z, s):
0068     dx = x[1]-x[0]
0069     dy = y[1]-y[0]
0070     fig, ax = plt.subplots(1, 1)
0071     ax.plot(z, np.sum(s, axis=(0,1))*dx*dy, 'ro-')
0072     ax.set_xlabel(r'$z$ [fm]', size=15)
0073     ax.set_title(r'Local-rest-frame entropy in $t-z$ coordinates', size=15)
0074     ax.set_ylabel(r'$ds/dV$ [Arb. units.]', size=15)
0075     plt.subplots_adjust(wspace=0.15, bottom=0.2)
0076     plt.semilogy()
0077     plt.show()
0078
0079 def main():
0080     if len(sys.argv) <= 2:
0081         help()
0082         exit()
0083     f = h5py.File(sys.argv[1], 'r')
0084     t0 = float(sys.argv[2])
0085     elist = ['event_{}'.format(index) for index in sys.argv[3:]] \
0086             if len(sys.argv)>=4 else list(f.keys())
0087     for eid in elist:
0088         x, y, z, s = convert(f[eid], t0)
0089         plot(x, y, z, s)
0090 if __name__ == "__main__":
0091     main()