delete specific count

isotropical in space

CHAOS.py

@@ -0,0 +1,250 @@
+#! /usr/bin/python3
+import numpy as np
+from Struktur import run, Rays, vector, Plane, Scintillator
+from Geometry import regular_hexagonal_prism
+import argparse
+import csv
+import sys
+
+parser = argparse.ArgumentParser(prog = 'CHAOS', description = 'This program tracks the propagation of light rays in a somehow shaped scintillator', epilog = 'Whatever')
+
+parser.add_argument('-s', '--size', default = 10, type = int, help = 'size of the simulation(number of rays)')
+
+parser.add_argument('--origin', type = str, help = 'origin of the bunch of rays')
+
+parser.add_argument('--trigger', default = False, type = int, help = 'origin of the bunch of rays')
+
+parser.add_argument('--height', default = 10, type = int, help = 'height of the scintillator')
+
+parser.add_argument('-e', '--edge', default = 51.962, type = float, help = 'edge length of the scintillator')
+
+parser.add_argument('--stop', default = 5000, type = float, help = 'cut off path length')
+
+parser.add_argument('-z', '--z_value', default = False, type = str, help = 'cut off path length')
+
+parser.add_argument('-i', '--refraction_index', default = 2.15, type = float, help = 'index of refraction of the BGO cristall')
+
+parser.add_argument('-a', '--all', default = False, type = bool, help = 'parameter for rastering all')
+
+parser.add_argument('-x', '--x_value', default = False, type = str, help = 'vertical slice for rasterizing (from -25 to 25)')
+
+parser.add_argument('-o', '--outfn', default = None, help = 'Output file name')
+
+args = parser.parse_args()
+
+
+def check_range_x(value):
+    ivalue = int(value)
+    if ivalue < -50 or ivalue > 51:  # Modify the range as needed
+        raise argparse.ArgumentTypeError(f"{value} is not in the valid range [-50,51]") 
+    return ivalue
+
+size = args.size
+height = args.height
+edge = args.edge
+stop = args.stop
+Ind_S = args.refraction_index
+all = args.all
+x_value = args.x_value
+trigger = args.trigger
+z_value = args.z_value
+check_range_x(x_value)
+outfn = args.outfn
+
+if args.origin:
+    try:
+        x, y, z = map(float, args.origin.split(','))
+        origin = vector(x,y,z)
+    except ValueError:
+        print("Error: Invalid input format. Please provide the origin in the format 'x,y,z'")
+
+
+if all:
+    pos = 0
+    origin = np.empty((145992,3,1))
+    x_values = int((edge/2)+(np.cos(np.radians(60))*edge))
+    for i in range(x_values*2):
+        if x_values < 25 and x_values > -25:
+            apothem = int(np.sin(np.radians(60))*edge)
+            o = apothem*2*(height*2+1)
+            origin[pos : pos + o, 0] = x_values
+            k = np.empty((apothem*2,1))
+            k[:,0] = np.array(np.arange(-apothem, apothem, 1))
+            for i in range(height*2+1):        
+                origin[pos + apothem*2*i:pos + (i+1)*apothem*2,1] = k
+                origin[pos + apothem*2*i:pos + (i+1)*apothem*2,2] = height - i
+            x_values = x_values - 1
+            pos = pos + o
+        else:
+            if x_values >= 25:
+                apothem = int(np.sin(np.radians(60))*edge - abs(np.tan(np.radians(60))*(x_values-edge/2)))
+            else:
+                apothem = int(np.sin(np.radians(60))*edge - abs(np.tan(np.radians(60))*(x_values+edge/2)))
+            o = apothem*2*(height*2+1)
+            origin[pos:pos + o, 0] = x_values
+            k = np.empty((apothem*2, 1))
+            k[:,0] = np.array(np.arange(-apothem, apothem, 1))
+            for i in range(height*2+1):
+                origin[pos + apothem*2*i:pos + (i+1)*apothem*2, 1] = k
+                origin[pos + apothem*2*i:pos + (i+1)*apothem*2, 2] = height - i
+            x_values = x_values - 1
+            pos = pos + o
+
+if z_value:
+    z = int(z_value)
+    origin = origin[(np.where(origin[:,2]==z)[0])]
+
+
+if trigger == 12:
+    origin = origin[(np.where(origin[:,0] <= 30)[0])]
+    origin = origin[(np.where(origin[:,0] >= 20)[0])]
+    origin = origin[(np.where(origin[:,1] <= -3)[0])]
+    origin = origin[(np.where(origin[:,1] >= -15)[0])]
+
+if trigger == 9:
+    origin = origin[(np.where(origin[:,0] <= 5)[0])]
+    origin = origin[(np.where(origin[:,0] >= -5)[0])]
+    origin = origin[(np.where(origin[:,1] <= 5)[0])]
+    origin = origin[(np.where(origin[:,1] >= -5)[0])]
+
+if trigger == 2:
+    origin = origin[(np.where(origin[:,0] <= 5)[0])]
+    origin = origin[(np.where(origin[:,0] >= -5)[0])]
+    origin = origin[(np.where(origin[:,1] <= 21)[0])]
+    origin = origin[(np.where(origin[:,1] >= 11)[0])]
+
+if trigger == 5:
+    origin = origin[(np.where(origin[:,0] <= 5)[0])]
+    origin = origin[(np.where(origin[:,0] >= -5)[0])]
+    origin = origin[(np.where(origin[:,1] <= 13)[0])]
+    origin = origin[(np.where(origin[:,1] >= 3)[0])]
+
+if trigger == 7:
+    origin = origin[(np.where(origin[:,0] <= 39)[0])]
+    origin = origin[(np.where(origin[:,0] >= 29)[0])]
+    origin = origin[(np.where(origin[:,1] <= 5)[0])]
+    origin = origin[(np.where(origin[:,1] >= -5)[0])]
+
+if trigger == 13:
+    origin = origin[(np.where(origin[:,0] <= 5)[0])]
+    origin = origin[(np.where(origin[:,0] >= -5)[0])]
+    origin = origin[(np.where(origin[:,1] >= -13)[0])]
+    origin = origin[(np.where(origin[:,1] <= -3)[0])]
+
+if trigger == 16:
+    origin = origin[(np.where(origin[:,0] <= 5)[0])]
+    origin = origin[(np.where(origin[:,0] >= -5)[0])]
+    origin = origin[(np.where(origin[:,1] >= -21)[0])]
+    origin = origin[(np.where(origin[:,1] <= -11)[0])]
+
+
+if x_value:
+    x = int(x_value)
+    if x < 25 and x > -25:
+        apothem = int(np.sin(np.radians(60))*edge)
+        origin = np.empty(((apothem*2-1)*(height*2),3,1))
+        origin[:,0] = x
+        k = np.empty((apothem*2-1,1))
+        k[:,0] = np.array(np.arange(-apothem+0.5, apothem-0.5, 1))
+        for i in range(height*2):      
+            origin[(apothem*2*i)-i:((i+1)*apothem*2)-i-1,1] = k
+            origin[(apothem*2*i)-i:((i+1)*apothem*2)-i-1,2] = height - 0.5 - i
+    else:
+        if x >= 25:
+            apothem = int(np.sin(np.radians(60))*edge - abs(np.tan(np.radians(60))*(x-edge/2)))
+        else:
+            apothem = int(np.sin(np.radians(60))*edge - abs(np.tan(np.radians(60))*(-edge/2-x)))
+        origin = np.empty(((apothem*2-1)*(height*2),3,1))
+        origin[:,0] = x
+        k = np.empty((apothem*2-1,1))
+        k[:,0] = np.array(np.arange(-apothem+0.5, apothem-0.5, 1))
+        for i in range(height*2):
+            origin[(apothem*2*i)-i:((i+1)*apothem*2)-i-1,1] = k
+            origin[(apothem*2*i)-i:((i+1)*apothem*2)-i-1,2] = height - 0.5 - i
+
+
+
+det=[Plane(vector(-10,-35,5),vector(-10,-35,-5),vector(-10,-30,0),[]),
+      Plane(vector(-10,-35,-5),vector(10,-45,-5),vector(0,-35,-5),[]),
+      Plane(vector(10,-35,5),vector(10,-35,0),vector(10,-45,-5),[]),
+      Plane(vector(-10,-35,5),vector(0,-35,5),vector(10,-45,5),[])]
+
+o = regular_hexagonal_prism([height*2, edge],[1,4])
+p = Scintillator(o, det)
+
+if outfn is None or outfn=='-':
+    csvfile = sys.stdout
+
+else: 
+    csvfile = open(outfn, "w" )
+
+if x_value:
+    x = int(x_value)
+    for i in range(int(len(origin)/(2*height))):
+        r = Rays(size, origin[2*height*i:2*height*(i+1)])
+        n = np.empty((2*height*size,1,3))
+        n[...,0] = r.origin[:,0]
+        n[...,1] = r.origin[:,1]
+        n[...,2] = r.origin[:,2]
+        n = n.reshape(-1, 3)
+        n = [','.join(map(str, row)) for row in n]
+        #calculate the incident th and phi angles
+        th, phi = Rays.angle(r)
+        sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
+        #edit hit_rec so it can be written in one row
+        hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
+        #create data list
+        data = list(zip(n, th, phi, sc, tr, path, hit_rec_processed))
+        #write output in csv_file
+        for row in data:
+            csvfile.write(','.join(map(str, row)) + '\n')
+
+else:
+    ind = len(origin) // 100
+    rest = len(origin) - ind * 100
+    for i in range(ind):
+        ori = origin[i*100:(i+1)*100]
+        r = Rays(size, ori)
+        th, phi = Rays.angle(r)
+        sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
+        #create origin vector for output file
+        n = np.empty((100*size,1,3))
+        if len(origin)>3:
+            for i in range(100):
+                n[size*i:size*(i+1),0] = ori[i].T
+        else: 
+            n[...,0] = origin[0]
+            n[...,1] = origin[1]
+            n[...,2] = origin[2]
+        n = n.reshape(-1, 3)
+        n = [','.join(map(str, row)) for row in n]
+        #edit hit_rec so it can be written in one row
+        hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
+        #create data list
+        data = list(zip(n, th, phi, sc, tr, path, hit_rec_processed))
+        #write output in csv_file
+        for row in data:
+            csvfile.write(','.join(map(str, row)) + '\n')
+    ori = origin[ind*100:]
+    r = Rays(size, ori)
+    th, phi = Rays.angle(r)
+    sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
+    #create origin vector for output file
+    n = np.empty((rest*size,1,3))
+    if len(origin)>3:
+        for i in range(rest):
+            n[size*i:size*(i+1),0] = ori[i].T
+    else: 
+        n[...,0] = origin[0]
+        n[...,1] = origin[1]
+        n[...,2] = origin[2]
+    n = n.reshape(-1, 3)
+    n = [','.join(map(str, row)) for row in n]
+    #edit hit_rec so it can be written in one row
+    hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
+    #create data list
+    data = list(zip(n, th, phi, sc, tr, path, hit_rec_processed))
+    #write output in csv_file
+    for row in data:
+        csvfile.write(','.join(map(str, row)) + '\n')
+

Geometry.py

@@ -0,0 +1,48 @@
+import math
+import numpy as np
+
+#to see the definitions of the plane numbers
+#and dimensions see geometry.jpg in this repository
+
+def cube(dimensions, det):
+    '''generates points for Planes 
+    dimensions for cube: [edge length]
+    det: python indices of planes with detectors on them (first plane has index 0)'''
+    edge=dimensions
+    return np.array(([[edge/2,edge/2,0],[edge/4,edge/2,0],[edge/2,edge/4,0]],
+                    [[edge/2,edge,edge/2],[edge/4,edge,edge/2],[edge/2,edge,edge/4]],
+                    [[edge/2,edge/2,edge],[edge/4,edge/2,edge],[edge/2,edge/4,edge]],
+                    [[edge/2,0,edge/2],[edge/4,0,edge/2],[edge/2,0,edge/4]],
+                    [[0,edge/2,edge/2],[0,edge/4,edge/2],[0,edge/2,edge/4]],
+                    [[edge,edge/2,edge/2],[edge,edge/4,edge/2],[edge,edge/2,edge/4]])),det
+
+def cuboid(dimensions, det):
+        '''generates points on planes
+        dimensions for cuboid: [length, heigth, depth]
+        det: python indices of planes with detectors on them (first plane has index 0)'''
+        length=dimensions[0]
+        height=dimensions[1]
+        depth=dimensions[2]
+        return np.array(([[depth/2,length/2,0],[depth/4,length/2,0],[depth/2,length/4,0]],
+                         [[depth/2,length,height/2],[depth/4,length,height/2],[depth/2,length,height/4]],
+                         [[depth/2,length/2,height],[depth/4,length/2,height],[depth/2,length/4,height]],
+                         [[depth/2,0,height/2],[depth/4,0,height/2],[depth/2,0,height/4]],
+                         [[depth,length/2,height/2],[depth,length/4,height/2],[depth,length/2,height/4]],
+                         [[0,length/2,height/2],[0,length/4,height/2],[0,length/2,height/4]])),det
+
+def regular_hexagonal_prism(dimensions, det):
+        '''generates points for planes
+        dimensions for hexagonal prism: [height, edge length]
+        det: python indices of planes with detectors on them (first plane has index 0)'''
+        height=dimensions[0]
+        edge=dimensions[1]
+        apothem=np.sin(np.radians(60))*edge
+        dia=np.cos(np.radians(60))*edge
+        return np.array(([[-edge/2, -apothem, 0],[ -dia/2-edge/2, -apothem/2, -height/4],[-dia/2-edge/2, -apothem/2, 0]],
+                         [[0, -apothem, -height/4],[-edge/4, -apothem, 0],[0, -apothem, 0]],
+                         [[dia/2+edge/2, -apothem/2, 0],[dia/2+edge/2, -apothem/2, -height/4],[edge/2, -apothem, 0]],
+                         [[dia/2+edge/2, apothem/2, height/2],[edge/2, apothem, 0],[dia/2+edge/2, apothem/2, height/4]],
+                         [[0, apothem, 0],[-edge/4, apothem, 0],[0, apothem, -height/4]],
+                         [[-dia/2-edge/2, apothem/2, 0],[-dia/2-edge/2, apothem/2, -height/4],[-edge/2, apothem, 0]],
+                         [[0, 0, height/2],[0, -apothem/2, height/2],[-edge/4, 0, height/2]],
+                         [[0, apothem/4, -height/2],[0, apothem/2 ,-height/2],[-edge/4, apothem/4, -height/2]])),det

Random_gen.py

@@ -0,0 +1,41 @@
+import random
+from Struktur import Plane,Ray
+
+def test_planes(test_size):
+    rand_planes=[]
+    rand_vec=[]
+    vec=[]
+    for i in range(test_size*9):
+        vec.append(random.randrange(-100,100,1))
+        if len(vec)==3:
+            rand_vec.append(vec)
+            vec=[]
+        if len(rand_vec)==3:
+            rand_planes.append(Plane(rand_vec[0],rand_vec[1],rand_vec[2]))
+            rand_vec=[]
+    return rand_planes
+
+def test_lines(test_size):
+    rand_lines=[]
+    rand_vec=[]
+    vec=[]
+    for i in range(test_size*6):
+        vec.append(random.randrange(-100,100,1))
+        if len(vec)==3:
+            rand_vec.append(vec)
+            vec=[]
+        if len(rand_vec)==2:
+            rand_lines.append(Ray(rand_vec[0],rand_vec[1]))
+            rand_vec=[]
+    return rand_lines
+
+p=test_planes(1)
+g=test_lines(1)
+print(p[0].dir1)
+print(p[0].dir2)
+print(p[0].o)
+print(g[0].origin)
+print(g[0].dir)
+print(Plane.cross(p[0],g[0]))
+print(p[0].a)
+print(p[0].norm)

Struktur.py

