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 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
|
||||
|
|
|
|||
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