diff --git a/optika/sensors/materials/_materials.py b/optika/sensors/materials/_materials.py
index 3860633..716ac7f 100644
--- a/optika/sensors/materials/_materials.py
+++ b/optika/sensors/materials/_materials.py
@@ -13,6 +13,7 @@
     _thickness_oxide,
     _thickness_implant,
     _thickness_substrate,
+    _width_pixel,
     _cce_backsurface,
 )
 from ._ramanathan_2020 import (
@@ -945,18 +946,22 @@ def electrons_measured_approx(
     return result
 
 
-_absorbance = absorbance
+_transmittance = transmittance
 
 
 def signal(
     photons_expected: u.Quantity | na.AbstractScalar,
     wavelength: u.Quantity | na.ScalarArray,
-    absorbance: None | float | na.AbstractScalar = None,
+    transmittance: None | float | na.AbstractScalar = None,
     absorption: None | u.Quantity | na.AbstractScalar = None,
     thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant,
+    thickness_depletion: u.Quantity | na.AbstractScalar = _thickness_substrate,
+    thickness_substrate: None | na.AbstractScalar = _thickness_substrate,
+    width_pixel: u.Quantity | na.AbstractScalar = _width_pixel,
     cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface,
     temperature: u.Quantity | na.ScalarArray = 300 * u.K,
-    method: Literal["exact", "approx", "expected"] = "exact",
+    method: Literal["monte-carlo", "expected"] = "monte-carlo",
+    axis_xy: None | tuple[str, str] = None,
     shape_random: None | dict[str, int] = None,
 ) -> na.AbstractScalar:
     r"""
@@ -974,9 +979,9 @@ def signal(
         The `expected` number of photons incident on the detector surface.
     wavelength
         The vacuum wavelength of the absorbed photons.
-    absorbance
-        The fraction of incident energy absorbed by the light-sensitive layer
-        of the detector computed using the average of :func:`absorbance`.
+    transmittance
+        The fraction of incident energy transmitted to the light-sensitive layer
+        of the detector.
         If :obj:`None` (the default), the result of :func:`absorbance`
         called with default values will be used.
     absorption
@@ -987,6 +992,16 @@ def signal(
     thickness_implant
         The thickness of the implant layer.
         Default is the value given in :cite:t:`Stern1994`.
+    thickness_depletion
+        The thickness of the depletion region, the region with significant electric
+        field.
+        If :obj:`None` (the default), this is set to the same value as
+        `thickness_substrate`.
+    thickness_substrate
+        The thickness of the entire light-sensitive region of the device.
+        The default
+    width_pixel
+        The size of a single pixel on the sensor.
     cce_backsurface
         The differential charge collection efficiency on the back surface
         of the sensor.
@@ -1001,6 +1016,10 @@ def signal(
         of photons is high.
         The `expected` method does not add any noise to the signal and just
         returns the expected number of electrons.
+    axis_xy
+        The two logical axes corresponding to the pixel grid of the sensor
+        along which electrons will diffuse.
+        If :obj:`None` (the default), there is no charge diffusion.
     shape_random
         Additional shape used to specify the number of samples to draw.
 
@@ -1048,8 +1067,8 @@ def signal(
         ax.set_ylabel(f"variance-to-mean ratio ({electrons.unit:latex_inline})");
     """
 
-    if absorbance is None:
-        absorbance = _absorbance(wavelength).average
+    if transmittance is None:
+        transmittance = _transmittance(wavelength).average
 
     if absorption is None:
         absorption = optika.chemicals.Chemical("Si").absorption(wavelength)
@@ -1059,6 +1078,7 @@ def signal(
             wavelength=wavelength,
             temperature=temperature,
         )
+        absorbance = transmittance * np.exp(-absorption * thickness_substrate)
         cce = charge_collection_efficiency(
             absorption=absorption,
             thickness_implant=thickness_implant,
@@ -1066,26 +1086,28 @@ def signal(
         )
         return iqy * absorbance * cce * photons_expected.to(u.ph)
 
-    photons_absorbed_expected = absorbance * photons_expected.to(u.ph)
-    photons_absorbed = na.random.poisson(
-        lam=photons_absorbed_expected,
-        shape_random=shape_random,
-    ).astype(int)
+    elif method == "monte-carlo":
 
-    kwargs = dict(
-        photons_absorbed=photons_absorbed,
-        wavelength=wavelength,
-        absorption=absorption,
-        thickness_implant=thickness_implant,
-        cce_backsurface=cce_backsurface,
-        temperature=temperature,
-        shape_random=shape_random,
-    )
+        photons_expected = transmittance * photons_expected.to(u.ph)
+        photons = na.random.poisson(
+            lam=photons_expected,
+            shape_random=shape_random,
+        ).astype(int)
+
+        return electrons_measured(
+            photons_transmitted=photons,
+            wavelength=wavelength,
+            absorption=absorption,
+            thickness_implant=thickness_implant,
+            thickness_depletion=thickness_depletion,
+            thickness_substrate=thickness_substrate,
+            width_pixel=width_pixel,
+            cce_backsurface=cce_backsurface,
+            temperature=temperature,
+            axis_xy=axis_xy,
+            shape_random=shape_random,
+        )
 
-    if method == "exact":
-        return electrons_measured(**kwargs)
-    elif method == "approx":
-        return electrons_measured_approx(**kwargs)
     else:  # pragma: nocover
         raise ValueError(f"Unrecognized method: {method}")
 
@@ -1258,22 +1280,31 @@ class AbstractSensorMaterial(
     @abc.abstractmethod
     def signal(
         self,
-        rays: optika.rays.RayVectorArray,
+        photons: u.Quantity | na.AbstractScalar,
+        wavelength: u.Quantity | na.AbstractScalar,
+        direction: na.AbstractCartesian3dVectorArray,
         normal: na.AbstractCartesian3dVectorArray,
+        axis_xy: None | tuple[str, str] = None,
         noise: bool = True,
-    ) -> optika.rays.RayVectorArray:
+    ) -> na.AbstractScalar:
         """
         Given a set of absorbed rays, compute the number of electrons
         measured by the sensor using :func:`signal`.
 
         Parameters
         ----------
-        rays
-            The rays absorbed by the light-sensitive silicon layer.
-            The :attr:`optika.rays.RayVectorArray.intensity` field should
-            either be in units of photons or energy.
+        photons
+            The number of photons incident on each pixel.
+        wavelength
+            An assumed grid of wavelengths for the incident photons.
+        direction
+            An assumed propagation direction for the incident photons.
         normal
             The vector perpendicular to the surface of the sensor.
+        axis_xy
+            The two logical axes corresponding to the pixel grid of the sensor.
+            If provided, charge diffusion will occur along these two axes.
+            If :obj:`None` (the default), no diffusion is performed.
         noise
             Whether to add noise to the result.
         """
@@ -1294,7 +1325,7 @@ def photons_incident(
         Parameters
         ----------
         electrons
-            The number of electrons measured by the sensor.
+            The number of electrons measured by each pixel.
         wavelength
             An assumed grid of wavelengths for the incident photons.
         direction
@@ -1303,24 +1334,6 @@ def photons_incident(
             The vector perpendicular to the surface of the sensor.
         """
 
-    @abc.abstractmethod
-    def charge_diffusion(
-        self,
-        rays: optika.rays.RayVectorArray,
-        normal: na.AbstractCartesian3dVectorArray,
-    ) -> optika.rays.RayVectorArray:
-        """
-        Given a set of incident rays, compute how much the position of each ray
-        is perturbed due to charge diffusion within the sensor.
-
-        Parameters
-        ----------
-        rays
-            The rays incident on the sensor surface.
-        normal
-            The vector perpendicular to the surface of the sensor.
-        """
-
 
 @dataclasses.dataclass(eq=False, repr=False)
 class IdealSensorMaterial(
@@ -1334,15 +1347,18 @@ class IdealSensorMaterial(
 
     def signal(
         self,
-        rays: optika.rays.RayVectorArray,
+        photons: u.Quantity | na.AbstractScalar,
+        wavelength: u.Quantity | na.AbstractScalar,
+        direction: na.AbstractCartesian3dVectorArray,
         normal: na.AbstractCartesian3dVectorArray,
-        noise: bool = False,
+        axis_xy: None | tuple[str, str] = None,
+        noise: bool = True,
     ) -> optika.rays.RayVectorArray:
-        intensity = rays.intensity
-        if not intensity.unit.is_equivalent(u.photon):
+
+        if not photons.unit.is_equivalent(u.photon):
             h = astropy.constants.h
             c = astropy.constants.c
-            intensity = intensity / (h * c / rays.wavelength) * u.photon
+            intensity = photons / (h * c / wavelength) * u.photon
 
         if noise:
             intensity = na.random.poisson(intensity.to(u.ph)).astype(int)
@@ -1350,9 +1366,7 @@ def signal(
         electrons = intensity * u.electron / u.photon
         electrons = electrons.to(u.electron)
 
-        result = dataclasses.replace(rays, intensity=electrons)
-
-        return result
+        return electrons
 
     def photons_incident(
         self,
@@ -1363,13 +1377,6 @@ def photons_incident(
     ) -> na.AbstractScalar:
         return electrons * u.photon / u.electron
 
-    def charge_diffusion(
-        self,
-        rays: optika.rays.RayVectorArray,
-        normal: na.AbstractCartesian3dVectorArray,
-    ) -> optika.rays.RayVectorArray:
-        return rays
-
 
 @dataclasses.dataclass(eq=False, repr=False)
 class AbstractSiliconSensorMaterial(
@@ -1499,8 +1506,7 @@ def depletion(self) -> AbstractDepletionModel:
 
     def width_charge_diffusion(
         self,
-        rays: optika.rays.RayVectorArray,
-        normal: na.AbstractCartesian3dVectorArray,
+        wavelength: u.Quantity | na.AbstractScalar,
     ) -> na.AbstractScalar:
         """
         The standard deviation of the charge diffusion kernel for this sensor.
@@ -1508,21 +1514,62 @@ def width_charge_diffusion(
 
         Parameters
         ----------
-        rays
-            The rays incident on the sensor surface.
-        normal
-            The vector perpendicular to the surface of the sensor.
+        wavelength
+            The wavelength of the incident light in vacuum.
         """
         return optika.sensors.charge_diffusion(
-            self._chemical.absorption(rays.wavelength),
+            self._chemical.absorption(wavelength),
             thickness_substrate=self.thickness_substrate,
             thickness_depletion=self.depletion.thickness,
         )
 
+    def transmittance(
+        self,
+        wavelength: u.Quantity | na.AbstractScalar,
+        direction: None | na.AbstractCartesian3dVectorArray = None,
+        index_refraction: float | na.AbstractScalar = 1,
+        normal: None | na.AbstractCartesian3dVectorArray = None,
+    ) -> optika.vectors.PolarizationVectorArray:
+        """
+        Compute the fraction of energy transmitted to the light-sensitive region
+        of the sensor.
+
+        Parameters
+        ----------
+        wavelength
+            The wavelength of the incident light in vacuum.
+        direction
+            The propagation direction of the incident light in the ambient medium.
+            If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed.
+        index_refraction
+            The complex index of refraction of the ambient medium.
+        normal
+            The vector perpendicular to the surface of the CCD sensor.
+        """
+        if direction is None:
+            direction = na.Cartesian3dVectorArray(0, 0, 1)
+
+        if normal is None:
+            normal = na.Cartesian3dVectorArray(0, 0, 1)
+
+        return transmittance(
+            wavelength=wavelength,
+            direction=-direction @ normal,
+            n=index_refraction,
+            thickness_oxide=self.thickness_oxide,
+            thickness_substrate=self.thickness_substrate,
+            chemical_oxide=self._chemical_oxide,
+            chemical_substrate=self._chemical,
+            roughness_oxide=self.roughness_oxide,
+            roughness_substrate=self.roughness_substrate,
+        )
+
     def absorbance(
         self,
-        rays: optika.rays.AbstractRayVectorArray,
-        normal: na.AbstractCartesian3dVectorArray,
+        wavelength: u.Quantity | na.AbstractScalar,
+        direction: None | na.AbstractCartesian3dVectorArray = None,
+        index_refraction: float | na.AbstractScalar = 1,
+        normal: None | na.AbstractCartesian3dVectorArray = None,
     ) -> optika.vectors.PolarizationVectorArray:
         """
         Compute the fraction of energy absorbed by the light-sensitive region
@@ -1530,15 +1577,26 @@ def absorbance(
 
         Parameters
         ----------
-        rays
-            The light rays incident on the CCD surface.
+        wavelength
+            The wavelength of the incident light in vacuum.
+        direction
+            The propagation direction of the incident light in the ambient medium.
+            If :obj:`None` (default), normal incidence (:math:`\hat{z}`) is assumed.
+        n
+            The complex index of refraction of the ambient medium.
         normal
             The vector perpendicular to the surface of the CCD sensor.
         """
+        if direction is None:
+            direction = na.Cartesian3dVectorArray(0, 0, 1)
+        
+        if normal is None:
+            normal = na.Cartesian3dVectorArray(0, 0, 1)
+        
         return absorbance(
-            wavelength=rays.wavelength,
-            direction=-rays.direction @ normal,
-            n=rays.index_refraction,
+            wavelength=wavelength,
+            direction=-direction @ normal,
+            n=index_refraction,
             thickness_oxide=self.thickness_oxide,
             thickness_substrate=self.thickness_substrate,
             chemical_oxide=self._chemical_oxide,
diff --git a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py
index af58fa3..795cda9 100644
--- a/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py
+++ b/optika/sensors/materials/_ramanathan_2020/_ramanathan_2020.py
@@ -8,6 +8,8 @@
 import optika
 from .._stern_1994 import (
     _thickness_implant,
+    _thickness_substrate,
+    _width_pixel,
     _cce_backsurface,
 )
 
@@ -464,12 +466,16 @@ def probability_of_n_pairs(
 
 
 def electrons_measured(
-    photons_absorbed: u.Quantity | na.AbstractScalar,
+    photons_transmitted: u.Quantity | na.AbstractScalar,
     wavelength: u.Quantity | na.ScalarArray,
     absorption: None | u.Quantity | na.AbstractScalar = None,
     thickness_implant: u.Quantity | na.AbstractScalar = _thickness_implant,
+    thickness_depletion: None | u.Quantity | na.AbstractScalar = None,
+    thickness_substrate: u.Quantity | na.AbstractScalar = _thickness_substrate,
+    width_pixel: u.Quantity | na.ScalarArray = _width_pixel,
     cce_backsurface: u.Quantity | na.AbstractScalar = _cce_backsurface,
     temperature: u.Quantity | na.ScalarArray = 300 * u.K,
+    axis_xy: None | tuple[str, str] = None,
     shape_random: None | dict[str, int] = None,
 ) -> na.AbstractScalar:
     r"""
@@ -482,20 +488,34 @@ def electrons_measured(
 
     Parameters
     ----------
-    photons_absorbed
-        The number of photons absorbed by the light-sensitive layer of the sensor.
+    photons_transmitted
+        The number of photons transmitted into the light-sensitive region
+        of the sensor.
     wavelength
         The vacuum wavelength of the absorbed photons.
     absorption
         The absorption coefficient of silicon at the given wavelength.
     thickness_implant
         The thickness of the implant layer, where partial-charge collection occurs.
+    thickness_depletion
+        The thickness of the depletion region, the region with significant electric
+        field.
+        If :obj:`None` (the default), this is set to the same value as
+        `thickness_substrate`.
+    thickness_substrate
+        The thickness of the entire light-sensitive region of the device.
+    width_pixel
+        The size of a single pixel on the sensor.
     cce_backsurface
         The differential charge collection efficiency on the back surface
         of the sensor.
     temperature
         The temperature of the silicon detector.
         Default is room temperature.
+    axis_xy
+        The two logical axes corresponding to the pixel grid of the sensor
+        along which electrons will diffuse.
+        If :obj:`None` (the default), there is no charge diffusion.
     shape_random
         Additional shape used to specify the number of samples to draw.
 
@@ -518,7 +538,7 @@ def electrons_measured(
 
         # Define the expected number of photons
         # for each experiment
-        photons_absorbed = (100 * u.photon).astype(int)
+        photons_transmitted = (100 * u.photon).astype(int)
 
         # Define the wavelength at which to sample the distribution
         wavelength = 5.9 * u.keV
@@ -526,7 +546,7 @@ def electrons_measured(
 
         # Compute the actual number of electrons measured for each experiment
         electrons = optika.sensors.electrons_measured(
-            photons_absorbed=photons_absorbed,
+            photons_transmitted=photons_transmitted,
             wavelength=wavelength,
             shape_random=dict(experiment=num_experiments),
         )
@@ -561,22 +581,41 @@ def electrons_measured(
     if absorption is None:
         absorption = optika.chemicals.Chemical("Si").absorption(wavelength)
 
+    if thickness_depletion is None:
+        thickness_depletion = thickness_substrate
+
     if shape_random is None:
         shape_random = dict()
 
     shape = na.broadcast_shapes(
-        na.shape(photons_absorbed),
+        na.shape(photons_transmitted),
         na.shape(wavelength),
         na.shape(absorption),
         na.shape(thickness_implant),
+        na.shape(thickness_depletion),
+        na.shape(thickness_substrate),
+        na.shape(width_pixel),
         na.shape(cce_backsurface),
         na.shape(temperature),
         shape_random,
     )
 
-    photons_absorbed = na.broadcast_to(photons_absorbed, shape)
+    if axis_xy is not None:
+        axis_x, axis_y = axis_xy
+        shape[axis_x] = shape.pop(axis_x)
+        shape[axis_y] = shape.pop(axis_y)
+    else:
+        axis_x = "__dummy_x__"
+        axis_y = "__dummy_y__"
+        shape[axis_x] = 1
+        shape[axis_y] = 1
+
+    photons_transmitted = na.broadcast_to(photons_transmitted, shape)
     absorption = na.broadcast_to(absorption, shape)
     thickness_implant = na.broadcast_to(thickness_implant, shape)
+    thickness_depletion = na.broadcast_to(thickness_depletion, shape)
+    thickness_substrate = na.broadcast_to(thickness_substrate, shape)
+    width_pixel = na.broadcast_to(width_pixel, shape)
     cce_backsurface = na.broadcast_to(cce_backsurface, shape)
 
     if not isinstance(cce_backsurface, u.Quantity):
@@ -600,10 +639,13 @@ def electrons_measured(
     fano_inf = fano_inf.broadcast_to(shape)
 
     result = _electrons_measured_quantity(
-        photons_absorbed=photons_absorbed.ndarray,
+        photons_transmitted=photons_transmitted.ndarray,
         wavelength=wavelength.ndarray,
         absorption=absorption.ndarray,
         thickness_implant=thickness_implant.ndarray,
+        thickness_depletion=thickness_depletion.ndarray,
+        thickness_substrate=thickness_substrate.ndarray,
+        width_pixel=width_pixel,
         cce_backsurface=cce_backsurface.ndarray,
         p_n=p_n.ndarray,
         n=n.ndarray,
@@ -618,10 +660,13 @@ def electrons_measured(
 
 
 def _electrons_measured_quantity(
-    photons_absorbed: u.Quantity,
+    photons_transmitted: u.Quantity,
     wavelength: u.Quantity,
     absorption: u.Quantity,
     thickness_implant: u.Quantity,
+    thickness_depletion: u.Quantity,
+    thickness_substrate: u.Quantity,
+    width_pixel: u.Quantity,
     cce_backsurface: u.Quantity,
     p_n: np.ndarray,
     n: np.ndarray,
@@ -630,32 +675,44 @@ def _electrons_measured_quantity(
 ) -> u.Quantity:
 
     shape = np.broadcast_shapes(
-        photons_absorbed.shape,
+        photons_transmitted.shape,
         absorption.shape,
         thickness_implant.shape,
+        thickness_depletion.shape,
+        thickness_substrate.shape,
         cce_backsurface.shape,
+        width_pixel.shape,
     )
 
+    num_x = shape[~1]
+    num_y = shape[~0]
+
     unit_length = u.mm
 
-    photons_absorbed = photons_absorbed.to_value(u.photon)
+    photons_transmitted = photons_transmitted.to_value(u.photon)
     energy = wavelength.to(u.eV, equivalencies=u.spectral())
     absorption = absorption.to_value(1 / unit_length)
     thickness_implant = thickness_implant.to_value(unit_length)
+    thickness_depletion = thickness_depletion.to_value(unit_length)
+    thickness_substrate = thickness_substrate.to_value(unit_length)
+    width_pixel = width_pixel.to_value(unit_length)
     cce_backsurface = cce_backsurface.to_value(u.dimensionless_unscaled)
     energy_pair_inf = energy_pair_inf.to_value(u.eV)
     fano_inf = fano_inf.to_value(u.electron / u.photon)
 
     result = _electrons_measured_numba(
-        photons_absorbed=photons_absorbed.reshape(-1),
-        energy=energy.reshape(-1),
-        absorption=absorption.reshape(-1),
-        thickness_implant=thickness_implant.reshape(-1),
-        cce_backsurface=cce_backsurface.reshape(-1),
-        p_n=p_n.reshape(-1, p_n.shape[~0]),
-        n=n.reshape(-1, n.shape[~0]),
-        energy_pair_inf=energy_pair_inf.reshape(-1),
-        fano_inf=fano_inf.reshape(-1),
+        photons_transmitted=photons_transmitted.reshape(-1, num_x, num_y),
+        energy=energy.reshape(-1, num_x, num_y),
+        absorption=absorption.reshape(-1, num_x, num_y),
+        thickness_implant=thickness_implant.reshape(-1, num_x, num_y),
+        thickness_depletion=thickness_depletion.reshape(-1, num_x, num_y),
+        thickness_substrate=thickness_substrate.reshape(-1, num_x, num_y),
+        width_pixel=width_pixel.reshape(-1, num_x, num_y),
+        cce_backsurface=cce_backsurface.reshape(-1, num_x, num_y),
+        p_n=p_n.reshape(-1, num_x, num_y, p_n.shape[~0]),
+        n=n.reshape(-1, num_x, num_y, n.shape[~0]),
+        energy_pair_inf=energy_pair_inf.reshape(-1, num_x, num_y),
+        fano_inf=fano_inf.reshape(-1, num_x, num_y),
     )
 
     result = result.reshape(shape)
@@ -667,73 +724,92 @@ def _electrons_measured_quantity(
 
 @numba.njit(fastmath=True, parallel=True)
 def _electrons_measured_numba(
-    photons_absorbed: np.ndarray,
+    photons_transmitted: np.ndarray,
     energy: np.ndarray,
     absorption: np.ndarray,
     thickness_implant: np.ndarray,
+    thickness_depletion: np.ndarray,
+    thickness_substrate: np.ndarray,
+    width_pixel: np.ndarray,
     cce_backsurface: np.ndarray,
     p_n: np.ndarray,
     n: np.ndarray,
     energy_pair_inf: np.ndarray,
     fano_inf: np.ndarray,
-) -> np.ndarray:  # pragma: nocover
+) -> np.ndarray:
+    
+    num_i, num_x, num_y, num_n = p_n.shape
 
-    num_i, num_n = p_n.shape
-
-    result = np.empty(num_i)
+    result = np.zeros((num_i, num_x, num_y))
 
     for i in numba.prange(num_i):
+        for x in range(num_x):
+            for y in range(num_y):
+                num_photon = int(photons_transmitted[i, x, y])
+                energy_i = energy[i, x, y]
+                a = absorption[i, x, y]
+                W = thickness_implant[i, x, y]
+                h_0 = cce_backsurface[i, x, y]
+                cmf_i = np.cumsum(p_n[i, x, y])
+                n_i = n[i, x, y]
+                energy_pair_inf_i = energy_pair_inf[i, x, y]
+                fano_inf_i = fano_inf[i, x, y]
+                z_substrate = thickness_substrate[i, x, y]
+                z_ff = z_substrate - thickness_depletion[i, x, y]
+
+                d = 1 / a
 
-        num_photon = int(photons_absorbed[i])
-        energy_i = energy[i]
-        a = absorption[i]
-        W = thickness_implant[i]
-        h_0 = cce_backsurface[i]
-        cmf_i = np.cumsum(p_n[i])
-        n_i = n[i]
-        energy_pair_inf_i = energy_pair_inf[i]
-        fano_inf_i = fano_inf[i]
+                mean_inf = energy_i / energy_pair_inf_i
+                std_inf = math.sqrt(fano_inf_i * mean_inf)
 
-        d = 1 / a
+                low_energy = energy_i <= 50
 
-        num_electron_i = 0
+                for j in range(num_photon):
+                    if low_energy:
+                        x_ij = random.uniform(0, 1)
 
-        mean_inf = energy_i / energy_pair_inf_i
-        std_inf = math.sqrt(fano_inf_i * mean_inf)
+                        for k_ij in range(num_n):
+                            if cmf_i[k_ij] > x_ij:
+                                break
 
-        low_energy = energy_i <= 50
+                        n_ij = n_i[k_ij]
 
-        for j in range(num_photon):
+                    else:
+                        n_ij = random.normalvariate(
+                            mu=mean_inf,
+                            sigma=std_inf,
+                        )
+                        n_ij = round(n_ij)
 
-            if low_energy:
+                    y_ij = random.uniform(0, 1)
+                    z_ij = -d * math.log(1 - y_ij)
 
-                x_ij = random.uniform(0, 1)
+                    if z_ij < W:
+                        h_ij = h_0 + (1 - h_0) * z_ij / W
+                    else:
+                        h_ij = 1
 
-                for k_ij in range(num_n):
-                    if cmf_i[k_ij] > x_ij:
-                        break
+                    m_ij = np.random.binomial(n=n_ij, p=h_ij)
 
-                n_ij = n_i[k_ij]
+                    u = random.uniform(-0.5, 0.5)
+                    v = random.uniform(-0.5, 0.5)
 
-            else:
-                n_ij = random.normalvariate(
-                    mu=mean_inf,
-                    sigma=std_inf,
-                )
-                n_ij = round(n_ij)
+                    for e in range(m_ij):
+                        if z_ij < z_ff:
+                            w = z_ff * math.sqrt(1 - z_ij / z_ff) / width_pixel
 
-            y_ij = random.uniform(0, 1)
-            z_ij = -d * math.log(1 - y_ij)
+                            p = random.gauss(u, w)
+                            q = random.gauss(v, w)
 
-            if z_ij < W:
-                h_ij = h_0 + (1 - h_0) * z_ij / W
-            else:
-                h_ij = 1
+                            p = round(p)
+                            q = round(q)
 
-            m_ij = np.random.binomial(n=n_ij, p=h_ij)
+                        elif z_ij > z_substrate:
+                            continue
 
-            num_electron_i = num_electron_i + m_ij
+                        else:
+                            p = q = 0
 
-        result[i] = num_electron_i
+                        result[i, (x + p) % num_x, (y + q) % num_y] += 1
 
     return result
diff --git a/optika/sensors/materials/_stern_1994.py b/optika/sensors/materials/_stern_1994.py
index 99f6615..ee336e1 100644
--- a/optika/sensors/materials/_stern_1994.py
+++ b/optika/sensors/materials/_stern_1994.py
@@ -3,4 +3,5 @@
 _thickness_oxide = 50 * u.AA
 _thickness_substrate = 7 * u.um
 _thickness_implant = 2317 * u.AA
+_width_pixel = 27 * u.um
 _cce_backsurface = 0.21