@@ -0,0 +1,398 @@
+from sys import argv, stdout, stderr
+import numpy as np
+import random
+from numpy.linalg import solve
+import time
+from Geometry import cube, cuboid, regular_hexagonal_prism
+
+#tolerance for calculation
+tolerance = 1e-9
+
+def matrix(vec1, vec2, vec3):
+    return np.array([vec1, vec2, vec3])
+
+def vector(x=0, y=0, z=0):
+    return np.array([[x], [y], [z]])
+
+def ref_angle(Ind_in,Ind_out):
+    return np.rad2deg(np.arcsin(Ind_out/Ind_in))
+
+
+def rotate_any(axis, vec, angle):
+    '''performs a rotation around any axis in the 3D space defined by a vector
+    using the Rodriguez formula'''
+    K = matrix( [0, -axis[0,2], axis[0,1]],
+                [axis[0,2], 0, axis[0,0]],
+                [-axis[0,1], axis[0,0], 0] )
+    A = np.identity(3) + np.sin(angle) * K + (1-np.cos(angle)) * np.dot(K,K)
+    A = A[np.newaxis, :, :]
+    A = np.tile(A, (len(vec), 1, 1))
+    h = np.sum(A*vec, axis=1)
+    h = h[:, :, np.newaxis]
+    return h
+
+    
+class Rays_custom:
+
+    def __init__(self, o, dir):
+        '''self.o is the starting point. self.d the direction vector of the line
+        o and p are two points on the line'''
+        self.origin = o
+        self.dir = dir
+        self.A = np.empty((len(dir), 3, 3))
+        self.A[...,0:1] = -self.dir
+
+    def angle(self):
+        '''gives back the angles in sphere coordinates theta and phi'''
+        th = np.arccos(self.dir[...,2,0] / np.sqrt(self.dir[...,0,0]**2+self.dir[...,1,0]**2+self.dir[...,2,0]**2))
+        phi = np.empty(len(self.dir))
+
+        s_1 = self.dir[...,0,0] > 0
+        s_2 = self.dir[...,0,0] == 0
+        s_3 = np.logical_and(self.dir[...,0,0] < 0, self.dir[...,1,0] > 0)
+        s_4 = np.logical_and(self.dir[...,0,0] < 0, self.dir[...,1,0] < 0)
+        phi_1 = np.arctan(self.dir[...,1,0][s_1] / self.dir[...,0,0][s_1])
+        phi_2 = np.pi/2 * np.sign(self.dir[...,1,0][s_2])
+        phi_3 = np.arctan(self.dir[...,1,0][s_3] / self.dir[...,0,0][s_3]) + np.pi
+        phi_4 = np.arctan(self.dir[...,1,0][s_4] / self.dir[...,0,0][s_4]) - np.pi
+       
+        phi[s_1] = phi_1
+        phi[s_2] = phi_2
+        phi[s_3] = phi_3
+        phi[s_4] = phi_4
+
+        th=np.rad2deg(th)
+        phi=np.rad2deg(phi)
+
+        return th, phi
+
+
+    def __str__(self):
+        return str(self.origin) + str(self.dir)
+    
+class Rays:
+
+    def __init__(self, size, origin_list):
+        '''creates a beam of rays with the same origin but equally distributed
+        direction on a sphere using the Fibonacci sphere algorithm
+        self.origin and self.dir are (n,3,1) arrays'''
+        self.size = size
+        self.origin = np.empty((self.size*len(origin_list), 3, 1))
+        for i in range(len(origin_list)):
+            self.origin[self.size*i:self.size*(i+1),0:3] = origin_list[i]
+        self.dir = np.empty((self.size*len(origin_list), 3, 1))
+        for i in range(len(origin_list)):
+            ga = (3 - np.sqrt(5)) * np.pi
+            th = ga * np.arange(self.size)
+            z = np.linspace(1 / self.size - 1, 1 - 1 / self.size, self.size)
+            radius = np.sqrt(1 - (z**2))
+            y = radius * np.sin(th)
+            x = radius * np.cos(th)
+            self.dir[self.size*i:self.size*(i+1),0,0] = x
+            self.dir[self.size*i:self.size*(i+1),1,0] = y
+            self.dir[self.size*i:self.size*(i+1),2,0] = z
+        self.A = np.empty((self.size*len(origin_list), 3, 3))
+        self.A[...,0:1] = -self.dir
+        #fig = plt.figure()
+        #ax = fig.add_subplot(111, projection='3d')
+        #ax.scatter(x, y, z)
+        #plt.show()
+
+    def angle(self):
+        '''gives back the angles in sphere coordinates theta and phi'''
+        th = np.arccos(self.dir[...,2,0] / np.sqrt(self.dir[...,0,0]**2+self.dir[...,1,0]**2+self.dir[...,2,0]**2))
+        phi = np.empty(len(self.dir))
+
+        s_1 = self.dir[...,0,0] > 0
+        s_2 = self.dir[...,0,0] == 0
+        s_3 = np.logical_and(self.dir[...,0,0] < 0, self.dir[...,1,0] > 0)
+        s_4 = np.logical_and(self.dir[...,0,0] < 0, self.dir[...,1,0] < 0)
+        phi_1 = np.arctan(self.dir[...,1,0][s_1] / self.dir[...,0,0][s_1])
+        phi_2 = np.pi/2 * np.sign(self.dir[...,1,0][s_2])
+        phi_3 = np.arctan(self.dir[...,1,0][s_3] / self.dir[...,0,0][s_3]) + np.pi
+        phi_4 = np.arctan(self.dir[...,1,0][s_4] / self.dir[...,0,0][s_4]) - np.pi
+       
+        phi[s_1] = phi_1
+        phi[s_2] = phi_2
+        phi[s_3] = phi_3
+        phi[s_4] = phi_4
+
+        th=np.rad2deg(th)
+        phi=np.rad2deg(phi)
+
+        return th, phi
+
+
+
+class Plane:
+
+    def __init__(self, o, v, w, det):
+        '''o, v, w are points on the plane'''
+        self.o = o
+        self.v = v
+        self.w = w
+        self.dir1 = o - v
+        self.dir2 = v - w
+        self.norm = np.cross(self.dir1, self.dir2, axis=0).T
+        self.norm = self.norm / np.linalg.norm(self.norm)
+        if np.any(det):
+            self.det = det
+        self.a = np.dot(self.norm,o)
+        self.tolerance = tolerance
+      
+    def rotate_matrix(self):
+        '''calculates the matrix for the rotate function'''
+        return matrix(  [0, -self.norm[2], self.norm[1]],
+                        [self.norm[2], 0, -self.norm[0]],
+                        [-self.norm[1], self.norm[0], 0] )
+    
+    const=1e-3
+    def cross(self, g):
+        '''takes an array of lines and a plane ,filters out all the linear systems
+        without solution and gives back an array with all cross_points and a 
+        boolean array signifiying their position in the original array'''
+        g.A[...,1:2] = self.o - self.v
+        g.A[...,2:3] = self.w - self.v
+        g.A[...,0:1] = -g.dir
+        slc = np.abs(np.linalg.det(g.A)) >= self.const #filters all rays that are not parrallel to plane
+        if not np.any(slc):
+            return [], slc
+        else:
+            A = g.A[slc]  #constructing A and b for the operation
+            b = (g.origin-self.v)[slc]  
+            t = np.linalg.solve(A, b)[:,0:1] #using numpy (linalg) to solve equation
+            return g.origin[slc] + t * g.dir[slc], slc #calculate resulting cross point
+    
+    def cross_angle(self, g, slc):
+        '''calculates the angle between the normal vector of the plane 
+        and the direction vector of the line'''
+        if np.any(slc):
+            ra = g.dir[slc]
+            n_len = np.linalg.norm(self.norm)
+            dir_len = np.linalg.norm(ra, axis=1)
+            m = np.empty((len(dir_len), 3, 1))
+            m[...,0] = self.norm
+            ra = ra.reshape(-1, 1, 3)
+            sc = (ra @ m).flatten()
+            sc_t = sc <= 0
+            sc = sc[sc_t]
+            dir_len = dir_len[sc_t]
+            arg = sc / (n_len * dir_len)[...,0]
+            res = np.rad2deg(np.arccos(arg))
+            res[np.where(res>90)]-=90
+            return res
+        else:
+            return None
+    
+    def cross_planes(self, pl2):
+        '''calculates the cross line of two planes
+        returns the origin and direction vector of the line'''
+        cv = np.cross(self.norm, pl2.norm)
+        if np.allclose(cv, 0):
+            return 'The planes are parallel'
+        else:
+            A = np.array([[self.norm[0], self.norm[1]], [pl2.norm[0], pl2.norm[1]]])
+            b = np.array([-self.a, -pl2.a])
+            cross_point = np.linalg.solve(A, b)
+            origin = cross_point[0] * self.norm
+        return origin, cv
+    
+    def reflect(self, g, slc):
+        '''calculates the reflection vector of a line g and a plane '''
+        rot_vec = -rotate_any(self.norm, g.dir[slc], np.pi)
+        return rot_vec
+        
+    
+    def scatter(self, g, Ind_S, Ind_A, slc):
+        '''calculates the new ray after being refracted outside the Scintillator, the
+        scattered at the tape and refracted back into the medium
+        Ind_S(global variable) is the refraction Inedx of the Scintillator, Ind_A the one of
+        the surrounding medium(Air)'''
+        u = np.random.uniform(0, 1, len(g.dir[slc]))
+        v = np.random.uniform(0, 1, len(g.dir[slc]))
+        phi = 2 * np.pi * u
+        sin_theta = Ind_A / Ind_S * np.sqrt(1 - v**2)
+        new_vec = np.empty((len(phi), 3, 1))
+        new_vec[:,0,0] = sin_theta * np.cos(phi)
+        new_vec[:,1,0] = sin_theta * np.sin(phi) 
+        new_vec[:,2,0] = np.sqrt(1 - (sin_theta)**2)
+        nr = np.empty((len(new_vec), 1, 3))
+        nr[...,0:3] = self.norm
+        b = (nr @ new_vec).flatten() <= tolerance
+        new_vec[b]=-new_vec[b]
+        return new_vec
+    
+ 
+    def line_in_plane(self, g):
+        '''checks wether a line g is in a plane '''
+        res = True
+        prod = np.dot(g.dir, self.norm)
+        if prod != 0 or self.point_in_plane(g.origin) != True:
+            res = False
+        return res
+              
+    def point_in_plane(self, point):
+        '''checks wether a point is in a plane  '''
+        distance = np.dot(self.norm, self.o - point)
+        if abs(distance) < self.tolerance:
+            return True
+        else:
+            return False
+
+    def __str__(self):
+        return str(self.norm) + str(self.a) + str(self.o) 
+
+
+class Scintillator:
+
+    def __init__(self, l, det):
+        self.scin=[]
+        pl=l[0]
+        det_ind=l[1]
+        self.det = det
+        for i in range(len(pl)):
+            self.scin.append(Plane(vector(pl[i,0,0],pl[i,0,1],pl[i,0,2]),vector(pl[i,1,0],pl[i,1,1],pl[i,1,2]),vector(pl[i,2,0],pl[i,2,1],pl[i,2,2]),[]))
+        for i in range(len(det_ind)):
+            self.scin[det_ind[i]].det = self.det
+        self.len = len(pl)
+
+    def real_hit(self, g, i, hit_rec, tr, sc, path):
+        '''finds all cross points that lie within the detector, and return the 
+        selection array(slc), the rays objects with all rays(g) and the real crosspoints(cp)'''
+        cp,slc = self.scin[i].cross(g)
+        if not np.any(slc):
+            return slc, g, cp, hit_rec
+        else:
+            slc_n = np.ones(len(cp), dtype=bool)
+            for j in range(self.len):
+                if j != i:
+                    nr = np.empty((len(cp), 1, 3))
+                    nr[...,0:3] = self.scin[j].norm
+                    b = ((nr @ (cp - self.scin[j].o)).flatten()) >= tolerance 
+                    slc_n = np.logical_and(slc_n, b)
+            nr[...,0:3] = self.scin[i].norm
+            c = (nr @ g.dir[slc]).flatten() <= 0
+            slc_n = np.logical_and(slc_n,c) 
+            cp = cp[slc_n]
+            slc[np.where(slc)] = slc_n
+            if i == 1:
+                hit_rec = Scintillator.det_hit(self, g, 1, cp, tr, sc, path, hit_rec, slc)
+            if i == 4 :
+                hit_rec = Scintillator.det_hit(self, g, 4, cp, tr, sc, path, hit_rec, slc)
+            return slc, g, cp, hit_rec
+        
+    def det_hit(self, g, i, cp, tr, sc, path, hit_rec, slc):
+        '''finds out wether or not a ray which hit a plane also hits the detector
+        writes hits in the hit_rec array'''
+        slc_dummy = np.ones(len(slc), dtype=bool)
+        slc_dummy = np.logical_and(slc_dummy, slc)
+        slc_1 = np.ones(len(cp), dtype=bool)
+        for j in range(len(self.scin[i].det)):
+            nr = np.empty((len(cp), 1, 3))
+            nr[...,0:3] = ((self.scin[i].det)[j]).norm
+            res = ((nr @ (cp-(((self.scin[i].det)[j]).o))).flatten()) >= tolerance
+            slc_1 = np.logical_and(slc_1, res)
+        slc_dummy[np.where(slc_dummy)] = slc_1
+        l = np.array((hit_rec[slc_dummy,0,0]).astype(int))
+        slc_l = l < 3
+        l = l[slc_l]
+        slc_dummy[np.where(slc_dummy)] = slc_l
+        if np.any(slc_dummy):
+            hit_rec[slc_dummy, :, 0] += 1
+            if i == 1:
+                hit_rec[slc_dummy, l, 1] = 1
+            else:
+                hit_rec[slc_dummy, l, 1] = 2
+            hit_rec[slc_dummy, l, 2] = tr[slc_dummy]
+            hit_rec[slc_dummy, l, 3] = sc[slc_dummy]
+            hit_rec[slc_dummy, l, 4] = path[slc_dummy]
+            hit_rec[slc_dummy, l, 5] = g.origin[slc_dummy,0].T
+            hit_rec[slc_dummy, l, 6] = g.origin[slc_dummy,1].T
+            hit_rec[slc_dummy, l, 7] = g.origin[slc_dummy,2].T
+        return hit_rec
+
+
+    def next_ray(self, g, i, cp, tr, sc, path, slc, Ind_S):
+        '''calculates the reflected and scattered rays that have cross points
+        with the plane i, alter the original ray beam and 
+        give back arrays the size of the beam ray with the amount of scattering/reflection
+        events and the path length of the rays'''
+        an = Plane.cross_angle(self.scin[i], g, slc)
+        if np.any(an):
+            r = an <= ref_angle(Ind_S,1)
+            s = an > ref_angle(Ind_S,1)
+            r_dummy = r[np.where(~r)]  #following lines are for reflected rays though angle of scattering
+            b = int(len(r_dummy) * ((Ind_S - 1) / (Ind_S + 1)**2))
+            u = np.array( np.random.choice(len(r_dummy), size=b, replace = False))
+            r_dummy[u] = True
+            r[np.where(~r)] = r_dummy
+            s[np.where(s)] = ~r_dummy
+            slc_reflect = np.zeros(len(slc), dtype=bool)
+            slc_reflect[slc] = True
+            if np.any(r):
+                slc_reflect[np.where(slc_reflect)] = r
+                tr[slc_reflect] += 1
+                ref = Plane.reflect(self.scin[i], g, slc_reflect)
+            slc_scatter = np.zeros(len(slc), dtype=bool)
+            slc_scatter[slc] = True
+            if np.any(s):
+                slc_scatter[np.where(slc_scatter)] = s
+                sc[slc_scatter] += 1
+                scat = Plane.scatter(self.scin[i], g, Ind_S, 1, slc_scatter)
+            path[slc] = path[slc] + (np.linalg.norm(g.origin[slc] - cp, axis=1))[:,0]
+            if np.any(s):
+                g.dir[slc_scatter] = scat
+                g.origin[slc_scatter] = cp[s]
+            if np.any(r):
+                g.dir[slc_reflect] = ref
+                g.origin[slc_reflect] = cp[r]
+        return g, tr, sc, path
+
+
+def run(pls, rys, stop, Ind_S):
+    '''main program for simulation, sc,tr,path are arrays including the num of reflection, 
+    scattering and the path length. Hit_rec includes for each ray: [num of detector hits,
+    det number,num of ref, num of scat,path_length] each for the moment the det hit happened.'''
+    start = time.time()
+    sc = np.zeros(len(rys.dir), dtype = np.int32)
+    tr = np.zeros(len(rys.dir), dtype = np.int32)
+    hit_rec = np.zeros((len(rys.dir), 3, 8), dtype = np.float32)  
+    path = np.zeros(len(rys.dir), dtype = np.float32)
+    while (not np.any(path > stop)):
+        for i in range(len(pls.scin)):
+            slc, rys, cp, hit_rec = Scintillator.real_hit(pls, rys, i, hit_rec, tr, sc, path)
+            rys, tr, sc, path = Scintillator.next_ray(pls, rys, i, cp, tr, sc, path, slc, Ind_S)
+    end=time.time()
+    return sc,tr,path,hit_rec
+
+# dir = np.ones((6,3,1))
+# print(dir)
+# origin = np.ones((6,3,1))
+# g = Rays_custom(origin, dir)
+# print(g.dir)
+# g=Rays(200,vector(0,0,0))
+# det=[Plane(vector(-10,-35,5),vector(-10,-35,-5),vector(-10,-30,0),[]),
+#        Plane(vector(-10,-35,-5),vector(10,-45,-5),vector(0,-35,-5),[]),
+#        Plane(vector(10,-35,5),vector(10,-35,0),vector(10,-45,-5),[]),
+#        Plane(vector(-10,-35,5),vector(0,-35,5),vector(10,-45,5),[])]
+# o = regular_hexagonal_prism([20, 51.962],[1,4])
+# p = Scintillator(o, det)
+
+# sc, tr, path, hit_rec = run(p,g,5000,2.15)
+# hit_rec = np.empty((200,3,5))
+# cp,slc = Plane.cross(p.scin[7],g)
+# tr = np.ones(200)
+# sc = np.ones(200)
+# path = np.empty(200)
+# slc, g, cp, hit_rec = Scintillator.real_hit(p,g,7,hit_rec,tr,sc,path)
+# print(slc)
+# g, tr, sc, path = Scintillator.next_ray(p,g,7,cp,tr,sc,path,slc)
+
+#p=Scintillator([Plane(vector(0,0,0),vector(0,-25.981,45),vector(20,0,0),[]),
+#                Plane(vector(0,0,0),vector(20,0,0),vector(0,51.962,0),det),
+#                Plane(vector(0,51.962,0),vector(20,51.962,0),vector(0,77.943,45),[]),
+#                Plane(vector(0,51.962,90),vector(0,77.943,45),vector(20,51.962,90),[]),
+#                Plane(vector(0,0,90),vector(0,51.962,90),vector(20,0,90),det),
+#                Plane(vector(0,0,90),vector(20,0,90),vector(0,-25.981,45),[]),
+#                Plane(vector(0,0,0),vector(0,51.962,0),vector(0,0,90),[]),
+#                Plane(vector(20,0,0),vector(20,0,90),vector(20,51.962,0),[])])

Visualize.gp

@@ -0,0 +1,29 @@
+#program to visualize the BGO, using the planes norm equation
+#Plane 1: 900*y+519.62*z=0
+#Plane 2: z=0
+#Plane 3: -900*y+519.62*z=-46765.8
+#Plane 4: -900*y-519.62*z=-93531.6
+#Plane 5: -1039.24*z=-93531.6
+#Plane 6: 900*y-519.62*z=-46765.8
+#Plane 7: x=0
+#Plane 8: -4676.58*x=-93531.6
+#[  0.   900.   519.62]0.0
+#[   0.      0.   1039.24]0.0
+#[   0.   -900.    519.62]-46765.8
+#[   0.   -900.   -519.62]-93531.59999999999
+#[    0.       0.   -1039.24]-93531.6
+#[   0.    900.   -519.62]-46765.8
+#[4676.58    0.      0.  ]0.0
+#[-4676.58     0.       0.  ]-93531.6
+
+
+#splot -(900/519.62)*y title "Plane 1", 0 title "Plane 2", (900/519.62)*y-(46765.8/519.62) title "Plane 3", -(900/519.62)*y+(93531.6/519.62) title "Plane 4", (93531.6/1039.24) title "Plane 5", (900/519.62)*y+(46765.8/519.62) title "Plane 6", '/Users/clarapittschellis/Desktop/Light/Light/journey.dat' using 1:2:3:4:5:6 with vectors lw 3
+splot '/Users/clarapittschellis/Desktop/Light/Light/journey.dat' using 1:2:3:4:5:6 with vectors lw 3
+
+set title "BGO"
+set xlabel "X"
+set ylabel "Y"
+set zlabel "Z"
+set xrange[-10:20]
+set yrange[-40:80]
+set zrange[0:100]

Visualize.py

@@ -0,0 +1,64 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from Struktur import Rays, Plane, Ray
+from mpl_toolkits.mplot3d import Axes3D
+
+#vektor = list(map(float, input("Bitte geben Sie die drei Zahlen für den Vektor ein (getrennt durch Leerzeichen): ").split()))
+
+def vis_vector(vektor):
+
+    #if len(vektor) != 3:
+        #print("Bitte geben Sie genau drei Zahlen ein.")
+    #else:
+        #print("Der erstellte Vektor lautet:", vektor)
+
+    ray=Rays(vektor, 1)
+    vec=ray.get_data()
+    ax = plt.figure().add_subplot(projection='3d')
+    u=[]
+    v=[]
+    w=[]
+    for i in range(len(vec)-1):
+        u.append(vec[i][0])
+        v.append(vec[i][1])
+        w.append(vec[i][2])
+    ax.quiver(vektor[0],vektor[1],vektor[2],u,v,w, length=3, normalize=True)
+    ax.set_xlim([-5, 5])
+    ax.set_ylim([-5, 5])
+    ax.set_zlim([-5, 5])
+ 
+
+
+def vis_plane(plane):
+
+    # Define a plane in the form ax + by + cz + d = 0
+    a=plane.norm[0]
+    b=plane.norm[1]
+    c=plane.norm[2]
+    d=plane.a
+
+    # Generate grid points for x and y
+    x = np.linspace(-10, 10, 100)
+    y = np.linspace(-10, 10, 100)
+    x, y = np.meshgrid(x, y)
+
+    # Calculate corresponding z values for the plane
+    z = (-a * x - b * y - d) / c
+
+    # Create a 3D plot
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+
+    # Plot the surface
+    ax.plot_surface(x, y, z, alpha=0.5)
+
+    # Set labels and title
+    ax.set_xlabel('X')
+    ax.set_ylabel('Y')
+    ax.set_zlabel('Z')
+    ax.set_title('Plane in 3D')
+
+
+
+
+

ahepam-het.gf

@@ -1,82 +0,0 @@
-#! /usr/bin/env ./geometryfactor.py
-# AHEPAM w/ HETA/HETB detectors.
-
-φ1=0
-φ2=π/12
-Δr=1
-Δθ=1°
-radius=65
-
-define=a=4
-define=a₁=${a}/sqrt(3)
-define=a₂=${a}*tan(π/12)
-define=r₁=2/sqrt(3)*${a}
-define=b=8.7
-define=b₁=${b}*cos(π/12)
-define=b₂=${b}*sin(π/12)
-define=c=18.25
-define=c₀=${c}/sqrt(2)
-define=c₁=${c}*cos(π/6)
-define=c₂=${c}*sin(π/6)
-define=c₃=${c}*cos(π/12)
-define=c₄=${c}*sin(π/12)
-define=d=30.0
-define=dd=50.0
-define=w=1.0
-define=v=(${d}-10)
-define=e=39.0
-
-#radius=sqrt((2*${d}+${w})**2 + 2*${b}**2)
-
-# center segment
-prism= 24, ${a}, ${w}/2
-
-# ring segment
-prism= 24, ${b}, ${w}/2
-
-# outer segment
-prism= 24, ${c}, ${w}/2
-
-copy=0
-move=0,0,${d}
-copy=1
-move=0,0,${d}
-copy=2
-move=0,0,${d}
-
-copy=0
-move=0,0,-${d}
-copy=1
-move=0,0,-${d}
-copy=2
-move=0,0,-${d}
-
-copy=0
-move=0,0, ${d}+${dd}
-copy=1
-move=0,0, ${d}+${dd}
-copy=2
-move=0,0, ${d}+${dd}
-
-copy=0
-move=0,0, -${d}-${dd}
-copy=1
-move=0,0, -${d}-${dd}
-copy=2
-move=0,0, -${d}-${dd}
-
-# BGO
-
-prism= 6, ${e}, ${v}/2
-move=0,0,${d}/2
-
-copy=-1
-move=0,0,-${d}
-
-# Cherenkov
-
-box= 31, 31, 20
-move= 0,0, ${d}+${dd}/2
-
-copy=-1
-move= 0,0, -2*${d}-${dd}

ahepam.gf

@@ -1,159 +0,0 @@
-# AHEPAM honey
-
-φ1=0
-φ2=π/6
-Δr=1
-Δθ=1°
-radius=65
-
-define=a=8.0
-define=a₁=${a}/sqrt(3)
-define=r₁=2/sqrt(3)*${a}
-define=b=${a}*sqrt(21/6/sqrt(3)/tan(π/12))
-define=b₁=${b}*cos(π/6)
-define=b₂=${b}*sin(π/6)
-define=Δc=2
-define=c=(2*${b}-${a}+${Δc})
-define=c₁=${c}*cos(π/6)
-define=c₂=${c}*sin(π/6)
-define=d=30.0
-define=w=1.0
-define=v=(${d}-10)
-define=e=(${c}-${a})
-
-#radius=sqrt((2*${d}+${w})**2 + 2*${b}**2)
-
-# center segment
-detector
-plane=    0, 0,-${w}/2;   1., 0, 0;   0, 1., 0
-plane=    0, 0, ${w}/2;   0, 1., 0;   1., 0, 0
-plane=   -${a}, 0, 0;   0, 1., 0;   0, 0,  1.
-plane=    ${a}, 0, 0;   0, 0,  1.;  0, 1., 0
-plane=  0, -${r₁}, 0;  0, 0,  1.;  ${a},  ${a₁}, 0
-plane=  0, -${r₁}, 0;  0, 0,  1.;  ${a}, -${a₁}, 0
-plane=  0,  ${r₁}, 0;  ${a},  ${a₁}, 0;  0, 0,  1.
-plane=  0,  ${r₁}, 0;  ${a}, -${a₁}, 0;  0, 0,  1.
-
-# hex segment
-detector
-plane=    0, 0,-${w}/2;   1., 0, 0;   0, 1., 0
-plane=    0, 0, ${w}/2;   0, 1., 0;   1., 0, 0
-plane=    ${a}, 0, 0;   0, 1., 0;   0, 0,  1.
-plane=    ${a}, -${a₁}, 0;  0, 0,  1.;  ${a}, -${a₁}, 0
-plane=    ${a},  ${a₁}, 0;  ${a},  ${a₁}, 0;  0, 0,  1.
-plane=    ${b}, 0, 0;   0, 0,  1.;  0, 1., 0
-plane=    ${b₁}, -${b₂}, 0;  0, 0,  1.;  ${b₂}, ${b₁}, 0
-plane=    ${b₁},  ${b₂}, 0;  ${b₂}, -${b₁}, 0;  0, 0,  1.
-
-copy=-1
-rz=π/3
-copy=-2
-rz=2*π/3
-copy=-3
-rz=π
-copy=-4
-rz=4*π/3
-copy=-5
-rz=5*π/3
-
-# outer segment
-detector
-plane=   0, 0,-${w}/2;   1., 0, 0;   0, 1., 0
-plane=   0, 0, ${w}/2;   0, 1., 0;   1., 0, 0
-plane=  -${c}, 0, 0;   0, 1., 0;   0, 0,  1.
-plane=   ${c}, 0, 0;   0, 0, 1.;   0, 1., 0
-plane=  0, -${c}, 0;   0, 0, 1.;   1., 0, 0
-plane=  0,  ${c}, 0;   1., 0, 0;   0, 0,  1.
-
-plane=  -${c₁},  ${c₂}, 0;   ${c₂}, ${c₁}, 0;   0, 0,  1.
-plane=  -${c₁}, -${c₂}, 0;  -${c₂}, ${c₁}, 0;   0, 0,  1.
-plane=   ${c₁},  ${c₂}, 0;   0, 0, 1.;  -${c₂}, ${c₁}, 0
-plane=   ${c₁}, -${c₂}, 0;   0, 0, 1.;   ${c₂}, ${c₁}, 0
-
-plane=  ${c₂}, -${c₁}, 0;   0, 0, 1.;   ${c₁},  ${c₂}, 0
-plane= -${c₂}, -${c₁}, 0;   0, 0, 1.;   ${c₁}, -${c₂}, 0
-plane=  ${c₂},  ${c₁}, 0;   ${c₁}, -${c₂}, 0;   0, 0,  1.
-plane= -${c₂},  ${c₁}, 0;   ${c₁},  ${c₂}, 0;   0, 0,  1.
-
-
-copy=-8
-move=0,0,${d}
-copy=-8
-move=0,0,${d}
-copy=-8
-move=0,0,${d}
-copy=-8
-move=0,0,${d}
-copy=-8
-move=0,0,${d}
-copy=-8
-move=0,0,${d}
-copy=-8
-move=0,0,${d}
-copy=-8
-move=0,0,${d}
-
-# define=Ar=π/6
-define=Ar=0
-copy=-16
-move=0,0,-${d}
-rz=${Ar}
-copy=-16
-move=0,0,-${d}
-rz=${Ar}
-copy=-16
-move=0,0,-${d}
-rz=${Ar}
-copy=-16
-move=0,0,-${d}
-rz=${Ar}
-copy=-16
-move=0,0,-${d}
-rz=${Ar}
-copy=-16
-move=0,0,-${d}
-rz=${Ar}
-copy=-16
-move=0,0,-${d}
-rz=${Ar}
-copy=-16
-move=0,0,-${d}
-rz=${Ar}
-
-copy=0
-move=0,0, 2*${d}
-copy=1
-move=0,0, 2*${d}
-copy=2
-move=0,0, 2*${d}
-copy=3
-move=0,0, 2*${d}
-copy=4
-move=0,0, 2*${d}
-copy=5
-move=0,0, 2*${d}
-copy=6
-move=0,0, 2*${d}
-
-copy=0
-move=0,0, -2*${d}
-copy=1
-move=0,0, -2*${d}
-copy=2
-move=0,0, -2*${d}
-copy=3
-move=0,0, -2*${d}
-copy=4
-move=0,0, -2*${d}
-copy=5
-move=0,0, -2*${d}
-copy=6
-move=0,0, -2*${d}
-
-# BGO
-copy=7
-scale= ${e}/${c}, ${e}/${c}, ${v}/${w}
-move=0,0,${d}/2
-
-copy=-1
-move=0,0,-${d}

