from astropy.io import ascii
import numpy as np
import matplotlib.pyplot as plt
ydelta_data=ascii.read("all_phase_space_at_12278_ns.dat",header_start=0)
#ydelta_data=ascii.read("all_phase_space_at_17868_ns.dat",header_start=0)
#plt.scatter('y','pz',data=ydelta_data, s=1)
#plt.xlabel('x')
#plt.ylabel('$\delta$')
H, xedges, yedges = np.histogram2d(ydelta_data['y'], ydelta_data['pz'], bins=(50,50))
H=H.T
plt.figure(1)
plt.imshow(H, interpolation='nearest', origin='lower',aspect='auto',extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])
H, xedges, yedges = np.histogram2d(ydelta_data['y'], ydelta_data['pz'], bins=(50,50))
H.ndim
H2=np.multiply(H[:,:],yedges[0:50]).sum(axis=0)
H2.ndim
H2.size
H, xedges, yedges = np.histogram2d(ydelta_data['pz'], ydelta_data['y'], bins=(50,50))
H2=np.multiply(H,xedges[0:50]).sum(axis=1)/H.sum(axis=1)
Ntot =H.sum(axis=1)
Ytot=np.multiply(H,xedges[0:50]).sum(axis=1)
Y2tot=np.multiply(H,xedges[0:50]**2).sum(axis=1)
err =(Y2tot*Ntot**2-Ytot**2*Ntot)**.5/Ntot**2
plt.figure(2)
plt.errorbar(yedges[0:50],H2,err,fmt='.')
plt.xlabel('y [m]')
plt.ylabel('$\delta$')
plt.title('config_72/12278_ns')
#plt.title('config_72/17868_ns')
coef=np.polyfit(yedges[0:50],H2,3)
print(coef)
def f(yedges): return coef[0]*yedges**3+coef[1]*yedges**2+coef[2]*yedges+coef[3]

plt.errorbar(yedges,f(yedges))
plt.show()