Newer
Older
import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
"""
Analysis class:
This class computes mean and error on the detector output angles.
The results can then be visualized with two plots.
"""
"""
The class is initialized with the angles from the detector results
and then automatically deletes 'None' entries and computes the minimum
range of angles and the error.
"""
...
def fill(self, bounds):
self.xBounds = pd.Series({'High' : bounds[:,0,0],
'Low' : bounds[:,0,1]})
self.yBounds = pd.Series({'High' : bounds[:,1,0],
'Low' : bounds[:,1,1]})
self.results = pd.DataFrame({'w' : np.zeros(2),
'Vx' : np.zeros(2),
'Vy' : np.zeros(2),
'Vz' : np.zeros(2),
'R' : np.zeros(2),
'y0' : np.zeros(2),
'z0' : np.zeros(2)})
#self.Bounds(self.xBounds['High'], self.xBounds['Low'], call=0)
#self.Bounds(self.yBounds['High'], self.yBounds['Low'], call=1)
lorenzennio
committed
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)
print(self.zPoints)
#plt.plot(self.zPoints['Mean'], self.xPoints['Mean'])
#plt.show()
#self.xFit()
#self.yFit()
#print(self.results)
Fill the angle data into the class variables.
#print(self.xBounds['High'])
self.xBounds['High'] = np.array([x for x in self.xBounds['High'] if x is not None])
self.xBounds['Low'] = np.array([x for x in self.xBounds['Low'] if x is not None])
self.yBounds['High'] = np.array([x for x in self.yBounds['High'] if x is not None])
self.yBounds['Low'] = np.array([x for x in self.yBounds['Low'] if x is not None])
def xEOM(self, B, x):
"""
B= [w, Vx, Vy, Vz]
"""
z = (1/B[0])*(-B[2]*(1-np.cos(x*B[0]/B[1])) + B[3]*np.sin(x*B[0]/B[1]))
def yEOM(self, B, z):
"""
B = [R, y0, z0]
"""
y = B[1] + np.sqrt(B[0]**2 - (z-B[2])**2)
lorenzennio
committed
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)
lorenzennio
committed
xData = odr.RealData(self.zPoints['Mean'],
self.xPoints['Mean'],
sx = self.zPoints['Error'],
sy = self.xPoints['Error'])
lorenzennio
committed
xODR = odr.ODR(xData, xModel, beta0=[1, 100, 100, 100])
xOUT = xODR.run()
xOUT.pprint()
lorenzennio
committed
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'],
self.zPoints['Mean']
#sigma=self.xPoints['Error'],
#absolute_sigma=True)
)
xPerr = np.sqrt(np.diag(xPcov))
self.results['w'][0] = xPopt[0]
self.results['w'][1] = xPerr[0]
self.results['Vx'][0] = xPopt[1]
self.results['Vx'][1] = xPerr[1]
self.results['Vy'][0] = xPopt[2]
self.results['Vy'][1] = xPerr[2]
self.results['Vz'][0] = xPopt[3]
self.results['Vz'][1] = xPerr[3]
"""
def yFit(self):
"""
yPopt, yPcov = curve_fit(self.yEOM,
self.zPoints['Mean'],
self.yPoints['Mean'],
#sigma=self.yPoints['Error'],
p0=[10^3, 1, 1],
absolute_sigma=True)
yPerr = np.sqrt(np.diag(xPcov))
lorenzennio
committed
yModel = odr.Model(self.pol)
yData = odr.RealData(self.zPoints['Mean'],
self.yPoints['Mean'],
sx = self.zPoints['Error'],
sy = self.yPoints['Error'])
lorenzennio
committed
yODR = odr.ODR(yData, yModel, beta0=[1e3, 100, 100,100])
yOUT = yODR.run()
yOUT.pprint()
lorenzennio
committed
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]
self.results['y0'][0] = yPopt[1]
self.results['y0'][1] = yPerr[1]
self.results['z0'][0] = yPopt[2]
self.results['z0'][1] = yPerr[2]
"""
def yEOM(self, t, w, yVel, zVel):
y = zVel + (1./w)*(-zVel*np.cos(w*t) + yVel*np.sin(w*t))
return y
def zEOM(self, t, w, yVel, zVel):
z = -yVel + (1./w)*(yVel*np.cos(w*t) + zVel*np.sin(w*t))
return z
"""
def Bounds(self, high, low, call=0):
"""
Algorithm to compute the minimum possible range of angles.
The upper and lower angles must be given as a argument and also
the 'call' variable secifies the plane you are working in
(x-z or y-z plane).
The algorithim loops over the angles from the last detector layer
to the first, and replaces the upper and lower angle bounds, if
a certain layer restricts the range further.
The mean is calculated as the arithmetic mean of the resulting
angles and the error is the difference from the mean to the angle
bounds.
"""
highRev = high[::-1]
lowRev = low[::-1]
highBound = highRev[0]
lowBound = lowRev[0]
for a,b in zip(highRev[1:], lowRev[1:]):
if a < highBound and a > lowBound:
highBound = a
if b < highBound and b > lowBound:
lowBound = b
error = (highBound - lowBound)/2
mean = highBound - error
self.results['Mean'][call] = mean
self.results['Error'][call] = error