ahepam.makefile

@@ -1,92 +0,0 @@
-#! /usr/bin/make -f
-
-RUN=ahepam-1mm_1°
-CUTS=H A B AB BC AC  AB1 AnBn  AB1~C AnBn~C  BC1 A1C BnCn AnCm  BC1~D A1C~D BnCn~D AnCm~D
-
-all: $(patsubst %, $(RUN)_%.ssv, $(CUTS))
-
-# Ray C D B E A G F
-# 4   8 8 8 7 7 1 1
-
-# A 36 31  C H H H H H H
-# B 21 16  C H H H H H H R
-# C  5  0  C H H H H H H R
-# D 13  8  C H H H H H H R
-# E 29 24  C H H H H H H
-# G 43 38  BGO
-# F 44 39  BGO
-
-# exactly one ' ' between path lengths
-# two space between ray path length
-# four columns 1..4 for the ray 
-
-# R segments are never part of a coincidence,
-#   but always part of an anticoincidence.
-
-%_H.ssv: %.ssv
-	egrep -v ' ( 0){40}$$' $< > $@
-
-%_A.ssv: %_H.ssv
-	egrep -v '  ([-+.e0-9]+ ){31}(0 ){7}' $< > $@
-
-%_AB.ssv: %_A.ssv
-	egrep -v '  ([-+.e0-9]+ ){16}(0 ){7}' $< > $@
-
-%_B.ssv: %_H.ssv
-	egrep -v '  ([-+.e0-9]+ ){16}(0 ){7}' $< > $@
-
-%_BC.ssv: %_B.ssv
-	egrep -v '  (0 ){7}' $< > $@
-
-%_AC.ssv: %_A.ssv
-	egrep -v '  (0 ){7}' $< > $@
-
-# Stopping in BGO1
-
-%_AB1.ssv: %_AB.ssv
-	egrep -v '  ([-+.e0-9]+ ){16}0 ' $< > $@
-
-%_AnBn.ssv: %_AB.ssv
-	egrep    '  ([-+.e0-9]+ ){16}0 ' $< \
-      | awk '$$37&&$$22||$$38&&$$23||$$39&&$$24||$$40&&$$25||$$41&&$$26||$$42&&$$27' \
-      > $@
-
-%_AB1~C.ssv: %_AB1.ssv
-	-egrep '  (0 ){8}' $< > $@
-
-%_AnBn~C.ssv: %_AnBn.ssv
-	-egrep '  (0 ){8}' $< > $@
-
-# Stopping in BGO2
-
-%_BC1.ssv: %_BC.ssv
-	egrep -v '  0 ' $< > $@
-
-%_A1C.ssv: %_BC.ssv
-	egrep    '  0 ' $< \
-      | egrep -v '  ([-+.e0-9]+ ){31}0 ' \
-      > $@
-
-%_BC~A1~C1.ssv: %_BC.ssv
-	egrep '  0 ([-+.e0-9]+ ){30}0 ' $< > $@
-
-%_BnCn.ssv: %_BC.ssv
-	egrep '  0 ([-+.e0-9]+ ){30}0 ' $< \
-      | awk '$$6&&$$22||$$7&&$$23||$$8&&$$24||$$9&&$$25||$$10&&$$26||$$11&&$$27' \
-      > $@
-
-%_AnCm.ssv: %_BC.ssv
-	egrep '  0 ([-+.e0-9]+ ){30}0 ' $< \
-      | awk   '$$6&&$$38||$$7&&$$39||$$8&&$$40||$$9&&$$41||$$10&&$$42||$$11&&$$37 \
-             ||$$6&&$$42||$$7&&$$37||$$8&&$$38||$$9&&$$39||$$10&&$$40||$$11&&$$41' \
-      | awk '!($$6&&$$22||$$7&&$$23||$$8&&$$24||$$9&&$$25||$$10&&$$26||$$11&&$$27)' \
-      > $@
-
-%_BC1~D.ssv: %_BC1.ssv
-	-egrep '  ([-+.e0-9]+ ){8}(0 ){8}' $< > $@
-%_A1C~D.ssv: %_A1C.ssv
-	-egrep '  ([-+.e0-9]+ ){8}(0 ){8}' $< > $@
-%_BnCn~D.ssv: %_BnCn.ssv
-	-egrep '  ([-+.e0-9]+ ){8}(0 ){8}' $< > $@
-%_AnCm~D.ssv: %_AnCm.ssv
-	-egrep '  ([-+.e0-9]+ ){8}(0 ){8}' $< > $@

ahepam24.gf

@@ -1,237 +0,0 @@
-#! /usr/bin/env ./geometryfactor.py
-# AHEPAM honey
-
-φ1=0
-φ2=π/12
-Δr=1
-Δθ=1°
-radius=65
-
-define=a=8.0
-define=a₁=${a}/sqrt(3)
-define=a₂=${a}*tan(π/12)
-define=r₁=2/sqrt(3)*${a}
-define=b=${a}*sqrt(21/6/sqrt(3)/tan(π/12))
-define=b₁=${b}*cos(π/12)
-define=b₂=${b}*sin(π/12)
-define=Δc=5
-define=c=(2*${b}-${a}+${Δc})
-define=c₀=${c}/sqrt(2)
-define=c₁=${c}*cos(π/6)
-define=c₂=${c}*sin(π/6)
-define=c₃=${c}*cos(π/12)
-define=c₄=${c}*sin(π/12)
-define=d=30.0
-define=w=1.0
-define=v=(${d}-10)
-#define=e=(${c}-${a})
-define=e=${c}
-
-#radius=sqrt((2*${d}+${w})**2 + 2*${b}**2)
-
-# center segment
-detector
-plane=    0, 0,-${w}/2;   1., 0, 0;   0, 1., 0
-plane=    0, 0, ${w}/2;   0, 1., 0;   1., 0, 0
-plane=   -${a}, 0, 0;   0, 1., 0;   0, 0,  1.
-plane=    ${a}, 0, 0;   0, 0,  1.;  0, 1., 0
-plane=   0,  ${a}, 0;   1., 0, 0;   0, 0,  1.
-plane=   0, -${a}, 0;   0, 0,  1.;  1., 0, 0
-plane=  0, -${r₁}, 0;  0, 0,  1.;  ${a},  ${a₁}, 0
-plane=  0, -${r₁}, 0;  0, 0,  1.;  ${a}, -${a₁}, 0
-plane=  0,  ${r₁}, 0;  ${a},  ${a₁}, 0;  0, 0,  1.
-plane=  0,  ${r₁}, 0;  ${a}, -${a₁}, 0;  0, 0,  1.
-
-plane=  -${r₁}, 0, 0;   ${a₁}, ${a}, 0;  0, 0,  1.
-plane=  -${r₁}, 0, 0;  -${a₁}, ${a}, 0;  0, 0,  1.
-plane=   ${r₁}, 0, 0;  0, 0,  1.;  ${a₁},  ${a}, 0
-plane=   ${r₁}, 0, 0;  0, 0,  1.;  -${a₁}, ${a}, 0
-
-# hex segment
-detector
-plane=    0, 0,-${w}/2;   1., 0, 0;   0, 1., 0
-plane=    0, 0, ${w}/2;   0, 1., 0;   1., 0, 0
-plane=    ${a}, 0, 0;   0, 1., 0;   0, 0,  1.
-plane=    ${b}, 0, 0;   0, 0,  1.;  0, 1., 0
-plane=    ${a}, -${a₂}, 0;  0, 0,  1.;  ${a}, -${a₂}, 0
-plane=    ${a},  ${a₂}, 0;  ${a},  ${a₂}, 0;  0, 0,  1.
-plane=    ${b₁}, -${b₂}, 0;  0, 0,  1.;  ${b₂}, ${b₁}, 0
-plane=    ${b₁},  ${b₂}, 0;  ${b₂}, -${b₁}, 0;  0, 0,  1.
-
-copy=-1
-rz=π/6
-copy=-2
-rz=π/3
-copy=-3
-rz=π/2
-copy=-4
-rz=2*π/3
-copy=-5
-rz=5*π/6
-copy=-6
-rz=π
-copy=-7
-rz=7*π/6
-copy=-8
-rz=4*π/3
-copy=-9
-rz=3*π/2
-copy=-10
-rz=5*π/3
-copy=-11
-rz=11*π/6
-
-# outer segment
-detector
-plane=   0, 0,-${w}/2;   1., 0, 0;   0, 1., 0
-plane=   0, 0, ${w}/2;   0, 1., 0;   1., 0, 0
-plane=  -${c}, 0, 0;   0, 1., 0;   0, 0,  1.
-plane=   ${c}, 0, 0;   0, 0, 1.;   0, 1., 0
-plane=  0, -${c}, 0;   0, 0, 1.;   1., 0, 0
-plane=  0,  ${c}, 0;   1., 0, 0;   0, 0,  1.
-
-plane=  -${c₀},  ${c₀}, 0;   1., 1., 0;   0, 0,  1.
-plane=  -${c₀}, -${c₀}, 0;  -1., 1., 0;   0, 0,  1.
-plane=   ${c₀},  ${c₀}, 0;   0, 0,  1.;  -1., 1., 0
-plane=   ${c₀}, -${c₀}, 0;   0, 0,  1.;   1., 1., 0
-
-plane=  -${c₁},  ${c₂}, 0;   ${c₂}, ${c₁}, 0;   0, 0,  1.
-plane=  -${c₁}, -${c₂}, 0;  -${c₂}, ${c₁}, 0;   0, 0,  1.
-plane=   ${c₁},  ${c₂}, 0;   0, 0, 1.;  -${c₂}, ${c₁}, 0
-plane=   ${c₁}, -${c₂}, 0;   0, 0, 1.;   ${c₂}, ${c₁}, 0
-
-plane=  -${c₃},  ${c₄}, 0;   ${c₄}, ${c₃}, 0;   0, 0,  1.
-plane=  -${c₃}, -${c₄}, 0;  -${c₄}, ${c₃}, 0;   0, 0,  1.
-plane=   ${c₃},  ${c₄}, 0;   0, 0, 1.;  -${c₄}, ${c₃}, 0
-plane=   ${c₃}, -${c₄}, 0;   0, 0, 1.;   ${c₄}, ${c₃}, 0
-
-plane=  ${c₂}, -${c₁}, 0;   0, 0, 1.;   ${c₁},  ${c₂}, 0
-plane= -${c₂}, -${c₁}, 0;   0, 0, 1.;   ${c₁}, -${c₂}, 0
-plane=  ${c₂},  ${c₁}, 0;   ${c₁}, -${c₂}, 0;   0, 0,  1.
-plane= -${c₂},  ${c₁}, 0;   ${c₁},  ${c₂}, 0;   0, 0,  1.
-
-plane=  ${c₄}, -${c₃}, 0;   0, 0, 1.;   ${c₃},  ${c₄}, 0
-plane= -${c₄}, -${c₃}, 0;   0, 0, 1.;   ${c₃}, -${c₄}, 0
-plane=  ${c₄},  ${c₃}, 0;   ${c₃}, -${c₄}, 0;   0, 0,  1.
-plane= -${c₄},  ${c₃}, 0;   ${c₃},  ${c₄}, 0;   0, 0,  1.
-
-
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-copy=-14
-move=0,0,${d}
-
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-copy=-28
-move=0,0,-${d}
-
-copy=0
-move=0,0, 2*${d}
-copy=1
-move=0,0, 2*${d}
-copy=2
-move=0,0, 2*${d}
-copy=3
-move=0,0, 2*${d}
-copy=4
-move=0,0, 2*${d}
-copy=5
-move=0,0, 2*${d}
-copy=6
-move=0,0, 2*${d}
-copy=7
-move=0,0, 2*${d}
-copy=8
-move=0,0, 2*${d}
-copy=9
-move=0,0, 2*${d}
-copy=10
-move=0,0, 2*${d}
-copy=11
-move=0,0, 2*${d}
-copy=12
-move=0,0, 2*${d}
-
-copy=0
-move=0,0, -2*${d}
-copy=1
-move=0,0, -2*${d}
-copy=2
-move=0,0, -2*${d}
-copy=3
-move=0,0, -2*${d}
-copy=4
-move=0,0, -2*${d}
-copy=5
-move=0,0, -2*${d}
-copy=6
-move=0,0, -2*${d}
-copy=7
-move=0,0, -2*${d}
-copy=8
-move=0,0, -2*${d}
-copy=9
-move=0,0, -2*${d}
-copy=10
-move=0,0, -2*${d}
-copy=11
-move=0,0, -2*${d}
-copy=12
-move=0,0, -2*${d}
-
-# BGO
-copy=13
-scale= ${e}/${c}, ${e}/${c}, ${v}/${w}
-move=0,0,${d}/2
-
-copy=-1
-move=0,0,-${d}

ahepam24.makefile

@@ -1,347 +0,0 @@
-#! /usr/bin/make -f
-
-RUN=ahepam-1mm_1°
-CUTS=H A B AB BC AC BD \
-     AB1 AnBn AnBm AnBl  AB1~C AnBn~C AnBm~C AnBl~C \
-     BC1 A1C BnCn BnCm AnCl  BC1~D A1C~D BnCn~D BnCm~D AnCl~D \
-     BcDc BcDn BnD0 BnD1 BnD2 BnD3 BnD4 BnD5 BnD6
-
-all: $(patsubst %, $(RUN)_%.ssv, $(CUTS))
-
-# Ray C  D  B  E  A   G  F
-# 4   14 14 14 13 13  1  1
-# 1    5 19 33 47 60 73 74
-# x    0 14 28 42 55 68 69
-
-# A 60 55  C H H H H H H H H H H H H
-# B 33 28  C H H H H H H H H H H H H R
-# C  5  0  C H H H H H H H H H H H H R
-# D 19 14  C H H H H H H H H H H H H R
-# E 47 42  C H H H H H H H H H H H H
-# G 73 68  BGO
-# F 74 69  BGO
-
-# exactly one ' ' between path lengths
-# two space between ray path length
-# four columns 1..4 for the ray 
-
-# R segments are never part of a coincidence,
-#   but always part of an anticoincidence.
-
-%_H.ssv: %.ssv
-	egrep -v ' ( 0){70}$$' $< > $@
-
-%_A.ssv: %_H.ssv
-	egrep -v '  ([-+.e0-9]+ ){55}(0 ){13}' $< > $@
-
-%_AB.ssv: %_A.ssv
-	egrep -v '  ([-+.e0-9]+ ){28}(0 ){13}' $< > $@
-
-%_B.ssv: %_H.ssv
-	egrep -v '  ([-+.e0-9]+ ){28}(0 ){13}' $< > $@
-
-%_BC.ssv: %_B.ssv
-	egrep -v '  (0 ){13}' $< > $@
-
-%_AC.ssv: %_A.ssv
-	egrep -v '  (0 ){13}' $< > $@
-
-# total penetrating, including outer ring.
-%_BD.ssv: %_H.ssv
-	egrep -v '  ([-+.e0-9]+ ){28}(0 ){14}' $< \
-      | egrep -v '  ([-+.e0-9]+ ){14}(0 ){14}' $< \
-      > $@
-
-# Stopping in BGO1
-
-%_AB1.ssv: %_AB.ssv
-	egrep -v '  ([-+.e0-9]+ ){28}0 ' $< > $@
-
-%_AnBn.ssv: %_AB.ssv
-	egrep    '  ([-+.e0-9]+ ){28}0 ' $< \
-      | awk ' $$61 && $$34 \
-           || $$62 && $$35 \
-           || $$63 && $$36 \
-           || $$64 && $$37 \
-           || $$65 && $$38 \
-           || $$66 && $$39 \
-           || $$67 && $$40 \
-           || $$68 && $$41 \
-           || $$69 && $$42 \
-           || $$70 && $$43 \
-           || $$71 && $$44 \
-           || $$72 && $$45' \
-      > $@
-
-%_AnBm.ssv: %_AB.ssv
-	egrep    '  ([-+.e0-9]+ ){28}0 ' $< \
-      | awk ' $$61 && $$34 \
-           || $$62 && $$35 \
-           || $$63 && $$36 \
-           || $$64 && $$37 \
-           || $$65 && $$38 \
-           || $$66 && $$39 \
-           || $$67 && $$40 \
-           || $$68 && $$41 \
-           || $$69 && $$42 \
-           || $$70 && $$43 \
-           || $$71 && $$44 \
-           || $$72 && $$45 \
-           {next};  \
-	      $$61 && $$45 \
-           || $$62 && $$34 \
-           || $$63 && $$35 \
-           || $$64 && $$36 \
-           || $$65 && $$37 \
-           || $$66 && $$38 \
-           || $$67 && $$39 \
-           || $$68 && $$40 \
-           || $$69 && $$41 \
-           || $$70 && $$42 \
-           || $$71 && $$43 \
-           || $$72 && $$44 \
-           || $$61 && $$35 \
-           || $$62 && $$36 \
-           || $$63 && $$37 \
-           || $$64 && $$38 \
-           || $$65 && $$39 \
-           || $$66 && $$40 \
-           || $$67 && $$41 \
-           || $$68 && $$42 \
-           || $$69 && $$43 \
-           || $$70 && $$44 \
-           || $$71 && $$45 \
-           || $$72 && $$34' \
-      > $@
-
-%_AnBl.ssv: %_AB.ssv
-	egrep    '  ([-+.e0-9]+ ){28}0 ' $< \
-      | awk ' $$61 && $$34 \
-           || $$62 && $$35 \
-           || $$63 && $$36 \
-           || $$64 && $$37 \
-           || $$65 && $$38 \
-           || $$66 && $$39 \
-           || $$67 && $$40 \
-           || $$68 && $$41 \
-           || $$69 && $$42 \
-           || $$70 && $$43 \
-           || $$71 && $$44 \
-           || $$72 && $$45 \
-	   || $$61 && $$45 \
-           || $$62 && $$34 \
-           || $$63 && $$35 \
-           || $$64 && $$36 \
-           || $$65 && $$37 \
-           || $$66 && $$38 \
-           || $$67 && $$39 \
-           || $$68 && $$40 \
-           || $$69 && $$41 \
-           || $$70 && $$42 \
-           || $$71 && $$43 \
-           || $$72 && $$44 \
-           || $$61 && $$35 \
-           || $$62 && $$36 \
-           || $$63 && $$37 \
-           || $$64 && $$38 \
-           || $$65 && $$39 \
-           || $$66 && $$40 \
-           || $$67 && $$41 \
-           || $$68 && $$42 \
-           || $$69 && $$43 \
-           || $$70 && $$44 \
-           || $$71 && $$45 \
-           || $$72 && $$34 \
-           {next};  \
-              $$61 && $$44 \
-           || $$62 && $$45 \
-           || $$63 && $$34 \
-           || $$64 && $$35 \
-           || $$65 && $$36 \
-           || $$66 && $$37 \
-           || $$67 && $$38 \
-           || $$68 && $$39 \
-           || $$69 && $$40 \
-           || $$70 && $$41 \
-           || $$71 && $$42 \
-           || $$72 && $$43 \
-           || $$61 && $$36 \
-           || $$62 && $$37 \
-           || $$63 && $$38 \
-           || $$64 && $$39 \
-           || $$65 && $$40 \
-           || $$66 && $$41 \
-           || $$67 && $$42 \
-           || $$68 && $$43 \
-           || $$69 && $$44 \
-           || $$70 && $$45 \
-           || $$71 && $$34 \
-           || $$72 && $$35' \
-      > $@
-
-%~C.ssv: %.ssv
-	-egrep '  (0 ){14}' $< > $@
-
-# Stopping in BGO2
-
-%_BC1.ssv: %_BC.ssv
-	egrep -v '  0 ' $< > $@
-
-%_A1C.ssv: %_BC.ssv
-	egrep    '  0 ' $< \
-      | egrep -v '  ([-+.e0-9]+ ){55}0 ' \
-      > $@
-
-%_BnCn.ssv: %_BC.ssv
-	egrep '  0 ([-+.e0-9]+ ){54}0 ' $< \
-      | awk ' $$6  && $$34 \
-           || $$7  && $$35 \
-           || $$8  && $$36 \
-           || $$9  && $$37 \
-           || $$10 && $$38 \
-           || $$11 && $$39 \
-           || $$12 && $$40 \
-           || $$13 && $$41 \
-           || $$14 && $$42 \
-           || $$15 && $$43 \
-           || $$16 && $$44 \
-           || $$17 && $$45' \
-      > $@
-
-%_BnCm.ssv: %_BC.ssv
-	egrep '  0 ([-+.e0-9]+ ){54}0 ' $< \
-      | awk ' $$6  && $$34 \
-           || $$7  && $$35 \
-           || $$8  && $$36 \
-           || $$9  && $$37 \
-           || $$10 && $$38 \
-           || $$11 && $$39 \
-           || $$12 && $$40 \
-           || $$13 && $$41 \
-           || $$14 && $$42 \
-           || $$15 && $$43 \
-           || $$16 && $$44 \
-           || $$17 && $$45 \
-	   {next}; \
-	      $$6  && $$45 \
-           || $$7  && $$34 \
-           || $$8  && $$35 \
-           || $$9  && $$36 \
-           || $$10 && $$37 \
-           || $$11 && $$38 \
-           || $$12 && $$39 \
-           || $$13 && $$40 \
-           || $$14 && $$41 \
-           || $$15 && $$42 \
-           || $$16 && $$43 \
-           || $$17 && $$44 \
-           || $$6  && $$35 \
-           || $$7  && $$36 \
-           || $$8  && $$37 \
-           || $$9  && $$38 \
-           || $$10 && $$39 \
-           || $$11 && $$40 \
-           || $$12 && $$41 \
-           || $$13 && $$42 \
-           || $$14 && $$43 \
-           || $$15 && $$44 \
-           || $$16 && $$45 \
-           || $$17 && $$34 ' \
-      > $@
-
-%_AnCl.ssv: %_AC.ssv
-	egrep '  0 ([-+.e0-9]+ ){54}0 ' $< \
-      | awk ' $$6  && $$34 \
-           || $$7  && $$35 \
-           || $$8  && $$36 \
-           || $$9  && $$37 \
-           || $$10 && $$38 \
-           || $$11 && $$39 \
-           || $$12 && $$40 \
-           || $$13 && $$41 \
-           || $$14 && $$42 \
-           || $$15 && $$43 \
-           || $$16 && $$44 \
-           || $$17 && $$45 \
-	   || $$6  && $$45 \
-           || $$7  && $$34 \
-           || $$8  && $$35 \
-           || $$9  && $$36 \
-           || $$10 && $$37 \
-           || $$11 && $$38 \
-           || $$12 && $$39 \
-           || $$13 && $$40 \
-           || $$14 && $$41 \
-           || $$15 && $$42 \
-           || $$16 && $$43 \
-           || $$17 && $$44 \
-           || $$6  && $$35 \
-           || $$7  && $$36 \
-           || $$8  && $$37 \
-           || $$9  && $$38 \
-           || $$10 && $$39 \
-           || $$11 && $$40 \
-           || $$12 && $$41 \
-           || $$13 && $$42 \
-           || $$14 && $$43 \
-           || $$15 && $$44 \
-           || $$16 && $$45 \
-           || $$17 && $$34 \
-	   {next}; \
-	      $$6  && $$63 \
-           || $$7  && $$64 \
-           || $$8  && $$65 \
-           || $$9  && $$66 \
-           || $$10 && $$67 \
-           || $$11 && $$68 \
-           || $$12 && $$69 \
-           || $$13 && $$70 \
-           || $$14 && $$71 \
-           || $$15 && $$72 \
-           || $$16 && $$61 \
-           || $$17 && $$62 \
-           || $$6  && $$71 \
-           || $$7  && $$72 \
-           || $$8  && $$61 \
-           || $$9  && $$62 \
-           || $$10 && $$63 \
-           || $$11 && $$64 \
-           || $$12 && $$65 \
-           || $$13 && $$66 \
-           || $$14 && $$67 \
-           || $$15 && $$68 \
-           || $$16 && $$69 \
-           || $$17 && $$70 ' \
-      > $@
-
-%~D.ssv: %.ssv
-	-egrep '  ([-+.e0-9]+ ){14}(0 ){14}' $< > $@
-
-%_BcDc.ssv: %_BC.ssv
-	egrep -v '  ([-+.e0-9]+ ){14}(0 )' $< \
-      | egrep -v '  ([-+.e0-9]+ ){28}(0 )' \
-      > $@ 
-%_BcDn.ssv: %_BC.ssv
-	egrep -v '  ([-+.e0-9]+ ){15}(0 ){12}' $< \
-      | egrep -v '  ([-+.e0-9]+ ){28}(0 )' \
-      > $@ 
-
-dN=0
-%_BnD$(dN).ssv: %_BC.ssv
-	awk '{for(i=0;i<12;i++) { \
-              if ($$(20+i) && $$(34+(i+N)%12)) { print; next } \
-              if ($$(20+i) && $$(34+(i+12-N)%12)) { print; next } \
-              }}' N=$(dN) $< > $@
-
-%_BnD1.ssv: %_BC.ssv$x
-	$(MAKE) -f $(word 1,$(MAKEFILE_LIST)) -o %.ssv dN=1 x=x $@
-%_BnD2.ssv: %_BC.ssv$x
-	$(MAKE) -f $(word 1,$(MAKEFILE_LIST)) -o %.ssv dN=2 x=x $@
-%_BnD3.ssv: %_BC.ssv$x
-	$(MAKE) -f $(word 1,$(MAKEFILE_LIST)) -o %.ssv dN=3 x=x $@
-%_BnD4.ssv: %_BC.ssv$x
-	$(MAKE) -f $(word 1,$(MAKEFILE_LIST)) -o %.ssv dN=4 x=x $@
-%_BnD5.ssv: %_BC.ssv$x
-	$(MAKE) -f $(word 1,$(MAKEFILE_LIST)) -o %.ssv dN=5 x=x $@
-%_BnD6.ssv: %_BC.ssv$x
-	$(MAKE) -f $(word 1,$(MAKEFILE_LIST)) -o %.ssv dN=6 x=x $@

