Merge branch 'geometry_restructure' into 'master'

Geometry restructure

See merge request !37

km3buu/geometry.py

@@ -23,6 +23,18 @@ class DetectorVolume(ABC):
     def __init__(self):
         self._volume = -1.0
         self._coord_origin = (0., 0., 0.)
+        self._solid_angle = 1
+
+    @abstractmethod
+    def random_dir(self):
+        """
+        Generate a random direction for the interaction
+
+        Return
+        ------
+        tuple [rad, 1] (phi, cos(theta))
+        """
+        pass
 
     @abstractmethod
     def random_pos(self):
@@ -47,6 +59,17 @@ class DetectorVolume(ABC):
         """
         pass
 
+    @property
+    def solid_angle(self):
+        """
+        Solid angle used for the direction sampling
+
+        Returns
+        -------
+        float [1]
+        """
+        return self._solid_angle
+
     @property
     def volume(self):
         """
@@ -70,6 +93,26 @@ class DetectorVolume(ABC):
         return self._coord_origin
 
 
+class NoVolume(DetectorVolume):
+    """
+    Dummy volume to write out data w/o geometry 
+    """
+
+    def __init__(self):
+        self._solid_angle = 1
+        self._volume = 1
+        self._coord_origin = (.0, .0, .0)
+
+    def header_entries(self, nevents=0):
+        return dict()
+
+    def random_pos(self):
+        return (.0, .0, .0)
+
+    def random_dir(self):
+        return (0, 1)
+
+
 class CanVolume(DetectorVolume):
     """
     Cylindrical detector geometry
@@ -97,6 +140,7 @@ class CanVolume(DetectorVolume):
         self._zmax = zmax
         self._volume = self._calc_volume()
         self._detector_center = detector_center
+        self._solid_angle = 4 * np.pi
 
     def _calc_volume(self):
         return np.pi * (self._zmax - self._zmin) * np.power(self._radius, 2)
@@ -109,6 +153,11 @@ class CanVolume(DetectorVolume):
         pos_z = np.random.uniform(self._zmin, self._zmax)
         return (pos_x, pos_y, pos_z)
 
+    def random_dir(self):
+        phi = np.random.uniform(0, 2 * np.pi)
+        cos_theta = np.random.uniform(-1, 1)
+        return (phi, cos_theta)
+
     def header_entries(self, nevents=0):
         retdct = dict()
         key = "genvol"
@@ -123,7 +172,7 @@ class CanVolume(DetectorVolume):
 class SphericalVolume(DetectorVolume):
     """
     Spherical detector geometry
-    
+
     Parameters
     ----------
     radius: float [m]
@@ -138,6 +187,7 @@ class SphericalVolume(DetectorVolume):
         self._radius = radius
         self._coord_origin = coord_origin
         self._volume = self._calc_volume()
+        self._solid_angle = 4 * np.pi
 
     def _calc_volume(self):
         return 4 / 3 * np.pi * np.power(self._radius, 3)
@@ -152,6 +202,11 @@ class SphericalVolume(DetectorVolume):
         pos = (pos_x, pos_y, pos_z)
         return tuple(np.add(self._coord_origin, pos))
 
+    def random_dir(self):
+        phi = np.random.uniform(0, 2 * np.pi)
+        cos_theta = np.random.uniform(-1, 1)
+        return (phi, cos_theta)
+
     def header_entries(self, nevents=0):
         retdct = dict()
         key = "sphere"

km3buu/output.py

@@ -710,8 +710,10 @@ def write_detector_file(gibuu_output,
     targets_per_volume = target_density / target[
         0].atomic_weight / constants.atomic_mass
 
-    w2 = gibuu_output.w2weights(geometry.volume, targets_per_volume, 4 * np.pi)
-    global_generation_weight = gibuu_output.global_generation_weight(4 * np.pi)
+    w2 = gibuu_output.w2weights(geometry.volume, targets_per_volume,
+                                geometry.solid_angle)
+    global_generation_weight = gibuu_output.global_generation_weight(
+        geometry.solid_angle)
     mean_xsec_func = gibuu_output.mean_xsec
 
     header_dct = EMPTY_KM3NET_HEADER_DICT.copy()
@@ -787,8 +789,7 @@ def write_detector_file(gibuu_output,
             # Vertex Position
             vtx_pos = np.array(geometry.random_pos())
             # Direction
-            phi = np.random.uniform(0, 2 * np.pi)
-            cos_theta = np.random.uniform(-1, 1)
+            phi, cos_theta = geometry.random_dir()
             sin_theta = np.sqrt(1 - cos_theta**2)
 
             dir_x = np.cos(phi) * sin_theta