add --equivalent

find equivalent planes in all detektor volumes and replace duplicates

BGO.py

@@ -1,114 +0,0 @@
-import json
-from Struktur import Plane, Ray, Rays
-import numpy as np
-
-#load the data from the json files
-with open('BGO_planes.json') as json_file:
-    data1 = json.load(json_file)
-planes_json=data1['Planes']
-with open('Detector_planes.json') as json_file:
-    data6 = json.load(json_file)
-planes_det_json=data6['Planes Det']
-
-
-def extract_planes(planes_json):
-    plane_list=[]
-    for key in planes_json:
-        plane_list.append(Plane(planes_json[key][0],planes_json[key][1],planes_json[key][2]))
-    return plane_list
-
-#extracts the example_points from the json file 
-#and returns them
-def extract_points(points_json):
-    points=[]
-    for key in points_json:
-        points.append(points_json[key])
-    return points
-
-
-def ray_way(ray,scinti_planes,scinti_shapes,point_list,ind):
-    hit_inc=False
-    num=0
-    new_ray=0
-    for i in range(len(scinti_shapes)):
-        cross_point=Plane.cross(scinti_planes[i],ray)
-        if cross_point!= None and i!=ind:
-            if Shape.hit(scinti_shapes[i],point_list[i],cross_point,len(scinti_shapes[i].g_list))==True:
-                p=np.dot(scinti_planes[i].norm,ray.dir)
-                if p<0:
-                    hit_inc=True
-                    num=i
-                    hit_point=cross_point
-    if hit_inc!=True:
-        ray.dir=-ray.dir
-        for i in range(len(scinti_shapes)):
-            cross_point=Plane.cross(scinti_planes[i],ray)
-            if cross_point!= None and i!=ind:
-                if Shape.hit(scinti_shapes[i],point_list[i],cross_point,len(scinti_shapes[i].g_list))==True:
-                    p=np.dot(scinti_planes[i].norm,ray.dir)
-                    if p<0:
-                        hit_inc=True
-                        num=i
-                        hit_point=cross_point
-    angle=Plane.cross_angle(scinti_planes[num],ray)
-    if angle<27.7177475:
-        new_dir=Plane.scatter(scinti_planes[num],ray)
-        new_ray=Ray(hit_point,hit_point+new_dir)
-        print('scatter')
-    else:
-        new_dir=Plane.reflect(scinti_planes[num],ray)
-        new_ray=Ray(hit_point,hit_point+new_dir)
-        print('reflect')
-    return [new_ray,num,hit_point]
-
-def destroy(det_shapes,det_planes,det_points,det_borders,ray):
-    res=False
-    for i in range(len(det_planes)):
-        cross_point=Plane.cross(det_planes[i],ray)
-        if cross_point!=None:
-            if Shape.hit(det_shapes[i],det_points[i],cross_point,len(det_borders[i]))==True:
-                res=i+1
-    if res!=False:
-        return res
-    else:
-        return ray
-    
-
-def ray_journey(ray,scinti_planes,scinti_shapes,scinti_points,det_shapes,det_planes,det_points,det_borders):
-    ray_list=[]
-    hit_ind=9
-    res=False
-    while res==False:
-        stop=destroy(det_shapes,det_planes,det_points,det_borders,ray)
-        if stop==1 or stop==2:
-            res=True
-            print('Hit Detector'+str(stop))
-        else:
-            if len(ray_list)<100:
-                obj=ray_way(ray,scinti_planes,scinti_shapes,scinti_points,hit_ind)
-                hit_ind=obj[1]
-                print(obj[0])
-                print(obj[1]+1)
-                ray_list.append(obj[0])
-                ray=obj[0]
-            else:
-                res=True
-    num=len(ray_list)
-    stretch=0
-    for i in range(num-1):
-        l=ray_list[i+1].origin-ray_list[i].origin
-        stretch=np.linalg.norm(l)*0.001
-    print(stretch)
-    return ray_list
-    
-def create_dat(filename, list):
-    with open(filename, 'w') as file:
-        for i in range(len(list)):
-            if isinstance(list[i],Ray):
-                a,b,c=list[i].origin
-                d,e,f=list[i].dir
-                file.write(f"{a}  {b}  {c}  {d}  {e}  {f}\n")
-    
-
-
-print(extract_planes(planes_json))

BGO_planes.json

