Compare commits
No commits in common. "master" and "master" have entirely different histories.
5 changed files with 15 additions and 10511 deletions
6480
ahepam24.gpt
6480
ahepam24.gpt
File diff suppressed because it is too large
Load diff
56
chaos.gf
56
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
|
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -4,7 +4,7 @@ from sys import argv, stdout, stderr
|
||||||
from math import pi as π, floor
|
from math import pi as π, floor
|
||||||
from math import atan, atan2
|
from math import atan, atan2
|
||||||
import numpy
|
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 numpy.linalg import solve
|
||||||
from getopt import getopt
|
from getopt import getopt
|
||||||
|
|
||||||
|
|
@ -61,13 +61,7 @@ class ray:
|
||||||
class rays:
|
class rays:
|
||||||
"Fast reimplementation of a beam of rays"
|
"Fast reimplementation of a beam of rays"
|
||||||
|
|
||||||
def __init__(self, r, Δr, spiral=False):
|
def __init__(self, r, Δr):
|
||||||
if spiral:
|
|
||||||
self.spiral(r, Δr)
|
|
||||||
else:
|
|
||||||
self.rings(r, Δr)
|
|
||||||
|
|
||||||
def rings(self, r, Δr):
|
|
||||||
nr = int(r/Δr + 1.001)
|
nr = int(r/Δr + 1.001)
|
||||||
self.n = 0
|
self.n = 0
|
||||||
rψ = 0
|
rψ = 0
|
||||||
|
|
@ -94,16 +88,6 @@ class rays:
|
||||||
self.xy[i,:,0] = (r*cos(ψ), r*sin(ψ), 0)
|
self.xy[i,:,0] = (r*cos(ψ), r*sin(ψ), 0)
|
||||||
i = i+1
|
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, θ, φ):
|
def bend(self, θ, φ):
|
||||||
self.θ = θ
|
self.θ = θ
|
||||||
self.φ = φ
|
self.φ = φ
|
||||||
|
|
@ -121,7 +105,6 @@ class plane:
|
||||||
self.v = v
|
self.v = v
|
||||||
self.w = w
|
self.w = w
|
||||||
self.normalize()
|
self.normalize()
|
||||||
self.cached_r = None
|
|
||||||
|
|
||||||
def normalize(self):
|
def normalize(self):
|
||||||
self.n = cross(self.v, self.w, axis=0).T
|
self.n = cross(self.v, self.w, axis=0).T
|
||||||
|
|
@ -144,29 +127,15 @@ class plane:
|
||||||
def inside(self, p):
|
def inside(self, p):
|
||||||
return (self.n @ (p-self.o))[0] >= -0.00001
|
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):
|
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
|
self.M[:,2:3] = r.v
|
||||||
try:
|
try:
|
||||||
c=solve(self.M, r.o-self.o)[...,2,:]
|
c=solve(self.M, r.o-self.o)[...,2,:]
|
||||||
self.solution = (c, r.o - c*r.v)
|
|
||||||
except Exception as e:
|
except Exception as e:
|
||||||
if verbose:
|
if verbose:
|
||||||
stderr.write("solve: %s: %s\n" % (str(self.M), repr(e)))
|
stderr.write("solve: %s: %s\n" % (str(self.M), repr(e)))
|
||||||
self.solution = (-1, None)
|
return -1, None
|
||||||
return self.solution
|
return c, r.o - c*r.v
|
||||||
|
|
||||||
def plane_intersection(self, other, limit=1e10, delta=1e-5):
|
def plane_intersection(self, other, limit=1e10, delta=1e-5):
|
||||||
p = [ self.intersection(r)[1] for r in other.rays]
|
p = [ self.intersection(r)[1] for r in other.rays]
|
||||||
|
|
@ -279,7 +248,7 @@ class detector(list):
|
||||||
if edjes is None:
|
if edjes is None:
|
||||||
edjes = self.edges()
|
edjes = self.edges()
|
||||||
for e in edjes:
|
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):
|
def planes(self):
|
||||||
# is this useful?
|
# 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)),
|
plane(vector(z= z), vector(y=1), vector(x=1)),
|
||||||
]
|
]
|
||||||
if n2 is None:
|
if n2 is None:
|
||||||
n2 = n-1
|
n2 = n
|
||||||
else:
|
else:
|
||||||
c = cos(2*π/n * n1)
|
c = cos(2*π/n * n1)
|
||||||
s = sin(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.)
|
µ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*π):
|
def run(world, r, Δr, Δθ, θ1=0, θ2=π/2, φ1=0, φ2=2*π):
|
||||||
fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
|
fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
|
||||||
nr = int(r/Δr + 1.001)
|
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]
|
pl = [d.pathlength(ry) for d in world]
|
||||||
stdout.write(fmt % ((θ, φ, r, ψ)+tuple(pl)))
|
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"
|
fmt = "%g %g %g %g "+len(world)*" %g"+"\n"
|
||||||
beam = rays(r, Δr)
|
beam = rays(r, Δr)
|
||||||
b = None
|
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])
|
ll = tuple(bb+[l[i] for l in pl])
|
||||||
stdout.write(fmt % ll)
|
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 = [
|
test_rays = [
|
||||||
("x-axis", ray().ov(vector(),vector(x=1))),
|
("x-axis", ray().ov(vector(),vector(x=1))),
|
||||||
("y-axis", ray().ov(vector(),vector(y=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",
|
long_opt = ["help",
|
||||||
"φ1=", "φ2=", "θ1=", "θ2=", "φ₁=", "φ₂=", "θ₁=", "θ₂=",
|
"φ1=", "φ2=", "θ1=", "θ2=", "φ₁=", "φ₂=", "θ₁=", "θ₂=",
|
||||||
"Δθ=", "µMustang", "test", "size", "run", "corners", "gnuplot",
|
"Δθ=", "µMustang", "test", "size", "run", "corners", "gnuplot",
|
||||||
"radius=", "resolution=", "Δr=", "equivalent",
|
"radius=", "resolution=", "Δr=" ,
|
||||||
"detector", "copy=", "select=", "box=", "prism=",
|
"detector", "copy=", "select=", "box=", "prism=",
|
||||||
"plane=", "ray=",
|
"plane=", "ray=",
|
||||||
"rx=", "ry=", "rz=", "move=", "scale=",
|
"rx=", "ry=", "rz=", "move=", "scale=",
|
||||||
|
|
@ -474,13 +406,12 @@ world = [detector()]
|
||||||
θ1=0
|
θ1=0
|
||||||
θ2=π/2
|
θ2=π/2
|
||||||
φ1=0
|
φ1=0
|
||||||
φ2=2*π
|
φ2=2*π/8
|
||||||
do_size = False
|
do_size = False
|
||||||
do_test = False
|
do_test = False
|
||||||
do_run = False
|
do_run = False
|
||||||
do_corners = False
|
do_corners = False
|
||||||
do_gnuplot = False
|
do_gnuplot = False
|
||||||
do_equiv = False
|
|
||||||
rot = rotate(0,0)
|
rot = rotate(0,0)
|
||||||
ori = vector()
|
ori = vector()
|
||||||
|
|
||||||
|
|
@ -543,8 +474,6 @@ for o,v in options:
|
||||||
do_gnuplot=True
|
do_gnuplot=True
|
||||||
if o in "-E --corners":
|
if o in "-E --corners":
|
||||||
do_corners=True
|
do_corners=True
|
||||||
if o in "--equivalent":
|
|
||||||
do_equiv=True
|
|
||||||
if o in "--rx":
|
if o in "--rx":
|
||||||
world[-1] = world[-1].transform(rotate(0, aval(v)))
|
world[-1] = world[-1].transform(rotate(0, aval(v)))
|
||||||
if o in "--ry":
|
if o in "--ry":
|
||||||
|
|
@ -585,17 +514,13 @@ for o,v in options:
|
||||||
if k in vv[1]:
|
if k in vv[1]:
|
||||||
raise ValueError("recursive macro %s" % v)
|
raise ValueError("recursive macro %s" % v)
|
||||||
macros[k]=vv[1]
|
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]:
|
if not world[-1]:
|
||||||
world[-1:]=[]
|
world[-1:]=[]
|
||||||
if not world:
|
if not world:
|
||||||
world=[µM]
|
world=[µM]
|
||||||
|
|
||||||
if do_equiv:
|
|
||||||
find_equivalent_planes(world)
|
|
||||||
|
|
||||||
if not r or do_corners:
|
if not r or do_corners:
|
||||||
max_r, c = corners(world)
|
max_r, c = corners(world)
|
||||||
if not r:
|
if not r:
|
||||||
|
|
@ -606,10 +531,11 @@ if do_gnuplot:
|
||||||
d.gnuplot()
|
d.gnuplot()
|
||||||
|
|
||||||
if do_size:
|
if do_size:
|
||||||
no = int(π * (r/Δr)**2 + 0.5)
|
nr = int(r/Δr+1.001)
|
||||||
nv = int((φ2 - φ1) / Δθ**2 * (cos(θ1) - cos(θ2)) + 0.5)
|
no = π*nr**2 - (π-0.5)*nr
|
||||||
sr = nv * Δθ**2
|
sr = (φ2-φ1)*(cos(θ1)-cos(θ2+Δθ))
|
||||||
A = π * r**2
|
A = π*((nr-0.5)*Δr)**2
|
||||||
|
nv = sr/Δθ**2
|
||||||
nn = int(nv)*int(no)
|
nn = int(nv)*int(no)
|
||||||
stderr.write("""Run size estimate:
|
stderr.write("""Run size estimate:
|
||||||
detectors: %d planes: %d
|
detectors: %d planes: %d
|
||||||
|
|
|
||||||
110
seth.gf
110
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
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue