1 files changed, 87 insertions, 0 deletions
diff --git a/detectors/sno.py b/detectors/sno.py
new file mode 100644
index 0000000..5439739
--- /dev/null
+++ b/detectors/sno.py
@@ -0,0 +1,87 @@
+import numpy as np
+import numpy.linalg
+import math
+from make import rotate_extrude
+from geometry import *
+from optics import *
+from solids import build_8inch_pmt
+from transform import rotate, make_rotation_matrix
+import sno_pmt_locations
+
+def sno_vessel(sphere_radius, neck_radius, neck_top, angle_step=math.pi/20):
+    '''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
+       angle_step: angular step size (radius) between points on 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.arange(-math.pi/2, max_angle, angle_step)
+    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
+
+angle_step = math.pi/20
+av_outside_profile = sno_vessel(sphere_radius=6.0604, neck_radius=0.79375, 
+                                neck_top=10.50, angle_step=angle_step)
+# 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, angle_step=angle_step)
+
+av_outside_mesh = rotate_extrude(av_outside_profile[0], av_outside_profile[1],
+                                 angle_step)
+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],
+                                 angle_step)
+av_inside_mesh.vertices = rotate(av_inside_mesh.vertices, np.pi/2, (-1,0,0))
+
+def build_sno():
+    pmt = build_8inch_pmt()
+    geo = Geometry()
+    geo.add_solid(Solid(av_outside_mesh, glass, water, color=0x99FFFFFF))
+    geo.add_solid(Solid(av_inside_mesh, water, glass, color=0x990000FF))
+    geo.pmtids = []
+    for x,y,z in zip(sno_pmt_locations.x, sno_pmt_locations.y,
+                     sno_pmt_locations.z):
+        direction = numpy.array((-x,-y,-z))
+        direction /= numpy.linalg.norm(direction)
+        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)
+        displacement = (x/1000.0, y/1000.0, z/1000.0)
+        geo.pmtids.append(geo.add_solid(pmt, rotation, displacement))
+        
+    return geo