import numpy as np
import matplotlib.pyplot as plt

# load data form SOHO/EPHIN 17.05.2012 (solar event)

#detA,Etot=np.loadtxt("epi12138.pl2",usecols=(5,10),unpack=True)
detA,Etot=np.loadtxt("artificial_data.pl2",usecols=(5,10),unpack=True)

# lets do a collection of scatter plots
plt.figure()
myalpha=0.01

# first plot will be a typical dE/dx-E plot
plt.subplot(2,2,1)
plt.scatter(Etot,detA,marker='o',color='k',rasterized=True,alpha=myalpha)
plt.xlim(0,100)
plt.ylim(0,20)
plt.xlabel("total kinetic energy / MeV")
plt.ylabel("energy loss in detector A / MeV")

# second plot will be the same as the first one, but with log scales
plt.subplot(2,2,2)
plt.scatter(Etot,detA,marker='o',color='k',rasterized=True,alpha=myalpha)
plt.xlim(0.1,100)
plt.ylim(0.1,20)
plt.xscale("log")
plt.yscale("log")
plt.xlabel("total kinetic energy / MeV")
plt.ylabel("energy loss in detector A / MeV")

# third plot will be a rotation of the first ones
plt.subplot(2,2,3)
plt.scatter(Etot/detA,Etot*detA,marker='o',color='k',rasterized=True,alpha=myalpha)
plt.xlim(9e-1,3e2)
plt.ylim(1e-3,1e3)
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r"E$_{tot}$ / E$_A$")
plt.ylabel(r"E$_{tot}$ * E$_A$ / MeV$^2$")

# lets now make a histogram over E*dE
plt.subplot(2,2,4)
mybins=np.logspace(-3,3,1000)
hist,binedges=np.histogram(Etot*detA,bins=mybins)

bincenter=[10**((np.log10(binedges[q])+np.log10(binedges[q+1]))/2.) for q in range(len(binedges)-1)]
plt.step(bincenter,hist,c='k')
plt.xscale("log")
plt.yscale("log")
plt.ylabel("number of entries")
plt.xlabel(r"E$_{tot}$ * E$_A$ / MeV$^2$")

#plt.show()
plt.savefig("artificial_data.pdf")