@@ -1,9 +0,0 @@
-{"Planes":{ "Plane 1":[[0,0,0],[0,-25.981,45],[20,0,0]],
-            "Plane 2":[[0,0,0],[20,0,0],[0,51.962,0]],
-            "Plane 3":[[0,51.962,0],[20,51.962,0],[0,77.943,45]],
-            "Plane 4":[[0,51.962,90],[0,77.943,45],[20,51.962,90]],
-            "Plane 5":[[0,0,90],[0,51.962,90],[20,0,90]],
-            "Plane 6":[[0,0,90],[20,0,90],[0,-25.981,45]],
-            "Plane 7":[[0,0,0],[0,51.962,0],[0,0,90]],
-            "Plane 8":[[20,0,0],[20,0,90],[20,51.962,0]]}}
-

BGO_point_in_plane.json

@@ -1,9 +0,0 @@
-{"Example_point":{  "Plane 1":[10,-10,17.321],
-                    "Plane 2":[10,20,0],
-                    "Plane 3":[10,61.962,17.321],
-                    "Plane 4":[10,61.962,62.321],
-                    "Plane 5":[10,20,90],
-                    "Plane 6":[10,-10,62.321],
-                    "Plane 7":[0,20,45],
-                    "Plane 8":[20,20,45]}}
-

Detector_planes.json

@@ -1,2 +0,0 @@
-{"Planes Det":{ "Det 1":[[0,0,0],[20,0,0],[0,51.962,0]],
-                "Det 2":[[0,0,90],[20,0,90],[0,51.962,90]]}}

Detector_point_in_plane.json

@@ -1,3 +0,0 @@
-{"Example_point Det":{  "Det 1":[10,20,0],
-                        "Det 2":[10,20,90]}}
-                   

Random_gen.py

@@ -1,41 +0,0 @@
-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