det_signal.py

@@ -0,0 +1,138 @@
+import csv
+import sys
+import matplotlib.pyplot as plt
+import argparse
+import numpy as np
+import random
+
+def read_csv_1(file_path):
+    file = open(file_path, newline='')
+    unique_points = set()
+    for line in file:
+        row = line.strip().split(',')
+        if len(row) == 32:
+            point = (float(row[0]), float(row[1]), float(row[2]))
+            unique_points.add(point)
+    unique_points = list(unique_points)
+    unique_points = np.array(unique_points)
+    point_list = np.zeros((len(unique_points),7))
+    point_list[:, 0:3] = unique_points
+    point_list[:,3] = 0
+    point_list[:,4] = 0
+    return point_list
+
+def read_csv_2(file_path, slope_sc, slope_stop, Ind_A, Ind_S, point_list):
+    file = open(file_path, newline='')
+    ind = 0 
+    for line in file:
+        row = line.strip().split(',')
+        if len(row) == 32:
+            point = (float(row[0]), float(row[1]), float(row[2]))
+            matches = np.all(point_list[:, 0:3] == point, axis=1)
+            position = np.where(matches)[0]
+            if int(float(row[9])) == 1:
+                point_list[position,5] += 1
+            elif int(float(row[9])) == 2:
+                point_list[position,6] += 1
+            an = []
+            ta = np.arcsin(Ind_A/Ind_S)
+            for i in range(3):
+                x = float(row[13+8*i])
+                y = float(row[14+8*i])
+                z = float(row[15+8*i])
+                v=np.array((x,y,z))
+                norm = np.array((0,1,0))
+                an.append(np.arccos(abs(np.dot(v,norm))/(np.linalg.norm(norm)*np.linalg.norm(v))))
+            if int(float(row[8])) == 1 :
+                u = random.random()
+                v = random.random()
+                sc = np.exp(-slope_sc*(int(float(row[11]))))
+                p = np.exp(-slope_stop*(int(float(row[12]))))
+                if u <= sc and v <= p and an[0] < ta:
+                    if int(float(row[9])) == 1:
+                        point_list[position,3] += 1
+                    else:
+                        point_list[position,4] += 1
+            elif int(float(row[8])) == 2:
+                u = random.random()
+                u2 = random.random()
+                v = random.random()
+                v2 = random.random()
+                sc = np.exp(-slope_sc*(int(float(row[11]))))
+                sc2 = np.exp(-slope_sc*int(float(row[19])))
+                p = np.exp(-slope_stop*(int(float(row[12]))))
+                p2 = np.exp(-slope_stop*(int(float(row[20]))))
+                if u <= sc and v <= p and an[0] < ta:
+                    if int(float(row[9])) == 1:
+                        point_list[position,3] += 1
+                    else:
+                        point_list[position,4] += 1
+                elif u2 <= sc2 and v2 <= p2 and an[1] < ta:
+                    if int(float(row[17])) == 1:
+                        point_list[position,3] += 1
+                    else:
+                        point_list[position,4] += 1
+            elif int(float(row[8])) == 3:
+                u = random.random()
+                u2 = random.random()
+                u3 = random.random()
+                v = random.random()
+                v2 = random.random()
+                v3 = random.random()
+                sc = np.exp(-slope_sc*(int(float(row[11]))))
+                sc2 = np.exp(-slope_sc*int(float(row[19])))
+                sc3 = np.exp(-slope_sc*int(float(row[27])))
+                p = np.exp(-slope_stop*(int(float(row[12]))))
+                p2 = np.exp(-slope_stop*(int(float(row[20]))))
+                p3 = np.exp(-slope_stop*(int(float(row[28]))))
+                if u <= sc and v <= p and an[0] < ta:
+                    if int(float(row[9])) == 1:
+                        point_list[position,3] += 1
+                    else:
+                        point_list[position,4] += 1 
+                elif u2 <= sc2 and v2 <= p2 and an[1] < ta:
+                    if int(float(row[17])) == 1:
+                        point_list[position,3] += 1
+                    else:
+                        point_list[position,4] += 1 
+                elif u3 <= sc3 and v3 <= p3 and an[2] < ta:
+                    if int(float(row[25])) == 1:
+                        point_list[position,3] += 1
+                    else:
+                        point_list[position,4] += 1  
+        ind += 1         
+    return point_list
+
+
+def main(csv_file, slope_sc, stop_slope, Ind_A, Ind_S):
+    slope_sc = float(slope_sc)
+    stop_slope = float(stop_slope)
+    point_list = read_csv_1(csv_file)
+    result = read_csv_2(csv_file, slope_sc, stop_slope, Ind_A, Ind_S, point_list) 
+    csvfile_new = sys.stdout
+    for row in result:
+        csvfile_new.write(','.join(map(str, row))+ '\n')
+    
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(prog = 'filter_not')
+
+    parser.add_argument('-c', '--csv_file', help="Path to the CSV file")
+
+    parser.add_argument('-sc', '--sc_slope', default = 0.00005, help = 'parameter for the exponential decay function that determines after how many scatter events the ray is lost')
+
+    parser.add_argument('-s', '--stop_slope', default = 0.0000001, help = 'parameter for the exponential decay function that determines after what path length the ray is lost')
+
+    parser.add_argument('--index_glue', default = 1.5, help = 'refraction index of the glue')
+    
+    parser.add_argument('--index_scin', default = 2.15 , help = 'refraction index of the scintillator')
+
+    args = parser.parse_args()
+    
+    # Call the main function with parsed arguments
+    main(args.csv_file, args.sc_slope, args.stop_slope, args.index_glue, args.index_scin)
+
+
+
+
+

dutoit-rpi.gf

@@ -1,21 +0,0 @@
-#! /usr/bin/env ./geometryfactor.py
-# RPiRENA over Africa 2019
-
-φ1=0
-φ2=π/2
-Δr=0.2
-Δθ=1°
-radius=12.5
-
-# length units are mm, distances are center-center
-
-# dimensions of the SSDs
-define=wx=11.0
-define=wy=21.0
-define=wz=0.3
-define=d=5.0
-
-box=${wx}/2, ${wy}/2, ${wz}/2
-move=0, 0, -${d}/2
-copy=-1
-move=0, 0, ${d}

filter_not.py

@@ -0,0 +1,95 @@
+import csv
+import sys
+import numpy as np
+import time
+import argparse
+from Struktur import run, Rays, vector, Plane, Scintillator, Rays_custom
+from Geometry import regular_hexagonal_prism  
+
+def read_csv(file_path, notorious_rays_origin, th, phi,):
+    f = open(file_path, mode='r', encoding='ascii', errors='ignore')
+    ind_not = 0
+    for line in f:
+        row = line.strip().split(',')
+        if int(float(row[8])) == 0:
+            notorious_rays_origin[ind_not,0] = float(row[0])
+            notorious_rays_origin[ind_not,1] = float(row[1])
+            notorious_rays_origin[ind_not,2] = float(row[2])
+            th[ind_not] = float(row[3])
+            phi[ind_not] = float(row[4])
+            ind_not += 1
+    notorious_rays_origin = notorious_rays_origin[:ind_not,:]
+    th = th[:ind_not]
+    phi = phi[:ind_not]
+    return notorious_rays_origin, th, phi
+
+def main(csv_file, stop, Ind_S, height, edge, track):
+    notorious_rays_origin = np.empty((90000000,3))
+    stop = int(stop)
+    th = np.empty((90000000))
+    phi = np.empty((90000000))
+    origin, th, phi = read_csv(csv_file, notorious_rays_origin, th, phi) 
+    origin = np.expand_dims(origin, axis=2)
+    new_dir = np.empty((len(th),3,1))
+    new_dir[:,0]= (np.sin(np.radians(th)) * np.cos(np.radians(phi))).reshape(-1,1)
+    new_dir[:,1] = (np.sin(np.radians(th)) * np.sin(np.radians(phi))).reshape(-1,1)
+    new_dir[:,2] = (np.cos(np.radians(th))).reshape(-1,1)
+    csvfile_new = sys.stdout
+    det=[Plane(vector(-10,-35,5),vector(-10,-35,-5),vector(-10,-30,0),[]),
+        Plane(vector(-10,-35,-5),vector(10,-45,-5),vector(0,-35,-5),[]),
+        Plane(vector(10,-35,5),vector(10,-35,0),vector(10,-45,-5),[]),
+        Plane(vector(-10,-35,5),vector(0,-35,5),vector(10,-45,5),[])]
+    o = regular_hexagonal_prism([height*2, edge],[1,4])
+    p = Scintillator(o, det)
+    iter = (len(origin) // (10**6))
+    rest = (len(origin) % (10**6))
+    for i in range(iter):
+        origin_s = origin[i*(10**6):(i+1)*(10**6)]
+        new_o = np.empty((len(origin_s),3,1))
+        new_o[:,0:3] = origin_s[:,:,0:3]
+        new_o = np.transpose(new_o, (0,2,1))
+        new_o = new_o.reshape(-1, 3) 
+        new_o = [','.join(map(str, row)) for row in new_o]
+        r = Rays_custom(origin[i*(10**6):(i+1)*(10**6)], new_dir[i*(10**6):(i+1)*(10**6)])
+        th, phi = Rays_custom.angle(r)
+        sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
+        hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
+        data = list(zip(new_o, th, phi, sc, tr, path, hit_rec_processed))
+        for row in data:
+            csvfile_new.write(','.join(map(str, row)) + '\n')
+    if rest != 0:
+        origin_s = origin[iter*(10**6):]
+        new_o = np.empty((len(origin_s),3,1))
+        new_o[:,0:3] = origin_s[:,:,0:3]
+        new_o = np.transpose(new_o, (0,2,1))
+        new_o = new_o.reshape(-1, 3) 
+        new_o = [','.join(map(str, row)) for row in new_o]
+        r = Rays_custom(origin[iter*(10**6):], new_dir[iter*(10**6):])
+        th, phi = Rays_custom.angle(r)
+        sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
+        hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
+        data = list(zip(new_o, th, phi, sc, tr, path, hit_rec_processed))
+        for row in data:
+            csvfile_new.write(','.join(map(str, row)) + '\n')
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(prog = 'filter_not')
+
+    parser.add_argument('-c', '--csv_file', help="Path to the CSV file")
+
+    parser.add_argument('--Ind_S', default = 2.15, type = int, help = 'refraction index of scintillator')
+
+    parser.add_argument('-s', '--stop', default = 100000, help = 'path length after which programme quits')
+
+    parser.add_argument('--height', default = 10, type = int, help = 'height of the scintillator')
+
+    parser.add_argument('-e', '--edge', default = 51.962, type = float, help = 'edge length of the scintillator')
+
+    parser.add_argument('-t', '--track', help = 'set value to True if every cross_point with a plane should be in the output')
+    
+    args = parser.parse_args()
+    
+    # Call the main function with parsed arguments
+    main(args.csv_file, args.stop, args.Ind_S, args.height, args.edge, args.track)
+
+

gamma.awk

@@ -1,20 +0,0 @@
-#! /usr/bin/gawk -f
-
-BEGIN {
-    gamma = 1
-    thres = 10
-}
-
-!/^#/ && NF==3 {
-    w = cos($1)
-    if (w>0) w = $3 * exp(log(w)*gamma)
-    else w = 0
-    l = $2
-    L[l] += w
-    if (l>=thres) N += w
-}
-
-END {
-    print "#", N
-    for (l in L) print l, L[l] | "sort -n"
-}

geometryfactor-µM-Zγ.gpt

@@ -1,194 +0,0 @@
-#!/usr/local/bin/gnuplot -persist
-#
-#    
-#    	G N U P L O T
-#    	Version 5.5 patchlevel 0    last modified 2020-12-30
-#    
-#    	Copyright (C) 1986-1993, 1998, 2004, 2007-2020
-#    	Thomas Williams, Colin Kelley and many others
-#    
-#    	gnuplot home:     http://www.gnuplot.info
-#    	mailing list:     gnuplot-beta@lists.sourceforge.net
-#    	faq, bugs, etc:   type "help FAQ"
-#    	immediate help:   type "help"  (plot window: hit 'h')
-# set terminal qt 0 font "Sans,9"
-# set output
-unset clip points
-set clip one
-unset clip two
-unset clip radial
-set errorbars front 1.000000 
-set border 31 front lt black linewidth 1.000 dashtype solid
-set cornerpoles
-set zdata 
-set ydata 
-set xdata 
-set y2data 
-set x2data 
-set boxwidth
-set boxdepth 0
-set style fill  empty border
-set style rectangle back fc  bgnd fillstyle   solid 1.00 border lt -1
-set style circle radius graph 0.02 
-set style ellipse size graph 0.05, 0.03 angle 0 units xy
-set dummy x, y
-set format x "% h" 
-set format y "% h" 
-set format x2 "% h" 
-set format y2 "% h" 
-set format z "% h" 
-set format cb "% h" 
-set format r "% h" 
-set ttics format "% h"
-set timefmt "%d/%m/%y,%H:%M"
-set angles radians
-set tics back
-set grid nopolar
-set grid xtics nomxtics noytics nomytics noztics nomztics nortics nomrtics \
- nox2tics nomx2tics y2tics nomy2tics nocbtics nomcbtics
-set grid layerdefault   lt 0 linecolor 0 linewidth 0.500,  lt 0 linecolor 0 linewidth 0.500
-unset raxis
-set theta counterclockwise right
-set style parallel front  lt black linewidth 2.000 dashtype solid
-set key notitle
-set key fixed right top vertical Right noreverse enhanced autotitle nobox
-set key noinvert samplen 4 spacing 1 width 0 height 0 
-set key maxcolumns 0 maxrows 0
-set key noopaque
-unset label
-unset arrow
-unset style line
-unset style arrow
-set style histogram clustered gap 2 title textcolor lt -1
-unset object
-unset walls
-set style textbox  transparent margins  1.0,  1.0 border  lt -1 linewidth  1.0
-set offsets 0, 0, 0, 0
-set pointsize 1
-set pointintervalbox 1
-set encoding default
-unset polar
-unset parametric
-unset spiderplot
-unset decimalsign
-unset micro
-unset minussign
-set view 60, 30, 1, 1
-set view azimuth 0
-set rgbmax 255
-set samples 100, 100
-set isosamples 10, 10
-set surface 
-unset contour
-set cntrlabel  format '%8.3g' font '' start 5 interval 20
-set mapping cartesian
-set datafile separator whitespace
-set datafile nocolumnheaders
-unset hidden3d
-set cntrparam order 4
-set cntrparam linear
-set cntrparam levels 5
-set cntrparam levels auto
-set cntrparam firstlinetype 0 unsorted
-set cntrparam points 5
-set size ratio 0 1,1
-set origin 0,0
-set style data points
-set style function lines
-unset xzeroaxis
-unset yzeroaxis
-unset zzeroaxis
-unset x2zeroaxis
-unset y2zeroaxis
-set xyplane relative 0.5
-set tics scale  1, 0.5, 1, 1, 1
-set mxtics default
-set mytics default
-set mztics default
-set mx2tics default
-set my2tics default
-set mcbtics default
-set mrtics default
-set nomttics
-set xtics border in scale 1,0.5 mirror norotate  autojustify
-set xtics  norangelimit autofreq 
-set ytics border in scale 1,0.5 nomirror norotate  autojustify
-set ytics  norangelimit autofreq 
-set ztics border in scale 1,0.5 nomirror norotate  autojustify
-set ztics  norangelimit autofreq 
-unset x2tics
-set y2tics border in scale 1,0.5 nomirror norotate  autojustify
-set y2tics  norangelimit autofreq 
-set cbtics border in scale 1,0.5 mirror norotate  autojustify
-set cbtics  norangelimit autofreq 
-set rtics axis in scale 1,0.5 nomirror norotate  autojustify
-set rtics  norangelimit autofreq 
-unset ttics
-set title "\302\265Mustang count rate ratio by flux-\316\263 and sensor orientation" 
-set title  font "" textcolor lt -1 norotate
-set timestamp bottom 
-set timestamp "" 
-set timestamp  font "" textcolor lt -1 norotate
-set trange [ * : * ] noreverse nowriteback
-set urange [ * : * ] noreverse nowriteback
-set vrange [ * : * ] noreverse nowriteback
-set xlabel "\316\263" 
-set xlabel  font "" textcolor lt -1 norotate
-set x2label "" 
-set x2label  font "" textcolor lt -1 norotate
-set xrange [ * : * ] noreverse writeback
-set x2range [ * : * ] noreverse writeback
-set ylabel "zenith-\316\223 [cm\302\262 sr]" 
-set ylabel  font "" textcolor lt -1 rotate
-set y2label "" 
-set y2label  font "" textcolor lt -1 rotate
-set yrange [ 0.00000 : 350.000 ] noreverse writeback
-set y2range [ * : 1.00000 ] noreverse writeback
-set zlabel "" 
-set zlabel  font "" textcolor lt -1 norotate
-set zrange [ * : * ] noreverse writeback
-set cblabel "" 
-set cblabel  font "" textcolor lt -1 rotate
-set cbrange [ * : * ] noreverse writeback
-set rlabel "" 
-set rlabel  font "" textcolor lt -1 norotate
-set rrange [ * : * ] noreverse writeback
-unset logscale
-unset jitter
-set zero 1e-08
-set lmargin  -1
-set bmargin  -1
-set rmargin  -1
-set tmargin  -1
-set locale "de_DE.UTF-8"
-set pm3d explicit at s
-set pm3d scansautomatic
-set pm3d interpolate 1,1 flush begin noftriangles noborder corners2color mean
-set pm3d clip z 
-set pm3d nolighting
-set palette positive nops_allcF maxcolors 0 gamma 1.5 color model RGB 
-set palette rgbformulae 7, 5, 15
-set colorbox default
-set colorbox vertical origin screen 0.9, 0.2 size screen 0.05, 0.6 front  noinvert border -1 cbtics 0
-set style boxplot candles range  1.50 outliers pt 7 separation 1 labels auto unsorted
-set loadpath 
-set fontpath
-set psdir
-set fit brief errorvariables nocovariancevariables errorscaling prescale nowrap v5
-turbo_red(z)  =turbo_r[1]+z*(turbo_r[2]+z*(turbo_r[3]+z*(turbo_r[4]+z*(turbo_r[5]+z*turbo_r[6]))))
-turbo_green(z)=turbo_g[1]+z*(turbo_g[2]+z*(turbo_g[3]+z*(turbo_g[4]+z*(turbo_g[5]+z*turbo_g[6]))))
-turbo_blue(z) =turbo_b[1]+z*(turbo_b[2]+z*(turbo_b[3]+z*(turbo_b[4]+z*(turbo_b[5]+z*turbo_b[6]))))
-GNUTERM = "qt"
-I = {0.0, 1.0}
-VoxelDistance = 0.0
-GridDistance = 0.0
-cois = " C-3+3 C-2+3 C-1+3 C+0+3 C+1+3 C+2+3 C+3+3 C-3+2 C-2+2 C-1+2 C+0+2 C+1+2 C+2+2 C+3+2 C-3+1 C-2+1 C-1+1 C+0+1 C+1+1 C+2+1 C+3+1 C-3+0 C-2+0 C-1+0 C+0+0 C+1+0 C+2+0 C+3+0 C-3-1 C-2-1 C-1-1 C+0-1 C+1-1 C+2-1 C+3-1 C-3-2 C-2-2 C-1-2 C+0-2 C+1-2 C+2-2 C+3-2 C-3-3 C-2-3 C-1-3 C+0-3 C+1-3 C+2-3 C+3-3"
-x = 0.0
-array turbo_r[6] = [0.13572138,4.6153926,-42.66032258,132.13108234,-152.94239396,59.28637943]
-array turbo_g[6] = [0.09140261,2.19418839,4.84296658,-14.18503333,4.27729857,2.82956604]
-array turbo_b[6] = [0.1066733,12.64194608,-60.58204836,110.36276771,-89.90310912,27.34824973]
-array gammas[4] = [0,1,2,3]
-plotext = "pdf"
-## Last datafile plotted: "<paste data/geometryfactor-µM.γ data/geometryfactor-µM-rx.γ"
-plot "data/geometryfactor-µM.γ" u 1:($3*8/52592.7/100) tit "horizontal" w lp pt 7, "data/geometryfactor-µM-rx.γ" u 1:($3*4/52592.7/100) tit "vertical" w lp pt 7, "<paste data/geometryfactor-µM.γ data/geometryfactor-µM-rx.γ" u 1:($6/$3/2) axis x1y2 tit "vertical/horizontal" w lp pt 7
-#    EOF

