diff options
Diffstat (limited to 'detectors/sno.py')
| -rw-r--r-- | detectors/sno.py | 71 |
1 files changed, 58 insertions, 13 deletions
diff --git a/detectors/sno.py b/detectors/sno.py index 5439739..156dcbe 100644 --- a/detectors/sno.py +++ b/detectors/sno.py @@ -2,11 +2,12 @@ import numpy as np import numpy.linalg import math from make import rotate_extrude +from stl import mesh_from_stl from geometry import * from optics import * -from solids import build_8inch_pmt +from solids import build_8inch_pmt_with_lc from transform import rotate, make_rotation_matrix -import sno_pmt_locations +import os 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 @@ -66,22 +67,66 @@ 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() +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() - geo.add_solid(Solid(av_outside_mesh, glass, water, color=0x99FFFFFF)) - geo.add_solid(Solid(av_inside_mesh, water, glass, color=0x990000FF)) + + 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 = [] - 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) + + 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) - displacement = (x/1000.0, y/1000.0, z/1000.0) - geo.pmtids.append(geo.add_solid(pmt, rotation, displacement)) + + # 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 |