add my simulation setup for my experiment

chaos.gf

@@ -1,56 +0,0 @@
-#! /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

@@ -4,7 +4,7 @@ from sys import argv, stdout, stderr
 from math import pi as π, floor
 from math import atan, atan2
 import numpy
-from numpy import sqrt, sin, cos, arccos, tan, array, arange, empty, cross, newaxis
+from numpy import sqrt, sin, cos, tan, array, arange, empty, cross, newaxis
 from numpy.linalg import solve
 from getopt import getopt
 
@@ -61,13 +61,7 @@ class ray:
 class rays:
     "Fast reimplementation of a beam of rays"
     
-    def __init__(self, r, Δr, spiral=False):
-        if spiral:
-            self.spiral(r, Δr)
-        else:
-            self.rings(r, Δr)
-
-    def rings(self, r, Δr):
+    def __init__(self, r, Δr):
         nr = int(r/Δr + 1.001)
         self.n = 0
         rψ = 0
@@ -94,16 +88,6 @@ 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.φ = φ
@@ -121,7 +105,6 @@ 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
@@ -144,29 +127,15 @@ 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)))
-            self.solution = (-1, None)
-        return self.solution
+            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]
@@ -279,7 +248,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][:,0])+tuple(e[1][:,0])))
+            stdout.write("%g %g %g\n%g %g %g\n\n\n" % (tuple(e[0])+tuple(e[1])))
 
     def planes(self):
         # is this useful?
@@ -307,7 +276,7 @@ def prism(n=6, a=1.0, z=1.0, φ=0.0, n1=0, n2=None):
         plane(vector(z= z), vector(y=1), vector(x=1)),
     ]
     if n2 is None:
-        n2 = n-1
+        n2 = n
     else:
         c = cos(2*π/n * n1)
         s = sin(2*π/n * n1)
@@ -330,25 +299,6 @@ def prism(n=6, a=1.0, z=1.0, φ=0.0, n1=0, n2=None):
 
 µ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)
@@ -374,7 +324,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_rings(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
+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
@@ -394,24 +344,6 @@ def run_rings(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))),
@@ -440,7 +372,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=", "equivalent",
+            "radius=", "resolution=", "Δr=" ,
             "detector", "copy=", "select=", "box=", "prism=",
             "plane=", "ray=",
             "rx=", "ry=", "rz=", "move=", "scale=",
@@ -474,13 +406,12 @@ world = [detector()]
 θ1=0
 θ2=π/2
 φ1=0
-φ2=2*π
+φ2=2*π/8
 do_size = False
 do_test = False
 do_run = False
 do_corners = False
 do_gnuplot = False
-do_equiv = False
 rot = rotate(0,0)
 ori = vector()
 
@@ -543,8 +474,6 @@ 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":
@@ -585,17 +514,13 @@ for o,v in options:
             if k in vv[1]:
                 raise ValueError("recursive macro %s" % v)
             macros[k]=vv[1]
-            if verbose:
-                stderr.write("define macro %s=%s\n" % (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 do_equiv:
-    find_equivalent_planes(world)
-
 if not r or do_corners:
     max_r, c = corners(world)
     if not r:
@@ -606,10 +531,11 @@ if do_gnuplot:
         d.gnuplot()
 
 if do_size:
-    no = int(π * (r/Δr)**2 + 0.5)
-    nv = int((φ2 - φ1) / Δθ**2 * (cos(θ1) - cos(θ2)) + 0.5)
-    sr = nv * Δθ**2
-    A = π * r**2
+    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

seth.gf

@@ -1,110 +0,0 @@
-#! /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

toms_erste_welt.gf

@@ -0,0 +1,147 @@
+#! /usr/bin/env ./geometryfactor.py
+# Toms position dependency. My quick first model
+
+# all units in mm
+#BGO but Big 
+define=W=20
+# side length
+define=D=51.962
+#apothem
+define=apo=${D}*sqrt(3/4)
+
+# number of sides, apothem of hexagon, thickness but here divided by 2
+prism= 6, ${apo},${W}/2
+
+# define SSDs variables. photosensitive area is 1cmx1cm and thickness is 0.007cm
+define=wx=5
+define=wy=5
+define=wz=0.07
+
+# define 
+define=a=18.85
+define=b=36.65
+define=c=19.5
+define=d=37.5
+define=e=18.85
+define=f=27.5
+
+# sum thickness of teflon and millipor
+define=dc=0.5
+
+# box in the middle
+box= ${wx},${wx},${wz}
+move=0,0,(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# box on the sides
+box= ${wx},${wx},${wz}
+move=${b},0,(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# box on the sides
+box= ${wx},${wx},${wz}
+move=-${b},0,(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# box on the corners
+box= ${wx},${wx},${wz}
+move=0,${d},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# box on the corners
+box= ${wx},${wx},${wz}
+move=0,-${d},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# boxes between the edges
+box= ${wx},${wx},${wz}
+move= ${a},0, (${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# boxes between the edges
+define=dm=17
+box= ${wx},${wx},${wz}
+move= -${a},0,(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# boxes between the edges
+define=dm=17
+box= ${wx},${wx},${wz}
+move= 0,${c},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# boxes between the edges
+define=dm=17
+box= ${wx},${wx},${wz}
+move= 0,-${c},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# x,y,z: x:photodiodes, y: corner, z:height
+
+#boxes in on the sides of the edges
+box= ${wx},${wx},${wz}
+move=${b},${c},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+#boxes in on the sides of the edges
+box= ${wx},${wx},${wz}
+move=${b},-${c},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+#boxes in on the sides of the edges
+box= ${wx},${wx},${wz}
+move=-${b},${c},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+#boxes in on the sides of the edges
+box= ${wx},${wx},${wz}
+move=-${b},-${c},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# sidies
+define=dm=17
+box= ${wx},${wx},${wz}
+move= ${e},${f},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# sidie
+define=dm=17
+box= ${wx},${wx},${wz}
+move= -${e},${f},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# sidies
+define=dm=17
+box= ${wx},${wx},${wz}
+move= ${e},-${f},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+# sidie
+define=dm=17
+box= ${wx},${wx},${wz}
+move= -${e},-${f},(${W}+${dc})/2
+copy=-1
+move=0,0, -${W}-${dc}
+
+
+
+
+
+
+