geometryfactor-µM-rx_θΔ.gpt

@@ -1,192 +0,0 @@
-#!/usr/local/bin/gnuplot -persist
-#
-#    
-#    	G N U P L O T
-#    	Version 5.5 patchlevel 0    last modified 2020-12-30
-#    
-#    	Copyright (C) 1986-1993, 1998, 2004, 2007-2020
-#    	Thomas Williams, Colin Kelley and many others
-#    
-#    	gnuplot home:     http://www.gnuplot.info
-#    	mailing list:     gnuplot-beta@lists.sourceforge.net
-#    	faq, bugs, etc:   type "help FAQ"
-#    	immediate help:   type "help"  (plot window: hit 'h')
-# set terminal qt 0 font "Sans,9"
-# set output
-unset clip points
-set clip one
-unset clip two
-unset clip radial
-set errorbars front 1.000000 
-set border 31 front lt black linewidth 1.000 dashtype solid
-set cornerpoles
-set zdata 
-set ydata 
-set xdata 
-set y2data 
-set x2data 
-set boxwidth
-set boxdepth 0
-set style fill  empty border
-set style rectangle back fc  bgnd fillstyle   solid 1.00 border lt -1
-set style circle radius graph 0.02 
-set style ellipse size graph 0.05, 0.03 angle 0 units xy
-set dummy x, y
-set format x "% h" 
-set format y "% h" 
-set format x2 "% h" 
-set format y2 "% h" 
-set format z "% h" 
-set format cb "% h" 
-set format r "% h" 
-set ttics format "% h"
-set timefmt "%d/%m/%y,%H:%M"
-set angles radians
-set tics back
-unset grid
-unset raxis
-set theta counterclockwise right
-set style parallel front  lt black linewidth 2.000 dashtype solid
-set key notitle
-set key fixed right top vertical Right noreverse enhanced autotitle nobox
-set key noinvert samplen 4 spacing 1 width 0 height 0 
-set key maxcolumns 0 maxrows 0
-set key noopaque
-unset label
-unset arrow
-unset style line
-unset style arrow
-set style histogram clustered gap 2 title textcolor lt -1
-unset object
-unset walls
-set style textbox  transparent margins  1.0,  1.0 border  lt -1 linewidth  1.0
-set offsets 0, 0, 0, 0
-set pointsize 1
-set pointintervalbox 1
-set encoding default
-unset polar
-unset parametric
-unset spiderplot
-unset decimalsign
-unset micro
-unset minussign
-set view 60, 30, 1, 1
-set view azimuth 0
-set rgbmax 255
-set samples 100, 100
-set isosamples 10, 10
-set surface 
-unset contour
-set cntrlabel  format '%8.3g' font '' start 5 interval 20
-set mapping cartesian
-set datafile separator whitespace
-set datafile nocolumnheaders
-unset hidden3d
-set cntrparam order 4
-set cntrparam linear
-set cntrparam levels 5
-set cntrparam levels auto
-set cntrparam firstlinetype 0 unsorted
-set cntrparam points 5
-set size ratio 0 1,1
-set origin 0,0
-set style data points
-set style function lines
-unset xzeroaxis
-unset yzeroaxis
-unset zzeroaxis
-unset x2zeroaxis
-unset y2zeroaxis
-set xyplane relative 0.5
-set tics scale  1, 0.5, 1, 1, 1
-set mxtics default
-set mytics default
-set mztics default
-set mx2tics default
-set my2tics default
-set mcbtics default
-set mrtics default
-set nomttics
-set xtics border in scale 1,0.5 mirror norotate  autojustify
-set xtics  norangelimit 
-set xtics   (0.00000, "\342\205\233\317\200" 0.392699, "\302\274\317\200" 0.785398, "\342\205\234\317\200" 1.17810, "\302\275\317\200" 1.57080)
-set ytics border in scale 1,0.5 mirror norotate  autojustify
-set ytics  norangelimit autofreq 
-set ztics border in scale 1,0.5 nomirror norotate  autojustify
-set ztics  norangelimit autofreq 
-unset x2tics
-unset y2tics
-set cbtics border in scale 1,0.5 mirror norotate  autojustify
-set cbtics  norangelimit logscale autofreq 
-set rtics axis in scale 1,0.5 nomirror norotate  autojustify
-set rtics  norangelimit autofreq 
-unset ttics
-set title "geometryfactor-\302\265M-rx" 
-set title  font "" textcolor lt -1 norotate
-set timestamp bottom 
-set timestamp "" 
-set timestamp  font "" textcolor lt -1 norotate
-set trange [ * : * ] noreverse nowriteback
-set urange [ * : * ] noreverse nowriteback
-set vrange [ * : * ] noreverse nowriteback
-set xlabel "\316\270" 
-set xlabel  font "" textcolor lt -1 norotate
-set x2label "" 
-set x2label  font "" textcolor lt -1 norotate
-set xrange [ 0.00000 : 1.57080 ] noreverse writeback
-set x2range [ * : * ] noreverse writeback
-set ylabel "\316\224 [mm]" 
-set ylabel  font "" textcolor lt -1 rotate
-set y2label "" 
-set y2label  font "" textcolor lt -1 rotate
-set yrange [ 0.00000 : 120.000 ] noreverse writeback
-set y2range [ * : * ] noreverse writeback
-set zlabel "" 
-set zlabel  font "" textcolor lt -1 norotate
-set zrange [ * : * ] noreverse writeback
-set cblabel "\316\223 [(mm\302\262 sr)/(mm \302\260)]" 
-set cblabel  font "" textcolor lt -1 rotate
-set cbrange [ * : * ] noreverse writeback
-set rlabel "" 
-set rlabel  font "" textcolor lt -1 norotate
-set rrange [ * : * ] noreverse writeback
-unset logscale
-set logscale cb 10
-unset jitter
-set zero 1e-08
-set lmargin  -1
-set bmargin  -1
-set rmargin  -1
-set tmargin  -1
-set locale "de_DE.UTF-8"
-set pm3d explicit at s
-set pm3d scansautomatic
-set pm3d interpolate 1,1 flush begin noftriangles noborder corners2color mean
-set pm3d clip z 
-set pm3d nolighting
-set palette positive nops_allcF maxcolors 0 gamma 1.5 color model RGB 
-set palette functions turbo_red(gray), turbo_green(gray), turbo_blue(gray)
-set colorbox default
-set colorbox vertical origin screen 0.9, 0.2 size screen 0.05, 0.6 front  noinvert border -1 cbtics -1
-set style boxplot candles range  1.50 outliers pt 7 separation 1 labels auto unsorted
-set loadpath 
-set fontpath
-set psdir
-set fit brief errorvariables nocovariancevariables errorscaling prescale nowrap v5
-turbo_red(z)  =turbo_r[1]+z*(turbo_r[2]+z*(turbo_r[3]+z*(turbo_r[4]+z*(turbo_r[5]+z*turbo_r[6]))))
-turbo_green(z)=turbo_g[1]+z*(turbo_g[2]+z*(turbo_g[3]+z*(turbo_g[4]+z*(turbo_g[5]+z*turbo_g[6]))))
-turbo_blue(z) =turbo_b[1]+z*(turbo_b[2]+z*(turbo_b[3]+z*(turbo_b[4]+z*(turbo_b[5]+z*turbo_b[6]))))
-GNUTERM = "qt"
-I = {0.0, 1.0}
-VoxelDistance = 0.0
-GridDistance = 0.0
-cois = " C-3+3 C-2+3 C-1+3 C+0+3 C+1+3 C+2+3 C+3+3 C-3+2 C-2+2 C-1+2 C+0+2 C+1+2 C+2+2 C+3+2 C-3+1 C-2+1 C-1+1 C+0+1 C+1+1 C+2+1 C+3+1 C-3+0 C-2+0 C-1+0 C+0+0 C+1+0 C+2+0 C+3+0 C-3-1 C-2-1 C-1-1 C+0-1 C+1-1 C+2-1 C+3-1 C-3-2 C-2-2 C-1-2 C+0-2 C+1-2 C+2-2 C+3-2 C-3-3 C-2-3 C-1-3 C+0-3 C+1-3 C+2-3 C+3-3"
-x = 0.0
-array turbo_r[6] = [0.13572138,4.6153926,-42.66032258,132.13108234,-152.94239396,59.28637943]
-array turbo_g[6] = [0.09140261,2.19418839,4.84296658,-14.18503333,4.27729857,2.82956604]
-array turbo_b[6] = [0.1066733,12.64194608,-60.58204836,110.36276771,-89.90310912,27.34824973]
-gammas = "1 2 3"
-plotext = "pdf"
-## Last datafile plotted: "/data/blaulicht/stephan/src/geometryfactor-µM-rx.θΔ"
-plot "/data/blaulicht/stephan/src/geometryfactor-µM-rx.θΔ" u 1:2:($3/52592.7*8)  with image
-#    EOF

geometryfactor-µM_θΔ.gpt

@@ -1,192 +0,0 @@
-#!/usr/local/bin/gnuplot -persist
-#
-#    
-#    	G N U P L O T
-#    	Version 5.5 patchlevel 0    last modified 2020-12-30
-#    
-#    	Copyright (C) 1986-1993, 1998, 2004, 2007-2020
-#    	Thomas Williams, Colin Kelley and many others
-#    
-#    	gnuplot home:     http://www.gnuplot.info
-#    	mailing list:     gnuplot-beta@lists.sourceforge.net
-#    	faq, bugs, etc:   type "help FAQ"
-#    	immediate help:   type "help"  (plot window: hit 'h')
-# set terminal qt 0 font "Sans,9"
-# set output
-unset clip points
-set clip one
-unset clip two
-unset clip radial
-set errorbars front 1.000000 
-set border 31 front lt black linewidth 1.000 dashtype solid
-set cornerpoles
-set zdata 
-set ydata 
-set xdata 
-set y2data 
-set x2data 
-set boxwidth
-set boxdepth 0
-set style fill  empty border
-set style rectangle back fc  bgnd fillstyle   solid 1.00 border lt -1
-set style circle radius graph 0.02 
-set style ellipse size graph 0.05, 0.03 angle 0 units xy
-set dummy x, y
-set format x "% h" 
-set format y "% h" 
-set format x2 "% h" 
-set format y2 "% h" 
-set format z "% h" 
-set format cb "% h" 
-set format r "% h" 
-set ttics format "% h"
-set timefmt "%d/%m/%y,%H:%M"
-set angles radians
-set tics back
-unset grid
-unset raxis
-set theta counterclockwise right
-set style parallel front  lt black linewidth 2.000 dashtype solid
-set key notitle
-set key fixed right top vertical Right noreverse enhanced autotitle nobox
-set key noinvert samplen 4 spacing 1 width 0 height 0 
-set key maxcolumns 0 maxrows 0
-set key noopaque
-unset label
-unset arrow
-unset style line
-unset style arrow
-set style histogram clustered gap 2 title textcolor lt -1
-unset object
-unset walls
-set style textbox  transparent margins  1.0,  1.0 border  lt -1 linewidth  1.0
-set offsets 0, 0, 0, 0
-set pointsize 1
-set pointintervalbox 1
-set encoding default
-unset polar
-unset parametric
-unset spiderplot
-unset decimalsign
-unset micro
-unset minussign
-set view 60, 30, 1, 1
-set view azimuth 0
-set rgbmax 255
-set samples 100, 100
-set isosamples 10, 10
-set surface 
-unset contour
-set cntrlabel  format '%8.3g' font '' start 5 interval 20
-set mapping cartesian
-set datafile separator whitespace
-set datafile nocolumnheaders
-unset hidden3d
-set cntrparam order 4
-set cntrparam linear
-set cntrparam levels 5
-set cntrparam levels auto
-set cntrparam firstlinetype 0 unsorted
-set cntrparam points 5
-set size ratio 0 1,1
-set origin 0,0
-set style data points
-set style function lines
-unset xzeroaxis
-unset yzeroaxis
-unset zzeroaxis
-unset x2zeroaxis
-unset y2zeroaxis
-set xyplane relative 0.5
-set tics scale  1, 0.5, 1, 1, 1
-set mxtics default
-set mytics default
-set mztics default
-set mx2tics default
-set my2tics default
-set mcbtics default
-set mrtics default
-set nomttics
-set xtics border in scale 1,0.5 mirror norotate  autojustify
-set xtics  norangelimit 
-set xtics   (0.00000, "\342\205\233\317\200" 0.392699, "\302\274\317\200" 0.785398, "\342\205\234\317\200" 1.17810, "\302\275\317\200" 1.57080)
-set ytics border in scale 1,0.5 mirror norotate  autojustify
-set ytics  norangelimit autofreq 
-set ztics border in scale 1,0.5 nomirror norotate  autojustify
-set ztics  norangelimit autofreq 
-unset x2tics
-unset y2tics
-set cbtics border in scale 1,0.5 mirror norotate  autojustify
-set cbtics  norangelimit logscale autofreq 
-set rtics axis in scale 1,0.5 nomirror norotate  autojustify
-set rtics  norangelimit autofreq 
-unset ttics
-set title "geometryfactor-\302\265M" 
-set title  font "" textcolor lt -1 norotate
-set timestamp bottom 
-set timestamp "" 
-set timestamp  font "" textcolor lt -1 norotate
-set trange [ * : * ] noreverse nowriteback
-set urange [ * : * ] noreverse nowriteback
-set vrange [ * : * ] noreverse nowriteback
-set xlabel "\316\270" 
-set xlabel  font "" textcolor lt -1 norotate
-set x2label "" 
-set x2label  font "" textcolor lt -1 norotate
-set xrange [ 0.00000 : 1.57080 ] noreverse writeback
-set x2range [ * : * ] noreverse writeback
-set ylabel "\316\224 [mm]" 
-set ylabel  font "" textcolor lt -1 rotate
-set y2label "" 
-set y2label  font "" textcolor lt -1 rotate
-set yrange [ 0.00000 : 120.000 ] noreverse writeback
-set y2range [ * : * ] noreverse writeback
-set zlabel "" 
-set zlabel  font "" textcolor lt -1 norotate
-set zrange [ * : * ] noreverse writeback
-set cblabel "\316\223 [(mm\302\262 sr)/(mm \302\260)]" 
-set cblabel  font "" textcolor lt -1 rotate
-set cbrange [ * : * ] noreverse writeback
-set rlabel "" 
-set rlabel  font "" textcolor lt -1 norotate
-set rrange [ * : * ] noreverse writeback
-unset logscale
-set logscale cb 10
-unset jitter
-set zero 1e-08
-set lmargin  -1
-set bmargin  -1
-set rmargin  -1
-set tmargin  -1
-set locale "de_DE.UTF-8"
-set pm3d explicit at s
-set pm3d scansautomatic
-set pm3d interpolate 1,1 flush begin noftriangles noborder corners2color mean
-set pm3d clip z 
-set pm3d nolighting
-set palette positive nops_allcF maxcolors 0 gamma 1.5 color model RGB 
-set palette functions turbo_red(gray), turbo_green(gray), turbo_blue(gray)
-set colorbox default
-set colorbox vertical origin screen 0.9, 0.2 size screen 0.05, 0.6 front  noinvert border -1 cbtics 0
-set style boxplot candles range  1.50 outliers pt 7 separation 1 labels auto unsorted
-set loadpath 
-set fontpath
-set psdir
-set fit brief errorvariables nocovariancevariables errorscaling prescale nowrap v5
-turbo_red(z)  =turbo_r[1]+z*(turbo_r[2]+z*(turbo_r[3]+z*(turbo_r[4]+z*(turbo_r[5]+z*turbo_r[6]))))
-turbo_green(z)=turbo_g[1]+z*(turbo_g[2]+z*(turbo_g[3]+z*(turbo_g[4]+z*(turbo_g[5]+z*turbo_g[6]))))
-turbo_blue(z) =turbo_b[1]+z*(turbo_b[2]+z*(turbo_b[3]+z*(turbo_b[4]+z*(turbo_b[5]+z*turbo_b[6]))))
-GNUTERM = "qt"
-I = {0.0, 1.0}
-VoxelDistance = 0.0
-GridDistance = 0.0
-cois = " C-3+3 C-2+3 C-1+3 C+0+3 C+1+3 C+2+3 C+3+3 C-3+2 C-2+2 C-1+2 C+0+2 C+1+2 C+2+2 C+3+2 C-3+1 C-2+1 C-1+1 C+0+1 C+1+1 C+2+1 C+3+1 C-3+0 C-2+0 C-1+0 C+0+0 C+1+0 C+2+0 C+3+0 C-3-1 C-2-1 C-1-1 C+0-1 C+1-1 C+2-1 C+3-1 C-3-2 C-2-2 C-1-2 C+0-2 C+1-2 C+2-2 C+3-2 C-3-3 C-2-3 C-1-3 C+0-3 C+1-3 C+2-3 C+3-3"
-x = 0.0
-array turbo_r[6] = [0.13572138,4.6153926,-42.66032258,132.13108234,-152.94239396,59.28637943]
-array turbo_g[6] = [0.09140261,2.19418839,4.84296658,-14.18503333,4.27729857,2.82956604]
-array turbo_b[6] = [0.1066733,12.64194608,-60.58204836,110.36276771,-89.90310912,27.34824973]
-gammas = "1 2 3"
-plotext = "pdf"
-## Last datafile plotted: "/data/blaulicht/stephan/src/geometryfactor-µM.θΔ"
-plot "/data/blaulicht/stephan/src/geometryfactor-µM.θΔ" u 1:2:($3/52592.7*16)  with image
-#    EOF

geometryfactor.awk

@@ -1,109 +0,0 @@
-#! /usr/bin/gawk -i
-
-function mult() {
-    nHITS=0
-    for (i=5; i<=NF; i++) 
-	if ($i > 0) {
-	    TLEN[nHITS] = $i
-	    HITS[nHITS++] = i-5
-	}
-    for (i=0; i<nHITS-1; i++) 
-	for (j=i+1; j<nHITS; j++)
-	    if (TLEN[i] < TLEN[j]) {
-		a=HITS[j]
-		HITS[j]=HITS[i]
-		HITS[i]=a
-		a=TLEN[j]
-		TLEN[j]=TLEN[i]
-		TLEN[i]=a
-	    }
-    return nHITS
-}
-
-BEGIN {
-    pi = 4*atan2(1,1)
-    deg = pi/180
-    res_xy = 0.5*deg
-}
-
-function hasDZ() {
-    return mult()>=2 && and(xor(HITS[0], HITS[1]), 16)
-}
-
-function Y(c) { return and(c,4)/2  + and(c,1)   - 1.5 }
-function X(c) { return and(c,8)/4  + and(c,2)/2 - 1.5 }
-function Z(c) { return and(c,16) }
-    
-#
-# Mirror φ=0…π/2 to the full circle.
-# print -z first
-#
-
-function doCoinc1(f, phi, c0, c1) {
-    if (Z(c1)) @f($1, phi, c0, c1, TLEN[0], TLEN[1])
-    else       @f($1, phi, c1, c0, TLEN[1], TLEN[0])
-}
-
-function doCoinc(f) {
-    doCoinc1(f, $2,        HITS[0],         HITS[1]    )
-    if (!$1) return
-    doCoinc1(f, pi-$2, xor(HITS[0],10), xor(HITS[1],10))
-    doCoinc1(f,   -$2, xor(HITS[0], 5), xor(HITS[1], 5))
-    doCoinc1(f, pi+$2, xor(HITS[0],15), xor(HITS[1],15))
-}
-
-function printCoinc(theta, phi, c0, c1, l0, l1) {
-    print $1, phi, X(c1)-X(c0), Y(c1)-Y(c0), c0, c1, l0, l1
-}
-
-BEGIN {
-    for (i=0; i<32; i++)
-	for (j=0; j<32; j++)
-	    HKEYxy[i,j] = sprintf("C%+d%+d", X(j)-X(i), Y(j)-Y(i))
-}
-
-function histxyCoinc(theta, phi, c0, c1, l0, l1) {
-    x = int(sin(theta)*cos(phi)/res_xy+1000.5)-1000
-    y = int(sin(theta)*sin(phi)/res_xy+1000.5)-1000
-    k = HKEYxy[c0,c1]
-    if (!iHIST2D) {
-	min_x = x
-	max_x = x
-	min_y = y
-	max_y = y
-    } else {
-	if (x>max_x) max_x = x
-	if (x<min_x) min_x = x
-	if (y>max_y) max_y = y
-	if (y<min_y) min_y = y
-    }
-    iHIST2D++
-    HIST2DKEY[k]++
-    HIST2D[k,x,y]++
-}
-
-END { if (iHIST2D) emit_HIST2D() }
-
-function emit_HIST2D() {
-    n=0;
-    for (k in HIST2DKEY) kk[n++] = k
-    for (i=0; i<n-1; i++) 
-	for (j=i+1; j<n; j++)
-	    if (kk[i] > kk[j]) {
-		k=kk[j]
-		kk[j]=kk[i]
-		kk[i]=k
-	    }
-    printf "x y "
-    for (i=0; i<n; i++) printf " %s", kk[i]
-    printf "\n"
-    for (x=min_x; x<=max_x; x++) {
-	for (y=min_y; y<=max_y; y++) {
-	    printf "%.3f %.3f ", x*res_xy, y*res_xy
-	    for (i=0; i<n; i++)
-		printf " %d", HIST2D[kk[i],x,y]+0
-	    printf "\n"
-	}
-	printf "\n"
-    }
-}

geometryfactor.md

