Detector.py

import numpy as np

class Layer:
    def __init__(self, Position, BField):
        self.Position = Position
        self.BField = BField
        
    
    #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
    
    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)

    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
        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))
        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))
        

        if 0 <= x < 50:
            xHigh = int(x/0.0025+1)*0.0025   #upper edge of hit grid
            xLow = int(x/0.0025)*0.0025      #lower edge of hit grid
            
        elif -50 < x < 0:                   #slight different if x<0 or x>0
            xHigh = int(x/0.0025)*0.0025 
            xLow = int(x/0.0025-1)*0.0025    
            
        else:
            return ((None, None), (None, None), (self.Position,	self.Position))			 # if the Layer is not hit, return "None"
            
        if 0 <= y < 50:
            yHigh = int(y/0.0025+1)*0.0025
            yLow = int(y/0.0025)*0.0025
        
        elif -50 < y < 0:
            yHigh = int(y/0.0025)*0.0025
            yLow = int(y/0.0025-1)*0.0025
        
        else:
            return ((None, None), (None, None), (self.Position,	self.Position))
        
        return ((xHigh, xLow), (yHigh, yLow), (self.Position+0.005,	self.Position-0.005))



class Detector:
    
    def __init__ (self, Magnet=0.522):
        """5 Layer of Tracker"""
        self.BField = Magnet #Tesla
        self.Layer1 = Layer (100, self.BField)
        self.Layer2 = Layer (110, self.BField)
        self.Layer3 = Layer (120, self.BField)
        self.Layer4 = Layer (130, self.BField)
        self.Layer5 = Layer (140, 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 or "None" if one Layer is not hitted"""
        
        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)
            if result[-1][0][0] == None:    #If any layer is not hit, exit the detection
                return None
            
        
        return result