@@ -1,269 +0,0 @@
-from sys import argv, stdout, stderr
-import numpy as np
-import random
-from numpy.linalg import solve
-
-#size of the simulation
-size=3
-
-def vtuple(v):
-    return tuple(v[:,0])
-
-def matrix(vec1,vec2,vec3):
-    return np.array([vec1,vec2,vec3])
-
-def norm_vec(vec):
-    return vec/np.linalg.norm(vec)
-
-def vector(x=0., y=0., z=0.):
-    return np.array([[x],[y],[z]])
-
-def length(v):
-    x=0
-    for i in range(len(v)):
-        x=x+v[i]**2
-    res=np.sqrt(x)
-    return res
-
-def scale(x=1., y=1., z=1.):
-    return np.array(((x,0,0),
-                  (0,y,0),
-                  (0,0,z)))
-
-def rotate(axis, φ):
-    m = scale()
-    c = np.cos(φ)
-    s = -np.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
-
-#performs a rotation around any axis in the 3D space defined by a vector
-#using the Rodriguez formula
-def rotate_any(axis, vec, angle):
-    K=matrix([0,-axis[2],axis[1]],
-            [axis[2],0,axis[0]],
-            [-axis[1],axis[0],0])
-    A=np.identity(3) + np.sin(angle) * K + (1-np.cos(angle)) * np.dot(K,K)
-    return np.dot(A,vec)
-    
-class Ray:
-    #self.o is the starting point. self.d the direction vector of the line
-    #o and p are two points on the line
-    def __init__(self, o, p):
-        self.origin=o
-        self.dir=np.subtract(o,p) #calculates o-p
-
-    def __str__(self):
-        return str(self.origin)+str(self.dir)
-
-class rays:
-    "Fast reimplementation of a beam of rays"
-    
-    def __init__(self, r, Δr, origin):
-        self.origin=origin
-        nr = int(r/Δr + 1.001)
-        self.n = 0
-        rpsi = 0
-        for ir in range(nr):
-            nnpsi = 2*np.pi*ir + 0.5 + rpsi
-            npsi = int(nnpsi + 0.5)
-            rpsi = nnpsi - npsi
-            self.n += npsi
-        self.xy = np.empty((self.n,3,1))
-        self.r = np.empty((self.n,))
-        self.psi = np.empty((self.n,))
-        i=0
-        rpsi = 0
-        for ir in range(nr):
-            r = ir*Δr
-            nnpsi = 2*np.pi*ir + 0.5 + rpsi
-            npsi = int(nnpsi + 0.5)
-            rpsi = nnpsi - npsi
-            dpsi = 2*np.pi/npsi
-            for ipsi in range(npsi):
-                psi = ipsi * dpsi + dpsi/2
-                self.r[i] = r
-                self.psi[i] = psi
-                self.xy[i,:,0] = (r*np.cos(psi), r*np.sin(psi), 0)
-                i = i+1
-
-    def bend(self, th, phi):
-        self.th = th
-        self.phi = phi
-        self.v = rotate(2,phi) @ rotate(1,-th)
-        self.dir = self.v @ vector(z=1)
-        self.origin = (self.v @ self.xy)[...,0].T
-        #print(self.origin)
-    
-class Plane:
-
-    #o, v, w are points on the plane
-    def __init__(self, o, v, w):
-        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.a=np.dot(self.norm,o)
-        self.tolerance=1e-9
-        self.A=np.empty((135,3,3))
-        self.A[...,1:2]=o-v
-        self.A[...,2:3]=w-v
-    
-    #calculates the matrix for the rotate function
-    def rotate_matrix(self):
-        return matrix([0,-self.norm[2],self.norm[1]],
-                        [self.norm[2],0,-self.norm[0]],
-                        [-self.norm[1],self.norm[0],0])
-    
-    #claculates the cross point of a plane and a line and returns it
-    #also checks wether the line and plane are parallel
-    def cross(self,g):
-        res=True
-        self.A[...,0:1]=-g.xy
-        print(self.A)
-        try:
-            t = np.linalg.solve(self.A, g.origin-self.v)
-        except np.linalg.LinAlgError:  
-            return False,0
-        return g.origin+t*g.xy
-    
-    #calculates the angle between the normal vector of the plane 
-    #and the direction vector of the line
-    def cross_angle(self,g):
-        n_len=length(self.norm)
-        dir_len=length(g.dir)
-        sc=np.dot(g.dir, self.norm)
-        arg=sc/(n_len*dir_len)
-        res=np.rad2deg(np.arccos(arg))
-        if sc<0:
-            res=360-res-180
-        return res
-    
-    #calculates the cross line of two planes
-    #returns the origin and direction vector of the line
-    def cross_planes(self,pl2):
-        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
-    
-    #calculates the reflection vector of a line g and a plane 
-    def reflect(self,g):
-        rot_vec=-rotate_any(self.norm,g.dir,np.pi)
-        if Plane.cross_angle(self,g)==0:
-            rot_vec=-g.dir
-        if np.dot(rot_vec,self.norm)<0:
-            rot_vec=-rot_vec
-        return rot_vec
-        
-    
-    #calcultes the refracted vector of a line g, a plane and the refraction index 
-    #of the scintillator n1 and the surrounding medium n2
-    def refract(self,g,n1,n2):
-        norm_dir=norm_vec(g.dir)
-        refrac_ind=n2/n1
-        angle=Plane.cross_angle(self,g)
-        l=length(g.dir)
-        vec_parallel=np.dot(norm_vec(self.norm),norm_dir)
-        refrac_angle=np.arcsin(refrac_ind * np.sin(angle))
-        vec_orth=[0,0,0]
-        for i in range(3):
-            vec_orth[i]=np.sin(refrac_angle)*norm_dir[i]
-        new_vec_norm=vec_parallel+vec_orth
-        new_vec=-l*new_vec_norm
-        return new_vec
-        
-    def scatter(self,g):
-        u = np.random.uniform()
-        v = np.random.uniform()
-        phi = 2 * np.pi * u
-        theta = np.arccos(2 * v - 1)
-        new_vec = np.array([(np.sin(theta) * np.cos(phi))*40,(np.sin(theta) * np.sin(phi))*40,(np.cos(theta))*40])
-        if np.dot(new_vec, self.norm) < 0:
-            new_vec = -new_vec
-        return new_vec
- 
-    def line_in_plane(self,g):
-        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):
-        distance=np.dot(self.norm,self.o - point)
-        if abs(distance)<self.tolerance:
-            return True
-        else:
-            return False
-    #Attention: the tolerance might need to be higher for certain examples
-    #source of error!
-
-    def __str__(self):
-        return str(self.norm)+str(self.a)+str(self.o).str(self.dir_vec1)
-    
-class Scintillator(list):
-
-    def real_hit(self,g,i):
-        cp=self[i].cross(g)
-        hit=True
-        j=0
-        while j<len(self) and hit==True:
-            if np.dot(self[j].norm,cp-self[j].o)<0:
-                hit=False
-        return hit
-    
-
-#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)),
-#               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)),
-#               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))])
-
-pl=Plane(vector(0,0,0),vector(20,0,0),vector(0,51.962,0))
-r=rays(3,0.5,vector(10,10,10))
-print(Plane.cross(pl,r))
-
-
-def create_dat(filename, list):
-    with open(filename, 'w') as file:
-        for i in range(len(list)):
-            d=list[i][0][0]
-            e=list[i][1][0]
-            f=list[i][2][0]
-            a,b,c=0,0,0
-            file.write(f"{a}  {b}  {c} {d} {e} {f}\n")
-    
-
-    
-
-
-
-  
-
-        
-
-
-
-
-
-
-
- 