@@ -1,124 +0,0 @@
-Integrate geometry factors of convex polyhedrons
-================================================
-The world can contain multiple detectors.  Detectors are defined as
-sets of planes.  The pathlengths of the intersection of rays with the
-detectors are computed.
-
-The rays are generated in a range of angles θ₁…θ₂ from the *z*-axis in
-steps Δθ, and a range of angles φ₁…φ₂ around the *z*-axis with density
-1/Δθ².
-
-For each direction a beam of radius *r* is generated, with parallel
-rays on circles around the ray through the origin separated by Δ*r*
-and density 1/Δ*r*².
-
-Data model
-----------
-A plane is defined by a point and two vectors, `plane(o,v,w)`.
-The cross product `n=v×w` points inside the detector.
-
-A ray is defined by a point and a vector `ray().ov(o,v)`.  The
-constructor of a ray() takes four arguments `ray(θ, φ, r, ψ)` and
-computes the point and vector.  (*r*,ψ) define a point in the
-*x*,*y*-plane, that is rotated by θ around the *y*-axis and then by φ
-around the *z*-axis.  The ray goes through the rotated point
-perpendicular to the rotated *x*,*y*-plane.
-
-Running the script
-------------------
-The script `geometryfactor.py` is run with commandline options.  Any
-non-options are files containing more options.
-
-Each line in a file is treated as a further commandline argument,
-except for empty lines and lines starting with `#`.  Lines starting
-with a `-` are used as they are. Lines starting with a `@` are used
-with the `@` stripped off.  All other lines get `--` prefixed.
-
-Action options
---------------
-These actions are executed after all other options were applied, in
-this order:
-
-- `-E` `--corners`: find all corners of all detectors.
-- `-S` `--size`: compute the size of the integration run.
-- `-T` `--test`: run the test rays through the world.
-- `-R` `--run`: perform the integration.
-- `-G` `--gnuplot`: output data to splot the world with gnuplot.
-
-Beam options
-------------
-Option values are evaluated as python expressions. `°` is translated
-as `*deg`.  `math.pi` is imported as `π`.  Macros are expanded,
-recursively.
-
-- `-a` `--Δθ=`Δθ: angular resolution (1°)
-- `-d` `--Δr=`Δr: spacial resolution (1)
-- `--θ1=`θ₁: θ starting angle (0)
-- `--θ2=`θ₂: θ stopping angle (π/2)
-- `--φ1=`φ₁: φ starting angle (0)
-- `--φ2=`φ₂: φ stopping angle (π/4)
-- `-r` `--radius=`*r*: beam radius (0)
-
-Defaults are given in parenthesis. `--radius=0` implies `--corners` and
-sets the beam radius to the radius of the world.
-
-World options
--------------
-Points or vectors are space or comma separated expressions evaluated
-as explained in section _beam options_.  Multiple vectors are
-separated by `;`.
-
-### Creating detectors
-
-- `-N` `--new`: Append a new empty detector to the world.
-- `-B` `--box=`*x*,*y*,*z*: Append a new rectangular box detector to
-  the world. The dimensions are the half-widths of the box, i.e., the
-  coordinates of the corners.
-- `C` `--copy=`*n*: Append a copy of detector number *n* to the
-  world. (python list index)
-- `--select=`*n*: Move detector number *n* to the end of the world.
-- `-µ` `--µMustang`: Append a µMustang Detector to the end of the
-  world. 
-
-Empty detectors are not retained, except at the end of the world
-during genesis.  When no detectors are given, `-µ` is implied.
-
-### Detector manipulation
-
-Detector manipulation always affects the last detector of the world.
-
-- `-p` `--plane=`*o*;*v*;*w*: Add a plane to the detector.
-- `--move=`*v*: Move the detector by the vector *v*.
-- `--scale=`*v*: Scale the detector by factors given for each axis.
-- `--rx=`φ: Rotate the detector by the angle φ around the *x*-axis,
-- `--ry=`φ: *y*-axis,
-- `--rz=`φ: or *z*-axis.
-
-Test beam
----------
-Three rays along the coordinate axes are defined for the
-execution of `--test`.  Further rays may be given by
-
-- `-x` `--ray=`*name*;*o*;*v*: with point *o* and vector *v*, or
-- `-x` `--ray=`*name*;θ,φ,r,ψ: with direction and offset.
-
-Macros
-------
-Macros are simple textual substitutions in expressions evaluated as
-part of option arguments.  Strings of the form `${`*macro*`}` are
-replaced by the definition of *macro*.
-
-- `-D` `--define=`*macro*`=`*text*:  Define *macro* unless it is
-  already defined.
-- `-U` `--undefine=`*macro*: Undefine *macro*.
-- `--replace=`*macro*`=`*text*: Define or redefine *macro*.
-
-Output
-------
-The output of `--corners`, `--size`, and `--test` goes to _stderr_.
-
-`--run` produces one line of output on _stdout_ for each ray sent
-through the world.  The columns are: θ, φ, r, ψ, and the intersection
-length for each detector in the world.
-
-`--gnuplot` output goes to _stdout_.

geometryfactor.py

