1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
import numpy as np
import numpy.linalg
import math
from chroma.make import rotate_extrude
from chroma.stl import mesh_from_stl
from chroma.geometry import *
from chroma.optics import *
from chroma.solids import build_8inch_pmt_with_lc
from chroma.transform import rotate, make_rotation_matrix
import os
def sno_vessel(sphere_radius, neck_radius, neck_top, nsteps=40):
'''Compute the 2D coordinates of the profile of one side of a
SNO-style acrylic vessel. The center of the sphere is at (0,0), with
the neck extending along the positive y direction.
sphere_radius: Radius of spherical part of profile
neck_radius: Radius of neck part
neck_top: y coordinate of top of neck
nsteps: number of points around the circumference of the sphere
Returns: Tuple of x and y coordinate numpy arrays.
'''
if neck_radius >= sphere_radius:
raise ValueError('neck_radius must be less than sphere_radius')
intersect_height = (sphere_radius**2 - neck_radius**2)**0.5
max_angle = math.atan2(intersect_height, neck_radius)
if neck_top < intersect_height:
raise ValueError('neck_top must be greater than the y-value where the sphere and cylinder intersect')
# Start with point at bottom
angles = np.linspace(-math.pi/2, max_angle, nsteps)
x = list(np.cos(angles) * sphere_radius)
y = list(np.sin(angles) * sphere_radius)
x[0] = 0.0 # Round-off error might make cos(-pi/2) not exactly zero
# Neck intersection point
x.append(neck_radius)
y.append(intersect_height)
# Top of neck
x.append(neck_radius)
y.append(neck_top)
# Top of neck on axis
x.append(0.0)
y.append(neck_top)
return x, y
##### SNO Parts
nsteps = 40
av_outside_profile = sno_vessel(sphere_radius=6.0604, neck_radius=0.79375,
neck_top=10.50, nsteps=nsteps)
# For simplicity, cap the top of the AV with acrylic
av_inside_profile = sno_vessel(sphere_radius=6.0053, neck_radius=0.72898,
neck_top=10.00, nsteps=nsteps)
av_outside_mesh = rotate_extrude(av_outside_profile[0], av_outside_profile[1],
nsteps)
av_outside_mesh.vertices = rotate(av_outside_mesh.vertices, np.pi/2, (-1,0,0))
av_inside_mesh = rotate_extrude(av_inside_profile[0], av_inside_profile[1],
nsteps)
av_inside_mesh.vertices = rotate(av_inside_mesh.vertices, np.pi/2, (-1,0,0))
dir = os.path.split(os.path.realpath(__file__))[0]
def build_sno(real_av=False):
import h5py
pmtinfo = h5py.File(dir+'/sno_pmt_info.hdf5')
pmt = build_8inch_pmt_with_lc()
geo = Geometry()
if real_av:
real_av_mesh = mesh_from_stl(dir+'/sno_av_ascii.stl.bz2')
real_av_mesh.vertices *= 0.0254 # inch -> meter
geo.add_solid(Solid(real_av_mesh, glass, water, color=0xBBAAAAFF))
else:
geo.add_solid(Solid(av_outside_mesh, acrylic_sno, water,
color=0xBBFFFFFF))
geo.add_solid(Solid(av_inside_mesh, labppo_scintillator, acrylic_sno,
color=0xBB0000FF))
geo.pmtids = []
snoman_id_dict = {}
pmt_offset = 122.0 - 2 * 0.01# max bucket height - 2 * facegap (mm)
for position, direction, idnum, pmt_type in zip(pmtinfo['pos'], pmtinfo['dir'],
pmtinfo['snoman_id'],
pmtinfo['pmt_type']):
# PMT type codes:
# 1 = normal; 2 = OWL; 3 = Low Gain; 4 = BUTT; 5 = Neck
# 6 = Calib channel; 10 = spare; 11 = invalid
if pmt_type != 1: # All the other PMTs have nonsense directions
continue
# Flip and renormalize
direction *= -1.0/numpy.linalg.norm(direction)
# Orient PMT that starts facing Y axis
y_axis = np.array((0.0,1.0,0.0))
axis = np.cross(direction, y_axis)
angle = np.arccos(np.dot(y_axis, direction))
rotation = make_rotation_matrix(angle, axis)
# Values in database are for front face of concentrator.
# need to shift so position is for equator of PMT
if pmt_type == 1:
displacement = position - direction * pmt_offset
displacement /= 1000.0 # mm -> m
chroma_id = geo.add_solid(pmt, rotation, displacement)
snoman_id_dict[chroma_id] = idnum
geo.pmtids.append(chroma_id)
# Convert dict to numpy array for fast array remapping later
chroma_ids = np.array(snoman_id_dict.keys())
snoman_ids = np.array(snoman_id_dict.values())
snoman_id_map = np.zeros(max(chroma_ids)+1) - 1 # Fill with -1 everywhere
snoman_id_map[chroma_ids] = snoman_ids
geo.snoman_id_map = snoman_id_map
return geo
|
