Merged and started testing with new fitting

Detector.py

@@ -9,11 +9,11 @@ class Layer:
     #Newton-method to calculate the hit time    
     def f(self, x, vector):
         Omega = vector.q * self.BField / vector.m
-        return -vector.vz/Omega * np.sin(Omega * x)+vector.vy/Omega * (np.cos(Omega * x)-1)-self.Position
+        return vector.vz/Omega * np.sin(Omega * x)+vector.vy/Omega * (np.cos(Omega * x)-1)-self.Position
     
     def df(self, x, vector):
         Omega = vector.q * self.BField / vector.m
-        return -vector.vz * np.cos(Omega* x)+vector.vy * np.sin(Omega * x)
+        return vector.vz * np.cos(Omega* x)+vector.vy * np.sin(Omega * x)
 
     def dx(self, x, vector):
         return abs(0-self.f(x, vector))
@@ -66,7 +66,7 @@ class Layer:
         else:
             return ((None, None), (None, None), (self.Position,	self.Position))
         
-        return ((xHigh, xLow), (yHigh, yLow), (self.Position,	self.Position))
+        return ((xHigh, xLow), (yHigh, yLow), (self.Position+0.005,	self.Position-0.005))
 
 
 
@@ -85,7 +85,7 @@ class Detector:
         """Calculate for a given particle class the hit grid. The class needs the membervalues: vx, vy, vz, q, m - Velocity in 3 dimensions, charge, mass.
         Returns the known passed volumes or "None" if one Layer is not hitted"""
         
-        result = np.array([[(0,0),(0,0),(0,0)]])
+        result = np.array([[(0+0.0001,0-0.0001),(0+0.0001,0-0.0001),(0+0.0001,0-0.0001)]])
         
         for Layer in [self.Layer1, self.Layer2, self.Layer3, self.Layer4, self.Layer5]:
             result = np.append(result, [Layer.detect(particle)], axis = 0)
@@ -93,4 +93,4 @@ class Detector:
                 return None
             
         
-        return result
+        return result
\ No newline at end of file

analysis.py

 import numpy as np
 import pandas as pd
-import matplotlib as plt
+import matplotlib.pyplot as plt
 from scipy.optimize import curve_fit
+from scipy import odr
 
 
 class Analysis:
@@ -11,12 +12,15 @@ class Analysis:
         The results can then be visualized with two plots.
     """
 
-    def __init__(self, bounds):
+    def __init__(self):
         """
         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],
@@ -41,8 +45,16 @@ class Analysis:
                                     '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.fit()
-        print(self.results)
+        
+        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)
         #TODO: Plot
 
     def rmNone(self):
@@ -55,15 +67,30 @@ class Analysis:
         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, x, w, Vx, Vy, Vz):
-        z = (1/w)*(-Vy*(1-np.cos(x*w/Vx)) + Vz*np.sin(x*w/Vx))
+    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]))
         return z
 
-    def yEOM(self, z, R, y0, z0):
-        y = y0 + np.sqrt(R**2 - (z-z0)**2)
+    def yEOM(self, B, z):
+        """
+        B = [R, y0, z0]
+        """
+        y = B[1] + np.sqrt(B[0]**2 - (z-B[2])**2)
         return y
 
-    def fit(self):
+    def xFit(self):
+        xModel  = odr.Model(self.xEOM)
+        xData   = odr.RealData(self.xPoints['Mean'],
+                               self.zPoints['Mean'],
+                               sx = self.xPoints['Error'],
+                               sy = self.zPoints['Error'])
+        xODR = odr.ODR(xData, xModel, beta0=[1, 1, 1, 1])
+        xOUT = xODR.run()
+        xOUT.pprint()
+        """
         xPopt, xPcov = curve_fit(self.xEOM,
                                  self.xPoints['Mean'],
                                  self.zPoints['Mean']
@@ -79,20 +106,35 @@ class Analysis:
         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'],
+                                 #sigma=self.yPoints['Error'],
+                                 p0=[10^3, 1, 1],
                                  absolute_sigma=True)
         yPerr = np.sqrt(np.diag(xPcov))
+        """
+        #ODR
+        
+        yModel  = odr.Model(self.yEOM)
+        yData   = odr.RealData(self.zPoints['Mean'],
+                               self.yPoints['Mean'],
+                               sx = self.zPoints['Error'],
+                               sy = self.yPoints['Error'])
+        yODR = odr.ODR(yData, yModel, beta0=[1e3, 1, 1])
+        yOUT = yODR.run()
+        yOUT.pprint()
+        """
         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))

test.py

+from analysis import Analysis
+from Detector import Detector
+
+class Particle:
+    def __init__(self):
+        self.vx = 1
+        self.vy = 1
+        self.vz = 300
+        self.m = 0.522
+        self.q = 1
+
+
+DET = Detector()
+DETRes = DET.detect(Particle())
+
+ANA = Analysis(DETRes)
+

testAnalysis.py

 import unittest
 import numpy as np
 from analysis import Analysis
+import matplotlib.pyplot as plt
 
 class testAnalysis(unittest.TestCase):
     """
     A simple test to check the angle means ans error computation.
     """
     def test(self):
-        ang = np.array([6,5,6,5,8,4,8,4,7,1,7,1,10,3,10,3,None,None,None,None])
-        ang = np.reshape(ang, (5,2,2))
-        A = Analysis(ang)
-        self.assertEqual(A.results['Mean'][0], 5.5)
-        self.assertEqual(A.results['Error'][0], 0.5)
+        #ang = np.array([6,5,6,5,8,4,8,4,7,1,7,1,10,3,10,3,None,None,None,None])
+        #ang = np.reshape(ang, (5,2,2))
+        A = Analysis()
+        #self.assertEqual(A.results['Mean'][0], 5.5)
+        #self.assertEqual(A.results['Error'][0], 0.5)
+        
+
+        z = np.arange(30)
+        y = A.yEOM([100,1,1], z)
+        #plt.plot(z,y)
+        #plt.show()
+        
+        Data = []
+        for i in range(0,len(z)):
+            err = np.random.random(size=(4))
+            print(err)
+            Data.append([[0,0], [y[i]+err[0],y[i]-err[1]], [z[i]+err[2], z[i]-err[3]]])
+            #np.concatenate((Data,
+            #        [[0,0], [y[i]+0.1,y[i]-0.1], [z[i]+0.1, z[i]-0.1]]), axis=0)
+        Data = np.array(Data)
+
+        A.fill(Data)
+        A.yFit()
+        #print(y)
+        #print(z)
+        #print(Data)