Visualize.gp

@@ -1,29 +0,0 @@
-#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

@@ -1,64 +0,0 @@
-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')
-
-
-
-
-

chaos-jr.gf

@@ -0,0 +1,26 @@
+#! /usr/bin/env ./geometryfactor.py
+# CHAOS Jr.  HET without A-detectors
+
+φ1=0
+φ2=π/6
+Δr=0.1
+Δθ=0.5°
+
+#BGO
+define=W=20
+define=D=40
+prism= 6, ${D}/2*sqrt(3/4), ${W}/2
+
+# SSDs
+define=a=24.3
+define=w=0.3
+define=d=36.5
+define=n=36
+
+prism= ${n}, ${d}/2, ${w}/2
+move= 0,0,-${a}/2
+copy=-1
+move= 0,0,${a}
+
+θ₂=atan2(${d},${a}+${w})
+radius=sqrt(${d}**2 + (${a}+${w})**2)/2

chaos.gf

@@ -0,0 +1,56 @@
+#! /usr/bin/env ./geometryfactor.py
+# CHAOS Jr.  HET without A-detectors
+
+φ1=0
+φ2=π/6
+Δr=0.1
+Δθ=0.5°
+
+#BGO
+define=W=20
+define=D=2*51.962
+prism= 6, ${D}/2*sqrt(3/4), ${W}/2
+
+#Aerogel
+#box = 62,62,40
+#move = 0,0,${W}/2+10
+prism= 4,31,20
+move= 0,0,40
+
+# SSDs
+define=w=0.3
+define=d=36.5
+define=n=36
+define=d1=17.4
+
+prism= ${n}, ${d}/2, ${w}/2
+#A
+move= 0,0,${W}/2+55
+#C
+copy=-1
+move= 0,0,-50
+#E
+copy=-1
+move= 0,0,-30
+
+#E123
+copy=-1
+move=-18.5,-10.53,-4
+copy=-1
+move=2*18.5,0,-4
+copy=-1
+move=-18.25,31.6,-4
+
+#A1
+prism= ${n}, ${d1}/2, ${w}/2
+move= 0,0,${W}/2+55
+#C1
+copy=-1
+move= 0,0,-50
+
+#?
+define=a=24.7
+θ₂=atan2(${d},${a}+${w})
+radius=sqrt(${d}**2 + (${a}+${w})**2)/2
+
+

geometryfactor.py

@@ -1,9 +1,10 @@
 #! /usr/bin/python3
 
 from sys import argv, stdout, stderr
-from math import pi as π
+from math import pi as π, floor
+from math import atan, atan2
 import numpy
-from numpy import sqrt, sin, cos, tan, array, arange, empty, cross, newaxis
+from numpy import sqrt, sin, cos, arccos, tan, array, arange, empty, cross, newaxis
 from numpy.linalg import solve
 from getopt import getopt
 
@@ -60,7 +61,13 @@ class ray:
 class rays:
     "Fast reimplementation of a beam of rays"
     
-    def __init__(self, r, Δr):
+    def __init__(self, r, Δr, spiral=False):
+        if spiral:
+            self.spiral(r, Δr)
+        else:
+            self.rings(r, Δr)
+
+    def rings(self, r, Δr):
         nr = int(r/Δr + 1.001)
         self.n = 0
         rψ = 0
@@ -87,6 +94,16 @@ class rays:
                 self.xy[i,:,0] = (r*cos(ψ), r*sin(ψ), 0)
                 i = i+1
 
