from astropy.io import ascii
import numpy as np
import matplotlib.pyplot as plt
data=ascii.read("muon_all_end.dat",header_start=1, fast_reader={'exponent_style': 'D'})
data=ascii.read("muon_all_end.dat")
print(data)
data.columns
data.index_column('s_start')
from matplotlib import pyplot as plt 
plt.scatter(data['x'], data['y'])
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('muon.pdf', format='pdf')
print(data)
data.columns
plt.scatter(data['x'], data['px=Px/p0c'])
plt.scatter(data['y'], data['pz=dP/p0c=(P-P0)/p0c'])
plt.scatter(data['y'], data['pz=dP/p0c=(P-P0)/p0c'],s=1)     #reduce point size
plt.scatter(data['x'], data['px=Px/p0c'],s=1)
plt.scatter(data['x'], data['px=Px/p0c']*2,s=1)     #do something to second column
counts, bins = np.histogram(data['x'])
plt.stairs(counts, bins)
counts, bins = np.histogram(data['x'],bins=20)
plt.stairs(counts, bins)
data['pz']=((1+data['pz=dP/p0c=(P-P0)/p0c'])**2-(data['px=Px/p0c'])**2-(data['px=Py/p0c'])**2)**0.5  #add s olumn with the true pz
data['sdotp']=(data['px=Px/p0c']*data['sx']+data['px=Py/p0c']*data['sy']+data['pz']*data['sz'])/(1+data['pz=dP/p0c=(P-P0)/p0c'])    #add column with s.p
data['phi']=np.arccos(data['sdotp'])
H, xedges, yedges = np.histogram2d(data['x'], data['phi'], bins=(100,100))
H=H.T    #some kind of translation
plt.imshow(H, interpolation='nearest', origin='lower',aspect='auto',
        extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])