1 files changed, 58 insertions, 0 deletions
diff --git a/detectors/lbne.py b/detectors/lbne.py
new file mode 100644
index 0000000..9f74975
--- /dev/null
+++ b/detectors/lbne.py
@@ -0,0 +1,58 @@
+import os
+import sys
+import numpy as np
+
+dir = os.path.split(os.path.realpath(__file__))[0]
+sys.path.append(dir + '/..')
+
+import models
+from geometry import *
+from materials import *
+from solids import *
+from transform import rotate, make_rotation_matrix
+from itertools import product
+import make
+
+class LBNE(Geometry):
+    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, radius, height+height/(pmts_per_string-1))
+        cylinder_mesh.vertices = rotate(cylinder_mesh.vertices, np.pi/2, (-1,0,0))
+        cylinder_solid = Solid(cylinder_mesh, lightwater_sno, vacuum, black_surface, 0x0000ff)
+        self.add_solid(cylinder_solid)
+
+        self.pmtids = []
+
+        # 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((0,-radius,-height/2+i*height/(pmts_per_string-1)), j*2*np.pi/nstrings, (0,0,1))
+                self.pmtids.append(self.add_solid(pmt_mesh, rotation, displacement))
+
+        # 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, (1,0,0))
+                displacement = (x,y,+height/2+height/(pmts_per_string-1)/2)
+                self.pmtids.append(self.add_solid(pmt_mesh, rotation, displacement))
+
+        # 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, (1,0,0))
+                displacement = (x,y,-height/2-height/(pmts_per_string-1)/2)
+                self.pmtids.append(self.add_solid(pmt_mesh, rotation, displacement))