diff --git a/pyproject.toml b/pyproject.toml
index 4ba54b80..3b39c0ae 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -27,7 +27,7 @@ dependencies = [
   "Pygments",
   "siphash24; python_version < '3.13'",
   "pyopencl",
-  "vtk>=9.4",
+  "vtk>=9.4,<9.5",
   "mayavi-dev",
   "pyqt5",
 ]
@@ -67,3 +67,14 @@ ignore = [
 [tool.ruff.format]
 quote-style = "double"
 indent-style = "space"
+
+[tool.coverage.run]
+source = ["pytissueoptics"]
+omit = [
+    "*/tests/*",
+    "*/test*.py",
+    "*/examples/*",
+]
+
+[tool.coverage.report]
+show_missing = true
diff --git a/pytissueoptics/examples/rayscattering/ex06.py b/pytissueoptics/examples/rayscattering/ex06.py
new file mode 100644
index 00000000..9f1a74ae
--- /dev/null
+++ b/pytissueoptics/examples/rayscattering/ex06.py
@@ -0,0 +1,54 @@
+import env  # noqa: F401
+
+from pytissueoptics import *  # noqa: F403
+
+TITLE = "Detectors"
+
+DESCRIPTION = """ Sampling volume simulation with a directional source and two circle detectors. """
+
+
+def exampleCode():
+    import numpy as np
+
+    N = 1000000 if hardwareAccelerationIsAvailable() else 1000
+
+    material = ScatteringMaterial(mu_s=10, mu_a=1, g=0.98, n=1.0)
+
+    cube = Cuboid(a=3, b=3, c=3, position=Vector(0, 0, 0), material=material, label="cube")
+
+    detectorA = Circle(
+        radius=0.25,
+        orientation=Vector(0, 0, -1),
+        position=Vector(0, 1, 1.501),
+        label="detectorA",
+    ).asDetector(halfAngle=np.pi / 4)
+
+    detectorB = Circle(
+        radius=0.25,
+        orientation=Vector(0, 0, -1),
+        position=Vector(0, -1, 1.501),
+        label="detectorB",
+    ).asDetector(halfAngle=np.pi / 4)
+
+    scene = ScatteringScene([cube, detectorA, detectorB])
+
+    logger = EnergyLogger(scene)
+    source = DirectionalSource(position=Vector(0, 0, -2), direction=Vector(0, 0, 1), N=N, diameter=0.5)
+
+    source.propagate(scene, logger)
+
+    viewer = Viewer(scene, source, logger)
+    viewer.reportStats()
+
+    # Either filter the whole logger
+    logger.filter(detectedBy=["detectorA", "detectorB"])
+    viewer.show2D(View2DProjectionX())
+    viewer.show3D()
+
+    # ... or keep the logger intact and filter individual views
+    # viewer.show2D(View2DProjectionX(detectedBy="detectorA"))
+    # viewer.show3D(pointCloudStyle=PointCloudStyle(detectedBy=["detectorA", "detectorB"]))
+
+
+if __name__ == "__main__":
+    exampleCode()
diff --git a/pytissueoptics/examples/scene/example2.py b/pytissueoptics/examples/scene/example2.py
index b58415ee..bb3d3117 100644
--- a/pytissueoptics/examples/scene/example2.py
+++ b/pytissueoptics/examples/scene/example2.py
@@ -9,9 +9,9 @@
 def exampleCode():
     from pytissueoptics.scene import Cuboid, Vector, get3DViewer
 
-    cuboid1 = Cuboid(1, 1, 1, position=Vector(2, 0, 0))
-    cuboid2 = Cuboid(2, 1, 1, position=Vector(0, 2, 0))
-    cuboid3 = Cuboid(3, 1, 1, position=Vector(0, 0, 2))
+    cuboid1 = Cuboid(1, 1, 1, position=Vector(2, 0, 0), label="1")
+    cuboid2 = Cuboid(2, 1, 1, position=Vector(0, 2, 0), label="2")
+    cuboid3 = Cuboid(3, 1, 1, position=Vector(0, 0, 2), label="3")
 
     viewer = get3DViewer()
     viewer.add(cuboid1, cuboid2, cuboid3, representation="wireframe", lineWidth=5)
diff --git a/pytissueoptics/rayscattering/display/viewer.py b/pytissueoptics/rayscattering/display/viewer.py
index adf529cb..ba3b398f 100644
--- a/pytissueoptics/rayscattering/display/viewer.py
+++ b/pytissueoptics/rayscattering/display/viewer.py
@@ -66,6 +66,7 @@ def __init__(
         showSurfacePointsLeaving: bool = True,
         showSurfacePointsEntering: bool = False,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
         showPointsAsSpheres: bool = False,
         pointSize: float = 0.15,
         scaleWithValue: bool = True,
@@ -82,6 +83,7 @@ def __init__(
         self.showSurfacePointsLeaving = showSurfacePointsLeaving
         self.showSurfacePointsEntering = showSurfacePointsEntering
         self.energyType = energyType
+        self.detectedBy = detectedBy
         self.showPointsAsSpheres = showPointsAsSpheres
 
         self.pointSize = pointSize
@@ -103,7 +105,6 @@ def __init__(self, scene: ScatteringScene, source: Source, logger: EnergyLogger)
         self._logger = logger
 
         self._viewer3D = None
-        self._pointCloudFactory = PointCloudFactory(logger)
         self._profileFactory = ProfileFactory(scene, logger)
 
     def listViews(self):
@@ -169,7 +170,7 @@ def show3DVolumeSlicer(
                 f"Consider using a larger binSize or tighter limits."
             )
 
-        points = self._pointCloudFactory.getPointCloudOfSolids(energyType=energyType).solidPoints
+        points = PointCloudFactory(self._logger).getPointCloudOfSolids(energyType=energyType).solidPoints
         try:
             hist, _ = np.histogramdd(points[:, 1:], bins=bins, weights=points[:, 0], range=limits)
         except MemoryError:
@@ -215,8 +216,12 @@ def reportStats(self, solidLabel: str = None, saveToFile: str = None, verbose=Tr
         stats.report(solidLabel=solidLabel, saveToFile=saveToFile, verbose=verbose)
 
     def _addPointCloud(self, style: PointCloudStyle):
-        pointCloud = self._pointCloudFactory.getPointCloud(
-            solidLabel=style.solidLabel, surfaceLabel=style.surfaceLabel, energyType=style.energyType
+        logger = self._logger if style.detectedBy is None else self._logger.getFiltered(style.detectedBy)
+
+        pointCloud = PointCloudFactory(logger).getPointCloud(
+            solidLabel=style.solidLabel,
+            surfaceLabel=style.surfaceLabel,
+            energyType=style.energyType,
         )
 
         self._drawPointCloudOfSolids(pointCloud, style)
@@ -241,7 +246,7 @@ def _drawPointCloudOfSurfaces(self, pointCloud: PointCloud, style: PointCloudSty
         if pointCloud.surfacePoints is None:
             return
 
-        surfacePoints = np.empty((0, 4))
+        surfacePoints = np.empty((0, 5))
 
         if style.showSurfacePointsLeaving:
             surfacePoints = np.vstack((surfacePoints, pointCloud.leavingSurfacePoints))
diff --git a/pytissueoptics/rayscattering/display/views/defaultViews.py b/pytissueoptics/rayscattering/display/views/defaultViews.py
index a18092eb..18f0ebfa 100644
--- a/pytissueoptics/rayscattering/display/views/defaultViews.py
+++ b/pytissueoptics/rayscattering/display/views/defaultViews.py
@@ -21,6 +21,7 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             projectionDirection,
@@ -29,6 +30,7 @@ def __init__(
             limits=limits,
             binSize=binSize,
             energyType=energyType,
+            detectedBy=detectedBy,
         )
 
     def _filter(self, dataPoints: np.ndarray) -> np.ndarray:
@@ -42,9 +44,15 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
-            *DEFAULT_X_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType
+            *DEFAULT_X_VIEW_DIRECTIONS,
+            solidLabel=solidLabel,
+            limits=limits,
+            binSize=binSize,
+            energyType=energyType,
+            detectedBy=detectedBy,
         )
 
 
@@ -55,9 +63,15 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
-            *DEFAULT_Y_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType
+            *DEFAULT_Y_VIEW_DIRECTIONS,
+            solidLabel=solidLabel,
+            limits=limits,
+            binSize=binSize,
+            energyType=energyType,
+            detectedBy=detectedBy,
         )
 
 
@@ -68,9 +82,15 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
-            *DEFAULT_Z_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType
+            *DEFAULT_Z_VIEW_DIRECTIONS,
+            solidLabel=solidLabel,
+            limits=limits,
+            binSize=binSize,
+            energyType=energyType,
+            detectedBy=detectedBy,
         )
 
 
@@ -84,6 +104,7 @@ def __init__(
         surfaceEnergyLeaving: bool = True,
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             projectionDirection,
@@ -93,6 +114,7 @@ def __init__(
             surfaceEnergyLeaving=surfaceEnergyLeaving,
             limits=limits,
             binSize=binSize,
+            detectedBy=detectedBy,
         )
 
     def _filter(self, dataPoints: np.ndarray) -> np.ndarray:
@@ -116,6 +138,7 @@ def __init__(
         surfaceEnergyLeaving: bool = True,
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             *DEFAULT_X_VIEW_DIRECTIONS,
@@ -124,6 +147,7 @@ def __init__(
             surfaceEnergyLeaving=surfaceEnergyLeaving,
             limits=limits,
             binSize=binSize,
+            detectedBy=detectedBy,
         )
 
 
@@ -135,6 +159,7 @@ def __init__(
         surfaceEnergyLeaving: bool = True,
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             *DEFAULT_Y_VIEW_DIRECTIONS,
@@ -143,6 +168,7 @@ def __init__(
             surfaceEnergyLeaving=surfaceEnergyLeaving,
             limits=limits,
             binSize=binSize,
+            detectedBy=detectedBy,
         )
 
 
@@ -154,6 +180,7 @@ def __init__(
         surfaceEnergyLeaving: bool = True,
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             *DEFAULT_Z_VIEW_DIRECTIONS,
@@ -162,6 +189,7 @@ def __init__(
             surfaceEnergyLeaving=surfaceEnergyLeaving,
             limits=limits,
             binSize=binSize,
+            detectedBy=detectedBy,
         )
 
 
@@ -176,6 +204,7 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             projectionDirection,
@@ -186,6 +215,7 @@ def __init__(
             limits=limits,
             binSize=binSize,
             energyType=energyType,
+            detectedBy=detectedBy,
         )
 
     def setContext(self, limits3D: List[Tuple[float, float]], binSize3D: Tuple[float, float, float]):
@@ -210,6 +240,7 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             *DEFAULT_X_VIEW_DIRECTIONS,
@@ -219,6 +250,7 @@ def __init__(
             limits=limits,
             binSize=binSize,
             energyType=energyType,
+            detectedBy=detectedBy,
         )
 
 
@@ -231,6 +263,7 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             *DEFAULT_Y_VIEW_DIRECTIONS,
@@ -240,6 +273,7 @@ def __init__(
             limits=limits,
             binSize=binSize,
             energyType=energyType,
+            detectedBy=detectedBy,
         )
 
 
@@ -252,6 +286,7 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         super().__init__(
             *DEFAULT_Z_VIEW_DIRECTIONS,
@@ -261,4 +296,5 @@ def __init__(
             limits=limits,
             binSize=binSize,
             energyType=energyType,
+            detectedBy=detectedBy,
         )
diff --git a/pytissueoptics/rayscattering/display/views/view2D.py b/pytissueoptics/rayscattering/display/views/view2D.py
index e61e6de2..025d5c8a 100644
--- a/pytissueoptics/rayscattering/display/views/view2D.py
+++ b/pytissueoptics/rayscattering/display/views/view2D.py
@@ -50,6 +50,7 @@ def __init__(
         limits: Tuple[Tuple[float, float], Tuple[float, float]] = None,
         binSize: Union[float, Tuple[int, int]] = None,
         energyType=EnergyType.DEPOSITION,
+        detectedBy: Union[str, List[str]] = None,
     ):
         """
         The 2D view plane is obtained by looking towards the 'projectionDirection'. The 'horizontalDirection'
@@ -78,6 +79,7 @@ def __init__(
         self._position = position
         self._thickness = thickness
         self._energyType = energyType
+        self._detectedBy = detectedBy
 
         limits = [sorted(limit) for limit in limits] if limits else [None, None]
         self._limitsU, self._limitsV = limits
@@ -103,6 +105,7 @@ def setContext(self, limits3D: List[Tuple[float, float]], binSize3D: Tuple[float
 
         limits = [self._limitsU, self._limitsV]
         self._binsU, self._binsV = [int((limit[1] - limit[0]) / bin_) for limit, bin_ in zip(limits, self._binSize)]
+        self._binsU, self._binsV = max(1, self._binsU), max(1, self._binsV)
 
         if self._verticalIsNegative:
             self._limitsV = self._limitsV[::-1]
@@ -121,6 +124,10 @@ def setContext(self, limits3D: List[Tuple[float, float]], binSize3D: Tuple[float
     def energyType(self) -> EnergyType:
         return self._energyType
 
+    @property
+    def detectedBy(self) -> Union[str, List[str], None]:
+        return self._detectedBy
+
     def extractData(self, dataPoints: np.ndarray):
         """
         Used internally by Logger2D to store 3D datapoints into this 2D view.
@@ -245,6 +252,8 @@ def isEqualTo(self, other: "View2D") -> bool:
             return False
         if self._energyType != other._energyType:
             return False
+        if self._detectedBy != other._detectedBy:
+            return False
         return True
 
     def isContainedBy(self, other: "View2D") -> bool:
@@ -260,6 +269,8 @@ def isContainedBy(self, other: "View2D") -> bool:
             return False
         if self._thickness != other._thickness:
             return False
+        if self._detectedBy != other._detectedBy:
+            return False
 
         # TODO: change/remove the following once the algorithm can extract a view contained inside a bigger view.
         isTransposed = self.axisU != other.axisU
diff --git a/pytissueoptics/rayscattering/energyLogging/energyLogger.py b/pytissueoptics/rayscattering/energyLogging/energyLogger.py
index 98fda268..db43b650 100644
--- a/pytissueoptics/rayscattering/energyLogging/energyLogger.py
+++ b/pytissueoptics/rayscattering/energyLogging/energyLogger.py
@@ -1,7 +1,9 @@
+from __future__ import annotations
+
 import json
 import os
 import pickle
-from typing import List, Optional, TextIO, Union
+from typing import Dict, List, Optional, TextIO, Union
 
 import numpy as np
 
@@ -13,7 +15,7 @@
 from pytissueoptics.scene.logger.listArrayContainer import ListArrayContainer
 from pytissueoptics.scene.logger.logger import DataType, InteractionData, InteractionKey, Logger
 
-from ..opencl.CLScene import NO_SOLID_LABEL
+from ..opencl.CLScene import WORLD_SOLID_LABEL
 from .energyType import EnergyType
 
 
@@ -76,7 +78,7 @@ def addView(self, view: View2D) -> bool:
             return True
 
         if self.has3D:
-            self._compileViews([view])
+            self._compileViews([view], detectedBy=view.detectedBy)
             self._views.append(view)
             return True
 
@@ -233,11 +235,24 @@ def logDataPointArray(self, array: np.ndarray, key: InteractionKey):
             self._compileViews(self._views)
             self._delete3DData()
 
-    def logDataPoint(self, value: float, position: Vector, key: InteractionKey):
-        self.logDataPointArray(np.array([[value, *position.array]]), key)
+    def logDataPoint(self, value: float, position: Vector, key: InteractionKey, ID: Optional[int] = None):
+        dataPoint = [value, *position.array]
+        if ID is not None:
+            dataPoint.append(ID)
+        self.logDataPointArray(np.array([dataPoint]), key)
+
+    def _compileViews(self, views: List[View2D], detectedBy: Union[str, List[str]] = None):
+        if detectedBy is None:
+            dataPerInteraction = self._data
+            if any(view.detectedBy for view in views):
+                utils.warn(
+                    "Ignoring the detectedBy property of a view. Can only use detectedBy when adding a view after the "
+                    "simulation with keep3D=True."
+                )
+        else:
+            dataPerInteraction = self.getFiltered(detectedBy)._data
 
-    def _compileViews(self, views: List[View2D]):
-        for key, data in self._data.items():
+        for key, data in dataPerInteraction.items():
             datapointsContainer: Optional[ListArrayContainer] = data.dataPoints
             if datapointsContainer is None or len(datapointsContainer) == 0:
                 continue
@@ -304,6 +319,71 @@ def getDataPoints(self, key: InteractionKey, energyType=EnergyType.DEPOSITION) -
 
         return self._getData(DataType.DATA_POINT, key)
 
+    def filter(self, detectedBy: Union[str, List[str]]) -> None:
+        """Keeps only the data points from photons detected by one of the specified detector(s)."""
+        if not self._keep3D:
+            utils.warn("Cannot filter a logger that has discarded the 3D data.")
+            return
+
+        filteredPhotonIDs = self._getDetectedPhotonIDs(detectedBy)
+        self._data = self._getDataForPhotons(filteredPhotonIDs)
+        self._outdatedViews = set(self._views)
+
+    def getFiltered(self, detectedBy: Union[str, List[str]]) -> "EnergyLogger":
+        """
+        Returns a new logger with only data from photons detected by one of the specified detector(s).
+        """
+        filteredLogger = EnergyLogger(self._scene, views=[])
+        filteredPhotonIDs = self._getDetectedPhotonIDs(detectedBy)
+        filteredLogger._data = self._getDataForPhotons(filteredPhotonIDs)
+        return filteredLogger
+
+    def _getDetectedPhotonIDs(self, detectedBy: Union[str, List[str]]) -> np.ndarray:
+        """Helper to get photon IDs detected by one of the specified detector(s)."""
+        detector_labels = [detectedBy] if isinstance(detectedBy, str) else detectedBy
+        detector_keys = [InteractionKey(label) for label in detector_labels]
+        photonIDs = self._getPhotonIDs(detector_keys)
+        if len(photonIDs) == 0:
+            utils.warn(f"No photons detected by: {detectedBy}")
+        return photonIDs
+
+    def _getPhotonIDs(self, key: Union[InteractionKey, List[InteractionKey], None] = None) -> np.ndarray:
+        """Get all unique photon IDs that interacted with one of the given interaction key(s)."""
+        if isinstance(key, list):
+            all_photon_ids = []
+            for k in key:
+                data = self.getRawDataPoints(k)
+                if data is not None and data.shape[1] >= 5:
+                    all_photon_ids.append(data[:, 4].astype(np.uint32))
+
+            if len(all_photon_ids) == 0:
+                return np.array([], dtype=np.uint32)
+            combined = np.concatenate(all_photon_ids)
+            return np.unique(combined)
+
+        data = self.getRawDataPoints(key)
+        if data is None or data.shape[1] < 5:
+            return np.array([], dtype=np.uint32)
+        return np.unique(data[:, 4].astype(np.uint32))
+
+    def _getDataForPhotons(self, photonIDs: np.ndarray) -> Dict[InteractionKey, InteractionData]:
+        keyToData: Dict[InteractionKey, InteractionData] = {}
+        photonIDs = np.asarray(photonIDs, dtype=np.uint32)
+        for key, interactionData in self._data.items():
+            points: Optional[ListArrayContainer] = interactionData.dataPoints
+            if points is None:
+                continue
+            data = points.getData()
+            if data.shape[1] < 5:
+                continue
+            mask = np.isin(data[:, 4].astype(np.uint32), photonIDs)
+            filteredData = data[mask]
+            if filteredData.size > 0:
+                container = ListArrayContainer()
+                container.append(filteredData)
+                keyToData[key] = InteractionData(dataPoints=container)
+        return keyToData
+
     def _fluenceTransform(self, key: InteractionKey, data: Optional[np.ndarray]) -> Optional[np.ndarray]:
         # Converts volumetric data to fluence rate when needed.
         if not key.volumetric or data is None:
@@ -312,19 +392,19 @@ def _fluenceTransform(self, key: InteractionKey, data: Optional[np.ndarray]) ->
         data[:, 0] = data[:, 0] / self._scene.getMaterial(key.solidLabel).mu_a
         return data
 
-    def export(self, exportPath: str):
+    def export(self, exportName: str):
         """
         Export the raw 3D data points to a CSV file, along with the scene information to a JSON file.
 
-        The data file <exportPath>.csv will be comma-delimited and will contain the following columns:
-        - energy, x, y, z, solid_index, surface_index
+        The data file <exportName>.csv will be comma-delimited and will contain the following columns:
+        - energy, x, y, z, photon_index, solid_index, surface_index
 
         Two types of interactions are logged: scattering and surface crossings. In the first case, the energy will be
         the delta energy deposited at the point and the surface index will be -1. In the second case, the energy
         will be the total photon energy when crossing the surface, either as positive if leaving the surface
         (along the normal) or as negative if entering the surface.
 
-        The scene information will be saved in a JSON file named <exportPath>.json, which includes details for each solid
+        The scene information will be saved in a JSON file named <exportName>.json, which includes details for each solid
         index and surface index, such as their labels, materials, and geometry. The world information is also exported
         as solid index -1.
         """
@@ -341,16 +421,39 @@ def export(self, exportPath: str):
         solidLabels.sort()
 
         print("Exporting raw data to file...")
-        filepath = f"{exportPath}.csv"
+        filepath = f"{exportName}.csv"
         with open(filepath, "w") as file:
-            file.write("energy,x,y,z,solid_index,surface_index\n")
-            self._writeKeyData(file, InteractionKey(NO_SOLID_LABEL), -1, -1)
+            file.write("energy,x,y,z,photon_index,solid_index,surface_index\n")
+            self._writeKeyData(file, InteractionKey(WORLD_SOLID_LABEL), -1, -1)
             for i, solidLabel in enumerate(solidLabels):
                 self._writeKeyData(file, InteractionKey(solidLabel), i, -1)
                 for j, surfaceLabel in enumerate(self._scene.getSurfaceLabels(solidLabel)):
                     self._writeKeyData(file, InteractionKey(solidLabel, surfaceLabel), i, j)
         print(f"Exported data points to {filepath}")
 
+        self._exportSceneInfo(f"{exportName}.json", solidLabels)
+
+    def _writeKeyData(self, file: TextIO, key: InteractionKey, solidIndex: int, surfaceIndex: int):
+        if key not in self._data or self._data[key].dataPoints is None:
+            return
+
+        dataArray = self._data[key].dataPoints.getData()
+        n_rows = dataArray.shape[0]
+
+        output = np.empty((n_rows, 7), dtype=np.float64)
+        output[:, :4] = dataArray[:, :4]
+        output[:, 4] = dataArray[:, 4].astype(np.uint32)
+        output[:, 5] = solidIndex
+        output[:, 6] = surfaceIndex
+
+        np.savetxt(
+            file,
+            output,
+            delimiter=",",
+            fmt=["%.8e", "%.8e", "%.8e", "%.8e", "%d", "%d", "%d"],
+        )
+
+    def _exportSceneInfo(self, filepath: str, solidLabels: List[str]):
         sceneInfo = {}
         material = self._scene.getWorldEnvironment().material
         sceneInfo["-1"] = {"label": "world", "material": material.__dict__ if material else None}
@@ -374,15 +477,6 @@ def export(self, exportPath: str):
                 "surfaces": surfaces,
             }
 
-        sceneFilepath = f"{exportPath}.json"
-        with open(sceneFilepath, "w") as file:
+        with open(filepath, "w") as file:
             json.dump(sceneInfo, file, indent=4)
-        print(f"Exported scene information to {sceneFilepath}")
-
-    def _writeKeyData(self, file: TextIO, key: InteractionKey, solidIndex: int, surfaceIndex: int):
-        if key not in self._data or self._data[key].dataPoints is None:
-            return
-        dataArray = self._data[key].dataPoints.getData().astype(str)
-        dataArray = np.hstack((dataArray, np.full((dataArray.shape[0], 2), str(solidIndex))))
-        dataArray[:, 5] = str(surfaceIndex)
-        file.write("\n".join([",".join(row) for row in dataArray]) + "\n")
+        print(f"Exported scene information to {filepath}")
diff --git a/pytissueoptics/rayscattering/opencl/CLPhotons.py b/pytissueoptics/rayscattering/opencl/CLPhotons.py
index 7e793901..43ed16d4 100644
--- a/pytissueoptics/rayscattering/opencl/CLPhotons.py
+++ b/pytissueoptics/rayscattering/opencl/CLPhotons.py
@@ -54,6 +54,7 @@ def propagate(self, IPP: float, verbose: bool = False):
             self._directions[params.maxPhotonsPerBatch :],
             materialID=scene.getMaterialID(self._initialMaterial),
             solidID=scene.getSolidID(self._initialSolid),
+            startID=params.maxPhotonsPerBatch,
         )
         photonPool.make(program.device)
         seeds = SeedCL(params.maxPhotonsPerBatch)
diff --git a/pytissueoptics/rayscattering/opencl/CLScene.py b/pytissueoptics/rayscattering/opencl/CLScene.py
index 22793a3d..01de2db2 100644
--- a/pytissueoptics/rayscattering/opencl/CLScene.py
+++ b/pytissueoptics/rayscattering/opencl/CLScene.py
@@ -12,10 +12,10 @@
 from pytissueoptics.rayscattering.scatteringScene import ScatteringScene
 
 NO_LOG_ID = 0
-NO_SOLID_ID = -1
+WORLD_SOLID_ID = -1
 NO_SURFACE_ID = -1
 FIRST_SOLID_ID = 1
-NO_SOLID_LABEL = "world"
+WORLD_SOLID_LABEL = "world"
 
 
 class CLScene:
@@ -40,32 +40,35 @@ def __init__(self, scene: ScatteringScene, nWorkUnits: int):
         self.vertices = VertexCL(self._vertices)
 
     def getMaterialID(self, material):
+        if material is None:
+            # Detector case. Set dummy value (not used).
+            return 0
         return self._sceneMaterials.index(material)
 
     def getSolidID(self, solid):
         if solid is None:
-            return NO_SOLID_ID
+            return WORLD_SOLID_ID
         return self._solidLabels.index(solid.getLabel()) + FIRST_SOLID_ID
 
     def getSolidLabel(self, solidID):
-        if solidID == NO_SOLID_ID:
-            return NO_SOLID_LABEL
+        if solidID == WORLD_SOLID_ID:
+            return WORLD_SOLID_LABEL
         return self._solidLabels[solidID - FIRST_SOLID_ID]
 
     def getSolidIDs(self) -> List[int]:
         solidIDs = list(self._surfaceLabels.keys())
-        solidIDs.insert(0, NO_SOLID_ID)
+        solidIDs.insert(0, WORLD_SOLID_ID)
         return solidIDs
 
     def getSurfaceIDs(self, solidID):
-        if solidID == NO_SOLID_ID:
+        if solidID == WORLD_SOLID_ID:
             return [NO_SURFACE_ID]
         surfaceIDs = list(self._surfaceLabels[solidID].keys())
         surfaceIDs.insert(0, NO_SURFACE_ID)
         return surfaceIDs
 
     def getSurfaceLabel(self, solidID, surfaceID):
-        if solidID == NO_SOLID_ID:
+        if solidID == WORLD_SOLID_ID:
             return None
         if surfaceID == NO_SURFACE_ID:
             return None
@@ -100,6 +103,8 @@ def _compileSurface(self, polygonRef, firstPolygonID, lastPolygonID):
         insideSolidID = self.getSolidID(insideEnvironment.solid)
         outsideSolidID = self.getSolidID(outsideEnvironment.solid)
         toSmooth = polygonRef.toSmooth
+        isDetector = insideEnvironment.solid.isDetector if insideEnvironment.solid else False
+        detectorCosine = insideEnvironment.solid.detectorAcceptanceCosine if isDetector else 0.0
 
         self._surfacesInfo.append(
             SurfaceCLInfo(
@@ -110,6 +115,8 @@ def _compileSurface(self, polygonRef, firstPolygonID, lastPolygonID):
                 insideSolidID,
                 outsideSolidID,
                 toSmooth,
+                isDetector,
+                detectorCosine,
             )
         )
 
diff --git a/pytissueoptics/rayscattering/opencl/buffers/dataPointCL.py b/pytissueoptics/rayscattering/opencl/buffers/dataPointCL.py
index 5a3cfeb4..cc153956 100644
--- a/pytissueoptics/rayscattering/opencl/buffers/dataPointCL.py
+++ b/pytissueoptics/rayscattering/opencl/buffers/dataPointCL.py
@@ -12,6 +12,7 @@ class DataPointCL(CLObject):
             ("x", cl.cltypes.float),
             ("y", cl.cltypes.float),
             ("z", cl.cltypes.float),
+            ("photonID", cl.cltypes.uint),
             ("solidID", cl.cltypes.int),
             ("surfaceID", cl.cltypes.int),
         ]
diff --git a/pytissueoptics/rayscattering/opencl/buffers/photonCL.py b/pytissueoptics/rayscattering/opencl/buffers/photonCL.py
index dc6775c7..9a74782a 100644
--- a/pytissueoptics/rayscattering/opencl/buffers/photonCL.py
+++ b/pytissueoptics/rayscattering/opencl/buffers/photonCL.py
@@ -3,6 +3,8 @@
 
 from .CLObject import CLObject, cl
 
+NULL_SOLID_ID = 0
+
 
 class PhotonCL(CLObject):
     STRUCT_NAME = "Photon"
@@ -14,16 +16,21 @@ class PhotonCL(CLObject):
             ("weight", cl.cltypes.float),
             ("materialID", cl.cltypes.uint),
             ("solidID", cl.cltypes.int),
+            ("lastIntersectedDetectorID", cl.cltypes.int),
+            ("ID", cl.cltypes.uint),
         ]
     )
 
-    def __init__(self, positions: np.ndarray, directions: np.ndarray, materialID: int, solidID: int, weight=1.0):
+    def __init__(
+        self, positions: np.ndarray, directions: np.ndarray, materialID: int, solidID: int, weight=1.0, startID=0
+    ):
         self._positions = positions
         self._directions = directions
         self._N = positions.shape[0]
         self._materialID = materialID
         self._solidID = solidID
         self._weight = weight
+        self._startID = startID
 
         super().__init__()
 
@@ -36,4 +43,6 @@ def _getInitialHostBuffer(self) -> np.ndarray:
         buffer["weight"] = self._weight
         buffer["materialID"] = self._materialID
         buffer["solidID"] = self._solidID
+        buffer["lastIntersectedDetectorID"] = NULL_SOLID_ID
+        buffer["ID"] = np.arange(self._startID, self._startID + self._N, dtype=np.uint32)
         return buffer
diff --git a/pytissueoptics/rayscattering/opencl/buffers/surfaceCL.py b/pytissueoptics/rayscattering/opencl/buffers/surfaceCL.py
index 265aef05..c1346fc1 100644
--- a/pytissueoptics/rayscattering/opencl/buffers/surfaceCL.py
+++ b/pytissueoptics/rayscattering/opencl/buffers/surfaceCL.py
@@ -14,6 +14,8 @@
         ("insideSolidID", int),
         ("outsideSolidID", int),
         ("toSmooth", bool),
+        ("isDetector", bool),
+        ("detectorCosine", float),
     ],
 )
 
@@ -29,6 +31,8 @@ class SurfaceCL(CLObject):
             ("insideSolidID", cl.cltypes.int),
             ("outsideSolidID", cl.cltypes.int),
             ("toSmooth", cl.cltypes.uint),
+            ("isDetector", cl.cltypes.uint),
+            ("detectorCosine", cl.cltypes.float),
         ]
     )
 
@@ -47,4 +51,6 @@ def _getInitialHostBuffer(self) -> np.ndarray:
             buffer[i]["insideSolidID"] = np.int32(surfaceInfo.insideSolidID)
             buffer[i]["outsideSolidID"] = np.int32(surfaceInfo.outsideSolidID)
             buffer[i]["toSmooth"] = np.uint32(surfaceInfo.toSmooth)
+            buffer[i]["isDetector"] = np.uint32(surfaceInfo.isDetector)
+            buffer[i]["detectorCosine"] = np.float32(surfaceInfo.detectorCosine)
         return buffer
diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index 2d62b5fd..fecfa67e 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -115,13 +115,18 @@ GemsBoxIntersection _getBBoxIntersection(Ray ray, float3 minCornerVector, float3
     return intersection;
 }
 
-void _findBBoxIntersectingSolids(Ray ray, Scene *scene, uint gid, uint photonSolidID) {
+void _findBBoxIntersectingSolids(Ray ray, Scene *scene, uint gid, uint photonSolidID, uint ignoreSolidID) {
 
     for (uint i = 0; i < scene->nSolids; i++) {
         uint boxGID = gid * scene->nSolids + i;
         uint solidID = i + 1;
         scene->solidCandidates[boxGID].solidID = solidID;
 
+        if (solidID == ignoreSolidID) {
+            scene->solidCandidates[boxGID].distance = -1;
+            continue;
+        }
+
         if (solidID == photonSolidID) {
             scene->solidCandidates[boxGID].distance = 0;
             continue;
@@ -375,12 +380,12 @@ void _composeIntersection(Intersection *intersection, Ray *ray, Scene *scene) {
     intersection->distanceLeft = ray->length - intersection->distance;
 }
 
-Intersection findIntersection(Ray ray, Scene *scene, uint gid, uint photonSolidID) {
+Intersection findIntersection(Ray ray, Scene *scene, uint gid, uint photonSolidID, uint ignoreSolidID) {
     /*
     OpenCL implementation of the Python module SimpleIntersectionFinder
     See the Python module documentation for more details.
     */
-    _findBBoxIntersectingSolids(ray, scene, gid, photonSolidID);
+    _findBBoxIntersectingSolids(ray, scene, gid, photonSolidID, ignoreSolidID);
     _sortSolidCandidates(scene, gid);
 
     Intersection closestIntersection;
@@ -420,7 +425,7 @@ __kernel void findIntersections(__global Ray *rays, uint nSolids, __global Solid
         __global Triangle *triangles, __global Vertex *vertices, __global SolidCandidate *solidCandidates, __global Intersection *intersections) {
     uint gid = get_global_id(0);
     Scene scene = {nSolids, solids, surfaces, triangles, vertices, solidCandidates};
-    intersections[gid] = findIntersection(rays[gid], &scene, gid, -1);
+    intersections[gid] = findIntersection(rays[gid], &scene, gid, -1, 0);
 }
 
 
diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c
index 1a471e0e..d4e296bb 100644
--- a/pytissueoptics/rayscattering/opencl/src/propagation.c
+++ b/pytissueoptics/rayscattering/opencl/src/propagation.c
@@ -4,7 +4,8 @@
 #include "intersection.c"
 #include "fresnel.c"
 
-__constant int NO_SOLID_ID = -1;
+__constant int NULL_SOLID_ID = 0;
+__constant int WORLD_SOLID_ID = -1;
 __constant int NO_SURFACE_ID = -1;
 __constant float MIN_ANGLE = 0.0001f;
 
@@ -36,6 +37,7 @@ void interact(__global Photon *photons, __constant Material *materials, __global
     logger[logIndex].delta_weight = delta_weight;
     logger[logIndex].solidID = photons[photonID].solidID;
     logger[logIndex].surfaceID = NO_SURFACE_ID;
+    logger[logIndex].photonID = photons[photonID].ID;
 }
 
 void scatter(__global Photon *photons, __constant Material *materials, __global uint *seeds, __global DataPoint *logger,
@@ -79,6 +81,7 @@ void logIntersection(Intersection *intersection, __global Photon *photons, __glo
     logger[logID].z = photons[photonID].position.z;
     logger[logID].surfaceID = intersection->surfaceID;
     logger[logID].solidID = surfaces[intersection->surfaceID].insideSolidID;
+    logger[logID].photonID = photons[photonID].ID;
 
     bool isLeavingSurface = dot(photons[photonID].direction, intersection->normal) > 0;
     int sign = isLeavingSurface ? 1 : -1;
@@ -86,7 +89,7 @@ void logIntersection(Intersection *intersection, __global Photon *photons, __glo
     (*logIndex)++;
 
     int outsideSolidID = surfaces[intersection->surfaceID].outsideSolidID;
-    if (outsideSolidID == NO_SOLID_ID){
+    if (outsideSolidID == WORLD_SOLID_ID){
         return;
     }
     logID++;
@@ -96,9 +99,34 @@ void logIntersection(Intersection *intersection, __global Photon *photons, __glo
     logger[logID].surfaceID = intersection->surfaceID;
     logger[logID].solidID = outsideSolidID;
     logger[logID].delta_weight = -sign * photons[photonID].weight;
+    logger[logID].photonID = photons[photonID].ID;
     (*logIndex)++;
 }
 
+bool detectOrIgnore(Intersection *intersection, __global Photon *photons, __global Surface *surfaces,
+    __global DataPoint *logger, uint *logIndex, uint gid, uint photonID){
+    // If the incidence angle is within the numerical aperture, absorb photon.
+    float cosIncidence = -1 * dot(intersection->normal, photons[photonID].direction);
+    float cosDetector = surfaces[intersection->surfaceID].detectorCosine;
+
+    if (cosIncidence < cosDetector){
+        return false;  // Outside NA, ignore.
+    }
+
+    logger[*logIndex].x = photons[photonID].position.x;
+    logger[*logIndex].y = photons[photonID].position.y;
+    logger[*logIndex].z = photons[photonID].position.z;
+    logger[*logIndex].solidID = surfaces[intersection->surfaceID].insideSolidID;
+    logger[*logIndex].delta_weight = photons[photonID].weight;
+    logger[*logIndex].surfaceID = NO_SURFACE_ID;
+    logger[*logIndex].photonID = photons[photonID].ID;
+    (*logIndex)++;
+
+    // Absorb photon.
+    photons[photonID].weight = 0;
+    return true;
+}
+
 float reflectOrRefract(Intersection *intersection, __global Photon *photons, __constant Material *materials,
         __global Surface *surfaces, __global DataPoint *logger, uint *logIndex, __global uint *seeds, uint gid, uint photonID){
     FresnelIntersection fresnelIntersection = computeFresnelIntersection(photons[photonID].direction, intersection,
@@ -156,13 +184,27 @@ float propagateStep(float distance, __global Photon *photons, __constant Materia
     }
 
     Ray stepRay = {photons[photonID].position, photons[photonID].direction, distance};
-    Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID);
+    Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID, photons[photonID].lastIntersectedDetectorID);
+
+    photons[photonID].lastIntersectedDetectorID = NULL_SOLID_ID;  // Reset ignored detector ID.
 
     float distanceLeft = 0;
 
     if (intersection.exists){
         moveTo(intersection.position, photons, photonID);
-        distanceLeft = reflectOrRefract(&intersection, photons, materials, scene->surfaces, logger, logIndex, seeds, gid, photonID);
+        if (scene->surfaces[intersection.surfaceID].isDetector) {
+            if (detectOrIgnore(&intersection, photons, scene->surfaces, logger, logIndex, gid, photonID)) {;
+                return 0;  // Skip unnecessary vertex check if detected.
+            }
+
+            // Prevent re-intersecting with the same detector when passing through it.
+            photons[photonID].lastIntersectedDetectorID = scene->surfaces[intersection.surfaceID].insideSolidID;
+
+            // Skipping vertex check for now.
+            return intersection.distanceLeft;
+        } else {
+            distanceLeft = reflectOrRefract(&intersection, photons, materials, scene->surfaces, logger, logIndex, seeds, gid, photonID);
+        }
 
         // Check if intersection lies too close to a vertex.
         int closeToVertexID = -1;
diff --git a/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py b/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py
index fcf22250..7a2aeb57 100644
--- a/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py
+++ b/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py
@@ -3,18 +3,19 @@
 
 import numpy as np
 
-from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, NO_SOLID_LABEL, CLScene
+from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, WORLD_SOLID_LABEL, CLScene
 from pytissueoptics.scene.logger import InteractionKey, Logger
 
-SOLID_ID_COL = 4
-SURFACE_ID_COL = 5
+SOLID_ID_COL = 5
+SURFACE_ID_COL = 6
 
 
 class CLKeyLog:
-    """Parses a DataPointCL array of shape (N, 6) where each point is of the form
-    (weight, x, y, z, solidID, surfaceID) to extract a dictionary of InteractionKey
-    and their corresponding datapoint array of the form (weight, x, y, z). The
-    translation from IDs to their corresponding labels is done using the given CLScene."""
+    """Parses a DataPointCL array of shape (N, 7) where each point is of the form
+    (weight, x, y, z, photonID, solidID, surfaceID) to extract a dictionary of InteractionKey
+    and their corresponding datapoint array of the form (weight, x, y, z, photonID). The
+    translation from IDs to their corresponding labels is done using the given CLScene.
+    """
 
     def __init__(self, log: np.ndarray, sceneCL: CLScene):
         self._log = log
@@ -40,9 +41,9 @@ def _extractKeyLog(self):
 
     def _extractNoKeyLog(self):
         noInteractionIndices = np.where(self._log[:, SOLID_ID_COL] == NO_LOG_ID)[0]
-        self._log = self._log[:, :4]
+        self._log = self._log[:, :5]
         self._log = np.delete(self._log, noInteractionIndices, axis=0)
-        self._keyLog[InteractionKey(NO_SOLID_LABEL, None)] = self._log
+        self._keyLog[InteractionKey(WORLD_SOLID_LABEL, None)] = self._log
 
     def _sortLocal(self):
         """Sorts the log locally by solidID and surfaceID,
@@ -95,7 +96,7 @@ def _sort(log: np.ndarray, columns: List[int]):
 
     def _merge(self):
         """Merges the local batches into a single key log with unique interaction keys."""
-        self._log = self._log[:, :4]
+        self._log = self._log[:, :5]
 
         for i, batchKeyIndices in enumerate(self._keyIndices):
             batchStartIndex = i * self._batchSize
diff --git a/pytissueoptics/rayscattering/opencl/utils/CLParameters.py b/pytissueoptics/rayscattering/opencl/utils/CLParameters.py
index 00883007..21ba3a43 100644
--- a/pytissueoptics/rayscattering/opencl/utils/CLParameters.py
+++ b/pytissueoptics/rayscattering/opencl/utils/CLParameters.py
@@ -12,8 +12,8 @@ def __init__(self, N, AVG_IT_PER_PHOTON):
         nBatch = 1 / CONFIG.BATCH_LOAD_FACTOR
         avgPhotonsPerBatch = int(np.ceil(N / min(nBatch, CONFIG.N_WORK_UNITS)))
         self._maxLoggerMemory = self._calculateAverageBatchMemorySize(avgPhotonsPerBatch, AVG_IT_PER_PHOTON)
-        self._maxPhotonsPerBatch = min(2 * avgPhotonsPerBatch, N)
         self._workItemAmount = CONFIG.N_WORK_UNITS
+        self.maxPhotonsPerBatch = min(2 * avgPhotonsPerBatch, N)
 
         self._assertEnoughRAM()
 
diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index 523d5497..38c6a860 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -16,11 +16,12 @@
 
 
 class Photon:
-    def __init__(self, position: Vector, direction: Vector):
+    def __init__(self, position: Vector, direction: Vector, ID: int = 0):
         self._position = position
         self._direction = direction
         self._weight = 1
         self._environment: Environment = None
+        self._ID = ID
 
         self._er = self._direction.getAnyOrthogonal()
         self._hasContext = False
@@ -29,6 +30,8 @@ def __init__(self, position: Vector, direction: Vector):
         self._intersectionFinder: Optional[IntersectionFinder] = None
         self._logger: Optional[Logger] = None
 
+        self._lastIntersectedDetector: Optional[str] = None
+
     @property
     def isAlive(self) -> bool:
         return self._weight > 0
@@ -83,10 +86,23 @@ def step(self, distance=0) -> float:
                 distance = 0
 
         intersection = self._getIntersection(distance)
+        self._lastIntersectedDetector = None  # Reset ignored intersection label after each step.
 
         if intersection:
             self.moveTo(intersection.position)
-            distanceLeft = self.reflectOrRefract(intersection)
+            isDetector = intersection.insideEnvironment.solid.isDetector
+
+            if isDetector:
+                if self.detectOrIgnore(intersection):
+                    return 0  # Skip unnecessary vertex check if detected.
+
+                # Prevent re-intersecting with the same detector when passing through it.
+                self._lastIntersectedDetector = intersection.insideEnvironment.solid.getLabel()
+
+                # skipping vertex check if detector
+                return intersection.distanceLeft
+            else:
+                distanceLeft = self.reflectOrRefract(intersection)
 
             # Check if intersection lies too close to a vertex.
             for vertex in intersection.polygon.vertices:
@@ -118,7 +134,22 @@ def _getIntersection(self, distance) -> Optional[Intersection]:
             return None
 
         stepRay = Ray(self._position, self._direction, distance)
-        return self._intersectionFinder.findIntersection(stepRay, self.solidLabel)
+        return self._intersectionFinder.findIntersection(stepRay, self.solidLabel, self._lastIntersectedDetector)
+
+    def detectOrIgnore(self, intersection: Intersection) -> bool:
+        # If the incidence angle is within the numerical aperture, absorb photon.
+        cosIncidence = -intersection.normal.dot(self._direction)
+        cosDetector = intersection.insideEnvironment.solid.detectorAcceptanceCosine
+        if cosIncidence >= cosDetector:
+            self._detectAndLog(intersection.insideEnvironment.solidLabel)
+            return True
+        return False
+
+    def _detectAndLog(self, solidLabel: str):
+        if self._logger:
+            key = InteractionKey(solidLabel)
+            self._logger.logDataPoint(self._weight, self._position, key, self._ID)
+        self._weight = 0
 
     def reflectOrRefract(self, intersection: Intersection):
         fresnelIntersection = self._getFresnelIntersection(intersection)
@@ -212,16 +243,16 @@ def _logIntersection(self, intersection: Intersection):
         key = InteractionKey(solidLabelA, intersection.surfaceLabel)
         isLeavingSurface = self._direction.dot(intersection.normal) > 0
         sign = 1 if isLeavingSurface else -1
-        self._logger.logDataPoint(sign * self._weight, self._position, key)
+        self._logger.logDataPoint(sign * self._weight, self._position, key, self._ID)
 
         solidB = intersection.outsideEnvironment.solid
         if solidB is None:
             return
         solidLabelB = solidB.getLabel()
         key = InteractionKey(solidLabelB, intersection.surfaceLabel)
-        self._logger.logDataPoint(-sign * self._weight, self._position, key)
+        self._logger.logDataPoint(-sign * self._weight, self._position, key, self._ID)
 
     def _logWeightDecrease(self, delta):
         if self._logger:
             key = InteractionKey(self.solidLabel)
-            self._logger.logDataPoint(delta, self._position, key)
+            self._logger.logDataPoint(delta, self._position, key, self._ID)
diff --git a/pytissueoptics/rayscattering/scatteringScene.py b/pytissueoptics/rayscattering/scatteringScene.py
index 2cb05cbe..b068d90a 100644
--- a/pytissueoptics/rayscattering/scatteringScene.py
+++ b/pytissueoptics/rayscattering/scatteringScene.py
@@ -13,8 +13,11 @@ def __init__(self, solids: List[Solid], worldMaterial=ScatteringMaterial(), igno
         super().__init__(solids, worldMaterial=worldMaterial, ignoreIntersections=ignoreIntersections)
 
     def add(self, solid: Solid, position: Vector = None):
+        if solid.isFlat and not solid.isDetector:
+            raise Exception(f"Solid '{solid.getLabel()}' is flat. Flat solids must be used with asDetector().")
+
         polygonSample = solid.getPolygons()[0]
-        if not isinstance(polygonSample.insideEnvironment.material, ScatteringMaterial):
+        if not isinstance(polygonSample.insideEnvironment.material, ScatteringMaterial) and not solid.isDetector:
             raise Exception(
                 f"Solid '{solid.getLabel()}' has no ScatteringMaterial defined. "
                 f"This is required for any RayScatteringScene. "
diff --git a/pytissueoptics/rayscattering/source.py b/pytissueoptics/rayscattering/source.py
index 59c0c1f1..a1b0f289 100644
--- a/pytissueoptics/rayscattering/source.py
+++ b/pytissueoptics/rayscattering/source.py
@@ -138,7 +138,7 @@ def _loadPhotons(self):
     def _loadPhotonsCPU(self):
         positions, directions = self.getInitialPositionsAndDirections()
         for i in range(self._N):
-            self._photons.append(Photon(Vector(*positions[i]), Vector(*directions[i])))
+            self._photons.append(Photon(Vector(*positions[i]), Vector(*directions[i]), ID=i))
 
     def _loadPhotonsOpenCL(self):
         positions, directions = self.getInitialPositionsAndDirections()
diff --git a/pytissueoptics/rayscattering/statistics/statistics.py b/pytissueoptics/rayscattering/statistics/statistics.py
index 5416535b..4e59a874 100644
--- a/pytissueoptics/rayscattering/statistics/statistics.py
+++ b/pytissueoptics/rayscattering/statistics/statistics.py
@@ -1,3 +1,4 @@
+import math
 import os
 from dataclasses import dataclass
 from typing import Dict
@@ -7,7 +8,7 @@
 from pytissueoptics.rayscattering import utils
 from pytissueoptics.rayscattering.display.views.defaultViews import View2DProjection
 from pytissueoptics.rayscattering.energyLogging import EnergyLogger, PointCloud, PointCloudFactory
-from pytissueoptics.rayscattering.opencl.CLScene import NO_SOLID_LABEL
+from pytissueoptics.rayscattering.opencl.CLScene import WORLD_SOLID_LABEL
 
 
 @dataclass
@@ -49,11 +50,11 @@ def report(self, solidLabel: str = None, saveToFile: str = None, verbose=True):
 
     def _computeStats(self, solidLabel: str = None):
         solidLabels = [solidLabel]
-        if solidLabel is None or utils.labelsEqual(solidLabel, NO_SOLID_LABEL):
+        if solidLabel is None or utils.labelsEqual(solidLabel, WORLD_SOLID_LABEL):
             solidLabels = self._logger.getSeenSolidLabels()
 
         for solidLabel in solidLabels:
-            if solidLabel == NO_SOLID_LABEL:
+            if solidLabel == WORLD_SOLID_LABEL:
                 continue
             try:
                 absorbance = self.getAbsorbance(solidLabel)
@@ -69,7 +70,7 @@ def _computeStats(self, solidLabel: str = None):
 
     def _makeReport(self, solidLabel: str = None, reportString: str = ""):
         if solidLabel:
-            if solidLabel == NO_SOLID_LABEL:
+            if solidLabel == WORLD_SOLID_LABEL:
                 reportString += self._reportWorld(solidLabel)
             else:
                 reportString += self._reportSolid(solidLabel)
@@ -88,6 +89,11 @@ def _reportSolid(self, solidLabel: str):
         solidStats = self._solidStatsMap[solidLabel]
         reportString = "Report of solid '{}'\n".format(solidLabel)
 
+        if solidStats.absorbance == math.inf:
+            # Detectors don't log energy input.
+            reportString += "  Detected {:.2f}% of total power\n".format(solidStats.totalAbsorbance)
+            return reportString
+
         if solidStats.absorbance is None:
             reportString += "  Absorbance: N/A ({:.2f}% of total power)\n".format(solidStats.totalAbsorbance)
             reportString += "  Absorbance + Transmittance: N/A\n"
@@ -110,12 +116,12 @@ def getAbsorbance(self, solidLabel: str, useTotalEnergy=False) -> float:
             return self._getAbsorbanceFromViews(solidLabel, useTotalEnergy)
         points = self._getPointCloud(solidLabel).solidPoints
         energyInput = self.getEnergyInput(solidLabel) if not useTotalEnergy else self.getPhotonCount()
-        return 100 * self._sumEnergy(points) / energyInput
+        return 100 * self._sumEnergy(points) / energyInput if energyInput else math.inf
 
     def _getAbsorbanceFromViews(self, solidLabel: str, useTotalEnergy=False) -> float:
         energyInput = self.getEnergyInput(solidLabel) if not useTotalEnergy else self.getPhotonCount()
         absorbedEnergy = self._getAbsorbedEnergyFromViews(solidLabel)
-        return 100 * absorbedEnergy / energyInput
+        return 100 * absorbedEnergy / energyInput if energyInput else math.inf
 
     def _getAbsorbedEnergyFromViews(self, solidLabel: str) -> float:
         for view in self._logger.views:
@@ -207,7 +213,7 @@ def getTransmittance(self, solidLabel: str, surfaceLabel: str = None, useTotalEn
             points = self._getPointCloud(solidLabel, surfaceLabel).leavingSurfacePoints
 
         energyInput = self.getEnergyInput(solidLabel) if not useTotalEnergy else self.getPhotonCount()
-        return 100 * self._sumEnergy(points) / energyInput
+        return 100 * self._sumEnergy(points) / energyInput if energyInput else math.inf
 
     def _getTransmittanceFromViews(self, solidLabel: str, surfaceLabel: str = None, useTotalEnergy=False):
         if surfaceLabel is None:
@@ -216,7 +222,7 @@ def _getTransmittanceFromViews(self, solidLabel: str, surfaceLabel: str = None,
             energyLeaving = self._getSurfaceEnergyFromViews(solidLabel, surfaceLabel, leaving=True)
 
         energyInput = self.getEnergyInput(solidLabel) if not useTotalEnergy else self.getPhotonCount()
-        return 100 * energyLeaving / energyInput
+        return 100 * energyLeaving / energyInput if energyInput else math.inf
 
     @staticmethod
     def _sumEnergy(points: np.ndarray):
diff --git a/pytissueoptics/rayscattering/tests/display/testViewer.py b/pytissueoptics/rayscattering/tests/display/testViewer.py
index 468c4531..608276e2 100644
--- a/pytissueoptics/rayscattering/tests/display/testViewer.py
+++ b/pytissueoptics/rayscattering/tests/display/testViewer.py
@@ -8,7 +8,7 @@
 from pytissueoptics.rayscattering.display.profiles import Profile1D, ProfileFactory
 from pytissueoptics.rayscattering.display.viewer import PointCloudStyle, Viewer, Visibility
 from pytissueoptics.rayscattering.display.views import View2D
-from pytissueoptics.rayscattering.energyLogging import EnergyLogger, PointCloud, PointCloudFactory
+from pytissueoptics.rayscattering.energyLogging import EnergyLogger, PointCloud
 from pytissueoptics.rayscattering.scatteringScene import ScatteringScene
 from pytissueoptics.rayscattering.source import Source
 from pytissueoptics.scene.geometry import BoundingBox
@@ -112,14 +112,11 @@ def testWhenShow3DWithSource_shouldDisplaySource(self):
         self.mock3DViewer.show.assert_called_once()
 
     def testWhenShow3DWithDefaultPointCloud_shouldDisplayPointCloudOfSolidsAndSurfaceLeaving(self):
-        mockPointCloudFactory = mock(PointCloudFactory)
         aPointCloud = PointCloud(
-            solidPoints=np.array([[0.5, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0], [-1, 0, 0, 0]])
+            solidPoints=np.array([[0.5, 0, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0, 0], [-1, 0, 0, 0, 0]])
         )
-        when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud)
-        self.viewer._pointCloudFactory = mockPointCloudFactory
-
-        self.viewer.show3D(visibility=Visibility.POINT_CLOUD)
+        with self._mockPointCloud(aPointCloud):
+            self.viewer.show3D(visibility=Visibility.POINT_CLOUD)
 
         self.mock3DViewer.addDataPoints.assert_called()
         addedSolidPoints = self.mock3DViewer.addDataPoints.call_args_list[0][0][0]
@@ -130,25 +127,21 @@ def testWhenShow3DWithDefaultPointCloud_shouldDisplayPointCloudOfSolidsAndSurfac
         self.mock3DViewer.show.assert_called_once()
 
     def testGivenNoData_whenShow3DWithPointCloud_shouldNotDisplayPointCloud(self):
-        mockPointCloudFactory = mock(PointCloudFactory)
         aPointCloud = PointCloud()
-        when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud)
-        self.viewer._pointCloudFactory = mockPointCloudFactory
-
-        self.viewer.show3D(visibility=Visibility.POINT_CLOUD)
+        with self._mockPointCloud(aPointCloud):
+            self.viewer.show3D(visibility=Visibility.POINT_CLOUD)
 
         self.mock3DViewer.addDataPoints.assert_not_called()
         self.mock3DViewer.show.assert_called_once()
 
     def testWhenShow3DWithSurfacePointCloud_shouldOnlyDisplaySurfacePoints(self):
-        mockPointCloudFactory = mock(PointCloudFactory)
         aPointCloud = PointCloud(
-            solidPoints=np.array([[0.5, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0], [-1, 0, 0, 0]])
+            solidPoints=np.array([[0.5, 0, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0, 0], [-1, 0, 0, 0, 0]])
         )
-        when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud)
-        self.viewer._pointCloudFactory = mockPointCloudFactory
-
-        self.viewer.show3D(visibility=Visibility.POINT_CLOUD, pointCloudStyle=PointCloudStyle(showSolidPoints=False))
+        with self._mockPointCloud(aPointCloud):
+            self.viewer.show3D(
+                visibility=Visibility.POINT_CLOUD, pointCloudStyle=PointCloudStyle(showSolidPoints=False)
+            )
 
         self.mock3DViewer.addDataPoints.assert_called_once()
         self.assertTrue(
@@ -157,19 +150,16 @@ def testWhenShow3DWithSurfacePointCloud_shouldOnlyDisplaySurfacePoints(self):
         self.mock3DViewer.show.assert_called_once()
 
     def testWhenShow3DWithEnteringSurfacePointCloud_shouldOnlyDisplayEnteringSurfacePoints(self):
-        mockPointCloudFactory = mock(PointCloudFactory)
         aPointCloud = PointCloud(
-            solidPoints=np.array([[0.5, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0], [-1, 1, 1, 1]])
-        )
-        when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud)
-        self.viewer._pointCloudFactory = mockPointCloudFactory
-
-        self.viewer.show3D(
-            visibility=Visibility.POINT_CLOUD,
-            pointCloudStyle=PointCloudStyle(
-                showSolidPoints=False, showSurfacePointsLeaving=False, showSurfacePointsEntering=True
-            ),
+            solidPoints=np.array([[0.5, 0, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0, 0], [-1, 1, 1, 0, 1]])
         )
+        with self._mockPointCloud(aPointCloud):
+            self.viewer.show3D(
+                visibility=Visibility.POINT_CLOUD,
+                pointCloudStyle=PointCloudStyle(
+                    showSolidPoints=False, showSurfacePointsLeaving=False, showSurfacePointsEntering=True
+                ),
+            )
 
         self.mock3DViewer.addDataPoints.assert_called_once()
         self.assertTrue(
@@ -206,16 +196,13 @@ def testWhenShow3DWithViewsIndexList_shouldAdd2DImageOfTheseViewsInThe3DDisplay(
         self.mock3DViewer.show.assert_called_once()
 
     def testGiven3DLogger_whenShow3DDefault_shouldDisplayEverythingExceptViews(self):
-        mockPointCloudFactory = mock(PointCloudFactory)
         aPointCloud = PointCloud()
-        when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud)
-        self.viewer._pointCloudFactory = mockPointCloudFactory
 
-        self.viewer.show3D()
+        with self._mockPointCloud(aPointCloud):
+            self.viewer.show3D()
 
         verify(self.source, times=1).addToViewer(...)
         verify(self.scene, times=1).addToViewer(...)
-        verify(mockPointCloudFactory, times=1).getPointCloud(...)
         self.mock3DViewer.addImage.assert_not_called()
         self.mock3DViewer.show.assert_called_once()
 
@@ -223,15 +210,11 @@ def testGiven2DLogger_whenShow3DDefault_shouldDisplayEverythingExceptPointCloud(
         self._givenLoggerWithXSceneView()
         self.logger.has3D = False
 
-        mockPointCloudFactory = mock(PointCloudFactory)
-        when(mockPointCloudFactory).getPointCloud(...).thenReturn()
-        self.viewer._pointCloudFactory = mockPointCloudFactory
-
-        self.viewer.show3D()
+        with self._mockPointCloud(PointCloud()):
+            self.viewer.show3D()
 
         verify(self.source, times=1).addToViewer(...)
         verify(self.scene, times=1).addToViewer(...)
-        verify(mockPointCloudFactory, times=0).getPointCloud(...)
         self.mock3DViewer.addImage.assert_called()
         self.mock3DViewer.show.assert_called_once()
 
@@ -239,16 +222,12 @@ def testGiven2DLogger_whenShow3DWithDefault3DVisibility_shouldWarnAndDisplayDefa
         self._givenLoggerWithXSceneView()
         self.logger.has3D = False
 
-        mockPointCloudFactory = mock(PointCloudFactory)
-        when(mockPointCloudFactory).getPointCloud(...).thenReturn()
-        self.viewer._pointCloudFactory = mockPointCloudFactory
-
-        with self.assertWarns(UserWarning):
-            self.viewer.show3D(visibility=Visibility.DEFAULT_3D)
+        with self._mockPointCloud(PointCloud()):
+            with self.assertWarns(UserWarning):
+                self.viewer.show3D(visibility=Visibility.DEFAULT_3D)
 
         verify(self.source, times=1).addToViewer(...)
         verify(self.scene, times=1).addToViewer(...)
-        verify(mockPointCloudFactory, times=0).getPointCloud(...)
         self.mock3DViewer.addImage.assert_called()
         self.mock3DViewer.show.assert_called_once()
 
@@ -260,3 +239,10 @@ def _givenLoggerWithXSceneView(self):
         self.logger.views = [sceneView]
         when(self.logger).updateView(sceneView).thenReturn()
         when(self.scene).getBoundingBox().thenReturn(BoundingBox([-2, 2], [-2, 2], [0, 5]))
+
+    @staticmethod
+    def _mockPointCloud(pointCloud: PointCloud):
+        return patch(
+            "pytissueoptics.rayscattering.display.viewer.PointCloudFactory.getPointCloud",
+            return_value=pointCloud,
+        )
diff --git a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py
index 6b599cde..446d6e0d 100644
--- a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py
+++ b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py
@@ -18,7 +18,7 @@
 )
 from pytissueoptics.rayscattering.energyLogging import EnergyLogger
 from pytissueoptics.rayscattering.materials import ScatteringMaterial
-from pytissueoptics.rayscattering.opencl.CLScene import NO_SOLID_LABEL
+from pytissueoptics.rayscattering.opencl.CLScene import WORLD_SOLID_LABEL
 from pytissueoptics.rayscattering.samples import PhantomTissue
 from pytissueoptics.rayscattering.scatteringScene import ScatteringScene
 from pytissueoptics.scene.geometry import Vector
@@ -177,6 +177,17 @@ def testGiven2DLoggerWithData_whenShowWithCustomViewNotContainedByExistingViews_
 
         self.assertEqual(initialNumberOfViews, len(self.logger.views))
 
+    def testGiven3DLoggerWithData_whenAddViewWithDetectedBy_shouldInitializeViewWithFilteredData(self):
+        # Photon 0
+        self.logger.logDataPoint(0.9, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=0)
+        # Photon 1 interacted with "sphere".
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=1)
+        self.logger.logDataPoint(0.4, Vector(0.5, 0.5, 0.5), InteractionKey("sphere"), ID=1)
+
+        newView = View2DProjectionX(detectedBy="sphere")
+        self.logger.addView(newView)
+        self.assertAlmostEqual(0.1 + 0.4, newView.getSum())
+
     def testWhenSave_shouldSaveLoggerToFile(self):
         with tempfile.TemporaryDirectory() as tempDir:
             filePath = os.path.join(tempDir, "test.log")
@@ -296,10 +307,10 @@ def testWhenExport_shouldExport3DDataPointsToFile(self):
         self.logger = EnergyLogger(scene)
 
         # Log entering surface event, world scattering event and scattering event in both solids.
-        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("middleLayer"))
-        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("frontLayer", "interface1"))
-        self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"))
-        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(NO_SOLID_LABEL))
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("middleLayer"), ID=1)
+        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("frontLayer", "interface1"), ID=0)
+        self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"), ID=0)
+        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=0)
 
         with tempfile.TemporaryDirectory() as tempDir:
             filePath = os.path.join(tempDir, "test_sim")
@@ -309,12 +320,15 @@ def testWhenExport_shouldExport3DDataPointsToFile(self):
             with open(filePath + ".csv", "r") as f:
                 lines = f.readlines()
 
+            def parse_line(line):
+                return list(map(float, line.strip().split(",")))
+
             self.assertEqual(5, len(lines))
-            self.assertEqual("energy,x,y,z,solid_index,surface_index\n", lines[0])
-            self.assertEqual("0.2,0.0,0.0,0.0,-1,-1\n", lines[1])
-            self.assertEqual("-0.9,0.5,1.0,0.75,1,5\n", lines[2])
-            self.assertEqual("0.1,0.7,0.8,0.8,2,-1\n", lines[3])
-            self.assertEqual("0.4,0.0,5.0,0.0,3,-1\n", lines[4])
+            self.assertEqual("energy,x,y,z,photon_index,solid_index,surface_index\n", lines[0])
+            self.assertEqual(parse_line(lines[1]), [0.2, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0])
+            self.assertEqual(parse_line(lines[2]), [-0.9, 0.5, 1.0, 0.75, 0.0, 1.0, 5.0])
+            self.assertEqual(parse_line(lines[3]), [0.1, 0.7, 0.8, 0.8, 1.0, 2.0, -1.0])
+            self.assertEqual(parse_line(lines[4]), [0.4, 0.0, 5.0, 0.0, 0.0, 3.0, -1.0])
 
     def testWhenExport_shouldExportMetadataToFile(self):
         scene = PhantomTissue(worldMaterial=ScatteringMaterial(0.1, 0.1, 0.99))
@@ -349,3 +363,69 @@ def testWhenExport_shouldExportMetadataToFile(self):
         self.assertEqual(expectedLayerGeometry, sceneInfo["0"]["geometry"])
         self.assertEqual(16, len(sceneInfo["0"]["surfaces"]))
         self.assertEqual({"label": "interface0", "normal": [0.0, 0.0, -1.0]}, sceneInfo["0"]["surfaces"]["14"])
+
+    def testWhenFilterWithOneLabel_shouldOnlyKeepPhotonsThatInteractedWithTheGivenLabel(self):
+        # Photon 0
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=0)
+        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=0)
+        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=0)
+        # Photon 1 interacted with "sphere".
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=1)
+        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=1)
+        self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"), ID=1)
+        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=1)
+
+        self.logger.filter(detectedBy="sphere")
+
+        data = self.logger.getRawDataPoints()
+
+        # We should only have the data related to photon 1.
+        self.assertEqual(4, data.shape[0])
+        self.assertTrue(np.array_equal(np.unique(data[:, 4]), [1]))
+
+    def testWhenFilterWithMultipleLabels_shouldKeepPhotonsThatInteractedWithAnyOfTheGivenLabels(self):
+        # Photon 0
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=0)
+        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=0)
+        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=0)
+        # Photon 1 interacted with "sphere".
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=1)
+        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=1)
+        self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"), ID=1)
+        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=1)
+        # Photon 2 interacted with "otherSphere".
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=2)
+        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=2)
+        self.logger.logDataPoint(0.4, Vector(3, 5, 0), InteractionKey("otherSphere"), ID=2)
+        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=2)
+
+        self.logger.filter(detectedBy=["sphere", "otherSphere"])
+
+        data = self.logger.getRawDataPoints()
+
+        # We should only have the data related to photon 1 and photon 2.
+        self.assertEqual(8, data.shape[0])
+        self.assertTrue(np.array_equal(np.unique(data[:, 4]), [1, 2]))
+
+    def testWhenGetFiltered_shouldReturnANewFilteredLogger(self):
+        # Photon 0
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=0)
+        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=0)
+        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=0)
+        # Photon 1 interacted with "sphere".
+        self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=1)
+        self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=1)
+        self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"), ID=1)
+        self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=1)
+
+        filteredLogger = self.logger.getFiltered(detectedBy="sphere")
+
+        data = filteredLogger.getRawDataPoints()
+
+        # We should only have the data related to photon 1.
+        self.assertEqual(4, data.shape[0])
+        self.assertTrue(np.array_equal(np.unique(data[:, 4]), [1]))
+
+        # Should not modify the original logger.
+        originalData = self.logger.getRawDataPoints()
+        self.assertEqual(7, originalData.shape[0])
diff --git a/pytissueoptics/rayscattering/tests/opencl/src/testCLFresnel.py b/pytissueoptics/rayscattering/tests/opencl/src/testCLFresnel.py
index 0c34f98b..afcaef63 100644
--- a/pytissueoptics/rayscattering/tests/opencl/src/testCLFresnel.py
+++ b/pytissueoptics/rayscattering/tests/opencl/src/testCLFresnel.py
@@ -55,6 +55,8 @@ def _setUpWith(self, n1=1.0, n2=1.5, normal=Vector(0, 0, 1)):
                     self.INSIDE_SOLID_ID,
                     self.OUTSIDE_SOLID_ID,
                     False,
+                    False,
+                    0,
                 )
             ]
         )
diff --git a/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py b/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py
index 64c65846..92dccd66 100644
--- a/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py
+++ b/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py
@@ -21,7 +21,7 @@
     VertexCL,
 )
 from pytissueoptics.rayscattering.opencl.CLProgram import CLProgram
-from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, NO_SOLID_ID, NO_SURFACE_ID, CLScene
+from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, NO_SURFACE_ID, WORLD_SOLID_ID, CLScene
 from pytissueoptics.rayscattering.opencl.config.CLConfig import OPENCL_SOURCE_DIR
 from pytissueoptics.scene.geometry import Vertex
 
@@ -172,10 +172,10 @@ def testWhenLogIntersectionWithPhotonLeavingIntoWorld_shouldLogOnlyOneIntersecti
         self.INITIAL_DIRECTION.normalize()
         intersectionNormal = Vector(0, 0, 1)
         insideSolidID = self.INITIAL_SOLID_ID
-        outsideSolidID = NO_SOLID_ID
+        outsideSolidID = WORLD_SOLID_ID
 
         surfaceID = 0
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False, False, 0)])
         logger = DataPointCL(2)
 
         self._photonFunc("logIntersection", intersectionNormal, 0, surfaces, logger, 0)
@@ -196,13 +196,13 @@ def testWhenLogIntersectionWithPhotonEnteringFromWorld_shouldLogOnlyOneIntersect
     ):
         self.INITIAL_DIRECTION = Vector(0, 1, -1)
         self.INITIAL_DIRECTION.normalize()
-        self.INITIAL_SOLID_ID = NO_SOLID_ID
+        self.INITIAL_SOLID_ID = WORLD_SOLID_ID
         intersectionNormal = Vector(0, 0, 1)
         insideSolidID = 9
-        outsideSolidID = NO_SOLID_ID
+        outsideSolidID = WORLD_SOLID_ID
 
         surfaceID = 0
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False, False, 0)])
         logger = DataPointCL(2)
 
         self._photonFunc("logIntersection", intersectionNormal, surfaceID, surfaces, logger, 0)
@@ -225,7 +225,7 @@ def testWhenLogIntersectionBetweenTwoSolids_shouldLogIntersectionOnEachSolidWith
         insideSolidID = self.INITIAL_SOLID_ID
         outsideSolidID = self.INITIAL_SOLID_ID + 11
 
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False, False, 0)])
         logger = DataPointCL(2)
 
         self._photonFunc("logIntersection", intersectionNormal, 0, surfaces, logger, 0)
@@ -255,13 +255,13 @@ def testWhenInteract_shouldLogInteraction(self):
     def testWhenReflectOrRefractWithReflectingIntersection_shouldReflectPhoton(self):
         # => Photon is trying to enter solid
         intersectionNormal = Vector(0, 1, 0)
-        self.INITIAL_SOLID_ID = NO_SOLID_ID
+        self.INITIAL_SOLID_ID = WORLD_SOLID_ID
         self.INITIAL_DIRECTION = Vector(1, -1, 0)
         self.INITIAL_DIRECTION.normalize()
         self._mockFresnelIntersection(isReflected=True, incidencePlane=Vector(0, 0, 1), angleDeflection=np.pi / 2)
 
         logger = DataPointCL(2)
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=NO_SOLID_ID, toSmooth=False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, WORLD_SOLID_ID, False, False, 0)])
         photonResult = self._photonFunc(
             "reflectOrRefract",
             intersectionNormal,
@@ -279,7 +279,7 @@ def testWhenReflectOrRefractWithReflectingIntersection_shouldReflectPhoton(self)
         self._assertVectorAlmostEqual(expectedDirection, photonResult.direction, places=6)
         expectedPosition = self.INITIAL_POSITION
         self._assertVectorAlmostEqual(expectedPosition, photonResult.position)
-        self.assertEqual(NO_SOLID_ID, photonResult.solidID)
+        self.assertEqual(WORLD_SOLID_ID, photonResult.solidID)
 
         dataPoint = self._getDataPointResult(logger)
         self.assertEqual(NO_LOG_ID, dataPoint.solidID)
@@ -287,7 +287,7 @@ def testWhenReflectOrRefractWithReflectingIntersection_shouldReflectPhoton(self)
     def testWhenReflectOrRefractWithRefractingIntersection_shouldRefractPhoton(self):
         # => Photon enters solid
         intersectionNormal = Vector(0, 1, 0)
-        self.INITIAL_SOLID_ID = NO_SOLID_ID
+        self.INITIAL_SOLID_ID = WORLD_SOLID_ID
         self.INITIAL_DIRECTION = Vector(1, -1, 0)
         self.INITIAL_DIRECTION.normalize()
         insideSolidID = 9
@@ -301,7 +301,7 @@ def testWhenReflectOrRefractWithRefractingIntersection_shouldRefractPhoton(self)
         )
 
         logger = DataPointCL(2)
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID=NO_SOLID_ID, toSmooth=False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, WORLD_SOLID_ID, False, False, 0)])
         photonResult = self._photonFunc(
             "reflectOrRefract",
             intersectionNormal,
@@ -334,7 +334,7 @@ def testWhenStepToInfinityWithNoIntersection_shouldKillPhoton(self):
         self._mockFindIntersection(exists=False)
 
         logger = DataPointCL(2)
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)])
         photonResult = self._photonFunc(
             "propagateStep",
             stepDistance,
@@ -354,7 +354,7 @@ def testWhenStepWithNoIntersection_shouldMovePhotonAndScatter(self):
         self._mockFindIntersection(exists=False)
 
         logger = DataPointCL(2)
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)])
         photonResult = self._photonFunc(
             "propagateStep",
             stepDistance,
@@ -381,9 +381,7 @@ def testWhenStepWithIntersectionTooCloseToVertex_shouldMovePhotonAwayABitAndScat
         intersectionPosition = self.INITIAL_POSITION + self.INITIAL_DIRECTION * stepDistance
 
         logger = DataPointCL(2)
-        surfaces = SurfaceCL(
-            [SurfaceCLInfo(0, 0, 0, 0, insideSolidID=nextSolidID, outsideSolidID=self.INITIAL_SOLID_ID, toSmooth=False)]
-        )
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, nextSolidID, self.INITIAL_SOLID_ID, False, False, 0)])
         triangles = TriangleCL([TriangleCLInfo([0, 1, 2], normal)])
         # Put vertex slightly after the intersection point
         vertex = Vertex(intersectionPosition.x, intersectionPosition.y, intersectionPosition.z - 2e-7)
@@ -410,7 +408,7 @@ def testWhenStepWithNoDistance_shouldStepWithANewScatteringDistance(self):
         self._mockFindIntersection(exists=False)
 
         logger = DataPointCL(2)
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)])
         photonResult = self._photonFunc(
             "propagateStep",
             stepDistance,
@@ -432,7 +430,7 @@ def testWhenStepWithIntersectionReflecting_shouldMovePhotonToIntersection(self):
         self._mockFresnelIntersection(isReflected=True)
 
         logger = DataPointCL(2)
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)])
         triangles = TriangleCL([TriangleCLInfo([0, 1, 2], Vector(0, 0, 1))])
         vertices = VertexCL([Vertex(0, 0, 0)] * 3)
         photonResult = self._photonFunc(
@@ -457,7 +455,7 @@ def testWhenStepWithIntersectionRefracting_shouldMovePhotonToIntersection(self):
         self._mockFresnelIntersection(isReflected=False)
 
         logger = DataPointCL(2)
-        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)])
+        surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)])
         triangles = TriangleCL([TriangleCLInfo([0, 1, 2], Vector(0, 0, 1))])
         vertices = VertexCL([Vertex(0, 0, 0)] * 3)
         photonResult = self._photonFunc(
@@ -590,7 +588,7 @@ def _getPhotonResult(self, photonBuffer: PhotonCL):
 
     def _getDataPointResult(self, dataPointBuffer: DataPointCL, i=0):
         data = self.program.getData(dataPointBuffer)[i]
-        return DataPointResult(deltaWeight=data[0], position=Vector(*data[1:4]), solidID=data[4], surfaceID=data[5])
+        return DataPointResult(deltaWeight=data[0], position=Vector(*data[1:4]), solidID=data[5], surfaceID=data[6])
 
     def _assertVectorAlmostEqual(self, v1: Vector, v2: Vector, places=7):
         self.assertAlmostEqual(v1.x, v2.x, places=places)
@@ -644,9 +642,7 @@ def _mockFresnelIntersection(
 
     def _mockFindIntersection(self, exists=True, distance=8.0, normal=Vector(0, 0, 1), surfaceID=0, distanceLeft=2):
         expectedPosition = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance
-        intersectionCall = (
-            """Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID);"""
-        )
+        intersectionCall = """Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID, photons[photonID].lastIntersectedDetectorID);"""
         px, py, pz = expectedPosition.array
         nx, ny, nz = normal.array
         mockCall = """Intersection intersection;
diff --git a/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py b/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py
index 9e27708c..5f81db63 100644
--- a/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py
+++ b/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py
@@ -5,7 +5,13 @@
 
 from pytissueoptics import Cube, EnergyLogger, ScatteringMaterial, ScatteringScene, Sphere
 from pytissueoptics.rayscattering.opencl import OPENCL_OK
-from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, NO_SOLID_ID, NO_SOLID_LABEL, NO_SURFACE_ID, CLScene
+from pytissueoptics.rayscattering.opencl.CLScene import (
+    NO_LOG_ID,
+    NO_SURFACE_ID,
+    WORLD_SOLID_ID,
+    WORLD_SOLID_LABEL,
+    CLScene,
+)
 from pytissueoptics.rayscattering.opencl.utils import CLKeyLog
 from pytissueoptics.scene.logger import InteractionKey
 
@@ -27,13 +33,13 @@ def _createTestLog(self, sceneCL):
 
         log = np.array(
             [
-                [0, 0, 0, 0, cubeID, NO_SURFACE_ID],
-                [2, 0, 0, 0, sphereID, NO_SURFACE_ID],
-                [3, 0, 0, 0, NO_SOLID_ID, NO_SURFACE_ID],
-                [4, 9, 9, 9, NO_LOG_ID, 9],
-                [5, 0, 0, 0, cubeID, cubeSurfaceIDs[1]],
-                [6, 0, 0, 0, sphereID, sphereSurfaceIDs[1]],
-                [1, 0, 0, 0, cubeID, NO_SURFACE_ID],
+                [0, 0, 0, 0, 0, cubeID, NO_SURFACE_ID],
+                [2, 0, 0, 0, 0, sphereID, NO_SURFACE_ID],
+                [3, 0, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID],
+                [4, 9, 9, 9, 0, NO_LOG_ID, 9],
+                [5, 0, 0, 0, 0, cubeID, cubeSurfaceIDs[1]],
+                [6, 0, 0, 0, 0, sphereID, sphereSurfaceIDs[1]],
+                [1, 0, 0, 0, 0, cubeID, NO_SURFACE_ID],
             ]
         )
         return log
@@ -49,17 +55,17 @@ def testGivenCLKeyLog_whenTransferToSceneLogger_shouldLogDataWithInteractionKeys
 
         verify(sceneLogger, times=5).logDataPointArray(...)
 
-        expectedCubeData = arg_that(lambda arg: np.array_equal(arg, np.array([[0, 0, 0, 0], [1, 0, 0, 0]])))
+        expectedCubeData = arg_that(lambda arg: np.array_equal(arg, np.array([[0, 0, 0, 0, 0], [1, 0, 0, 0, 0]])))
         verify(sceneLogger).logDataPointArray(expectedCubeData, InteractionKey(self.cube.getLabel()))
 
         expectedValuesWithKeys = [
             (2, InteractionKey(self.sphere.getLabel())),
-            (3, InteractionKey(NO_SOLID_LABEL)),
+            (3, InteractionKey(WORLD_SOLID_LABEL)),
             (5, InteractionKey(self.cube.getLabel(), self.cube.surfaceLabels[0])),
             (6, InteractionKey(self.sphere.getLabel(), self.sphere.surfaceLabels[0])),
         ]
         for value, expectedKey in expectedValuesWithKeys:
-            expectedData = arg_that(lambda arg: np.array_equal(arg, np.array([[value, 0, 0, 0]])))
+            expectedData = arg_that(lambda arg: np.array_equal(arg, np.array([[value, 0, 0, 0, 0]])))
             verify(sceneLogger).logDataPointArray(expectedData, expectedKey)
 
     def testGivenCLKeyLogForInfiniteScene_whenTransferToSceneLogger_shouldLogDataWithInteractionKeys(self):
@@ -67,9 +73,9 @@ def testGivenCLKeyLogForInfiniteScene_whenTransferToSceneLogger_shouldLogDataWit
         sceneCL = CLScene(self.scene, nWorkUnits=10)
         log = np.array(
             [
-                [1, 0, 0, 0, NO_SOLID_ID, NO_SURFACE_ID],
-                [2, 0, 0, 0, NO_LOG_ID, 99],
-                [3, 0, 0, 0, NO_SOLID_ID, NO_SURFACE_ID],
+                [1, 0, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID],
+                [2, 0, 0, 0, 0, NO_LOG_ID, 99],
+                [3, 0, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID],
             ]
         )
         clKeyLog = CLKeyLog(log, sceneCL)
@@ -79,5 +85,5 @@ def testGivenCLKeyLogForInfiniteScene_whenTransferToSceneLogger_shouldLogDataWit
         clKeyLog.toSceneLogger(sceneLogger)
 
         verify(sceneLogger, times=1).logDataPointArray(...)
-        expectedWorldData = arg_that(lambda arg: np.array_equal(arg, np.array([[1, 0, 0, 0], [3, 0, 0, 0]])))
-        verify(sceneLogger).logDataPointArray(expectedWorldData, InteractionKey(NO_SOLID_LABEL))
+        expectedWorldData = arg_that(lambda arg: np.array_equal(arg, np.array([[1, 0, 0, 0, 0], [3, 0, 0, 0, 0]])))
+        verify(sceneLogger).logDataPointArray(expectedWorldData, InteractionKey(WORLD_SOLID_LABEL))
diff --git a/pytissueoptics/rayscattering/tests/testPhoton.py b/pytissueoptics/rayscattering/tests/testPhoton.py
index 592701ac..90f96817 100644
--- a/pytissueoptics/rayscattering/tests/testPhoton.py
+++ b/pytissueoptics/rayscattering/tests/testPhoton.py
@@ -31,8 +31,10 @@ def setUp(self):
         self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy())
 
         self.solidInside = mock(Solid)
+        self.solidInside.isDetector = False
         when(self.solidInside).getLabel().thenReturn(self.SOLID_INSIDE_LABEL)
         self.solidOutside = mock(Solid)
+        self.solidOutside.isDetector = False
         when(self.solidOutside).getLabel().thenReturn(self.SOLID_OUTSIDE_LABEL)
 
     def testShouldBeInTheGivenState(self):
@@ -310,7 +312,7 @@ def testGivenALogger_whenSteppingOutsideASolid_shouldLogPositiveWeightOnSurfaceI
 
         intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance
         interactionKey = InteractionKey(self.SOLID_INSIDE_LABEL, self.SURFACE_LABEL)
-        verify(logger).logDataPoint(self.photon.weight, intersectionPoint, interactionKey)
+        verify(logger).logDataPoint(self.photon.weight, intersectionPoint, interactionKey, 0)
 
     def testGivenALogger_whenSteppingOutsideASolid_shouldLogNegativeWeightOnOtherSolidSurfaceIntersection(self):
         distance = 8
@@ -318,7 +320,7 @@ def testGivenALogger_whenSteppingOutsideASolid_shouldLogNegativeWeightOnOtherSol
 
         intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance
         interactionKey = InteractionKey(self.SOLID_OUTSIDE_LABEL, self.SURFACE_LABEL)
-        verify(logger).logDataPoint(-self.photon.weight, intersectionPoint, interactionKey)
+        verify(logger).logDataPoint(-self.photon.weight, intersectionPoint, interactionKey, 0)
 
     def givenALogger_whenSteppingInsideASolidAt(self, distance):
         logger = self._createLogger()
@@ -338,7 +340,7 @@ def testGivenALogger_whenSteppingInsideASurface_shouldLogNegativeWeightAtInterse
 
         intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance
         interactionKey = InteractionKey(self.SOLID_INSIDE_LABEL, self.SURFACE_LABEL)
-        verify(logger).logDataPoint(-self.photon.weight, intersectionPoint, interactionKey)
+        verify(logger).logDataPoint(-self.photon.weight, intersectionPoint, interactionKey, 0)
 
     def testGivenALogger_whenSteppingInsideASurface_shouldLogPositiveWeightOnOtherSolidSurfaceIntersection(self):
         distance = 8
@@ -346,7 +348,7 @@ def testGivenALogger_whenSteppingInsideASurface_shouldLogPositiveWeightOnOtherSo
 
         intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance
         interactionKey = InteractionKey(self.SOLID_OUTSIDE_LABEL, self.SURFACE_LABEL)
-        verify(logger).logDataPoint(self.photon.weight, intersectionPoint, interactionKey)
+        verify(logger).logDataPoint(self.photon.weight, intersectionPoint, interactionKey, 0)
 
     def testGivenALogger_whenScatter_shouldLogWeightLossAtThisPositionWithSolidLabel(self):
         solidInside = mock(Solid)
@@ -357,7 +359,7 @@ def testGivenALogger_whenScatter_shouldLogWeightLossAtThisPositionWithSolidLabel
         self.photon.scatter()
 
         weightLoss = self.photon.material.getAlbedo()
-        verify(logger).logDataPoint(weightLoss, self.INITIAL_POSITION, InteractionKey(self.SOLID_INSIDE_LABEL))
+        verify(logger).logDataPoint(weightLoss, self.INITIAL_POSITION, InteractionKey(self.SOLID_INSIDE_LABEL), 0)
 
     def testGivenALogger_whenScatterInWorldMaterial_shouldLogWeightLossAtThisPositionWithWorldLabel(self):
         logger = self._createLogger()
@@ -366,7 +368,7 @@ def testGivenALogger_whenScatterInWorldMaterial_shouldLogWeightLossAtThisPositio
         self.photon.scatter()
 
         weightLoss = self.photon.material.getAlbedo()
-        verify(logger).logDataPoint(weightLoss, self.INITIAL_POSITION, InteractionKey(WORLD_LABEL))
+        verify(logger).logDataPoint(weightLoss, self.INITIAL_POSITION, InteractionKey(WORLD_LABEL), 0)
 
     def givenALoggerAndNoSolidOutside_whenSteppingInsideASolidAt(self, distance):
         self.solidOutside = None
@@ -421,6 +423,89 @@ def testWhenInteractWithWeightAtFloatLimit_shouldKillPhoton(self):
         self.photon.interact()
         self.assertFalse(self.photon.isAlive)
 
+    def testWhenIntersectingDetectorInsideNA_shouldKillPhoton(self):
+        # Mock detector intersection
+        halfAngle = math.pi / 8  # 22.5 degrees
+        detector = mock(Solid)
+        detector.isDetector = True
+        detector.detectorAcceptanceCosine = math.cos(halfAngle)
+        when(detector).getLabel().thenReturn("detector")
+        self.solidInside = self.solidOutside = detector
+
+        # Setup photon going towards detector within NA.
+        photonIncidence = halfAngle
+        self.INITIAL_DIRECTION = Vector(0, math.sin(photonIncidence), -math.cos(photonIncidence))
+        self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy())
+
+        # Setup intersection with detector.
+        distance = 8
+        intersectionFinder = self._createIntersectionFinder(distance)
+        self.photon.setContext(Environment(ScatteringMaterial()), intersectionFinder=intersectionFinder)
+
+        self.photon.step(distance + 2)
+
+        expectedIntersectionPosition = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance
+        self.assertVectorEqual(expectedIntersectionPosition, self.photon.position)
+        self.assertFalse(self.photon.isAlive)
+
+    def testWhenIntersectingDetectorOutsideNA_shouldNotKillPhoton(self):
+        # Mock detector intersection
+        halfAngle = math.pi / 8  # 22.5 degrees
+        detector = mock(Solid)
+        detector.isDetector = True
+        detector.detectorAcceptanceCosine = math.cos(halfAngle)
+        when(detector).getLabel().thenReturn("detector")
+        self.solidInside = self.solidOutside = detector
+
+        # Setup photon going towards detector just outside NA.
+        photonIncidence = halfAngle + math.pi / 100  # Just outside NA
+        self.INITIAL_DIRECTION = Vector(0, math.sin(photonIncidence), -math.cos(photonIncidence))
+        self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy())
+
+        # Setup intersection with detector.
+        distance = 8
+        intersectionFinder = self._createIntersectionFinder(distance)
+        self.photon.setContext(Environment(ScatteringMaterial()), intersectionFinder=intersectionFinder)
+
+        self.photon.step(distance + 2)
+
+        expectedIntersectionPosition = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance
+        self.assertVectorEqual(expectedIntersectionPosition, self.photon.position)
+        self.assertTrue(self.photon.isAlive)
+        self.assertEqual(self.photon.solidLabel, WORLD_LABEL)  # Still in world
+
+    def testGivenALogger_whenDetected_shouldLogWeightAtThisPositionWithDetectorLabel(self):
+        # Mock detector intersection
+        detectorLabel = "detector"
+        halfAngle = math.pi / 8  # 22.5 degrees
+        detector = mock(Solid)
+        detector.isDetector = True
+        detector.detectorAcceptanceCosine = math.cos(halfAngle)
+        when(detector).getLabel().thenReturn(detectorLabel)
+        self.solidInside = self.solidOutside = detector
+
+        # Setup photon going towards detector within NA.
+        photonIncidence = halfAngle  # Within NA
+        self.INITIAL_DIRECTION = Vector(0, math.sin(photonIncidence), -math.cos(photonIncidence))
+        self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy(), ID=42)
+        initial_weight = self.photon.weight
+
+        # Setup intersection with detector.
+        distance = 8
+        intersectionFinder = self._createIntersectionFinder(distance)
+        logger = self._createLogger()
+        self.photon.setContext(
+            Environment(ScatteringMaterial()),
+            intersectionFinder=intersectionFinder,
+            logger=logger,
+        )
+
+        self.photon.step(distance + 2)
+
+        self.assertFalse(self.photon.isAlive)
+        intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance
+        verify(logger).logDataPoint(initial_weight, intersectionPoint, InteractionKey(detectorLabel), 42)
+
     def assertVectorEqual(self, v1, v2):
         self.assertAlmostEqual(v1.x, v2.x, places=7)
         self.assertAlmostEqual(v1.y, v2.y, places=7)
diff --git a/pytissueoptics/rayscattering/tests/testScatteringScene.py b/pytissueoptics/rayscattering/tests/testScatteringScene.py
index e2a9a4e4..739eedbc 100644
--- a/pytissueoptics/rayscattering/tests/testScatteringScene.py
+++ b/pytissueoptics/rayscattering/tests/testScatteringScene.py
@@ -4,6 +4,7 @@
 
 from mockito import mock, verify, when
 
+from pytissueoptics import Circle
 from pytissueoptics.rayscattering.materials import ScatteringMaterial
 from pytissueoptics.rayscattering.scatteringScene import ScatteringScene
 from pytissueoptics.scene.solids import Cuboid
@@ -54,3 +55,12 @@ def testShouldHaveIPPEstimationUsingMeanAlbedoInInfiniteMedium(self):
         estimation = scene.getEstimatedIPP(weightThreshold)
         expectedEstimation = -math.log(weightThreshold) / meanAlbedo
         self.assertAlmostEqual(expectedEstimation, estimation, places=7)
+
+    def testCannotAddFlatSolidThatIsNotADetector(self):
+        flatSolid = Circle(radius=1)
+        with self.assertRaises(Exception):
+            ScatteringScene([flatSolid])
+
+    def testCanAddFlatDetectors(self):
+        flatSolid = Circle(radius=1).asDetector()
+        ScatteringScene([flatSolid])
diff --git a/pytissueoptics/scene/__init__.py b/pytissueoptics/scene/__init__.py
index c3edf9b6..7d9ea20c 100644
--- a/pytissueoptics/scene/__init__.py
+++ b/pytissueoptics/scene/__init__.py
@@ -4,6 +4,7 @@
 from .material import RefractiveMaterial
 from .scene import Scene
 from .solids import (
+    Circle,
     Cone,
     Cube,
     Cuboid,
@@ -11,6 +12,7 @@
     Ellipsoid,
     PlanoConcaveLens,
     PlanoConvexLens,
+    Rectangle,
     Sphere,
     SymmetricLens,
     ThickLens,
@@ -37,4 +39,6 @@
     "SymmetricLens",
     "PlanoConvexLens",
     "PlanoConcaveLens",
+    "Rectangle",
+    "Circle",
 ]
diff --git a/pytissueoptics/scene/geometry/utils.py b/pytissueoptics/scene/geometry/utils.py
index 0d61d4c7..6bdaf8ca 100644
--- a/pytissueoptics/scene/geometry/utils.py
+++ b/pytissueoptics/scene/geometry/utils.py
@@ -55,7 +55,13 @@ def getAxisAngleBetween(fromDirection: Vector, toDirection: Vector) -> Tuple[Vec
     dot = fromDirection.dot(toDirection)
     dot = max(min(dot, 1), -1)
     angle = np.arccos(dot)
-    axis = fromDirection.cross(toDirection)
+
+    # Special case: vectors are opposite (antiparallel)
+    if np.isclose(dot, -1):
+        axis = fromDirection.getAnyOrthogonal()
+    else:
+        axis = fromDirection.cross(toDirection)
+
     axis.normalize()
     axis = axis.array
     return Vector(*axis), angle
diff --git a/pytissueoptics/scene/intersection/intersectionFinder.py b/pytissueoptics/scene/intersection/intersectionFinder.py
index 2be19ec5..f66a98bf 100644
--- a/pytissueoptics/scene/intersection/intersectionFinder.py
+++ b/pytissueoptics/scene/intersection/intersectionFinder.py
@@ -35,11 +35,13 @@ def __init__(self, scene: Scene):
         self._polygonIntersect = MollerTrumboreIntersect()
         self._boxIntersect = GemsBoxIntersect()
 
-    def findIntersection(self, ray: Ray, currentSolidLabel: Optional[str]) -> Optional[Intersection]:
+    def findIntersection(
+        self, ray: Ray, currentSolidLabel: Optional[str], ignoreLabel: Optional[str] = None
+    ) -> Optional[Intersection]:
         raise NotImplementedError
 
     def _findClosestPolygonIntersection(
-        self, ray: Ray, polygons: List[Polygon], currentSolidLabel: str
+        self, ray: Ray, polygons: List[Polygon], currentSolidLabel: str, ignoreLabel: Optional[str] = None
     ) -> Optional[Intersection]:
         closestIntersection = Intersection(sys.maxsize)
         minSameSolidDistance = -sys.maxsize
@@ -52,6 +54,9 @@ def _findClosestPolygonIntersection(
             if insideLabel != currentSolidLabel and outsideLabel != currentSolidLabel:
                 continue
 
+            if ignoreLabel and insideLabel == ignoreLabel:
+                continue
+
             intersectionPoint = self._polygonIntersect.getIntersection(ray, polygon)
             if not intersectionPoint:
                 continue
@@ -107,7 +112,9 @@ def _composeIntersection(ray: Ray, intersection: Intersection) -> Optional[Inter
 
 
 class SimpleIntersectionFinder(IntersectionFinder):
-    def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Intersection]:
+    def findIntersection(
+        self, ray: Ray, currentSolidLabel: str, ignoreLabel: Optional[str] = None
+    ) -> Optional[Intersection]:
         """
         Find the closest intersection between a ray and the scene.
 
@@ -119,7 +126,7 @@ def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Interse
             If the solid bbox distance is greater than the closest intersection distance, then we can stop testing since
             they are ordered by distance. Note that bbox distance is zero if the ray starts inside.
         """
-        bboxIntersections = self._findBBoxIntersectingSolids(ray, currentSolidLabel)
+        bboxIntersections = self._findBBoxIntersectingSolids(ray, currentSolidLabel, ignoreLabel)
         bboxIntersections.sort(key=lambda x: x[0])
 
         closestDistance = sys.maxsize
@@ -134,12 +141,16 @@ def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Interse
 
         return self._composeIntersection(ray, closestIntersection)
 
-    def _findBBoxIntersectingSolids(self, ray: Ray, currentSolidLabel: str) -> Optional[List[Tuple[float, Solid]]]:
+    def _findBBoxIntersectingSolids(
+        self, ray: Ray, currentSolidLabel: str, ignoreLabel: Optional[str] = None
+    ) -> Optional[List[Tuple[float, Solid]]]:
         """We need to handle the special case where ray starts inside bbox. The Box Intersect will not compute
         the intersection for this case and will instead return ray.origin. When that happens, distance will be 0,
         and we continue to check for possibly other contained solids."""
         solidCandidates = []
         for solid in self._scene.solids:
+            if ignoreLabel and solid.getLabel() == ignoreLabel:
+                continue
             if solid.getLabel() == currentSolidLabel:
                 bboxIntersectionPoint = ray.origin
             else:
@@ -158,16 +169,22 @@ def __init__(self, scene: Scene, constructor=NoSplitThreeAxesConstructor(), maxD
             self._scene.getBoundingBox(), self._scene.getPolygons(), constructor, maxDepth, minLeafSize
         )
 
-    def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Intersection]:
+    def findIntersection(
+        self, ray: Ray, currentSolidLabel: str, ignoreLabel: Optional[str] = None
+    ) -> Optional[Intersection]:
         self._currentSolidLabel = currentSolidLabel
-        intersection = self._findIntersection(ray, self._partition.root)
+        intersection = self._findIntersection(ray, self._partition.root, ignoreLabel=ignoreLabel)
         return self._composeIntersection(ray, intersection)
 
-    def _findIntersection(self, ray: Ray, node: Node, closestDistance=sys.maxsize) -> Optional[Intersection]:
+    def _findIntersection(
+        self, ray: Ray, node: Node, closestDistance=sys.maxsize, ignoreLabel: Optional[str] = None
+    ) -> Optional[Intersection]:
         # todo: test if this node search is still compatible with the new core intersection logic
         #  which can now require testing a polygon slightly behind the ray origin.
         if node.isLeaf:
-            intersection = self._findClosestPolygonIntersection(ray, node.polygons, self._currentSolidLabel)
+            intersection = self._findClosestPolygonIntersection(
+                ray, node.polygons, self._currentSolidLabel, ignoreLabel
+            )
             return intersection
 
         if not self._nodeIsWorthExploring(ray, node, closestDistance):
@@ -175,7 +192,7 @@ def _findIntersection(self, ray: Ray, node: Node, closestDistance=sys.maxsize) -
 
         closestIntersection = None
         for child in node.children:
-            intersection = self._findIntersection(ray, child, closestDistance)
+            intersection = self._findIntersection(ray, child, closestDistance, ignoreLabel)
             if intersection is None:
                 continue
             if intersection.distance < closestDistance:
diff --git a/pytissueoptics/scene/logger/logger.py b/pytissueoptics/scene/logger/logger.py
index 2e813795..8f701bb0 100644
--- a/pytissueoptics/scene/logger/logger.py
+++ b/pytissueoptics/scene/logger/logger.py
@@ -78,8 +78,11 @@ def getStoredSurfaceLabels(self, solidLabel: str) -> List[str]:
     def logPoint(self, point: Vector, key: InteractionKey = None):
         self._appendData([point.x, point.y, point.z], DataType.POINT, key)
 
-    def logDataPoint(self, value: float, position: Vector, key: InteractionKey):
-        self._appendData([value, position.x, position.y, position.z], DataType.DATA_POINT, key)
+    def logDataPoint(self, value: float, position: Vector, key: InteractionKey, ID: Optional[int] = None):
+        dataPoint = [value, *position.array]
+        if ID is not None:
+            dataPoint.append(ID)
+        self._appendData(dataPoint, DataType.DATA_POINT, key)
 
     def logSegment(self, start: Vector, end: Vector, key: InteractionKey = None):
         self._appendData([start.x, start.y, start.z, end.x, end.y, end.z], DataType.SEGMENT, key)
@@ -90,8 +93,9 @@ def logPointArray(self, array: np.ndarray, key: InteractionKey = None):
         self._appendData(array, DataType.POINT, key)
 
     def logDataPointArray(self, array: np.ndarray, key: InteractionKey):
-        """'array' must be of shape (n, 4) where the second axis is (value, x, y, z)"""
-        assert array.shape[1] == 4 and array.ndim == 2, "Data point array must be of shape (n, 4)"
+        """'array' must be of shape (n, 4) or (n, 5) where the second axis is (value, x, y, z)
+        or (value, x, y, z, photonID). The photonID column is optional."""
+        assert array.shape[1] in [4, 5] and array.ndim == 2, "Data point array must be of shape (n, 4) or (n, 5)"
         self._appendData(array, DataType.DATA_POINT, key)
 
     def logSegmentArray(self, array: np.ndarray, key: InteractionKey = None):
@@ -127,8 +131,8 @@ def getPoints(self, key: InteractionKey = None) -> np.ndarray:
         return self._getData(DataType.POINT, key)
 
     def getRawDataPoints(self, key: InteractionKey = None) -> np.ndarray:
-        """All raw 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) where
-        the second axis is (value, x, y, z)."""
+        """All raw 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) or (n, 5) where
+        the second axis is (value, x, y, z) or (value, x, y, z, photonID) if photon IDs were logged."""
         return self._getData(DataType.DATA_POINT, key)
 
     def getSegments(self, key: InteractionKey = None) -> np.ndarray:
diff --git a/pytissueoptics/scene/scene/scene.py b/pytissueoptics/scene/scene/scene.py
index 7e7b4085..25ad9d1e 100644
--- a/pytissueoptics/scene/scene/scene.py
+++ b/pytissueoptics/scene/scene/scene.py
@@ -142,6 +142,8 @@ def getPolygons(self) -> List[Polygon]:
     def getMaterials(self) -> list:
         materials = [self._worldMaterial]
         for solid in self._solids:
+            if solid.isDetector:
+                continue
             surfaceLabels = solid.surfaceLabels
             for surfaceLabel in surfaceLabels:
                 material = solid.getPolygons(surfaceLabel)[0].insideEnvironment.material
diff --git a/pytissueoptics/scene/solids/__init__.py b/pytissueoptics/scene/solids/__init__.py
index 9e63b939..b4e464f8 100644
--- a/pytissueoptics/scene/solids/__init__.py
+++ b/pytissueoptics/scene/solids/__init__.py
@@ -1,9 +1,11 @@
+from .circle import Circle
 from .cone import Cone
 from .cube import Cube
 from .cuboid import Cuboid
 from .cylinder import Cylinder
 from .ellipsoid import Ellipsoid
 from .lens import PlanoConcaveLens, PlanoConvexLens, SymmetricLens, ThickLens
+from .rectangle import Rectangle
 from .solid import Solid
 from .solidFactory import SolidFactory
 from .sphere import Sphere
@@ -21,4 +23,6 @@
     "PlanoConcaveLens",
     "ThickLens",
     "SymmetricLens",
+    "Rectangle",
+    "Circle",
 ]
diff --git a/pytissueoptics/scene/solids/circle.py b/pytissueoptics/scene/solids/circle.py
new file mode 100644
index 00000000..c35b9b44
--- /dev/null
+++ b/pytissueoptics/scene/solids/circle.py
@@ -0,0 +1,44 @@
+from ..geometry import Vector, primitives
+from .cylinder import Cylinder
+
+
+class Circle(Cylinder):
+    def __init__(
+        self,
+        radius: float,
+        orientation: Vector = Vector(0, 0, 1),
+        position: Vector = Vector(0, 0, 0),
+        u: int = 32,
+        s: int = 1,
+        primitive: str = primitives.DEFAULT,
+        label="circle",
+    ):
+        super().__init__(
+            radius, length=0, u=u, v=1, s=s, position=position, label=label, primitive=primitive, smooth=False
+        )
+        self.orient(towards=orientation)
+
+    def _computeTriangleMesh(self):
+        # A circle is a cylinder of length 0 with a single (back) surface pointing towards the positive z axis before orientation.
+        lateralLayers = self._computeSectionVertices(lateralSteps=[0], radialSteps=[self._s])
+        backLayers = self._computeSectionVertices(lateralSteps=[0], radialSteps=list(range(self._s - 1, 0, -1)))
+        backLayers.insert(0, lateralLayers[-1])
+        backLayers.append([self._backCenter])
+        self._vertices.append(self._backCenter)
+        self._surfaces.add("surface", self._getSurfaceTriangles(backLayers))
+
+    def _geometryParams(self) -> dict:
+        return {
+            "radius": self._radius,
+            "position": self._position,
+            "orientation": self._orientation,
+            "u": self._u,
+            "s": self._s,
+        }
+
+    @property
+    def isFlat(self) -> bool:
+        return True
+
+    def __hash__(self):
+        return hash((self._radius, self._orientation, self._position))
diff --git a/pytissueoptics/scene/solids/cuboid.py b/pytissueoptics/scene/solids/cuboid.py
index 609efd84..29dfa204 100644
--- a/pytissueoptics/scene/solids/cuboid.py
+++ b/pytissueoptics/scene/solids/cuboid.py
@@ -154,3 +154,15 @@ def _getLayerPolygons(self, layerLabel: str) -> List[Polygon]:
 
     def _geometryParams(self) -> dict:
         return {"shape": self.shape}
+
+    def __hash__(self):
+        baseHash = super().__hash__()
+        if not self.isStack():
+            return baseHash
+
+        materials = set()
+        for surfaceLabel in self.surfaceLabels:
+            material = self.getPolygons(surfaceLabel)[0].insideEnvironment.material
+            materials.add(material)
+        materialsHash = hash(tuple(materials))
+        return hash((baseHash, materialsHash))
diff --git a/pytissueoptics/scene/solids/rectangle.py b/pytissueoptics/scene/solids/rectangle.py
new file mode 100644
index 00000000..736f151d
--- /dev/null
+++ b/pytissueoptics/scene/solids/rectangle.py
@@ -0,0 +1,41 @@
+from ..geometry import Quad, Triangle, Vector, Vertex, primitives
+from .solid import Solid
+
+
+class Rectangle(Solid):
+    def __init__(
+        self,
+        a: float,
+        b: float,
+        orientation=Vector(0, 0, 1),
+        position=Vector(0, 0, 0),
+        primitive: str = primitives.DEFAULT,
+        label="rectangle",
+    ):
+        self._a = a
+        self._b = b
+
+        vertices = [  # Mesh points towards (0, 0, 1) by default, as expected by the orient() method.
+            Vertex(-a / 2, -b / 2, 0),
+            Vertex(a / 2, -b / 2, 0),
+            Vertex(a / 2, b / 2, 0),
+            Vertex(-a / 2, b / 2, 0),
+        ]
+
+        super().__init__(vertices, position, label=label, primitive=primitive)
+        self.orient(towards=orientation)
+
+    def _computeTriangleMesh(self):
+        V = self._vertices
+        self._surfaces.add("surface", [Triangle(V[0], V[1], V[2]), Triangle(V[2], V[3], V[0])])
+
+    def _computeQuadMesh(self):
+        V = self._vertices
+        self._surfaces.add("surface", [Quad(V[0], V[1], V[2], V[3])])
+
+    def _geometryParams(self) -> dict:
+        return {"a": self._a, "b": self._b, "position": self._position, "orientation": self._orientation}
+
+    @property
+    def isFlat(self) -> bool:
+        return True
diff --git a/pytissueoptics/scene/solids/solid.py b/pytissueoptics/scene/solids/solid.py
index 007b9fd0..53874bdc 100644
--- a/pytissueoptics/scene/solids/solid.py
+++ b/pytissueoptics/scene/solids/solid.py
@@ -42,6 +42,7 @@ def __init__(
         self._bbox = None
         self._label = label
         self._layerLabels = {}
+        self._detectorAcceptanceCosine = None
 
         if not self._surfaces:
             self._computeMesh()
@@ -80,6 +81,36 @@ def primitive(self) -> str:
     def bbox(self) -> BoundingBox:
         return self._bbox
 
+    def asDetector(self, halfAngle: float = np.pi / 2) -> "Solid":
+        """Treat this solid as a detector with a given half angle in radians.
+
+        Detectors will fully absorb a photon when incident within the half angle, else the photon will go through
+        unaffected. Enabling this makes any previous material assigned to the solid irrelevant. Note that using a half
+        angle above pi/2 will detect rays coming from the back as well.
+        """
+        if halfAngle < 0 or halfAngle > np.pi:
+            raise ValueError("Detector half angle must be between 0 and pi radians.")
+
+        if self._material:
+            warnings.warn(
+                f"Solid '{self._label}' is now treated as a detector. Its material '{self._material}' will be ignored.",
+                RuntimeWarning,
+            )
+            self._material = None
+
+        self._detectorAcceptanceCosine = np.cos(halfAngle)
+        return self
+
+    @property
+    def isDetector(self) -> bool:
+        return self._detectorAcceptanceCosine is not None
+
+    @property
+    def detectorAcceptanceCosine(self) -> float:
+        if self._detectorAcceptanceCosine is None:
+            raise RuntimeError("This solid is not a detector. Set the half angle first with asDetector().")
+        return self._detectorAcceptanceCosine
+
     def getBoundingBox(self) -> BoundingBox:
         return self.bbox
 
@@ -348,4 +379,8 @@ def geometryExport(self) -> dict[str, Any]:
 
     def _geometryParams(self) -> dict:
         """To be implemented by Solid subclasses to detail other geometry parameters."""
-        return {}
+        raise NotImplementedError(f"Geometry parameters not implemented for Solids of type {type(self).__name__}")
+
+    @property
+    def isFlat(self) -> bool:
+        return False
diff --git a/pytissueoptics/scene/tests/geometry/testUtils.py b/pytissueoptics/scene/tests/geometry/testUtils.py
index ffa2c41b..f88d06cf 100644
--- a/pytissueoptics/scene/tests/geometry/testUtils.py
+++ b/pytissueoptics/scene/tests/geometry/testUtils.py
@@ -94,13 +94,12 @@ def testGivenTwoVectorsOnABaseAxisPlaneWithLessThan180DegreesBetween_shouldRetur
         testCases = self._generateAxisAlignedTestCases()
         for testCase in testCases:
             axis, angle = getAxisAngleBetween(testCase.v1, testCase.v2)
-            self.assertTrue(np.allclose(testCase.expectedAxis.array, axis.array))
             self.assertTrue(np.isclose(testCase.expectedAngle, angle))
+            self.assertTrue(np.isclose(1, axis.getNorm()))
 
-            axis, angle = getAxisAngleBetween(testCase.v2, testCase.v1)
-            testCase.expectedAxis *= -1
-            self.assertTrue(np.allclose(testCase.expectedAxis.array, axis.array))
-            self.assertTrue(np.isclose(testCase.expectedAngle, angle))
+            # For non-degenerate cases (angle not 0 or 180), check axis direction
+            if not np.isclose(testCase.expectedAngle % np.pi, 0):
+                self.assertTrue(np.allclose(testCase.expectedAxis.array, axis.array))
 
     def _generateAxisAlignedTestCases(self):
         testCases = []
@@ -126,3 +125,10 @@ def testGivenTwoVectorsWithSameOrientation_shouldReturnZeroAngle(self):
         v2 = Vector(1, 0, 0)
         _, angle = getAxisAngleBetween(v1, v2)
         self.assertTrue(np.isclose(0, angle))
+
+    def testGivenTwoOppositeVectors_shouldReturnAnyAxisAnd180Degrees(self):
+        v1 = Vector(1, 0, 0)
+        v2 = Vector(-1, 0, 0)
+        axis, angle = getAxisAngleBetween(v1, v2)
+        self.assertTrue(np.isclose(np.pi, angle))
+        self.assertTrue(np.isclose(1, axis.getNorm()))
diff --git a/pytissueoptics/scene/tests/intersection/testIntersectionFinder.py b/pytissueoptics/scene/tests/intersection/testIntersectionFinder.py
index 50e86bfe..c502a1be 100644
--- a/pytissueoptics/scene/tests/intersection/testIntersectionFinder.py
+++ b/pytissueoptics/scene/tests/intersection/testIntersectionFinder.py
@@ -227,6 +227,13 @@ def testGivenSmoothSolid_shouldNotSmoothTheIntersectionNormalIfItChangesTheSignO
         self.assertIsNotNone(intersection)
         self.assertNotEqual(smoothNormal, intersection.normal)
 
+    def testGivenIgnoreLabel_shouldNotIntersectWithIt(self):
+        ray = Ray(origin=Vector(0, 0.5, 0), direction=Vector(0, 0, 1))
+        solid = Cube(2, position=Vector(0, 0, 5), label="ignoreMe")
+
+        intersection = self.getIntersectionFinder([solid]).findIntersection(ray, WORLD_LABEL, ignoreLabel="ignoreMe")
+        self.assertIsNone(intersection)
+
     def assertVectorEqual(self, expected, actual):
         self.assertEqual(expected.x, actual.x)
         self.assertEqual(expected.y, actual.y)
diff --git a/pytissueoptics/scene/tests/solids/testCircle.py b/pytissueoptics/scene/tests/solids/testCircle.py
new file mode 100644
index 00000000..c3629e65
--- /dev/null
+++ b/pytissueoptics/scene/tests/solids/testCircle.py
@@ -0,0 +1,37 @@
+import unittest
+
+from pytissueoptics.scene.geometry import Vector
+from pytissueoptics.scene.solids import Circle
+
+
+class TestCircle(unittest.TestCase):
+    def testShouldBeFlat(self):
+        circle = Circle(radius=1)
+        self.assertTrue(circle.isFlat)
+
+    def testShouldHaveDefaultZOrientation(self):
+        circle = Circle(radius=1)
+        polygons = circle.getPolygons()
+        normals = [p.normal for p in polygons]
+        for n in normals:
+            self.assertEqual(n, Vector(0, 0, 1))
+
+    def testGivenAnOrientation_shouldOrientTheCircle(self):
+        orientation = Vector(1, 1, 0)
+        orientation.normalize()
+        circle = Circle(radius=1, orientation=orientation)
+        polygons = circle.getPolygons()
+        normals = [p.normal for p in polygons]
+        for n in normals:
+            self.assertEqual(n, orientation)
+
+    def testShouldHaveASingleSurface(self):
+        circle = Circle(radius=1)
+        assert len(circle.surfaceLabels) == 1
+
+    def testShouldHaveDifferentHashWhenOrientationChanges(self):
+        circle1 = Circle(radius=1, orientation=Vector(0, 0, 1))
+        circle2 = Circle(radius=1, orientation=Vector(1, 1, 0))
+        self.assertNotEqual(hash(circle1), hash(circle2))
+        circle3 = Circle(radius=1, orientation=Vector(0, 0, 1))
+        self.assertEqual(hash(circle1), hash(circle3))
diff --git a/pytissueoptics/scene/tests/solids/testRectangle.py b/pytissueoptics/scene/tests/solids/testRectangle.py
new file mode 100644
index 00000000..55a6ca11
--- /dev/null
+++ b/pytissueoptics/scene/tests/solids/testRectangle.py
@@ -0,0 +1,37 @@
+import unittest
+
+from pytissueoptics.scene.geometry import Vector
+from pytissueoptics.scene.solids import Rectangle
+
+
+class TestRectangle(unittest.TestCase):
+    def testShouldBeFlat(self):
+        rectangle = Rectangle(2, 2)
+        self.assertTrue(rectangle.isFlat)
+
+    def testShouldHaveDefaultZOrientation(self):
+        rectangle = Rectangle(2, 2)
+        polygons = rectangle.getPolygons()
+        normals = [p.normal for p in polygons]
+        for n in normals:
+            self.assertEqual(n, Vector(0, 0, 1))
+
+    def testGivenAnOrientation_shouldOrientTheRectangle(self):
+        orientation = Vector(1, 1, 0)
+        orientation.normalize()
+        rectangle = Rectangle(2, 2, orientation=orientation)
+        polygons = rectangle.getPolygons()
+        normals = [p.normal for p in polygons]
+        for n in normals:
+            self.assertEqual(n, orientation)
+
+    def testShouldHaveASingleSurface(self):
+        rectangle = Rectangle(2, 2)
+        assert len(rectangle.surfaceLabels) == 1
+
+    def testShouldHaveDifferentHashWhenOrientationChanges(self):
+        rectangle1 = Rectangle(2, 2, orientation=Vector(0, 0, 1))
+        rectangle2 = Rectangle(2, 2, orientation=Vector(1, 1, 0))
+        self.assertNotEqual(hash(rectangle1), hash(rectangle2))
+        rectangle3 = Rectangle(2, 2, orientation=Vector(0, 0, 1))
+        self.assertEqual(hash(rectangle1), hash(rectangle3))
diff --git a/pytissueoptics/scene/tests/solids/testSolid.py b/pytissueoptics/scene/tests/solids/testSolid.py
index 3dac4eb8..e5329e99 100644
--- a/pytissueoptics/scene/tests/solids/testSolid.py
+++ b/pytissueoptics/scene/tests/solids/testSolid.py
@@ -257,6 +257,24 @@ def testWhenCheckIfContainsAVertexPartiallyInside_shouldWarnAndReturnFalse(self)
         with self.assertWarns(RuntimeWarning):
             self.assertFalse(self.solid.contains(Vertex(2, 2, -0.75)))
 
+    def testShouldNotBeFlat(self):
+        self.assertFalse(self.solid.isFlat)
+
+    def shouldHaveNoDetectorProperties(self):
+        with self.assertRaises(RuntimeError):
+            _ = self.solid.detectorAcceptanceCosine
+
+    def testWhenConvertToDetector_shouldSetDetectorCosine(self):
+        halfAngle = math.pi / 8
+        self.solid.asDetector(halfAngle)
+        expectedAcceptanceCosine = math.cos(halfAngle)
+        self.assertEqual(expectedAcceptanceCosine, self.solid.detectorAcceptanceCosine)
+
+    def testWhenConvertToDetector_shouldRemoveMaterial(self):
+        self.assertIsNotNone(self.solid.getEnvironment().material)
+        self.solid.asDetector()
+        self.assertIsNone(self.solid.getEnvironment().material)
+
     @staticmethod
     def createPolygonMock() -> Polygon:
         polygon = mock(Polygon)