Detector.py

import numpy as np

class Layer:
    def __init__(self, Position, BField):
        self.Position = Position
        self.BField = BField
        self.GridSize = 0.0025 #in centimeter
        
    
    #Newton-method to calculate the hit time    
    def f(self, x, vector):				#EoM function of z-axis in Lorenzforce for magnetic fiel with B = Bx
        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
    
    def df(self, x, vector):				#first derivation of f(x,vector)
        Omega = vector.q * self.BField / vector.m
        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))
    
    def getT(self, vector, x0, e):
        delta = self.dx(x0, vector)
        i = 0
        while delta > e:
            x0 = x0 - self.f(x0, vector)/self.df(x0, vector)
            delta = self.dx(x0, vector)
            i += 1
            if i > 10000: return None		#diverge controll
        return x0
    #End of Newton-method
        
    def detect (self, vector):
        """Calculate for a given particle class the hit grid. The class needs the members: vx, vy, vz, q, m - Velocity in 3 dimensions, charge, mass.
        Returns the passed volume."""
        
        
        #First kollekt hit time, then calculate the position of particle
        t = self.getT(vector, 0 , 1e-10)
        if t == None:
            print ("ERROR - Newton-method do not convert")
            return ((None, None), (None, None), (self.Position,	self.Position), (0,0))
        Omega = vector.q * self.BField / vector.m
        x = vector.vx*t
        y = vector.vy/Omega * np.sin(Omega * t)+vector.vz/Omega * (1-np.cos(Omega * t))
        
	#controll if the particle hit the layer
        if 0 <= x < 50:
            xHigh = int(x/self.GridSize+1)*self.GridSize   #upper edge of hit grid
            xLow = int(x/self.GridSize)*self.GridSize      #lower edge of hit grid
            
        elif -50 < x < 0:               		    #slight different if x<0 or x>0
            xHigh = int(x/self.GridSize)*self.GridSize 
            xLow = int(x/self.GridSize)*self.GridSize    
            
        else:
            return ((None, None), (None, None), (self.Position,	self.Position), (t,t))			 # if the Layer is not hit, return "None"
            
        if 0 <= y < 50:
            yHigh = int(y/self.GridSize+1)*self.GridSize
            yLow = int(y/self.GridSize)*self.GridSize
        
        elif -50 < y < 0:
            yHigh = int(y/self.GridSize)*self.GridSize
            yLow = int(y/self.GridSize-1)*self.GridSize
        
        else:
            return ((None, None), (None, None), (self.Position,	self.Position), (t,t))
        
        return ((xHigh, xLow), (yHigh, yLow), (self.Position+0.0001,	self.Position-0.0001), (t,t))



class Detector:
    
    def __init__ (self, LayerArray = [100,110,120,130,140], Magnet=0.522):
        """Initialize the Layer and magnetic fiel of the detector"""

        self.BField = Magnet 	#Microesla
        self.Layers = []

	#Produce the Layers
        for zPos in LayerArray:
            self.Layers.append(Layer (zPos, self.BField))
        
    def detect (self, particle):
        """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 and time or "None" if one Layer is not hitted"""
        
	#Source
        result = np.array([[(0+0.0001,0-0.0001),(0+0.0001,0-0.0001),(0+0.0001,0-0.0001), (0.,0.)]])
        
        for Layer in self.Layers:
            result = np.append(result, [Layer.detect(particle)], axis = 0)
            if result[-1][0][0] == None:    #If any layer is not hit, exit the detection
                return None
            
        
        return result