+    def spiral(self, r, Δr):
+        self.n = int(π * (r/Δr)**2 + 0.5)
+        self.ψ = sqrt(4*π * arange(self.n))
+        self.r = Δr/(2*π) * self.ψ
+        self.ψ %= 2*π
+        self.xy = empty((self.n, 3, 1))
+        self.xy[:,0,0] = self.r * cos(self.ψ)
+        self.xy[:,1,0] = self.r * sin(self.ψ)
+        self.xy[:,2,0] = 0
+
     def bend(self, θ, φ):
         self.θ = θ
         self.φ = φ
@@ -104,6 +121,7 @@ class plane:
         self.v = v
         self.w = w
         self.normalize()
+        self.cached_r = None
 
     def normalize(self):
         self.n = cross(self.v, self.w, axis=0).T
@@ -126,15 +144,29 @@ class plane:
     def inside(self, p):
         return (self.n @ (p-self.o))[0] >= -0.00001
 
+    def equivalent(self, other):
+        if length(cross(self.n, other.n)) > 1e-10:
+            return False
+        if (self.n.T @ other.n)[0,0] <= 0:
+            return False
+        if length(self.n @ (other.o - self.o)) > 1e-10:
+            return False
+        return True
+
     def intersection(self, r):
+        if self.cached_r is r and self.cached_v is r.v:
+            return self.solution
+        self.cached_r = r
+        self.cached_v = r.v
         self.M[:,2:3] = r.v
         try:
             c=solve(self.M, r.o-self.o)[...,2,:]
+            self.solution = (c, r.o - c*r.v)
         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
+            self.solution = (-1, None)
+        return self.solution
 
     def plane_intersection(self, other, limit=1e10, delta=1e-5):
         p = [ self.intersection(r)[1] for r in other.rays]
@@ -247,7 +279,7 @@ class detector(list):
         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])))
+            stdout.write("%g %g %g\n%g %g %g\n\n\n" % (tuple(e[0][:,0])+tuple(e[1][:,0])))
 
     def planes(self):
         # is this useful?
@@ -269,12 +301,23 @@ def box(x=1.0, y=1.0, z=1.0):
         plane(vector(y= y), vector(x=1), vector(z=1)),
     ])
 
-def prism(n=6, a=1.0, z=1.0, φ=0.0):
+def prism(n=6, a=1.0, z=1.0, φ=0.0, n1=0, n2=None):
     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 + φ
