Skip to content
Snippets Groups Projects
Commit 4a7f903b authored by lorenzennio's avatar lorenzennio
Browse files

final corrections, everything works now

parent 08c8ccb8
No related branches found
No related tags found
No related merge requests found
Pipeline #4844 passed
...@@ -73,7 +73,7 @@ class Layer: ...@@ -73,7 +73,7 @@ class Layer:
class Detector: class Detector:
def __init__ (self, LayerArray = [100,110,120,130,140], Magnet=0.522): def __init__ (self, LayerArray = [100,110,120,130,140], Magnet=0.511):
"""Initialize the Layer and magnetic fiel of the detector""" """Initialize the Layer and magnetic fiel of the detector"""
self.BField = Magnet #Microesla self.BField = Magnet #Microesla
......
...@@ -14,17 +14,23 @@ class Particle: ...@@ -14,17 +14,23 @@ class Particle:
"""If we want two different masses (electron, proton), I have here """If we want two different masses (electron, proton), I have here
created them with equal probability""" created them with equal probability"""
self.e = np.random.uniform(1, 0) self.e = np.random.uniform(0,1)
if self.e > 0.5: if self.e >= 0.5:
self.m = 0.511 #MeV self.m = 0.511 #MeV
self.q = -1
else: else:
self.m = 938.272 #Mev self.m = 0.511 #Mev
self.q = 1
"""Velocity: random velocity between X and Y""" """Velocity: random velocity between X and Y"""
self.v = np.random.uniform(0.2, 0.6) #Natural units self.v = np.random.uniform(0.2, 0.6) #Natural units
#Alter the range of velocities #Alter the range of velocities
#Perhaps make vx, vy, vz #Perhaps make vx, vy, vz
self.vx = np.random.uniform(0, 80)
self.vy = np.random.uniform(0, 80)
self.vz = np.random.uniform(200, 500)
"""Charge, 1,-1""" """Charge, 1,-1"""
self.k = np.random.uniform(1, 0) self.k = np.random.uniform(1, 0)
......
...@@ -3,7 +3,7 @@ Script to run the full detector simulation. ...@@ -3,7 +3,7 @@ Script to run the full detector simulation.
""" """
#TODO: Import all classes #TODO: Import all classes
from Particles_software import Particles from Particles_software import Particle
from Detector import Layer, Detector from Detector import Layer, Detector
from analysis import Analysis from analysis import Analysis
......
import numpy as np import numpy as np
import pandas as pd import pandas as pd
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from scipy.optimize import curve_fit from mpl_toolkits.mplot3d import Axes3D
from scipy import odr from scipy import odr
...@@ -38,7 +38,8 @@ class Analysis: ...@@ -38,7 +38,8 @@ class Analysis:
self.results = pd.DataFrame({'w' : np.zeros(2), self.results = pd.DataFrame({'w' : np.zeros(2),
'vx' : np.zeros(2), 'vx' : np.zeros(2),
'vy' : np.zeros(2), 'vy' : np.zeros(2),
'vz' : np.zeros(2)}) 'vz' : np.zeros(2),
'm' : np.zeros(2)})
#Determine mean and std of input #Determine mean and std of input
...@@ -100,24 +101,22 @@ class Analysis: ...@@ -100,24 +101,22 @@ class Analysis:
#Extract correct results from fitting and fill to results container #Extract correct results from fitting and fill to results container
self.results['vx'] = np.array([xResults.beta[0], xResults.sd_beta[0]]) self.results['vx'] = np.array([xResults.beta[0], xResults.sd_beta[0]])
#w = (np.abs(yResults.beta[0]) + np.abs(zResults.beta[0]))/2.
#wStd = np.sqrt(yResults.sd_beta[0]**2 + zResults.sd_beta[0]**2)
w = wSgn*np.abs(zResults.beta[0]) w = wSgn*np.abs(zResults.beta[0])
wStd = zResults.sd_beta[0] wStd = zResults.sd_beta[0]
self.results['w'] = np.array([w , wStd]) self.results['w'] = np.array([w , wStd])
#vy = ySgn*(np.abs(yResults.beta[1]) + np.abs(zResults.beta[1]))/2.
#vyStd = np.sqrt(yResults.sd_beta[1]**2 + zResults.sd_beta[1]**2)
vy = vySgn*np.abs(zResults.beta[1]) vy = vySgn*np.abs(zResults.beta[1])
vyStd = zResults.sd_beta[1] vyStd = zResults.sd_beta[1]
self.results['vy'] = np.array([vy,vyStd]) self.results['vy'] = np.array([vy,vyStd])
#vz = zSgn*(np.abs(yResults.beta[2]) + np.abs(zResults.beta[2]))/2.
#vzStd = np.sqrt(yResults.sd_beta[2]**2 + zResults.sd_beta[2]**2)
vz = np.abs(zResults.beta[2]) vz = np.abs(zResults.beta[2])
vzStd = zResults.sd_beta[2] vzStd = zResults.sd_beta[2]
self.results['vz'] = np.array([vz,vzStd]) self.results['vz'] = np.array([vz,vzStd])
m = 0.511/np.abs(w)
mStd = m*wStd/np.abs(w)
self.results['m'] = np.array([m, mStd])
""" """
Equations of motion: Equations of motion:
...@@ -175,15 +174,6 @@ class Analysis: ...@@ -175,15 +174,6 @@ class Analysis:
print('The x fit was performend with chi2 = ', xOUT.res_var) print('The x fit was performend with chi2 = ', xOUT.res_var)
#print(xOUT.sd_beta) #print(xOUT.sd_beta)
return xOUT return xOUT
"""
print('x ', xOUT.res_var)
plt.plot(self.tPoints['Mean'],self.xEOM(xOUT.beta, self.tPoints['Mean']))
plt.plot(self.tPoints['Mean'], self.xPoints['Mean'], 'o')
plt.xlabel('t')
plt.ylabel('x')
plt.savefig('x.pdf', format='pdf')
plt.gcf().clear()
"""
def yFit(self): def yFit(self):
yModel = odr.Model(self.yEOM) yModel = odr.Model(self.yEOM)
...@@ -196,15 +186,6 @@ class Analysis: ...@@ -196,15 +186,6 @@ class Analysis:
print('The y fit was performend with chi2 = ', yOUT.res_var) print('The y fit was performend with chi2 = ', yOUT.res_var)
#yOUT.pprint() #yOUT.pprint()
return yOUT return yOUT
"""
print('y :', yOUT.res_var)
plt.plot(self.tPoints['Mean'], self.yEOM(yOUT.beta, self.tPoints['Mean']))
plt.plot(self.tPoints['Mean'], self.yPoints['Mean'], 'o')
plt.xlabel('t')
plt.ylabel('y')
plt.savefig('y.pdf', format='pdf')
plt.gcf().clear()
"""
def zFit(self): def zFit(self):
zModel = odr.Model(self.zEOM) zModel = odr.Model(self.zEOM)
...@@ -218,14 +199,7 @@ class Analysis: ...@@ -218,14 +199,7 @@ class Analysis:
#zOUT.pprint() #zOUT.pprint()
#print('z :', zOUT.res_var) #print('z :', zOUT.res_var)
return zOUT return zOUT
"""
plt.plot(self.tPoints['Mean'], self.zEOM(zOUT.beta, self.tPoints['Mean']))
plt.plot(self.tPoints['Mean'], self.zPoints['Mean'], 'o')
plt.xlabel('t')
plt.ylabel('z')
plt.savefig('z.pdf', format='pdf')
plt.gcf().clear()
"""
def output(self): def output(self):
""" """
A function to print results to terminal. A function to print results to terminal.
...@@ -251,7 +225,7 @@ class Analysis: ...@@ -251,7 +225,7 @@ class Analysis:
x_eom = np.array(self.xEOM([self.results["vx"][0]], times)) x_eom = np.array(self.xEOM([self.results["vx"][0]], times))
y_eom = np.array(self.yEOM([self.results["w"][0], self.results["vy"][0], self.results["vz"][0]], times), dtype=pd.Series) y_eom = np.array(self.yEOM([self.results["w"][0], self.results["vy"][0], self.results["vz"][0]], times), dtype=pd.Series)
z_eom = np.array(self.zEOM([self.results["w"][0], self.results["vy"][0], self.results["vz"][0]], times), dtype=pd.Series) z_eom = np.array(self.zEOM([self.results["w"][0], self.results["vy"][0], self.results["vz"][0]], times), dtype=pd.Series)
spaces = np.zeros(7) spaces = np.zeros(5)
fig = plt.figure() fig = plt.figure()
...@@ -295,4 +269,4 @@ class Analysis: ...@@ -295,4 +269,4 @@ class Analysis:
plt.subplots_adjust(left=0, right=0.99) plt.subplots_adjust(left=0, right=0.99)
fig.tight_layout() fig.tight_layout()
plt.savefig('3D.png', format='png') plt.savefig('3D.png', format='png')
\ No newline at end of file
...@@ -11,7 +11,7 @@ class Particle(): ...@@ -11,7 +11,7 @@ class Particle():
self.vx = 30 self.vx = 30
self.vy = 20 self.vy = 20
self.vz = 300 self.vz = 300
self.m = 0.522 self.m = 0.511
self.q = -1 self.q = -1
class Test(unittest.TestCase): class Test(unittest.TestCase):
......
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