visually tested optics code. added models of the inner and outer meshes for the 12" hamamatsu and sno pmts. ratdb.py is able to parse ratdb files. chromaticity.py provides a function to map wavelength -> rgb color. lbne detector model now includes an outer black cylinder and pmts with a glass layer and photocathode/reflective surfaces.

1 files changed, 63 insertions, 12 deletions

diff --git a/detectors/lbne.py b/detectors/lbne.py
index 3c28c58..a439682 100644
--- a/detectors/lbne.py
+++ b/detectors/lbne.py
@@ -11,6 +11,8 @@ from solid import Solid
 from geometry import Geometry
 from transform import rotate, make_rotation_matrix
 from itertools import product
+from materials import *
+from make import cylinder
 
 endcap_spacing = .485
 
@@ -24,35 +26,84 @@ class LBNE(Geometry):
     def __init__(self):
         super(LBNE, self).__init__()
 
-        pmt_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch.stl')
-        pmt_mesh.vertices /= 1000.0
+        # outer cylinder
+        cylinder_solid = Solid(0, cylinder(radius, height+height/(pmts_per_string-1)), material1=lightwater_sno, material2=vacuum, surface=black_surface, color=0x0000ff)
+        self.add_solid(cylinder_solid)
 
-        pmtid = 0
+        pmt_inner_mesh = \
+            mesh_from_stl(models.dir + '/hamamatsu_12inch_inner.stl')
+        pmt_outer_mesh = \
+            mesh_from_stl(models.dir + '/hamamatsu_12inch_outer.stl')
+
+        photocathode_triangles = np.mean(pmt_inner_mesh[:], axis=1)[:,1] > 0
+
+        inner_color = np.empty(len(pmt_inner_mesh.triangles), np.uint32)
+        inner_color[photocathode_triangles] = 0xff0000
+        inner_color[~photocathode_triangles] = 0x00ff00
+
+        outer_color = np.empty(len(pmt_outer_mesh.triangles), np.uint32)
+        outer_color[:] = 0xffffff
+
+        inner_surface = np.empty(len(pmt_inner_mesh.triangles), np.object)
+        inner_surface[photocathode_triangles] = black_surface
+        inner_surface[~photocathode_triangles] = shiny_surface
+
+        self.pmtids = [1]
 
         # construct the barrel
         for i in range(pmts_per_string):
             for j in range(nstrings):
                 rotation = make_rotation_matrix(j*2*np.pi/nstrings, (0,0,1))
-                displacement = rotate((-radius,0,-height/2+i*height/(pmts_per_string-1)), j*2*np.pi/nstrings, (0,0,1))
-                self.add_solid(Solid(pmtid, pmt_mesh, rotation, displacement))
-                pmtid += 1
+                displacement = rotate((0,-radius,-height/2+i*height/(pmts_per_string-1)), j*2*np.pi/nstrings, (0,0,1))
+                self.add_solid(Solid(self.pmtids[-1], pmt_inner_mesh,
+                                     rotation, displacement,
+                                     material1=vacuum,
+                                     material2=glass,
+                                     surface=inner_surface,
+                                     color=inner_color))
+                self.add_solid(Solid(self.pmtids[-1], pmt_outer_mesh,
+                                     rotation, displacement,
+                                     material1=glass,
+                                     material2=lightwater_sno,
+                                     color=outer_color))
+                self.pmtids += [self.pmtids[-1]+1]
 
         # construct the top endcap
         for x, y in np.array(tuple(product(\
                     np.arange(-radius+0.075, radius, endcap_spacing),
                     np.arange(-radius+0.075, radius, endcap_spacing)))):
             if np.sqrt(x**2 + y**2) <= radius:
-                rotation = make_rotation_matrix(-np.pi/2, (0,1,0))
+                rotation = make_rotation_matrix(+np.pi/2, (1,0,0))
                 displacement = (x,y,+height/2+height/(pmts_per_string-1)/2)
-                self.add_solid(Solid(pmtid, pmt_mesh, rotation, displacement))
-                pmtid += 1
+                self.add_solid(Solid(self.pmtids[-1], pmt_inner_mesh,
+                                     rotation, displacement,
+                                     material1=vacuum,
+                                     material2=glass,
+                                     surface=inner_surface,
+                                     color=inner_color))
+                self.add_solid(Solid(self.pmtids[-1], pmt_outer_mesh,
+                                     rotation, displacement,
+                                     material1=glass,
+                                     material2=lightwater_sno,
+                                     color=outer_color))
+                self.pmtids += [self.pmtids[-1]+1]
 
         # construct the bottom endcap
         for x, y in np.array(tuple(product(\
                     np.arange(-radius+0.075, radius, endcap_spacing),
                     np.arange(-radius+0.075, radius, endcap_spacing)))):
             if np.sqrt(x**2 + y**2) <= radius:
-                rotation = make_rotation_matrix(+np.pi/2, (0,1,0))
+                rotation = make_rotation_matrix(-np.pi/2, (1,0,0))
                 displacement = (x,y,-height/2-height/(pmts_per_string-1)/2)
-                self.add_solid(Solid(pmtid, pmt_mesh, rotation, displacement))
-                pmtid += 1
+                self.add_solid(Solid(self.pmtids[-1], pmt_inner_mesh,
+                                     rotation, displacement,
+                                     material1=vacuum,
+                                     material2=glass,
+                                     surface=inner_surface,
+                                     color=inner_color))
+                self.add_solid(Solid(self.pmtids[-1], pmt_outer_mesh,
+                                     rotation, displacement,
+                                     material1=glass,
+                                     material2=lightwater_sno,
+                                     color=outer_color))
+                self.pmtids += [self.pmtids[-1]+1]