+    if n2 is None:
+        n2 = n-1
+    else:
+        c = cos(2*π/n * n1)
+        s = sin(2*π/n * n1)
+        p.append(plane(vector(), vector(c,s), vector(z=-1)))
+        c = cos(2*π/n * n2)
+        s = sin(2*π/n * n2)
+        p.append(plane(vector(), vector(c,s), vector(z= 1)))
+ 
+    φ += 2*π/n * arange(int(floor(n1+0.4999)),
+                        int(floor(n2+1.5001)) )
     p.extend([
         plane(
             vector(a*c, a*s),
@@ -287,6 +330,25 @@ def prism(n=6, a=1.0, z=1.0, φ=0.0):
 
 µM = box(41., 41., 13.)
 
+def find_equivalent_planes(world):
+    done = []
+    for i,d in enumerate(world):
+        for j,p in enumerate(d):
+            if (i,j) in done:
+                continue
+            for ii,dd in enumerate(world):
+                if ii < i:
+                    continue
+                for jj,pp in enumerate(dd):
+                    if ii == i and jj <= j:
+                        continue
+                    if p.equivalent(pp):
+                        dd[jj] = p
+                        done.append((ii,jj))
+                        if do_size or verbose:
+                            print(f"equiv panes({len(done)}) {i}/{j} → {ii}/{jj}",
+                                  file=stderr)
+
 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)
@@ -312,7 +374,7 @@ def run(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
                     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*π):
+def run_rings(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
@@ -332,6 +394,24 @@ def run_beam(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
                 ll = tuple(bb+[l[i] for l in pl])
                 stdout.write(fmt % ll)
         
+def run_beam(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
+    if φ2 < 2*π - Δθ or φ1 > Δθ:
+        run_rings(world, r, Δr, Δθ, θ1, θ2, φ1, φ2)
+    fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
+    beam = rays(r, Δr, spiral=True)
+    N = 2*π / Δθ**2
+    t1 = N*(1-cos(θ1))
+    n = int(N*(cos(θ1) - cos(θ2)) + 0.5)
+    t = arange(n) + t1
+    for θ in arccos(1 - t/N):
+        φ = (2*π * θ/Δθ) % (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))),
@@ -360,7 +440,7 @@ 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=" ,
+            "radius=", "resolution=", "Δr=", "equivalent",
             "detector", "copy=", "select=", "box=", "prism=",
             "plane=", "ray=",
             "rx=", "ry=", "rz=", "move=", "scale=",
@@ -394,12 +474,13 @@ world = [detector()]
 θ1=0
 θ2=π/2
 φ1=0
-φ2=2*π/8
+φ2=2*π
 do_size = False
 do_test = False
 do_run = False
 do_corners = False
 do_gnuplot = False
+do_equiv = False
 rot = rotate(0,0)
 ori = vector()
 
@@ -462,6 +543,8 @@ for o,v in options:
         do_gnuplot=True
     if o in "-E --corners":
         do_corners=True
+    if o in "--equivalent":
+        do_equiv=True
     if o in "--rx":
         world[-1] = world[-1].transform(rotate(0, aval(v)))
     if o in "--ry":
@@ -502,13 +585,17 @@ for o,v in options:
             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 verbose:
+                stderr.write("define macro %s=%s\n" % (k,vv[1]))
 
 if not world[-1]:
     world[-1:]=[]
 if not world:
     world=[µM]
 
+if do_equiv:
+    find_equivalent_planes(world)
+
 if not r or do_corners:
     max_r, c = corners(world)
     if not r:
@@ -519,11 +606,10 @@ if do_gnuplot:
         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
+    no = int(π * (r/Δr)**2 + 0.5)
+    nv = int((φ2 - φ1) / Δθ**2 * (cos(θ1) - cos(θ2)) + 0.5)
+    sr = nv * Δθ**2
+    A = π * r**2
     nn = int(nv)*int(no)
     stderr.write("""Run size estimate:
 detectors: %d  planes: %d

journey.dat

@@ -1,65 +0,0 @@
-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

norm.dat

@@ -1,9 +0,0 @@
-10 -10 17.321 0 900 519.62
-#10 20 0 0 0 1039.24
-#10 61.962 17.321 0 -900 519.62
-#10 61.962 62.321 0 -900 -519.62
-#10 20 90 0 0 -1039.24
-#10 -10 62.321 0 900 -519.62
-#0 20 45 4676.58 0 0
-#20 20 45 -4676.58 0 0
-

seth.gf

@@ -0,0 +1,110 @@
+#! /usr/bin/env ./geometryfactor.py
+# CHAOS Jr.  HET without A-detectors
+
+Δr=1
+Δθ=5°
+radius= 66
+
+#BGO
+define=W=20
+define=D=2*51.962
+prism= 6, ${D}/2*sqrt(3/4), ${W}/2
+
+move= 0, 0, -15
+copy= -1
+move= 0, 0,  30
+
+box= 10, 5, 0.15
+move= 0, -49.7, -30.25
+copy= -1
+move= 0, 14.2, 0
+copy= -1
+move= 0, 14.2, 0
+copy= -1
+move= 0, 14.2, 0
+copy= -1
+move= 0, 14.2, 0
+copy= -1
+move= 0, 14.2, 0
+copy= -1
+move= 0, 14.2, 0
+copy= -1
+move= 0, 14.2, 0
+
+copy= -7
+move= -27, 0,0
+copy= -7
+move= -27, 0,0
+copy= -7
+move= -27, 0,0
+copy= -7
+move= -27, 0,0
+copy= -7
+move= -27, 0,0
+copy= -7
+move= -27, 0,0
+
+copy= -6
+move= 54
+copy= -6
+move= 54
+copy= -6
+move= 54
+copy= -6
+move= 54
+copy= -6
+move= 54
+copy= -6
+move= 54
+
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+copy= -20
+move= 0,0, 60.5
+
+# SSDs
+define=w=0.5
+define=d=36.5
+define=n=36
+define=d1=17.4
+
+prism= ${n}, ${d}/2, ${w}/2
+move= 0,0, -${w}/2
+prism= ${n}, ${d1}/2, ${w}/2
+move= 0,0, ${w}/2

test.py

@@ -1,34 +0,0 @@
-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))
-
-
-
-        
-
-
-

try.py

@@ -1,35 +0,0 @@
-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])))