Create a r7081_cut mesh that uses the cut PMT model.

Modify LBNE class to take a cut_pmt option to switch to cut PMT model.
Add corresponding lbne_cut, minilbne_cut, microlbne_cut detectors.

4 files changed, 65 insertions, 5 deletions

diff --git a/detectors/__init__.py b/detectors/__init__.py
index 285480d..90f8c22 100644
--- a/detectors/__init__.py
+++ b/detectors/__init__.py
@@ -7,6 +7,11 @@ nstrings = 324
 pmts_per_string = 102
 endcap_spacing = .485
 
-lbne = LBNE(radius, height, nstrings, pmts_per_string, endcap_spacing)
 minilbne = LBNE(radius/10, height/10, nstrings//10, pmts_per_string//10, endcap_spacing)
 microlbne = LBNE(radius/50, height/50, nstrings//50, pmts_per_string//50, endcap_spacing)
+
+minilbne_cut = LBNE(radius/10, height/10, nstrings//10, pmts_per_string//10, endcap_spacing,
+                cut_pmt=True)
+microlbne_cut = LBNE(radius/50, height/50, nstrings//50, pmts_per_string//50, endcap_spacing,
+                 cut_pmt=True)
+lbne_cut = LBNE(radius, height, nstrings, pmts_per_string, endcap_spacing, cut_pmt=True)
diff --git a/detectors/lbne.py b/detectors/lbne.py
index 04751e6..c98fa42 100644
--- a/detectors/lbne.py
+++ b/detectors/lbne.py
@@ -14,9 +14,15 @@ from itertools import product
 import make
 
 class LBNE(Geometry):
-    def __init__(self, radius, height, nstrings, pmts_per_string, endcap_spacing):
+    def __init__(self, radius, height, nstrings, pmts_per_string, endcap_spacing,
+                 cut_pmt=False):
         super(LBNE, self).__init__()
 
+        if cut_pmt:
+            pmt_mesh = r7081_cut
+        else:
+            pmt_mesh = r7081
+
         # outer cylinder
         cylinder_mesh = \
             make.cylinder(radius, height+height/(pmts_per_string-1))
@@ -30,7 +36,7 @@ class LBNE(Geometry):
             for j in range(nstrings):
                 rotation = make_rotation_matrix(j*2*np.pi/nstrings, (0,0,1))
                 displacement = rotate((0,-radius,-height/2+i*height/(pmts_per_string-1)), j*2*np.pi/nstrings, (0,0,1))
-                self.pmtids.append(self.add_solid(r7081, rotation, displacement))
+                self.pmtids.append(self.add_solid(pmt_mesh, rotation, displacement))
 
         # construct the top endcap
         for x, y in np.array(tuple(product(\
@@ -39,7 +45,7 @@ class LBNE(Geometry):
             if np.sqrt(x**2 + y**2) <= radius:
                 rotation = make_rotation_matrix(+np.pi/2, (1,0,0))
                 displacement = (x,y,+height/2+height/(pmts_per_string-1)/2)
-                self.pmtids.append(self.add_solid(r7081, rotation, displacement))
+                self.pmtids.append(self.add_solid(pmt_mesh, rotation, displacement))
 
         # construct the bottom endcap
         for x, y in np.array(tuple(product(\
@@ -48,4 +54,4 @@ class LBNE(Geometry):
             if np.sqrt(x**2 + y**2) <= radius:
                 rotation = make_rotation_matrix(-np.pi/2, (1,0,0))
                 displacement = (x,y,-height/2-height/(pmts_per_string-1)/2)
-                self.pmtids.append(self.add_solid(r7081, rotation, displacement))
+                self.pmtids.append(self.add_solid(pmt_mesh, rotation, displacement))
diff --git a/solids/__init__.py b/solids/__init__.py
index f01d97f..8fc66a4 100644
--- a/solids/__init__.py
+++ b/solids/__init__.py
@@ -1 +1,2 @@
 from r7081 import r7081
+from r7081_cut import r7081_cut
diff --git a/solids/r7081_cut.py b/solids/r7081_cut.py
new file mode 100644
index 0000000..d38d54b
--- /dev/null
+++ b/solids/r7081_cut.py
@@ -0,0 +1,48 @@
+import os
+import sys
+import numpy as np
+
+dir = os.path.split(os.path.realpath(__file__))[0]
+sys.path.append(dir + '/..')
+
+import models
+from stl import mesh_from_stl
+from geometry import *
+from materials import *
+
+r7081_outer_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_outer_cut.stl')
+r7081_inner_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_inner_cut.stl')
+
+photocathode_triangles = np.mean(r7081_inner_mesh[:], axis=1)[:,1] > 0
+
+inner_color = np.empty(len(r7081_inner_mesh.triangles), np.uint32)
+inner_color[photocathode_triangles] = 0xff0000
+inner_color[~photocathode_triangles] = 0x00ff00
+
+inner_surface = np.empty(len(r7081_inner_mesh.triangles), np.object)
+inner_surface[photocathode_triangles] = black_surface
+inner_surface[~photocathode_triangles] = shiny_surface
+
+r7081_inner_solid = Solid(r7081_inner_mesh, vacuum, glass, inner_surface, color=inner_color)
+r7081_outer_solid = Solid(r7081_outer_mesh, glass, lightwater_sno)
+
+r7081_cut = r7081_inner_solid + r7081_outer_solid
+
+if __name__ == '__main__':
+    from view import view
+    from copy import deepcopy
+
+    r7081_outer_mesh_cutaway = deepcopy(r7081_outer_mesh)
+    r7081_outer_mesh_cutaway.triangles = \
+        r7081_outer_mesh_cutaway.triangles[\
+        np.mean(r7081_outer_mesh_cutaway[:], axis=1)[:,0] > 0]
+
+    r7081_outer_solid_cutaway = Solid(r7081_outer_mesh_cutaway, glass, lightwater_sno)
+
+    r7081_cutaway = r7081_inner_solid + r7081_outer_solid_cutaway
+
+    geometry = Geometry()
+    geometry.add_solid(r7081_cutaway)
+    geometry.build(bits=8)
+
+    view(geometry, 'r7081_cutaway')