Skip to content
Snippets Groups Projects
Commit 7c47f61d authored by lorenzennio's avatar lorenzennio
Browse files

Fitting with polynomials, to z-x plane reversed - not sure if this is smart....

Fitting with polynomials, to z-x plane reversed - not sure if this is smart. and added simmple plots
parent 3145ffbb
No related branches found
No related tags found
No related merge requests found
......@@ -39,12 +39,12 @@ class Analysis:
#self.Bounds(self.xBounds['High'], self.xBounds['Low'], call=0)
#self.Bounds(self.yBounds['High'], self.yBounds['Low'], call=1)
self.xPoints = pd.DataFrame({'Mean' : (bounds[:,0,0] + bounds[:,0,1])/2.,
'Error' : (bounds[:,0,0] - bounds[:,0,1])/2.})
self.yPoints = pd.DataFrame({'Mean': (bounds[:,1,0] + bounds[:,1,1])/2.,
'Error': (bounds[:,1,0] - bounds[:,1,1])/2.})
self.zPoints = pd.DataFrame({'Mean': (bounds[:,2,0] + bounds[:,2,1])/2.,
'Error': (bounds[:,2,0] - bounds[:,2,1])/2.})
self.xPoints = pd.DataFrame({'Mean' : (bounds[1:,0,0] + bounds[1:,0,1])/2.,
'Error' : (bounds[1:,0,0] - bounds[1:,0,1])/2.})
self.yPoints = pd.DataFrame({'Mean': (bounds[1:,1,0] + bounds[1:,1,1])/2.,
'Error': (bounds[1:,1,0] - bounds[1:,1,1])/2.})
self.zPoints = pd.DataFrame({'Mean': (bounds[1:,2,0] + bounds[1:,2,1])/2.,
'Error': (bounds[1:,2,0] - bounds[1:,2,1])/2.})
print(self.xPoints)
print(self.yPoints)
......@@ -82,22 +82,31 @@ class Analysis:
return y
def xPol(self, B, x):
a = B[1]/B[0]
b = - (B[0]*B[0])/(2*B[1]**2)
c = - (B[3]*B[0]**2)/(6*B[1]**3)
return a*x + b*x**2 + c*x**3#B[0] + B[1]*x + B[2]*(x**2) + B[3]*(x**3)
def pol(self, B, x):
return B[0] + B[1]*x + B[2]*(x**2) + B[3]*(x**3)
def xFit(self):
xModel = odr.Model(self.pol)#self.xEOM)
xData = odr.RealData(self.xPoints['Mean'],
self.zPoints['Mean'],
sx = self.xPoints['Error'],
sy = self.zPoints['Error'])
xData = odr.RealData(self.zPoints['Mean'],
self.xPoints['Mean'],
sx = self.zPoints['Error'],
sy = self.xPoints['Error'])
xODR = odr.ODR(xData, xModel, beta0=[1, 100, 100, 100])
xOUT = xODR.run()
xOUT.pprint()
print(xOUT.beta)
plt.plot(self.pol(xOUT.beta, self.xPoints['Mean']), self.xPoints['Mean'])
plt.plot(self.zPoints['Mean']+5, self.xPoints['Mean'])
plt.show()
print('x ', xOUT.res_var)
plt.plot(self.zPoints['Mean'],self.pol(xOUT.beta, self.zPoints['Mean']))
plt.plot(self.zPoints['Mean'], self.xPoints['Mean'], 'o')
plt.xlabel('z')
plt.ylabel('x')
plt.savefig('x-z.pdf', format='pdf')
plt.gcf().clear()
"""
xPopt, xPcov = curve_fit(self.xEOM,
self.xPoints['Mean'],
......@@ -127,17 +136,21 @@ class Analysis:
"""
#ODR
yModel = odr.Model(self.yEOM)
yModel = odr.Model(self.pol)
yData = odr.RealData(self.zPoints['Mean'],
self.yPoints['Mean'],
sx = self.zPoints['Error'],
sy = self.yPoints['Error'])
yODR = odr.ODR(yData, yModel, beta0=[1e3, 100, 100])
yODR = odr.ODR(yData, yModel, beta0=[1e3, 100, 100,100])
yOUT = yODR.run()
yOUT.pprint()
plt.plot(self.zPoints['Mean'], self.yEOM(yOUT.beta, self.zPoints['Mean']))
plt.plot(self.zPoints['Mean'], self.yPoints['Mean'])
plt.show()
print('y :', yOUT.res_var)
plt.plot(self.zPoints['Mean'], self.pol(yOUT.beta, self.zPoints['Mean']))
plt.plot(self.zPoints['Mean'], self.yPoints['Mean'], 'o')
plt.xlabel('z')
plt.ylabel('y')
plt.savefig('y-z.pdf', format='pdf')
plt.gcf().clear()
"""
self.results['R'][0] = yPopt[0]
self.results['R'][1] = yPerr[0]
......
......@@ -3,8 +3,8 @@ from Detector import Detector
class Particle:
def __init__(self):
self.vx = 1
self.vy = 1
self.vx = 10
self.vy = 10
self.vz = 300
self.m = 0.522
self.q = 1
......@@ -15,5 +15,6 @@ DETRes = DET.detect(Particle())
ANA = Analysis()
ANA.fill(DETRes)
ANA.xFit()
ANA.yFit()
......@@ -15,8 +15,8 @@ class testAnalysis(unittest.TestCase):
#self.assertEqual(A.results['Mean'][0], 5.5)
#self.assertEqual(A.results['Error'][0], 0.5)
x = np.arange(0, 6)
z = A.xEOM([1,20,10,100], x)
x = np.arange(0, 10, 2)
z = A.xEOM([1,10,10,300], x)
#plt.plot(z,x)
#plt.show()
......@@ -43,8 +43,8 @@ class testAnalysis(unittest.TestCase):
#self.assertEqual(A.results['Error'][0], 0.5)
z = np.arange(0, 6)
y = A.yEOM([1000,100,100], z)
z = np.arange(90,150,10)
y = A.yEOM([1000,10,10,300], z)
#plt.plot(z,y)
#plt.show()
......
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