@@ -1,547 +0,0 @@
-#! /usr/bin/python3
-
-from sys import argv, stdout, stderr
-from math import pi as π
-import numpy
-from numpy import sqrt, sin, cos, tan, array, arange, empty, cross, newaxis
-from numpy.linalg import solve
-from getopt import getopt
-
-verbose=0
-
-def scale(x=1., y=1., z=1.):
-    return array(((x,0,0),
-                  (0,y,0),
-                  (0,0,z)))
-
-def rotate(axis, φ):
-    m = scale()
-    c = cos(φ)
-    s = -sin(φ)
-    for i in range(3):
-        if i != axis:
-            m[i,i]=c
-            for j in range(3):
-                if j != i and j != axis:
-                    m[i,j] = s
-                    s = -s
-    return m
-
-def vector(x=0., y=0., z=0.):
-    return array([[x],[y],[z]])
-
-def vtuple(v):
-    return tuple(v[:,0])
-
-def length(v):
-    return sqrt(v.T @ v)[...,0,0]
-
-class ray:
-    "Basic ray impelmentation"
-    
-    def __init__(self, θ=0, φ=0, r=0, ψ=0):
-        self.θ = θ
-        self.φ = φ
-        self.r = r
-        self.ψ = ψ
-        self.dir = rotate(2,φ) @ rotate(1,θ)
-        self.o = self.dir @ vector(r*cos(ψ), r*sin(ψ))
-        self.v = self.dir @ vector(z=1)
-
-    def ov(self, o, v):
-        self.o = o
-        self.v = v
-        return self
-        
-    def __repr__(self):
-        return "ray(o=(%g,%g,%g), v=(%g,%g,%g))" % (
-            vtuple(self.o)+vtuple(self.v))
-
-class rays:
-    "Fast reimplementation of a beam of rays"
-    
-    def __init__(self, r, Δr):
-        nr = int(r/Δr + 1.001)
-        self.n = 0
-        rψ = 0
-        for ir in range(nr):
-            nnψ = 2*π*ir + 0.5 + rψ
-            nψ = int(nnψ + 0.5)
-            rψ = nnψ - nψ
-            self.n += nψ
-        self.xy = empty((self.n,3,1))
-        self.r = empty((self.n,))
-        self.ψ = empty((self.n,))
-        i=0
-        rψ = 0
-        for ir in range(nr):
-            r = ir*Δr
-            nnψ = 2*π*ir + 0.5 + rψ
-            nψ = int(nnψ + 0.5)
-            rψ = nnψ - nψ
-            dψ = 2*π/nψ
-            for iψ in range(nψ):
-                ψ = iψ * dψ + dψ/2
-                self.r[i] = r
-                self.ψ[i] = ψ
-                self.xy[i,:,0] = (r*cos(ψ), r*sin(ψ), 0)
-                i = i+1
-
-    def bend(self, θ, φ):
-        self.θ = θ
-        self.φ = φ
-        self.dir = rotate(2,φ) @ rotate(1,-θ)
-        self.v = self.dir @ vector(z=1)
-        self.o = (self.dir @ self.xy)[...,0].T
-    
-class plane:
-    index=[0]
-    
-    def __init__(self, o, v, w):
-        self.idx = self.index[0]
-        self.index[0] += 1
-        self.o = o
-        self.v = v
-        self.w = w
-        self.normalize()
-
-    def normalize(self):
-        self.n = cross(self.v, self.w, axis=0).T
-        self.M = empty((3,3))
-        self.M[:,0:1]=self.v
-        self.M[:,1:2]=self.w
-        self.rays=(ray().ov(self.o, self.v),
-                   ray().ov(self.o, self.w),
-                   ray().ov(self.o+self.w, self.v),
-                   ray().ov(self.o+self.v, self.w),
-                   ray().ov(self.o, self.v+self.w),
-                   ray().ov(self.o, self.v-self.w), )
-
-    def transform(self, m):
-        return plane(m @ self.o, m @ self.v, m @ self.w)
-        
-    def move(self, v):
-        return plane(self.o+v, self.v, self.w)
-        
-    def inside(self, p):
-        return (self.n @ (p-self.o))[0] >= -0.00001
-
-    def intersection(self, r):
-        self.M[:,2:3] = r.v
-        try:
-            c=solve(self.M, r.o-self.o)[...,2,:]
-        except Exception as e:
-            if verbose:
-                stderr.write("solve: %s: %s\n" % (str(self.M), repr(e)))
-            return -1, None
-        return c, r.o - c*r.v
-
-    def plane_intersection(self, other, limit=1e10, delta=1e-5):
-        p = [ self.intersection(r)[1] for r in other.rays]
-        p += [ other.intersection(r)[1] for r in self.rays]
-        p = [(length(pp),pp) for pp in p if pp is not None]
-        p.sort(key=lambda x: x[0])
-        while p and p[-1][0] > limit:
-            p[-1:]=[]
-        while len(p) > 2:
-            if length(p[1][1]-p[0][1]) < delta*p[1][0]:
-                p[:1]=[]
-                continue
-            if p[1][0] < delta*p[-1][0]:
-                p[:1]=[]
-                continue
-            break
-        if len(p)<2:
-            return None
-        ry = ray().ov(p[0][1], p[1][1]-p[0][1])
-        return ry
-        
-    def planes_intersection(self, other, *more):
-        ry = self.plane_intersection(other)
-        if not ry:
-            return []
-        p = [pl.intersection(ry)[1] for pl in more]
-        return p
-    
-    def __repr__(self):
-        return "plane[%d](o=(%g,%g,%g), v=(%g,%g,%g), w=(%g,%g,%g))" % (
-            (self.idx,)+vtuple(self.o)+vtuple(self.v)+vtuple(self.w))
-
-class detector(list):
-
-    def hitplane(self, r, i):
-        c, p = self[i].intersection(r)
-        if p is None:
-            return False, 0
-        h = abs(c) < 1e10
-        for j, pl in enumerate(self):
-            if j != i:
-                inside = pl.inside(p)
-                h = numpy.logical_and(h, inside, out=h)
-        return h, p
-
-    def inside(self, p):
-        h = (p.T[...,newaxis,:] @ p.T[...,newaxis])[...,0,0] < 1e20
-        for pl in self:
-            inside = pl.inside(p)
-            h = numpy.logical_and(h, inside, out=h)
-        return h
-
-    def hit(self, r):
-        h=[]
-        for i in range(len(self)):
-            f,p = self.hitplane(r, i)
-            h.append((f, p))
-        return h
-
-    def pathlength(self, r):
-        h=[p for f,p in self.hit(r) if f]
-        if len(h)<2:
-            return 0
-        d = h[1]-h[0]
-        return sqrt((d.T @ d)[0,0])
-
-    def pathlength_beam(self, r):
-        h=self.hit(r)
-        c = [cc for cc,pp in h if cc is not None]
-        p = [pp for cc,pp in h if cc is not None]
-        p1 = numpy.select(c,p)
-        c.reverse()
-        p.reverse()
-        p2 = numpy.select(c,p)
-        d = p2 - p1
-        return sqrt((d.T[...,newaxis,:] @ d.T[...,newaxis])[...,0,0])
-    
-    def transform(self, m):
-        return detector([p.transform(m) for p in self])
-        
-    def move(self, v):
-        return detector([p.move(v) for p in self])
-
-    def corners(self):
-        c=[]
-        for i in range(len(self)-2):
-            for j in range(i+1,len(self)-1):
-                c.extend(self[i].planes_intersection(self[j], *self[j+1:]))
-        c = array([cc for cc in c if cc is not None])[...,0].T
-        c = c[:,self.inside(c)].T[...,newaxis]
-        return c
-
-    def edges(self):
-        ee=[]
-        for i in range(len(self)-1):
-            for j in range(i+1, len(self)):
-                ry = self[i].plane_intersection(self[j])
-                pp=[]
-                if ry is not None:
-                    for k in range(len(self)):
-                        if k!=i and k!=j:
-                            c,p = self[k].intersection(ry)
-                            if p is not None and length(p)<1e10 and self.inside(p):
-                                pp.append(p)
-                if len(pp) >= 2:
-                    ee.append(pp)
-        return ee
-
-    def gnuplot(self, edjes=None):
-        if edjes is None:
-            edjes = self.edges()
-        for e in edjes:
-            stdout.write("%g %g %g\n%g %g %g\n\n\n" % (tuple(e[0])+tuple(e[1])))
-
-    def planes(self):
-        # is this useful?
-        self.o = numpy.stack([p.o for p in self])
-        self.n = numpy.stack([p.n for p in self])
-        self.v = numpy.stack([p.v for p in self])
-        self.w = numpy.stack([p.w for p in self])
-        return self.o, self.n, self.v, self.w
-
-deg=π/180
-
-def box(x=1.0, y=1.0, z=1.0):
-    return detector([
-        plane(vector(z=-z), vector(x=1), vector(y=1)),
-        plane(vector(z= z), vector(y=1), vector(x=1)),
-        plane(vector(x=-x), vector(y=1), vector(z=1)),
-        plane(vector(x= x), vector(z=1), vector(y=1)),
-        plane(vector(y=-y), vector(z=1), vector(x=1)),
-        plane(vector(y= y), vector(x=1), vector(z=1)),
-    ])
-
-def prism(n=6, a=1.0, z=1.0, φ=0.0):
-    p = [
-        plane(vector(z=-z), vector(x=1), vector(y=1)),
-        plane(vector(z= z), vector(y=1), vector(x=1)),
-    ]
-    φ = 2*π * arange(int(n))/n + φ
-    p.extend([
-        plane(
-            vector(a*c, a*s),
-            vector(  s,  -c),
-            vector(z=1)
-        )
-        for c,s in zip(cos(φ), sin(φ))
-    ])
-    return detector(p)
-
-µM = box(41., 41., 13.)
-
-def run(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
-    fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
-    nr = int(r/Δr + 1.001)
-    rφ = 0;
-    for iθ in range(int((θ2-θ1)/Δθ + 1.001)):
-        θ = θ1 + iθ * Δθ
-        nnφ = (φ2-φ1)*(cos(θ)-cos(θ+Δθ))/Δθ**2 + rφ
-        nφ = max(1, int(nnφ+0.5))
-        rφ = nnφ - nφ
-        Δφ = (φ2-φ1)/nφ
-        for iφ in range(nφ):
-            φ = φ1 + iφ * Δφ + Δφ/2
-            rψ = 0
-            for ir in range(nr):
-                r = ir * Δr
-                nnψ = 2*π*ir + 0.5 + rψ
-                nψ = int(nnψ + 0.5)
-                rψ = nnψ - nψ
-                dψ = 2*π/nψ
-                for iψ in range(nψ):
-                    ψ = iψ * dψ + dψ/2
-                    ry = ray(θ, φ, r, ψ)
-                    pl = [d.pathlength(ry) for d in world]
-                    stdout.write(fmt % ((θ, φ, r, ψ)+tuple(pl)))
-
-def run_beam(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
-    fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
-    beam = rays(r, Δr)
-    b = None
-    rφ = 0;
-    for iθ in range(int((θ2-θ1)/Δθ + 1.001)):
-        θ = θ1 + iθ * Δθ
-        nnφ = (φ2-φ1)*(cos(θ)-cos(θ+Δθ))/Δθ**2 + rφ
-        nφ = max(1, int(nnφ+0.5))
-        rφ = nnφ - nφ
-        Δφ = (φ2-φ1)/nφ
-        for iφ in range(nφ):
-            φ = φ1 + iφ * Δφ + Δφ/2
-            beam.bend(θ, φ)
-            pl = [d.pathlength_beam(beam) for d in world]
-            for i in range(beam.n):
-                bb=[θ, φ, beam.r[i], beam.ψ[i]]
-                ll = tuple(bb+[l[i] for l in pl])
-                stdout.write(fmt % ll)
-        
-test_rays = [
-    ("x-axis", ray().ov(vector(),vector(x=1))),
-    ("y-axis", ray().ov(vector(),vector(y=1))),
-    ("z-axis", ray().ov(vector(),vector(z=1))),
-]
-
-def test(det, rays):
-    for name, ry in rays:
-        pl = [d.pathlength(ry) for d in det]
-        stderr.write(("test %s:"+len(pl)*" %g"+"\n") % ((name,)+tuple(pl)))
-
-def corners(world):
-    c = [d.corners() for d in world]
-    ll=[]
-    for i,d in enumerate(c):
-        stderr.write("Detector[%d]\n" %i)
-        for cc in d:
-            if cc is not None:
-                l = length(cc)
-                ll.append(l)
-                stderr.write(" |%g, %g, %g| = %g\n" % (cc[0,0],cc[1,0],cc[2,0],l))
-    max_r = max(ll)
-    return max_r, c
-        
-short_opt = "hSNTREGU:D:C:B:P:µMr:d:a:p:x:"
-long_opt = ["help",
-            "φ1=", "φ2=", "θ1=", "θ2=", "φ₁=", "φ₂=", "θ₁=", "θ₂=",
-            "Δθ=", "µMustang", "test", "size", "run", "corners", "gnuplot",
-            "radius=", "resolution=", "Δr=" ,
-            "detector", "copy=", "select=", "box=", "prism=",
-            "plane=", "ray=",
-            "rx=", "ry=", "rz=", "move=", "scale=",
-            "define=", "replace=", "undefine="]
-
-args = argv[1:]
-options = []
-while args:
-    oo, aa = getopt(args, short_opt, long_opt)
-    args = []
-    options.extend(oo)
-    if aa:
-        fn = aa[0]
-        with open(fn) as f:
-            for l in f:
-                l=l.strip()
-                if not l:
-                    continue
-                if l[0]=='-':
-                    args.append(l)
-                elif l[0]=='@':
-                    args.append(l[1:])
-                elif l[0] != "#":
-                    args.append("--"+l)
-        args.extend(aa[1:])
-
-Δr = 1.0
-Δθ = 1*deg
-r = 0
-world = [detector()]
-θ1=0
-θ2=π/2
-φ1=0
-φ2=2*π/8
-do_size = False
-do_test = False
-do_run = False
-do_corners = False
-do_gnuplot = False
-rot = rotate(0,0)
-ori = vector()
-
-macros={}
-
-def aval(v):
-    vv=""
-    while v!=vv and v.find("$") >= 0:
-        vv=v
-        for m in macros:
-            v = v.replace(m, macros[m])
-    v = v.replace("°", "*deg")
-    return eval(v)
-
-def avec(v):
-    v = v.strip().replace(",", " ")
-    return tuple([aval(vv) for vv in v.split()])
-
-for o,v in options:
-
-    if o in "-h --help":
-        with open("geometryfactor.md") as f:
-            stderr.write(f.read())
-        raise SystemExit(0)
-    
-    if o in "-µ -M --µMustang":
-        if world[-1]:
-            world.append(µM)
-        else:
-            world[-1] = µM
-    if o in "-p --plane":
-        world[-1].append(plane(*tuple([vector(*avec(vv)) for vv in v.split(";")])))
-    if o in "-x --ray":
-        vv = v.split(";")
-        if len(vv)==2:
-            test_rays.append((vv[0], ray(*avec(vv[1]))))
-        else:
-            test_rays.append((vv[0], ray().ov(vector(*avec(vv[1])), vector(*avec(vv[2])))))
-    if o in "-a --Δθ":
-        Δθ = aval(v)
-    if o in "-d --Δr --resolution":
-        Δr = aval(v)
-    if o in "--θ1 --θ₁":
-        θ1 = aval(v)
-    if o in "--θ2 --θ₂":
-        θ2 = aval(v)
-    if o in "--φ1 --φ₁":
-        φ1 = aval(v)
-    if o in "--φ2 --φ₂":
-        φ2 = aval(v)
-    if o in "-r --radius":
-        r = aval(v)
-    if o in "-S --size":
-        do_size=True
-    if o in "-T --test":
-        do_test=True
-    if o in "-R --run":
-        do_run=True
-    if o in "-G --gnuplot":
-        do_gnuplot=True
-    if o in "-E --corners":
-        do_corners=True
-    if o in "--rx":
-        world[-1] = world[-1].transform(rotate(0, aval(v)))
-    if o in "--ry":
-        world[-1] = world[-1].transform(rotate(1, aval(v)))
-    if o in "--rz":
-        world[-1] = world[-1].transform(rotate(2, aval(v)))
-    if o in "--move":
-        world[-1] = world[-1].move(vector(*avec(v)))
-    if o in "--scale":
-        world[-1] = world[-1].transform(scale(*avec(v)))
-    if o in "-N --detector":
-        if world[-1]:
-            world.append(detector())
-    if o in "-C --copy --select":
-        if not world[-1]:
-            world[-1:] = []
-        idx=int(v)
-        world.append(world[idx])
-        if o in "--select":
-            world[idx:idx+1]=[]
-    if o in "-B --box":
-        if not world[-1]:
-            world[-1:] = []
-        world.append(box(*avec(v)))
-    if o in "-P --prism":
-        if not world[-1]:
-            world[-1:] = []
-        world.append(prism(*avec(v)))
-        
-    if o in "-U --undefine -D --define --replace":
-        vv=v.split("=",1)
-        k = "${"+vv[0]+"}"
-        if o in "-U --undefine":
-            if k in macros:
-                stderr.write("undefine macro %s=%s\n" % (k,macros[k]))
-                del macros[k]
-        elif o in "--replace" or not k in macros:
-            if k in vv[1]:
-                raise ValueError("recursive macro %s" % v)
-            macros[k]=vv[1]
-            stderr.write("define macro %s=%s\n" % (k,vv[1]))
-
-if not world[-1]:
-    world[-1:]=[]
-if not world:
-    world=[µM]
-
-if not r or do_corners:
-    max_r, c = corners(world)
-    if not r:
-        r = max_r
-
-if do_gnuplot:
-    for d in world:
-        d.gnuplot()
-
-if do_size:
-    nr = int(r/Δr+1.001)
-    no = π*nr**2 - (π-0.5)*nr
-    sr = (φ2-φ1)*(cos(θ1)-cos(θ2+Δθ))
-    A = π*((nr-0.5)*Δr)**2
-    nv = sr/Δθ**2
-    nn = int(nv)*int(no)
-    stderr.write("""Run size estimate:
-detectors: %d  planes: %d
-r = %g
-Δr = %g
-Δθ = %g = %g°, θ=%g…%g φ=%g…%g
-ω = %g sr = %g π sr
-A = %g ⧠
-n = %.0f × %.0f = %.0f = %g /(sr ⧠)
-""" % (
-    len(world), sum([len(d) for d in world]),
-    r, Δr,
-    Δθ, Δθ/deg, θ1, θ2, φ1, φ2,
-    sr, sr/π, A,
-    nv, no, nn, nn/sr/A 
-))
-
-if do_test:
-    test(world, test_rays)
-if do_run:
-    run_beam(world, r, Δr, Δθ, θ1, θ2, φ1, φ2)

journey.dat

@@ -0,0 +1,65 @@
+0  0  0 -0.0 0.0 0.0
+0  0  0 0.43301270189221935 0.24999999999999997 0.0
+0  0  0 3.061616997868383e-17 0.5 0.0
+0  0  0 -0.43301270189221924 0.25000000000000017 0.0
+0  0  0 -0.4330127018922194 -0.2499999999999999 0.0
+0  0  0 -9.184850993605148e-17 -0.5 0.0
+0  0  0 0.4330127018922192 -0.2500000000000002 0.0
+0  0  0 0.970941817426052 0.23931566428755774 0.0
+0  0  0 0.7485107481711012 0.6631226582407952 0.0
+0  0  0 0.3546048870425358 0.9350162426854147 0.0
+0  0  0 -0.1205366802553229 0.992708874098054 0.0
+0  0  0 -0.5680647467311556 0.8229838658936566 0.0
+0  0  0 -0.8854560256532096 0.46472317204376906 0.0
+0  0  0 -1.0 5.66553889764798e-16 0.0
+0  0  0 -0.8854560256532101 -0.46472317204376806 0.0
+0  0  0 -0.5680647467311559 -0.8229838658936564 0.0
+0  0  0 -0.12053668025532356 -0.992708874098054 0.0
+0  0  0 0.35460488704253473 -0.9350162426854152 0.0
+0  0  0 0.7485107481711009 -0.6631226582407955 0.0
+0  0  0 0.9709418174260518 -0.23931566428755865 0.0
+0  0  0 1.4815325108927067 0.23465169756034632 0.0
+0  0  0 1.336509786282552 0.6809857496093201 0.0
+0  0  0 1.0606601717798214 1.0606601717798212 0.0
+0  0  0 0.6809857496093202 1.3365097862825517 0.0
+0  0  0 0.23465169756034637 1.4815325108927067 0.0
+0  0  0 -0.2346516975603462 1.4815325108927067 0.0
+0  0  0 -0.6809857496093201 1.336509786282552 0.0
+0  0  0 -1.0606601717798212 1.0606601717798214 0.0
+0  0  0 -1.3365097862825517 0.6809857496093203 0.0
+0  0  0 -1.4815325108927064 0.23465169756034648 0.0
+0  0  0 -1.4815325108927067 -0.23465169756034615 0.0
+0  0  0 -1.336509786282552 -0.68098574960932 0.0
+0  0  0 -1.0606601717798214 -1.0606601717798212 0.0
+0  0  0 -0.6809857496093203 -1.3365097862825517 0.0
+0  0  0 -0.2346516975603466 -1.4815325108927064 0.0
+0  0  0 0.234651697560346 -1.4815325108927067 0.0
+0  0  0 0.6809857496093199 -1.336509786282552 0.0
+0  0  0 1.0606601717798212 -1.0606601717798214 0.0
+0  0  0 1.3365097862825517 -0.6809857496093205 0.0
+0  0  0 1.4815325108927064 -0.23465169756034668 0.0
+0  0  0 1.9842294026289558 0.2506664671286085 0.0
+0  0  0 1.8595529717765027 0.736249105369356 0.0
+0  0  0 1.618033988749895 1.1755705045849463 0.0
+0  0  0 1.2748479794973793 1.5410264855515785 0.0
+0  0  0 0.8515585831301453 1.8096541049320392 0.0
+0  0  0 0.3747626291714494 1.9645745014573772 0.0
+0  0  0 -0.1255810390586268 1.9960534568565431 0.0
+0  0  0 -0.6180339887498951 1.902113032590307 0.0
+0  0  0 -1.0716535899579935 1.68865585100403 0.0
+0  0  0 -1.4579372548428233 1.369094211857377 0.0
+0  0  0 -1.7526133600877274 0.9635073482034304 0.0
+0  0  0 -1.9371663222572624 0.49737977432970876 0.0
+0  0  0 -2.0 -6.432490598706546e-16 0.0
+0  0  0 -1.9371663222572622 -0.49737977432971003 0.0
+0  0  0 -1.7526133600877265 -0.9635073482034315 0.0
+0  0  0 -1.4579372548428222 -1.3690942118573781 0.0
+0  0  0 -1.0716535899579926 -1.6886558510040306 0.0
+0  0  0 -0.6180339887498935 -1.9021130325903075 0.0
+0  0  0 -0.1255810390586264 -1.9960534568565431 0.0
+0  0  0 0.37476262917145026 -1.9645745014573772 0.0
+0  0  0 0.8515585831301453 -1.8096541049320392 0.0
+0  0  0 1.27484797949738 -1.541026485551578 0.0
+0  0  0 1.6180339887498958 -1.1755705045849452 0.0
+0  0  0 1.859552971776503 -0.7362491053693556 0.0
+0  0  0 1.9842294026289558 -0.2506664671286076 0.0

map_origin.py

@@ -0,0 +1,110 @@
+import sys
+import argparse
+import numpy as np
+import random
+
+
+def read_csv(file_path, slope_sc, Ind_A, Ind_S, slope_stop, x, y, z):
+    file = open(file_path, newline='')
+    det1 = 0
+    det2 = 0
+    for line in file:
+        row = line.strip().split(',')
+        if float(row[0]) == x and float(row[1]) == y and float(row[2]) == z:
+            an = []
+            ta = np.arcsin(Ind_A/Ind_S)
+            for i in range(3):
+                x = float(row[13+8*i])
+                y = float(row[14+8*i])
+                z = float(row[15+8*i])
+                v=np.array((x,y,z))
+                norm = np.array((0,1,0))
+                an.append(np.arccos(abs(np.dot(v,norm))/(np.linalg.norm(norm)*np.linalg.norm(v))))
+            if int(float(row[8])) == 1 :
+                u = random.random()
+                v = random.random()
+                sc = np.exp(-slope_sc*(int(float(row[11]))))
+                p = np.exp(-slope_stop*(int(float(row[12]))))
+                if u <= sc and v <= p and an[0] < ta:
+                    if int(float(row[9])) == 1:
+                        det1 += 1
+                    else:
+                        det2 += 1
+            elif int(float(row[8])) == 2:
+                u = random.random()
+                u2 = random.random()
+                sc = np.exp(-slope_sc*(int(float(row[11]))))
+                sc2 = np.exp(-slope_sc*int(float(row[19])))
+                v = random.random()
+                v2 = random.random()
+                p = np.exp(-slope_stop*(int(float(row[12]))))
+                p2 = np.exp(-slope_stop*(int(float(row[20]))))
+                if u <= sc and v <= p and an[0] < ta:
+                    if int(float(row[9])) == 1:
+                        det1 += 1
+                    else:
+                        det2 += 1
+                elif u2 <= sc2 and v2 <= p2 and an[1] < ta:
+                    if int(float(row[17])) == 1:
+                        det1 += 1
+                    else:
+                        det2 += 1
+            elif int(float(row[8])) == 3:
+                u = random.random()
+                u2 = random.random()
+                u3 = random.random()
+                sc = np.exp(-slope_sc*(int(float(row[11]))))
+                sc2 = np.exp(-slope_sc*int(float(row[19])))
+                sc3 = np.exp(-slope_sc*int(float(row[27])))
+                v = random.random()
+                v2 = random.random()
+                v3 = random.random()
+                p = np.exp(-slope_stop*(int(float(row[12]))))
+                p2 = np.exp(-slope_stop*(int(float(row[20]))))
+                p3 = np.exp(-slope_stop*(int(float(row[28]))))
+                if u <= sc and v <= p and an[0] < ta:
+                    if int(float(row[9])) == 1:
+                        det1 += 1
+                    else:
+                        det2 += 1 
+                elif u2 <= sc2 and v2 <= p2 and an[1] < ta:
+                    if int(float(row[17])) == 1:
+                        det1 += 1
+                    else:
+                        det2 += 1 
+                elif u3 <= sc3 and v3 <= p3 and an[2] < ta:
+                    if int(float(row[25])) == 1:
+                        det1 += 1
+                    else:
+                        det2 += 1          
+    return det1, det2
+
+
+def main(csv_file, slope_sc, stop_slope):
+    slope_sc = float(slope_sc)
+    stop_slope = float(stop_slope)
+    x, y, z = map(float, args.origin.split(','))
+    det1, det2 = read_csv(csv_file, slope_sc, 1, 2.15, stop_slope, x, y, z) 
+    output = f"{det1} {det2}\n"  
+    sys.stdout.write(output)
+    
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(prog = 'filter_not')
+
+    parser.add_argument('-c', '--csv_file', help="Path to the CSV file")
+
+    parser.add_argument('--stop_slope', default = 0.0000001, help = 'parameter for the exponential decay function that determines after what path length the ray is lost')
+
+    parser.add_argument('--sc_slope', default = 0.00000005, help = 'parameter for the exponential decay function that determines after how many scatter events the ray is lost')
+
+    parser.add_argument('-o', '--origin', help = 'origin to look at')
+    args = parser.parse_args()
+    
+    # Call the main function with parsed arguments
+    main(args.csv_file, args.sc_slope, args.stop_slope)
+
+
+
+
+

map_trigger.py

@@ -0,0 +1,111 @@
+import sys
+import argparse
+import numpy as np
+import random
+
+def read_csv(file_path, slope_sc, Ind_A, Ind_S, slope_stop):
+    file = open(file_path, newline='')
+    det1 = 0
+    det2 = 0
+    det1_first = 0
+    det2_first = 0
+    for line in file:
+        row = line.strip().split(',')
+        if int(float(row[9])) == 1:
+            det1_first += 1
+        elif int(float(row[9])) == 2:
+            det2_first += 1
+        an = []
+        ta = np.arcsin(Ind_A/Ind_S)
+        for i in range(3):
+            x = float(row[13+8*i])
+            y = float(row[14+8*i])
+            z = float(row[15+8*i])
+            v=np.array((x,y,z))
+            norm = np.array((0,1,0))
+            an.append(np.arccos(abs(np.dot(v,norm))/(np.linalg.norm(norm)*np.linalg.norm(v))))
+        if int(float(row[8])) == 1 :
+            u = random.random()
+            v = random.random()
+            sc = np.exp(-slope_sc*(int(float(row[11]))))
+            p = np.exp(-slope_stop*(int(float(row[12]))))
+            if u <= sc and v <= p and an[0] < ta:
+                if int(float(row[9])) == 1:
+                    det1 += 1
+                else:
+                    det2 += 1
+        elif int(float(row[8])) == 2:
+            u = random.random()
+            u2 = random.random()
+            sc = np.exp(-slope_sc*(int(float(row[11]))))
+            sc2 = np.exp(-slope_sc*int(float(row[19])))
+            v = random.random()
+            v2 = random.random()
+            p = np.exp(-slope_stop*(int(float(row[12]))))
+            p2 = np.exp(-slope_stop*(int(float(row[20]))))
+            if u <= sc and v <= p and an[0] < ta:
+                if int(float(row[9])) == 1:
+                    det1 += 1
+                else:
+                    det2 += 1
+            elif u2 <= sc2 and v2 <= p2 and an[1] < ta:
+                if int(float(row[17])) == 1:
+                    det1 += 1
+                else:
+                    det2 += 1
+        elif int(float(row[8])) == 3:
+            u = random.random()
+            u2 = random.random()
+            u3 = random.random()
+            sc = np.exp(-slope_sc*(int(float(row[11]))))
+            sc2 = np.exp(-slope_sc*int(float(row[19])))
+            sc3 = np.exp(-slope_sc*int(float(row[27])))
+            v = random.random()
+            v2 = random.random()
+            v3 = random.random()
+            p = np.exp(-slope_stop*(int(float(row[12]))))
+            p2 = np.exp(-slope_stop*(int(float(row[20]))))
+            p3 = np.exp(-slope_stop*(int(float(row[28]))))
+            if u <= sc and v <= p and an[0] < ta:
+                if int(float(row[9])) == 1:
+                    det1 += 1
+                else:
+                    det2 += 1 
+            elif u2 <= sc2 and v2 <= p2 and an[1] < ta:
+                if int(float(row[17])) == 1:
+                    det1 += 1
+                else:
+                    det2 += 1 
+            elif u3 <= sc3 and v3 <= p3 and an[2] < ta:
+                if int(float(row[25])) == 1:
+                    det1 += 1
+                else:
+                    det2 += 1      
+    return det1, det2, det1_first, det2_first
+
+def main(csv_file, slope_sc, stop_slope):
+    slope_sc = float(slope_sc)
+    stop_slope = float(stop_slope)
+    det1, det2, det1_first, det2_first= read_csv(csv_file, slope_sc, 1, 2.15, stop_slope) 
+    output = f"{det1} {det2} {det1_first} {det2_first} \n"  
+    sys.stdout.write(output)
+    
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(prog = 'filter_not')
+
+    parser.add_argument('-c', '--csv_file', help="Path to the CSV file")
+
+    parser.add_argument('--stop_slope', default = 0.0000001, help = 'parameter for the exponential decay function that determines after what path length the ray is lost')
+
+    parser.add_argument('--sc_slope', default = 0.00000005, help = 'parameter for the exponential decay function that determines after how many scatter events the ray is lost')
+
+    args = parser.parse_args()
+    
+    # Call the main function with parsed arguments
+    main(args.csv_file, args.sc_slope, args.stop_slope)
+
+
+
+
+

mustang.geometryfaktoren

@@ -1,256 +0,0 @@
-0 0 1341031.0 0.20290814892301096 0.05069796467951925 0.7972617590187603 0.8422581322330368
-0 1 1986080.0 0.4012516445497979 0.17904959352503397 1.2125982829906625 1.5510665850409142
-0 2 2040563.0 0.41529336939431527 0.1906563176744989 0.35548935921354874 0.20272691823122851
-0 3 3161775.0 0.5376438689389628 0.30547929261632606 0.7719350628713082 0.6825935382885188
-0 4 823384.0 0.7485357838038929 0.5686924262663288 1.4204535944616108 2.0299228484791847
-0 5 379835.0 0.936743837987913 0.8813951180553474 1.4671931586085913 2.158272764133479
-0 6 1422202.0 0.804091467508231 0.6540279786186627 1.14316421159702 1.3349474398490124
-0 7 690204.0 0.9671109603755198 0.9388755429776542 1.263596434088368 1.6114337793213878
-0 8 711040.0 0.75541376792732 0.5777127022180836 0.1545495044655409 0.036601216775153066
-0 9 1211133.0 0.8025648909121936 0.6503820078769514 0.4206968529282656 0.20591680701622722
-0 10 371584.0 0.9473615246954701 0.9013347804118663 0.10409242232948127 0.016434403026936266
-0 11 672570.0 0.9695401945575313 0.9434898225076767 0.2980007875472181 0.10353684749953741
-0 12 703322.0 0.9351701698111494 0.8781723538386255 0.787471032713443 0.6371009574940887
-0 13 408083.0 1.0388691224712987 1.0812803448590638 0.9652181380749233 0.9428738324934574
-0 14 466318.0 1.0401597892686825 1.0843108310750373 0.6031497703239409 0.3747726957166929
-0 15 313018.0 1.1052087943294517 1.223039903370744 0.7815175401034282 0.6190848774927195
-1 0 0.0 nan nan nan nan
-1 1 1341214.0 0.20289792039478732 0.050692150656630956 0.7972699518589962 0.842253407954774
-1 2 0.0 nan nan nan nan
-1 3 2041677.0 0.4153166543389886 0.19067669833027134 0.3554723553468729 0.20272169459781292
-1 4 1709275.0 0.4796552675096441 0.24637731911583272 1.2927698630391917 1.7170640967677737
-1 5 823397.0 0.7485763533478835 0.5687542968950995 1.420406879238328 2.029792562195492
-1 6 2767358.0 0.5923264896931107 0.3654121730924992 0.8860582811184025 0.8515371690683806
-1 7 1421661.0 0.8040734998224996 0.6540011409929922 1.1432010306350813 1.3350302739450746
-1 8 0.0 nan nan nan nan
-1 9 712147.0 0.7554206376170401 0.5777217408633354 0.15462228111105142 0.036640224013440414
-1 10 0.0 nan nan nan nan
-1 11 372210.0 0.9473623135123402 0.9013364176076847 0.10409344509752312 0.016435398657037498
-1 12 1113920.0 0.8262367461201479 0.688466188161221 0.510981785437805 0.2872548770884155
-1 13 702488.0 0.9351600041294312 0.8781537164231598 0.7874236963735674 0.6370324673592848
-1 14 642566.0 0.9811683563316167 0.9659900472917671 0.36722076460139086 0.14894016013587774
-1 15 466385.0 1.0401727567485883 1.0843373103573553 0.6031258160337903 0.3747522357601341
-2 0 0.0 nan nan nan nan
-2 1 0.0 nan nan nan nan
-2 2 1483887.0 0.21420903279744924 0.05654864969081774 0.7760866291640545 0.8067038414473829
-2 3 2170029.0 0.4125478688548589 0.1891182621223371 1.1977331447953943 1.5179035011458824
-2 4 0.0 nan nan nan nan
-2 5 0.0 nan nan nan nan
-2 6 867310.0 0.7623715248380875 0.5898226519088194 1.4165618098295538 2.0191579872124996
-2 7 390740.0 0.9503211647110371 0.907045866911589 1.4659307811082531 2.1546119575675986
-2 8 1803401.0 0.4781970594623001 0.24542674909226228 0.29078043947273247 0.13355497370804023
-2 9 2837825.0 0.5822060175961518 0.35383128930328717 0.6813836311835574 0.5341219822019206
-2 10 986867.0 0.7567279062350565 0.5824294157441062 0.15149317007109525 0.035163562602109816
-2 11 1676737.0 0.8023338885031722 0.6523569025525789 0.41516778500156915 0.20049790810800733
-2 12 1328583.0 0.8228601582532988 0.6840977338838917 1.0710464579549914 1.1730772844158073
-2 13 665660.0 0.9765751389116561 0.9571169469413244 1.2078516437325386 1.4731367969129678
-2 14 964915.0 0.9332167006125296 0.8759244958424913 0.7828853851681628 0.6294000750942009
-2 15 555370.0 1.0370948572416403 1.0783886784601142 0.9624207572790557 0.9371083540928041
-3 0 0.0 nan nan nan nan
-3 1 0.0 nan nan nan nan
-3 2 0.0 nan nan nan nan
-3 3 1485205.0 0.2142024499279579 0.056547952301451335 0.7760513358014539 0.8066949326617141
-3 4 0.0 nan nan nan nan
-3 5 0.0 nan nan nan nan
-3 6 1857507.0 0.49142819234033447 0.25856983581410053 1.2815077719822965 1.6898110332998868
-3 7 867286.0 0.7623600308068184 0.5898020869022559 1.4166371964014903 2.019369568538402
-3 8 0.0 nan nan nan nan
-3 9 1806590.0 0.47821051352570443 0.24543980143638666 0.29080430320399286 0.133586135307775
-3 10 0.0 nan nan nan nan
-3 11 988020.0 0.7567298850635752 0.5824319172235912 0.1515086517322584 0.03516553762298842
-3 12 2521393.0 0.6300973010096894 0.41036799508776395 0.7953013227413666 0.6883604948470025
-3 13 1328571.0 0.8228843269382178 0.6841370914362079 1.0709278899849832 1.1728196231717225
-3 14 1540476.0 0.8254724239811463 0.6893715726573858 0.5052265879823434 0.2806923785854828
-3 15 965758.0 0.93326028703934 0.8760002666179418 0.782886629834323 0.6293988297406045
-4 0 0.0 nan nan nan nan
-4 1 0.0 nan nan nan nan
-4 2 0.0 nan nan nan nan
-4 3 0.0 nan nan nan nan
-4 4 1341140.0 0.20290214683051372 0.05069742967272359 0.7973013486984952 0.8423233901510625
-4 5 1985816.0 0.40125667947058397 0.17904839349581256 1.2125802803851613 1.5510467952108336
-4 6 2041246.0 0.41530423954343 0.19066219328248782 0.3554605411333851 0.20269758488162645
-4 7 3161646.0 0.5376294985004041 0.30546698395273514 0.7719764775018184 0.6826719959096396
-4 8 0.0 nan nan nan nan
-4 9 0.0 nan nan nan nan
-4 10 0.0 nan nan nan nan
-4 11 0.0 nan nan nan nan
-4 12 712066.0 0.7554322158396256 0.5777386878976862 0.15457836283102516 0.036619204825268335
-4 13 1210980.0 0.8025423350816383 0.6503471982687781 0.4206924369214766 0.20592286263842452
-4 14 372338.0 0.9473631656248545 0.9013388819857729 0.10406663330417283 0.016430633847943955
-4 15 672542.0 0.9695160708930571 0.9434422044055087 0.29806066907826384 0.10357149736641522
-5 0 0.0 nan nan nan nan
-5 1 0.0 nan nan nan nan
-5 2 0.0 nan nan nan nan
-5 3 0.0 nan nan nan nan
-5 4 0.0 nan nan nan nan
-5 5 1340884.0 0.20288715357072348 0.05068884172493339 0.7973341890132793 0.8423767617685847
-5 6 0.0 nan nan nan nan
-5 7 2040561.0 0.41528021537458415 0.19064749547267573 0.35548413741357815 0.20271750023782714
-5 8 0.0 nan nan nan nan
-5 9 0.0 nan nan nan nan
-5 10 0.0 nan nan nan nan
-5 11 0.0 nan nan nan nan
-5 12 0.0 nan nan nan nan
-5 13 711144.0 0.755402040700679 0.5776975819861864 0.154594060636747 0.03662526032080497
-5 14 0.0 nan nan nan nan
-5 15 371550.0 0.9473392597320444 0.9012937741184909 0.10411087145062095 0.016440109895515337
-6 0 0.0 nan nan nan nan
-6 1 0.0 nan nan nan nan
-6 2 0.0 nan nan nan nan
-6 3 0.0 nan nan nan nan
-6 4 0.0 nan nan nan nan
-6 5 0.0 nan nan nan nan
-6 6 1485316.0 0.21421956705300368 0.05655426326747592 0.7760855656888564 0.8067465373150426
-6 7 2169948.0 0.4125482719439972 0.18912020954765354 1.1977532014177814 1.5179467620778604
-6 8 0.0 nan nan nan nan
-6 9 0.0 nan nan nan nan
-6 10 0.0 nan nan nan nan
-6 11 0.0 nan nan nan nan
-6 12 1806442.0 0.47822045655459805 0.24544342960384483 0.2908183190019047 0.13361337325813688
-6 13 2837787.0 0.5821946493802479 0.3538249226038536 0.6814016445517569 0.5341453281118242
-6 14 988472.0 0.7567417229223797 0.5824521133937383 0.15148738316695803 0.035164516666261934
-6 15 1677734.0 0.8023135775416316 0.6523257740521955 0.4152626702966833 0.2005939130625495
-7 0 0.0 nan nan nan nan
-7 1 0.0 nan nan nan nan
-7 2 0.0 nan nan nan nan
-7 3 0.0 nan nan nan nan
-7 4 0.0 nan nan nan nan
-7 5 0.0 nan nan nan nan
-7 6 0.0 nan nan nan nan
-7 7 1484050.0 0.21421811728119416 0.0565551089168354 0.7760663854049075 0.8066874578266343
-7 8 0.0 nan nan nan nan
-7 9 0.0 nan nan nan nan
-7 10 0.0 nan nan nan nan
-7 11 0.0 nan nan nan nan
-7 12 0.0 nan nan nan nan
-7 13 1803592.0 0.47819646775954405 0.24542573564122175 0.2908328890519232 0.1336149271277042
-7 14 0.0 nan nan nan nan
-7 15 986981.0 0.7567243477305784 0.5824260030189728 0.15150355127544696 0.035166829834109785
-8 0 0.0 nan nan nan nan
-8 1 0.0 nan nan nan nan
-8 2 0.0 nan nan nan nan
-8 3 0.0 nan nan nan nan
-8 4 0.0 nan nan nan nan
-8 5 0.0 nan nan nan nan
-8 6 0.0 nan nan nan nan
-8 7 0.0 nan nan nan nan
-8 8 1341207.0 0.20289757398315392 0.05069294433281922 0.7972905736177281 0.8423092189334862
-8 9 1986872.0 0.4012517402913564 0.17904450999094404 1.2125769603002217 1.5510335132064303
-8 10 2040548.0 0.41529339007921995 0.19065326881734398 0.3554650121565654 0.20268278245104598
-8 11 3161189.0 0.5376067812129309 0.30544078402204417 0.7719621851245546 0.6826427050176065
-8 12 824620.0 0.748569775348104 0.568747149265904 1.4204299055239333 2.0298614424815105
-8 13 380597.0 0.9367378001616958 0.8813822004005688 1.4672222979689982 2.1583602518600213
-8 14 1421566.0 0.8040713285525599 0.6539938814446126 1.1432159266011541 1.3350598462437635
-8 15 690109.0 0.9670969641796117 0.9388476500594812 1.2635758422894483 1.611378585878433
-9 0 0.0 nan nan nan nan
-9 1 0.0 nan nan nan nan
-9 2 0.0 nan nan nan nan
-9 3 0.0 nan nan nan nan
-9 4 0.0 nan nan nan nan
-9 5 0.0 nan nan nan nan
-9 6 0.0 nan nan nan nan
-9 7 0.0 nan nan nan nan
-9 8 0.0 nan nan nan nan
-9 9 1342598.0 0.20291434195937144 0.05070147997738545 0.797343559970149 0.8423970104306718
-9 10 0.0 nan nan nan nan
-9 11 2041411.0 0.41530121121296176 0.19065946589836497 0.3554877901486851 0.20272047565439882
-9 12 1712978.0 0.47963750146320894 0.24636413289530398 1.2927151926862892 1.7169784950863791
-9 13 824445.0 0.7485790525620942 0.5687612287092344 1.4203891835115015 2.029741822146791
-9 14 2767363.0 0.592292157254801 0.36537457583347166 0.8860514754527217 0.851532350431903
-9 15 1421914.0 0.8040530539660952 0.653969025364184 1.1432695545073235 1.335161573543511
-10 0 0.0 nan nan nan nan
-10 1 0.0 nan nan nan nan
-10 2 0.0 nan nan nan nan
-10 3 0.0 nan nan nan nan
-10 4 0.0 nan nan nan nan
-10 5 0.0 nan nan nan nan
-10 6 0.0 nan nan nan nan
-10 7 0.0 nan nan nan nan
-10 8 0.0 nan nan nan nan
-10 9 0.0 nan nan nan nan
-10 10 1483386.0 0.21419278309028092 0.05654190606941044 0.7761500649018663 0.8067964968284699
-10 11 2168677.0 0.4125304255109409 0.18910927127520574 1.1977652146587263 1.5179648070405773
-10 12 0.0 nan nan nan nan
-10 13 0.0 nan nan nan nan
-10 14 865967.0 0.7623644231608726 0.5898088423537944 1.4165788743327292 2.019203647747152
-10 15 389985.0 0.9503608076572463 0.9071237659294952 1.4659043921951258 2.1545306307738503
-11 0 0.0 nan nan nan nan
-11 1 0.0 nan nan nan nan
-11 2 0.0 nan nan nan nan
-11 3 0.0 nan nan nan nan
-11 4 0.0 nan nan nan nan
-11 5 0.0 nan nan nan nan
-11 6 0.0 nan nan nan nan
-11 7 0.0 nan nan nan nan
-11 8 0.0 nan nan nan nan
-11 9 0.0 nan nan nan nan
-11 10 0.0 nan nan nan nan
-11 11 1483979.0 0.21421628771677925 0.0565536755739456 0.7760195538692563 0.8066495028555718
-11 12 0.0 nan nan nan nan
-11 13 0.0 nan nan nan nan
-11 14 1853346.0 0.4914818400329044 0.258612364797809 1.2813482549225448 1.6894068331160028
-11 15 866088.0 0.7624076728708337 0.5898750974132595 1.4165606207998398 2.01915349134398
-12 0 0.0 nan nan nan nan
-12 1 0.0 nan nan nan nan
-12 2 0.0 nan nan nan nan
-12 3 0.0 nan nan nan nan
-12 4 0.0 nan nan nan nan
-12 5 0.0 nan nan nan nan
-12 6 0.0 nan nan nan nan
-12 7 0.0 nan nan nan nan
-12 8 0.0 nan nan nan nan
-12 9 0.0 nan nan nan nan
-12 10 0.0 nan nan nan nan
-12 11 0.0 nan nan nan nan
-12 12 1342586.0 0.2028993014260298 0.05069511543833706 0.7973447369081822 0.8423938727990237
-12 13 1987134.0 0.40128599995740827 0.17907448354808042 1.2126392341278718 1.5511592608521678
-12 14 2041179.0 0.4152797102618002 0.19064123386634735 0.3554754141607334 0.20269084587690894
-12 15 3161664.0 0.5376337548088288 0.30546968192310237 0.7719598353707505 0.6826368783514183
-13 0 0.0 nan nan nan nan
-13 1 0.0 nan nan nan nan
-13 2 0.0 nan nan nan nan
-13 3 0.0 nan nan nan nan
-13 4 0.0 nan nan nan nan
-13 5 0.0 nan nan nan nan
-13 6 0.0 nan nan nan nan
-13 7 0.0 nan nan nan nan
-13 8 0.0 nan nan nan nan
-13 9 0.0 nan nan nan nan
-13 10 0.0 nan nan nan nan
-13 11 0.0 nan nan nan nan
-13 12 0.0 nan nan nan nan
-13 13 1341292.0 0.2028914693706316 0.05068869666163884 0.7973064202947792 0.8423155210311275
-13 14 0.0 nan nan nan nan
-13 15 2040539.0 0.41527975318578586 0.1906449723853907 0.3554760992137822 0.20269854151013766
-14 0 0.0 nan nan nan nan
-14 1 0.0 nan nan nan nan
-14 2 0.0 nan nan nan nan
-14 3 0.0 nan nan nan nan
-14 4 0.0 nan nan nan nan
-14 5 0.0 nan nan nan nan
-14 6 0.0 nan nan nan nan
-14 7 0.0 nan nan nan nan
-14 8 0.0 nan nan nan nan
-14 9 0.0 nan nan nan nan
-14 10 0.0 nan nan nan nan
-14 11 0.0 nan nan nan nan
-14 12 0.0 nan nan nan nan
-14 13 0.0 nan nan nan nan
-14 14 1483833.0 0.21421043548688512 0.056551816039018675 0.7760746487063221 0.8067121722367818
-14 15 2168891.0 0.4125328102788482 0.18910686518441988 1.197757269282572 1.517945233552834
-15 0 0.0 nan nan nan nan
-15 1 0.0 nan nan nan nan
-15 2 0.0 nan nan nan nan
-15 3 0.0 nan nan nan nan
-15 4 0.0 nan nan nan nan
-15 5 0.0 nan nan nan nan
-15 6 0.0 nan nan nan nan
-15 7 0.0 nan nan nan nan
-15 8 0.0 nan nan nan nan
-15 9 0.0 nan nan nan nan
-15 10 0.0 nan nan nan nan
-15 11 0.0 nan nan nan nan
-15 12 0.0 nan nan nan nan
-15 13 0.0 nan nan nan nan
-15 14 0.0 nan nan nan nan
-15 15 1483586.0 0.21420590648786764 0.05654804359524735 0.7761049763490119 0.8067742008976206

mustang.gf

@@ -1,164 +0,0 @@
-#! /usr/bin/env ./geometryfactor.py
-# Mustang 
-φ1=0
-φ2=π/2
-Δr=10
-Δθ=0.5°
-
-# length units are mm, distances are center-center
-
-# dimensions of the Scintillators
-define=wx=500.
-define=wy=${wx}
-define=wz=50.
-define=αy=6°
-
-# position of the szintillators inside the boxs
-define=ox=523.
-define=oy=510.
-define=oz=0
-
-# position of the boxes
-define=px=1150.
-define=py=${px}
-define=pz=1155.
-
-#                        +y
-#          __    __               __    __
-#     __---  |  |  ---__  ^  __---  |  |  ---__
-#    |    21 |  |    23 | | |    29 |  |    31 |
-#    |       |  |       | | |       |  |       |
-#    |  5    |  |  7    | | | 13    |  | 15    |
-#    |     __|  |__     | | |     __|  |__     |
-#    |__---        ---__| | |__---        ---__|
-#          __    __       |       __    __
-#     __---  |  |  ---__  |  __---  |  |  ---__
-#    |    20 |  |    22 | | |    28 |  |    30 |
-#    |       |  |       | | |       |  |       |
-#    |  4    |  |  6    | | | 12    |  | 14    |
-#    |     __|  |__     | | |     __|  |__     |
-#    |__---        ---__| | |__---        ---__|
-#  -----------------------+----------------------> +x
-#          __    __       |       __    __
-#     __---  |  |  ---__  |  __---  |  |  ---__
-#    |    17 |  |    19 | | |    25 |  |    27 |
-#    |       |  |       | | |       |  |       |
-#    |  1    |  |  3    | | |  9    |  | 11    |
-#    |     __|  |__     | | |     __|  |__     |
-#    |__---        ---__| | |__---        ---__|
-#          __    __       |       __    __
-#     __---  |  |  ---__  |  __---  |  |  ---__
-#    |    16 |  |    18 | | |    24 |  |    26 |
-#    |   (+z)|  |       | | |       |  |       |
-#    |  0    |  |  2    | | |  8    |  | 10    |
-#    | (-z)  |  |__     | | |     __|  |__     |
-#    |__---        ---__| | |__---        ---__|
-#                         |
-
-
-box=${wx}/2 ${wy}/2 ${wz}/2
-ry=${αy}
-move=-${ox}/2-${px}/2,-${oy}/2-${py}/2,-${oz}/2-${pz}/2
-copy=-1
-move=0,${oy},0
-
-box=${wx}/2 ${wy}/2 ${wz}/2
-ry=-${αy}
-move=${ox}/2-${px}/2,-${oy}/2-${py}/2,-${oz}/2-${pz}/2
-copy=-1
-move=0,${oy},0
-
-copy=0
-move=0,${py},0
-copy=1
-move=0,${py},0
-copy=2
-move=0,${py},0
-copy=3
-move=0,${py},0
-
-copy=0
-move=${px},0,0
-copy=1
-move=${px},0,0
-copy=2
-move=${px},0,0
-copy=3
-move=${px},0,0
-copy=4
-move=${px},0,0
-copy=5
-move=${px},0,0
-copy=6
-move=${px},0,0
-copy=7
-move=${px},0,0
-
-copy=0
-move=0,0,${pz}
-copy=1
-move=0,0,${pz}
-copy=2
-move=0,0,${pz}
-copy=3
-move=0,0,${pz}
-copy=4
-move=0,0,${pz}
-copy=5
-move=0,0,${pz}
-copy=6
-move=0,0,${pz}
-copy=7
-move=0,0,${pz}
-copy=8
-move=0,0,${pz}
-copy=9
-move=0,0,${pz}
-copy=10
-move=0,0,${pz}
-copy=11
-move=0,0,${pz}
-copy=12
-move=0,0,${pz}
-copy=13
-move=0,0,${pz}
-copy=14
-move=0,0,${pz}
-copy=15
-move=0,0,${pz}
-
-ray=z00; -${ox}/2-${px}/2, -${oy}/2-${py}/2;    0,0,1
-ray=z01; -${ox}/2-${px}/2, +${oy}/2-${py}/2;    0,0,1
-ray=z02; +${ox}/2-${px}/2, -${oy}/2-${py}/2;    0,0,1
-ray=z03; +${ox}/2-${px}/2, +${oy}/2-${py}/2;    0,0,1
-ray=z10; -${ox}/2-${px}/2, -${oy}/2+${py}/2;    0,0,1
-ray=z11; -${ox}/2-${px}/2, +${oy}/2+${py}/2;    0,0,1
-ray=z12; +${ox}/2-${px}/2, -${oy}/2+${py}/2;    0,0,1
-ray=z13; +${ox}/2-${px}/2, +${oy}/2+${py}/2;    0,0,1
-ray=z20; -${ox}/2+${px}/2, -${oy}/2-${py}/2;    0,0,1
-ray=z21; -${ox}/2+${px}/2, +${oy}/2-${py}/2;    0,0,1
-ray=z22; +${ox}/2+${px}/2, -${oy}/2-${py}/2;    0,0,1
-ray=z23; +${ox}/2+${px}/2, +${oy}/2-${py}/2;    0,0,1
-ray=z30; -${ox}/2+${px}/2, -${oy}/2+${py}/2;    0,0,1
-ray=z31; -${ox}/2+${px}/2, +${oy}/2+${py}/2;    0,0,1
-ray=z32; +${ox}/2+${px}/2, -${oy}/2+${py}/2;    0,0,1
-ray=z33; +${ox}/2+${px}/2, +${oy}/2+${py}/2;    0,0,1
-
-ray=x00; 0, -${oy}/2-${py}/2, -${oz}/2-${pz}/2; 1,0,0
-ray=x01; 0, +${oy}/2-${py}/2, -${oz}/2-${pz}/2; 1,0,0
-ray=x02; 0, -${oy}/2+${py}/2, -${oz}/2-${pz}/2; 1,0,0
-ray=x03; 0, +${oy}/2+${py}/2, -${oz}/2-${pz}/2; 1,0,0
-ray=x10; 0, -${oy}/2-${py}/2, -${oz}/2+${pz}/2; 1,0,0
-ray=x11; 0, +${oy}/2-${py}/2, -${oz}/2+${pz}/2; 1,0,0
-ray=x12; 0, -${oy}/2+${py}/2, -${oz}/2+${pz}/2; 1,0,0
-ray=x13; 0, +${oy}/2+${py}/2, -${oz}/2+${pz}/2; 1,0,0
-
-ray=y00; -${ox}/2-${px}/2, 0, -${oz}/2-${pz}/2; 0,1,0
-ray=y01; +${ox}/2-${px}/2, 0, -${oz}/2-${pz}/2; 0,1,0
-ray=y02; -${ox}/2+${px}/2, 0, -${oz}/2-${pz}/2; 0,1,0
-ray=y03; +${ox}/2+${px}/2, 0, -${oz}/2-${pz}/2; 0,1,0
-ray=y10; -${ox}/2-${px}/2, 0, -${oz}/2+${pz}/2; 0,1,0
-ray=y11; +${ox}/2-${px}/2, 0, -${oz}/2+${pz}/2; 0,1,0
-ray=y12; -${ox}/2+${px}/2, 0, -${oz}/2+${pz}/2; 0,1,0
-ray=y13; +${ox}/2+${px}/2, 0, -${oz}/2+${pz}/2; 0,1,0
-

plot-z.awk

@@ -1,32 +0,0 @@
-#! /usr/bin/gawk -f
-
-BEGIN {
-    Z=0
-}
-
-/^[0-9]/&& NF>=4 {
-    theta = -$1+0
-    phi = $2+0
-    r = $3+0
-    psi = $4+0
-    x0 = r*cos(psi)
-    y0 = r*sin(psi)
-    # rotate by -theta around y
-    costh = cos(theta)
-    sinth = sin(theta)
-    xx1 =  x0*costh
-    xx2 =  x0*costh
-    z1 =  -x0*sinth
-    z2 =  -x0*sinth
-    # rotate by phi around z
-    cosphi = cos(phi)
-    sinphi = sin(phi)
-    x1 =  xx1*cosphi - y0*sinphi
-    x2 =  xx2*cosphi - y0*sinphi
-    y1 =  xx1*sinphi + y0*cosphi
-    y2 =  xx2*sinphi + y0*cosphi
-    print x1, y1, z1
-    print x2, y2, z2
-    print Nix
-    print Nix
-}

plot_isotrop.py

@@ -0,0 +1,56 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def Rays(size):
+        dir = np.empty((size, 3, 1))
+        ga = (3 - np.sqrt(5)) * np.pi
+        th = ga * np.arange(size)
+        z = np.linspace(1 / size - 1, 1 - 1 / size, size)
+        radius = np.sqrt(1 - (z**2))
+        y = radius * np.sin(th)
+        x = radius * np.cos(th)
+        dir[:,0,0] = x
+        dir[:,1,0] = y
+        dir[:,2,0] = z
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')
+        ax.scatter(x, y, z)
+        plt.title(f'Distribution of points in 3D')
+        ax.set_xlabel('x')
+        ax.set_ylabel('y')
+        ax.set_zlabel('z')
+        plt.savefig(f'create_rays.png')
+        plt.close()
+
+def scatter_iso(size):
+        u = np.random.uniform(0, 1, size)
+        v = np.random.uniform(0, 1, size)
+        phi = 2 * np.pi * u
+        sin_theta = np.sqrt(1 - v**2)
+        theta = np.arcsin(sin_theta)
+        cos_theta = np.cos(theta)
+        plt.figure(figsize=(10, 6))
+        plt.scatter(cos_theta, phi, s=1, color='blue')
+        plt.xlabel('cos(theta)')
+        plt.ylabel('phi')
+        plt.title('Plot of phi against cos(theta)')
+        plt.grid(True)
+        plt.show()
+        # x = sin_theta * np.cos(phi)
+        # y = sin_theta * np.sin(phi) 
+        # z = np.sqrt(1 - (sin_theta)**2)  
+        # fig = plt.figure()
+        # ax = fig.add_subplot(111, projection='3d')
+        # ax.scatter(x, y, z, s=1)
+        # ax.set_box_aspect([1,1,1])
+        # max_range = np.array([x.max()-x.min(), y.max()-y.min(), z.max()-z.min()]).max() / 2.0
+        # mid_x = (x.max() + x.min()) * 0.5
+        # mid_y = (y.max() + y.min()) * 0.5
+        # mid_z = (z.max() + z.min()) * 0.5
+        # ax.set_xlim(mid_x - max_range, mid_x + max_range)
+        # ax.set_ylim(mid_y - max_range, mid_y + max_range)
+        # ax.set_zlim(mid_z - max_range, mid_z + max_range)
+        # plt.show()
+
+print(scatter_iso(100000))

plot_rays.awk

@@ -1,33 +0,0 @@
-#! /usr/bin/gawk -f
-
-BEGIN {
-    R=100
-}
-
-/^[0-9]/&& NF>=4 {
-    theta = -$1+0
-    phi = $2+0
-    r = $3+0
-    psi = $4+0
-    x0 = r*cos(psi)
-    y0 = r*sin(psi)
-    z0 = sqrt(R**2-r**2)
-    # rotate by -theta around y
-    costh = cos(theta)
-    sinth = sin(theta)
-    xx1 =  x0*costh + z0*sinth
-    xx2 =  x0*costh - z0*sinth
-    z1 =  -x0*sinth + z0*costh
-    z2 =  -x0*sinth - z0*costh
-    # rotate by phi around z
-    cosphi = cos(phi)
-    sinphi = sin(phi)
-    x1 =  xx1*cosphi - y0*sinphi
-    x2 =  xx2*cosphi - y0*sinphi
-    y1 =  xx1*sinphi + y0*cosphi
-    y2 =  xx2*sinphi + y0*cosphi
-    print x1, y1, z1
-    print x2, y2, z2
-    print Nix
-    print Nix
-}

test.py

@@ -0,0 +1,34 @@
+import random
+from Random_gen import test_planes,test_lines
+from Struktur import Plane
+
+def test_cross(test_size):
+    rand_planes=test_planes(test_size)
+    rand_lines=test_lines(test_size)
+    cross_point=[]
+    for i in range(len(rand_planes)):
+        cross_point_fun=Plane.cross(rand_planes[i],rand_lines[i])
+        if cross_point!=str:
+            cross_point.append(cross_point_fun)
+        if Plane.point_in_plane(rand_planes[i],cross_point)==False:
+            raise Exception('Catastrophe!')
+    return cross_point
+
+def test_cross_angle(test_size):
+    rand_planes=test_planes(test_size)
+    rand_lines=test_lines(test_size)
+    cross_angle=[]
+    for i in range(len(rand_planes)):
+        angle=Plane.cross_angle(rand_planes[i],rand_lines[i])
+        cross_angle.append(angle)
+    return cross_angle
+
+ 
+print(test_cross_angle(60))
+
+
+
+        
+
+
+

theta_path.awk

@@ -1,30 +0,0 @@
-#! /usr/bin/gawk -f
-
-BEGIN {
-    dr = 1
-    r_max = 100
-    theta = 9999
-}
-
-$1 != theta {
-    if (n) emit()
-    theta = $1
-}
-
-{
-    i = int($5/dr+0.5)
-    N[i]++
-    n++
-    if (i) m++
-}
-
-function emit() {
-    nr = r_max/dr+1
-    print "#", n, m
-    for (i=0; i<nr; i++) print theta, i*dr, N[i]+0
-    print ""
-    delete N
-    n = 0
-}
-
-END { if (n) emit() }

try.py

@@ -0,0 +1,35 @@
+import numpy as np
+
+A=np.array([[[2,2,2],[4,4,4]],[[4,5,6],[7,8,9]]])
+b=np.array([6,6,6])
+c=np.array(6)
+#print(A.shape)
+#print(b.shape)
+#print(c.shape)
+#print(A*b)
+b=b[np.newaxis,np.newaxis,:]
+#print(b.shape)
+#print(b)
+b=np.tile(b,(2,2,1))
+#print(b.shape)
+#print(A*b)
+
+P=np.zeros((3,3,3))
+#print(P)
+n=np.array([[4],[5],[6]])
+n=np.squeeze(n.reshape(1,-1))
+print(n)
+P[:,1]=n.T
+print(P)
+P[...,2:3]=np.array([[7],[8],[9]])
+#print(P[...,0:1])
+n=np.array([[[2],[1],[3]],[[6],[4],[3]],[[8],[4],[0]]])
+#print(n)
+P[...,0:1]=n
+v=np.array([4,9,2])
+v=v[np.newaxis,:]
+v=np.tile(v,(3,1))
+#print(P)
+#print(v.shape)
+#print(P.shape)
+#print(np.linalg.solve(P,v-np.array([2,4,3])))