add --equivalent

find equivalent planes in all detektor volumes and replace duplicates

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

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