import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
# PLOTTING OPTIONS
fig_width_pt = 3*246.0  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0         # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt  # width in inches
fig_height = fig_width*golden_mean      # height in inches
fig_size =  [fig_width,fig_height]

params = { 'axes.labelsize': 24,
          'font.family': 'serif',
          'font.serif': 'Computer Modern Raman',
          'font.size': 24,
          'legend.fontsize': 20,
          'xtick.labelsize': 24,
          'ytick.labelsize': 24,
          'axes.grid' : True,
          'text.usetex': True,
          'savefig.dpi' : 100,
          'lines.markersize' : 14,
          'figure.figsize': fig_size}

mpl.rcParams.update(params)

# The index of A221 is 3
i = 3

r1 = np.loadtxt('/work/yifan.wang/ringdown/GW150914/pyring/reproduce/t6/221/GW150914_PROD1_Kerr_221_0M/Nested_sampler/posterior.dat')

r2=np.loadtxt('/work/yifan.wang/ringdown/GW150914/pyring/re-pyring-t0p2/t6/221/GW150914_PROD1_Kerr_221_0M/Nested_sampler/posterior.dat')
r8k = np.loadtxt('/work/yifan.wang/more_ringdown/runs/pyring/t1-1024hz/221-8k/GW150914_PROD1_Kerr_221_0M/Nested_sampler/posterior.dat')

r4k=np.loadtxt('/work/yifan.wang/more_ringdown/runs/pyring/t1-1024hz/221/GW150914_PROD1_Kerr_221_0M/Nested_sampler/posterior.dat')

bins = np.linspace(0,50,100)
plt.hist(r1[:,3],bins=bins,density=True,label='0.1s,nlive=2048,srate4096,band[20,2038]')
plt.hist(r2[:,3],bins=bins,density=True,alpha=0.5,label='0.2s,nlive=2048,srate4096,band[20,2038]')
plt.hist(r4k[:,3],bins=bins,density=True,alpha=0.5,label='0.2s,nlive=4000,srate2048,band[20,1014]')
plt.axhline(0.02,ls='--',color='black',label='y=0.02=1/prior')
plt.xlabel('$A_{221}/10^{-21}$')
plt.legend()
plt.title('PyRing Runs')
plt.savefig('/work/yifan.wang/public_html/GWevents/figure/pyringrun.png',bbox_inches='tight')

np.sqrt(16*np.pi/5)

3.1706618380848086

np.median(r4k[:,3]) / np.std(r4k[:,3])

3.3534632494610697