Skip to content
Snippets Groups Projects
Commit 06edd777 authored by Carl Gwilliam's avatar Carl Gwilliam
Browse files

Update to fluka reader

parent 7f457cf7
No related branches found
No related tags found
No related merge requests found
......@@ -82,25 +82,6 @@ class FlukaReader(EvgenAlg):
"Convert cos(theta) wrt x or y axis to theta wrt z axis"
return np.pi/2 - np.arccos(cosTheta)
# def theta_and_phi(self, thetaX, thetaY, twopi = True):
# "Converted thetaX and thetaY wrt Z to theta and phi"
#
# theta = thetaY
#
#
# phi = np.arctan(np.abs(thetaY)/np.abs(thetaX))
#
# if thetaX < 0 and thetaY > 0:
# phi = phi + np.pi/2.
# elif thetaX > 0 and thetaY < 0:
# phi = -phi
# elif thetaX < 0 and thetaY < 0:
# phi = -phi - np.pi/2.
#
# if twopi and phi < 0:
# phi = 2*np.pi + phi
return theta, phi
def pid(self, ftype):
"Convert fluka particle type to PID"
......@@ -291,9 +272,14 @@ class FlukaReader(EvgenAlg):
thetaX = self.angle(entry["cosX"])
thetaY = self.angle(entry["cosY"])
#theta, phi = self.theta_and_phi(thetaX, thetaY)
theta = thetaY
phi = np.arctan(entry["cosY"]/entry["cosX"])
# theta: just above z axis as phi deals with negative
theta = np.abs(thetaY)
# phi: 0 - 2pi
phi = np.arctan2(entry["cosY"], entry["cosX"])
#phi = np.arctan(entry["cosY"] / entry["cosX"])
if phi < 0: phi += 2*np.pi
if phi == 2*np.pi: phi = 0
px = p * np.sin(theta) * np.cos(phi)
py = p * np.sin(theta) * np.sin(phi)
......@@ -326,7 +312,7 @@ class FlukaReader(EvgenAlg):
plt.hist(self.after["E"], bins = ebins, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("energy.eps")
plt.savefig("energy.ong")
plt.figure()
plt.xlabel("Angle to beam in X dir")
......@@ -337,7 +323,7 @@ class FlukaReader(EvgenAlg):
plt.hist(thetaXout, bins = tbins, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("thetaX.eps")
plt.savefig("thetaX.png")
plt.figure()
plt.xlabel("Angle to beam in Y dir")
......@@ -347,7 +333,7 @@ class FlukaReader(EvgenAlg):
plt.hist(thetaYout, bins = tbins, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("thetaY.eps")
plt.savefig("thetaY.png")
plt.figure()
plt.xlabel("Dispacement in X dir")
......@@ -356,7 +342,7 @@ class FlukaReader(EvgenAlg):
plt.hist(self.after["x"], bins = xbins, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("x.eps")
plt.savefig("x.png")
plt.figure()
plt.xlabel("Dispacement in Y dir")
......@@ -364,7 +350,7 @@ class FlukaReader(EvgenAlg):
plt.hist(self.after["y"], bins = xbins, histtype='step', color = "r", fill = False, label = "after")
plt.gca().set_yscale('log')
plt.legend()
plt.savefig("y.eps")
plt.savefig("y.png")
return
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment