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