diff --git a/cli/simulate_pixels.py b/cli/simulate_pixels.py
index a6dde523..4d051ed6 100644
--- a/cli/simulate_pixels.py
+++ b/cli/simulate_pixels.py
@@ -1357,33 +1357,52 @@ def save_results(event_times, results, i_trig, i_mod=-1, light_only=False):
                     sorted_indices[idet] = np.argsort(light_inc[:,idet]['n_photons_det'])[::-1] # get the order in which to loop over tracks
                 ### END OF TEMPORARY FIX ###
 
+                # Create scintillation models with different component fractions for different particle types
+                # Heavy particles (neutrons, alphas, nuclei): higher fast component fraction (0.7)
+                # Light particles (electrons, muons, etc.): lower fast component fraction (0.3)
+                scint_model_heavy = np.zeros(n_light_ticks, dtype=np.float32)
+                scint_model_light = np.zeros(n_light_ticks, dtype=np.float32)
+
+                if light.USE_TRIPLE_EXPONENTIAL:
+                    # For triple exponential, vary fast_fraction while keeping intermediate_fraction constant
+                    light_sim.scintillation_array(scint_model_heavy, fast_fraction=0.7)
+                    light_sim.scintillation_array(scint_model_light, fast_fraction=0.3)
+                else:
+                    # For double exponential, vary singlet_fraction
+                    light_sim.scintillation_array(scint_model_heavy, singlet_fraction=0.7)
+                    light_sim.scintillation_array(scint_model_light, singlet_fraction=0.3)
+
+                # Prepare PDE correction factors for active optical channels
+                # PDE correction applies data/MC efficiency ratios to LUT visibility before convolutions
+                if light.ENABLE_PDE_CORRECTION:
+                    pde_correction = cp.array(light.OP_CHANNEL_PDE_CORRECTION[op_channel.get()])
+                else:
+                    pde_correction = cp.ones(op_channel.shape[0], dtype='f4')
+
                 TPB = (1,64)
                 BPG = (max(ceil(light_sample_inc.shape[0] / TPB[0]),1),
                         max(ceil(light_sample_inc.shape[1] / TPB[1]),1))
-                light_sim.sum_light_signals[BPG, TPB](
+
+                # Apply LAr scintillation convolution per segment with particle-dependent singlet fraction,
+                # then sum the convolved segments
+                light_sim.sum_light_signals_with_scintillation[BPG, TPB](
                     all_selected_tracks, track_light_voxel[batch_mask], selected_track_id,
                     light_inc, op_channel, lut, light_t_start, light_sample_inc, light_sample_inc_true_track_id,
-                    light_sample_inc_true_photons, sorted_indices, t0_profile_length)
+                    light_sample_inc_true_photons, sorted_indices, t0_profile_length, scint_model_heavy, scint_model_light, pde_correction)
                 RangePop()
                 if light_sample_inc_true_track_id.shape[-1] > 0 and cp.any(light_sample_inc_true_track_id[...,-1] != -1):
                     warnings.warn(f"Maximum number of true segments ({sim.MAX_MC_TRUTH_IDS}) reached in backtracking info, consider increasing MAX_MC_TRUTH_IDS (larndsim/consts/light.py)")
 
-                RangePush("sim_scintillation", 4)
-                light_sample_inc_scint = cp.zeros_like(light_sample_inc)
-                light_sample_inc_scint_true_track_id = cp.full_like(light_sample_inc_true_track_id, -1)
-                light_sample_inc_scint_true_photons = cp.zeros_like(light_sample_inc_true_photons)
-                scint_model = np.zeros(n_light_ticks, dtype=np.float32)
-                light_sim.scintillation_array(scint_model)
-                light_sim.calc_scintillation_effect[BPG, TPB](
-                    light_sample_inc, light_sample_inc_true_track_id, light_sample_inc_true_photons, light_sample_inc_scint,
-                    light_sample_inc_scint_true_track_id, light_sample_inc_scint_true_photons, scint_model)
-
+                # Apply Poisson fluctuations to the summed, scintillation-convolved signal
+                RangePush("sim_poisson_fluctuations", 4)
                 light_sample_inc_disc = cp.zeros_like(light_sample_inc)
                 rng_states = maybe_create_rng_states(int(np.prod(TPB) * np.prod(BPG)),
                                                         seed=rand_seed, rng_states=rng_states)
-                light_sim.calc_stat_fluctuations[BPG, TPB](light_sample_inc_scint, light_sample_inc_disc, rng_states)
+                light_sim.calc_stat_fluctuations[BPG, TPB](light_sample_inc, light_sample_inc_disc, rng_states)
                 RangePop()
 
+                # Apply SiPM response function to the summed signal
+                # (SPE response is invariant across particles, so it's applied after summing)
                 RangePush("sim_light_det_response", 4)
                 light_response = cp.zeros_like(light_sample_inc)
                 light_response_true_track_id = cp.full_like(light_sample_inc_true_track_id, -1)
@@ -1391,7 +1410,7 @@ def save_results(event_times, results, i_trig, i_mod=-1, light_only=False):
                 sipm_response = np.zeros(n_light_ticks, dtype=np.float32)
                 light_sim.sipm_response_array(sipm_response) #precalculate the sipm_response
                 light_sim.calc_light_detector_response[BPG, TPB](
-                    light_sample_inc_disc, light_sample_inc_scint_true_track_id, light_sample_inc_scint_true_photons,
+                    light_sample_inc_disc, light_sample_inc_true_track_id, light_sample_inc_true_photons,
                     light_response, light_response_true_track_id, light_response_true_photons, light_gain, sipm_response)
                 #light_response += cp.array(light_sim.gen_light_detector_noise(light_response.shape, light_noise[op_channel.get()]))
                 RangePop()
diff --git a/larndsim/consts/light.py b/larndsim/consts/light.py
index 7a672ae7..be925540 100644
--- a/larndsim/consts/light.py
+++ b/larndsim/consts/light.py
@@ -17,6 +17,21 @@
 OP_CHANNEL_TO_TPC = np.zeros(0)
 TPC_TO_OP_CHANNEL = np.zeros((0,0))
 
+#: PDE correction factors (data/MC efficiency ratios) applied to LUT visibility
+#: Shape: (n_tpc, n_detectors_per_tpc) - one correction per detector unit
+#: Each correction is applied uniformly to all channels in that detector
+#: (e.g., all 6 channels in an ACL, or both channels in an LCM)
+#: Internally expanded to per-channel array for kernel use
+#: Default value of 1.0 means no correction
+PDE_CORRECTION_2D = np.ones((0,0))  # 2D array: (n_tpc, n_detectors_per_tpc)
+OP_CHANNEL_PDE_CORRECTION = np.ones(0)  # 1D expanded version for kernel (n_op_channel,)
+#: Detector to channel mapping: list of (tpc_id, detector_id, [channel_ids]) for each detector unit
+#: Used to expand 2D PDE corrections to per-channel array
+#: Format: [(tpc, det, [ch0, ch1, ...]), (tpc, det, [ch0, ..., ch5]), ...]
+DETECTOR_CHANNEL_MAP = []
+#: Enable PDE correction feature
+ENABLE_PDE_CORRECTION = False
+
 #: Prescale factor analogous to ScintPreScale in LArSoft FIXME
 SCINT_PRESCALE = 1
 #: Ion + excitation work function in `MeV`
@@ -27,6 +42,10 @@
 #: Pre- and post-window for light simulation [microseconds]
 LIGHT_WINDOW = (1, 10) # us
 
+#: Use triple exponential scintillation model (default: False for double exponential)
+USE_TRIPLE_EXPONENTIAL = False
+
+#: Double exponential model parameters
 #: Fraction of total light emitted from singlet state
 SINGLET_FRACTION = 0.3
 #: Singlet decay time [microseconds]
@@ -34,6 +53,19 @@
 #: Triplet decay time [microseconds]
 TAU_T = 1.530
 
+#: Triple exponential model parameters (used when USE_TRIPLE_EXPONENTIAL = True)
+#: Fast component fraction
+FAST_FRACTION = 0.3
+#: Fast component decay time [microseconds]
+TAU_FAST = 0.001 # us
+#: Intermediate component fraction
+INTERMEDIATE_FRACTION = 0.0
+#: Intermediate component decay time [microseconds]
+TAU_INTERMEDIATE = 0.1 # us (placeholder, needs tuning)
+#: Slow component fraction (calculated as 1 - FAST_FRACTION - INTERMEDIATE_FRACTION)
+#: Slow component decay time [microseconds]
+TAU_SLOW = 1.530 # us
+
 #: Conversion from PE/microsecond to ADC
 DEFAULT_LIGHT_GAIN = -2.30 # ADC * us/PE
 LIGHT_GAIN = np.zeros((0,))
@@ -78,14 +110,24 @@ def set_light_properties(detprop_file):
     global OP_CHANNEL_EFFICIENCY
     global OP_CHANNEL_TO_TPC
     global TPC_TO_OP_CHANNEL
+    global PDE_CORRECTION_2D
+    global OP_CHANNEL_PDE_CORRECTION
+    global DETECTOR_CHANNEL_MAP
+    global ENABLE_PDE_CORRECTION
 
     global ENABLE_LUT_SMEARING
     global LIGHT_TICK_SIZE
     global LIGHT_WINDOW
 
+    global USE_TRIPLE_EXPONENTIAL
     global SINGLET_FRACTION
     global TAU_S
     global TAU_T
+    global FAST_FRACTION
+    global TAU_FAST
+    global INTERMEDIATE_FRACTION
+    global TAU_INTERMEDIATE
+    global TAU_SLOW
 
     global LIGHT_GAIN
     global SIPM_RESPONSE_MODEL
@@ -119,6 +161,70 @@ def set_light_properties(detprop_file):
         if OP_CHANNEL_EFFICIENCY.size == 1:
             OP_CHANNEL_EFFICIENCY = np.full(N_OP_CHANNEL, OP_CHANNEL_EFFICIENCY)
 
+        # Load PDE correction factors (data/MC efficiency ratios)
+        # Expected shape: (n_tpc, n_detectors_per_tpc) where detectors may have varying channel counts (2 or 6)
+        ENABLE_PDE_CORRECTION = bool(detprop.get('enable_pde_correction', ENABLE_PDE_CORRECTION))
+        pde_correction_file = str(detprop.get('pde_correction_file', ''))
+
+        # Build or load detector-to-channel mapping
+        # This handles mixed detector types (e.g., 2-channel LCMs and 6-channel ACLs)
+        detector_map_config = detprop.get('detector_channel_map', None)
+
+        if detector_map_config:
+            # Load explicit mapping from YAML: [(tpc, det, [channels]), ...]
+            DETECTOR_CHANNEL_MAP = [(entry['tpc'], entry['detector'], entry['channels'])
+                                     for entry in detector_map_config]
+        else:
+            # Infer mapping assuming uniform OP_CHANNEL_PER_TRIG channels per detector
+            # This is a fallback - explicit mapping is preferred for mixed detector types
+            DETECTOR_CHANNEL_MAP = []
+            for itpc in range(n_tpc):
+                tpc_channels = TPC_TO_OP_CHANNEL[itpc]
+                n_detectors_this_tpc = len(tpc_channels) // OP_CHANNEL_PER_TRIG
+                for idet in range(n_detectors_this_tpc):
+                    det_channels = tpc_channels[idet*OP_CHANNEL_PER_TRIG:(idet+1)*OP_CHANNEL_PER_TRIG]
+                    DETECTOR_CHANNEL_MAP.append((itpc, idet, list(det_channels)))
+
+        # Determine 2D shape from detector map
+        n_detectors_per_tpc = max([det for tpc, det, _ in DETECTOR_CHANNEL_MAP]) + 1 if DETECTOR_CHANNEL_MAP else 0
+
+        if ENABLE_PDE_CORRECTION and pde_correction_file:
+            # Load from file (numpy array with shape (n_tpc, n_detectors_per_tpc))
+            try:
+                # First try to load from current directory
+                PDE_CORRECTION_2D = np.load(pde_correction_file)
+            except FileNotFoundError:
+                # Then try from larnd-sim base directory
+                try:
+                    PDE_CORRECTION_2D = np.load(os.path.join(os.path.dirname(__file__), '../../') + pde_correction_file)
+                except FileNotFoundError:
+                    print("PDE correction file not found:", pde_correction_file, ", using default correction of 1.0")
+                    PDE_CORRECTION_2D = np.ones((n_tpc, n_detectors_per_tpc))
+        else:
+            # Load from YAML or use default
+            pde_corr_yaml = detprop.get('pde_correction_2d', None)
+            if pde_corr_yaml is None:
+                PDE_CORRECTION_2D = np.ones((n_tpc, n_detectors_per_tpc))
+            else:
+                pde_corr_yaml = np.array(pde_corr_yaml)
+                if pde_corr_yaml.size == 1:
+                    # Single value: apply to all detectors
+                    PDE_CORRECTION_2D = np.full((n_tpc, n_detectors_per_tpc), pde_corr_yaml)
+                else:
+                    PDE_CORRECTION_2D = pde_corr_yaml.reshape(n_tpc, n_detectors_per_tpc)
+
+        # Validate 2D shape
+        if PDE_CORRECTION_2D.shape[0] != n_tpc:
+            raise ValueError(f"PDE correction array has {PDE_CORRECTION_2D.shape[0]} TPCs, expected {n_tpc}")
+
+        # Expand 2D array to 1D per-channel array using detector-channel mapping
+        # This correctly handles mixed detector types (2-ch and 6-ch)
+        OP_CHANNEL_PDE_CORRECTION = np.ones(N_OP_CHANNEL)
+        for tpc_id, det_id, channels in DETECTOR_CHANNEL_MAP:
+            correction = PDE_CORRECTION_2D[tpc_id, det_id]
+            for ch in channels:
+                OP_CHANNEL_PDE_CORRECTION[ch] = correction
+
         try:
             tpc_to_op_channel = detprop['tpc_to_op_channel']
             OP_CHANNEL_TO_TPC = np.zeros((N_OP_CHANNEL,), int)
@@ -141,10 +247,20 @@ def set_light_properties(detprop_file):
         LIGHT_WINDOW = tuple(detprop.get('light_window', LIGHT_WINDOW))
         assert len(LIGHT_WINDOW) == 2
 
+        USE_TRIPLE_EXPONENTIAL = bool(detprop.get('use_triple_exponential', USE_TRIPLE_EXPONENTIAL))
+
+        # Double exponential model parameters
         SINGLET_FRACTION = float(detprop.get('singlet_fraction', SINGLET_FRACTION))
         TAU_S = float(detprop.get('tau_s', TAU_S))
         TAU_T = float(detprop.get('tau_t', TAU_T))
 
+        # Triple exponential model parameters
+        FAST_FRACTION = float(detprop.get('fast_fraction', FAST_FRACTION))
+        TAU_FAST = float(detprop.get('tau_fast', TAU_FAST))
+        INTERMEDIATE_FRACTION = float(detprop.get('intermediate_fraction', INTERMEDIATE_FRACTION))
+        TAU_INTERMEDIATE = float(detprop.get('tau_intermediate', TAU_INTERMEDIATE))
+        TAU_SLOW = float(detprop.get('tau_slow', TAU_SLOW))
+
         LIGHT_GAIN = np.array(detprop.get('light_gain', [DEFAULT_LIGHT_GAIN]))
         if LIGHT_GAIN.size == 1:
             LIGHT_GAIN = np.full(N_OP_CHANNEL, LIGHT_GAIN)
diff --git a/larndsim/light_sim.py b/larndsim/light_sim.py
index 7c166918..275deffe 100644
--- a/larndsim/light_sim.py
+++ b/larndsim/light_sim.py
@@ -56,9 +56,11 @@ def get_active_op_channel(light_incidence):
     return cp.empty((0,), dtype='i4')
     
 @cuda.jit
-def sum_light_signals(segments, segment_voxel, segment_track_id, light_inc, op_channel, lut, start_time, light_sample_inc, light_sample_inc_true_track_id, light_sample_inc_true_photons, sorted_indices, t0_profile_length):
+def sum_light_signals_with_scintillation(segments, segment_voxel, segment_track_id, light_inc, op_channel, lut, start_time, light_sample_inc, light_sample_inc_true_track_id, light_sample_inc_true_photons, sorted_indices, t0_profile_length, scint_model_heavy, scint_model_light, pde_correction):
     """
-    Sums the number of photons observed by each light detector at each time tick
+    Sums the number of photons observed by each light detector at each time tick,
+    applying LAr scintillation convolution per segment with particle-dependent singlet fraction
+    before summing. This allows for variation in pulse shape based on particle type.
 
     Args:
         segments(array): shape `(ntracks,)`, edep-sim tracks to simulate
@@ -68,10 +70,14 @@ def sum_light_signals(segments, segment_voxel, segment_track_id, light_inc, op_c
         op_channel(array): shape `(ntracks, ndet_active)`, optical channel index, will use lut[:,:,:,op_channel%lut.shape[3]] to look up timing information
         lut(array): shape `(nx,ny,nz,ndet_tpc)`, light look up table
         start_time(float): start time of light simulation in microseconds
-        light_sample_inc(array): output array, shape `(ndet, nticks)`, number of photons incident on each detector at each time tick (propogation delay only)
-        light_sample_inc_true_track_id(array): output array, shape `(ndet, nticks, maxtracks)`, true track ids on each detector at each time tick (propogation delay only)
+        light_sample_inc(array): output array, shape `(ndet, nticks)`, number of photons incident on each detector at each time tick (after scintillation convolution per segment)
+        light_sample_inc_true_track_id(array): output array, shape `(ndet, nticks, maxtracks)`, true track ids on each detector at each time tick
         light_sample_inc_true_photons(array): output array, shape `(ndet, nticks, maxtracks)`, number of photons incident on each detector at each time tick from each track
         sorted_indices(array): shape `(maxtracks,)`, indices of segments sorted by how much light they contribute
+        t0_profile_length(int): length of the LUT time profile
+        scint_model_heavy(array): shape `(nticks,)`, scintillation model for heavy particles (neutrons, alphas, nuclei) with singlet_fraction=0.7
+        scint_model_light(array): shape `(nticks,)`, scintillation model for light particles (electrons, muons, etc.) with singlet_fraction=0.3
+        pde_correction(array): shape `(ndet_active,)`, PDE correction factors (data/MC efficiency ratios) applied to LUT visibility before convolutions
     """
     idet,itick = cuda.grid(2)
 
@@ -89,75 +95,162 @@ def sum_light_signals(segments, segment_voxel, segment_track_id, light_inc, op_c
                     track_time = segments[itrk]['t0']
                     track_end_time = track_time + t0_profile_length * units.ns / units.mus # FIXME: assumes light LUT time profile bins are 1ns (might not be true in general)
 
-                    if track_end_time < start_tick_time or track_time > end_tick_time:
-                        continue
+                    # Determine which scintillation model to use based on particle type
+                    # PDG codes: neutron=2112, alpha=1000020040
+                    # Heavy particles (neutrons, alphas, nuclei) use singlet_fraction=0.7
+                    # Light particles (everything else) use singlet_fraction=0.3
+                    pdg_id = segments[itrk]['pdg_id']
+                    use_heavy_model = False
+                    if pdg_id == 2112:  # neutron
+                        use_heavy_model = True
+                    elif pdg_id > 1000000000:  # nuclei (including alphas)
+                        use_heavy_model = True
 
                     # use LUT time smearing
                     if light.ENABLE_LUT_SMEARING:
                         time_profile = lut[voxel[0],voxel[1],voxel[2],idet_lut]['time_dist'] # normalised
 
-                        # add photons to time tick
+                        # apply scintillation convolution for each profile time bin
                         for iprof in range(time_profile.shape[0]):
                             profile_time = track_time + iprof * units.ns / units.mus # FIXME: assumes light LUT time profile bins are 1ns (might not be true in general)
-                            if profile_time < end_tick_time and profile_time > start_tick_time:
-                                photons = light_inc['n_photons_det'][itrk,op_channel[idet]] * time_profile[iprof] / light.LIGHT_TICK_SIZE
-                                light_sample_inc[idet,itick] += photons
-
-                                if photons > sim.MC_TRUTH_THRESHOLD:
-                                    # get truth information for time tick
-                                    for itrue in range(light_sample_inc_true_track_id.shape[-1]):
-                                        if light_sample_inc_true_track_id[idet,itick,itrue] == -1 or light_sample_inc_true_track_id[idet,itick,itrue] == segment_track_id[itrk]:
-                                            light_sample_inc_true_track_id[idet,itick,itrue] = segment_track_id[itrk]
-                                            light_sample_inc_true_photons[idet,itick,itrue] += photons
-                                            break
+
+                            # Calculate which tick this profile time corresponds to (before convolution)
+                            profile_tick_float = (profile_time - start_time) / light.LIGHT_TICK_SIZE
+                            profile_tick = int(profile_tick_float)
+
+                            if profile_tick < 0 or profile_tick >= light_sample_inc.shape[1]:
+                                continue
+
+                            # Base photon rate from LUT for this profile bin
+                            base_photons = light_inc['n_photons_det'][itrk,op_channel[idet]] * time_profile[iprof] / light.LIGHT_TICK_SIZE
+
+                            # Apply PDE correction (data/MC efficiency ratio) to LUT visibility
+                            # This corrects for measured in-situ detector efficiency
+                            base_photons *= pde_correction[idet]
+
+                            if base_photons <= 0:
+                                continue
+
+                            # Apply scintillation convolution: this profile bin at profile_tick contributes to future ticks
+                            # The current tick (itick) receives contribution from profile bins that occurred (itick - profile_tick) ticks ago
+                            conv_offset = itick - profile_tick
+                            if conv_offset >= 0 and conv_offset < scint_model_heavy.shape[0]:
+                                if use_heavy_model:
+                                    tick_weight = scint_model_heavy[conv_offset]
+                                else:
+                                    tick_weight = scint_model_light[conv_offset]
+
+                                photons = base_photons * tick_weight
+
+                                if photons > 0:
+                                    light_sample_inc[idet,itick] += photons
+
+                                    if photons > sim.MC_TRUTH_THRESHOLD:
+                                        # get truth information for time tick
+                                        for itrue in range(light_sample_inc_true_track_id.shape[-1]):
+                                            if light_sample_inc_true_track_id[idet,itick,itrue] == -1 or light_sample_inc_true_track_id[idet,itick,itrue] == segment_track_id[itrk]:
+                                                light_sample_inc_true_track_id[idet,itick,itrue] = segment_track_id[itrk]
+                                                light_sample_inc_true_photons[idet,itick,itrue] += photons
+                                                break
                     # use average time only
                     else:
                         # calculate average delay time
-                        t0_avg = lut[voxel[0],voxel[1],voxel[2],idet_lut]['t0_avg'] * units.ns / units.mus # normalised averagein us
+                        t0_avg = lut[voxel[0],voxel[1],voxel[2],idet_lut]['t0_avg'] * units.ns / units.mus # normalised average in us
 
-                        # add photons to time tick
+                        # Calculate which tick this corresponds to (before convolution)
                         profile_time = track_time + t0_avg
-                        if profile_time < end_tick_time and profile_time > start_tick_time:
-                            photons = light_inc['n_photons_det'][itrk,op_channel[idet]] / light.LIGHT_TICK_SIZE
-                            light_sample_inc[idet,itick] += photons
-                            if photons > sim.MC_TRUTH_THRESHOLD:
-                                # get truth information for time tick
-                                for itrue in range(light_sample_inc_true_track_id.shape[-1]):
-                                    if light_sample_inc_true_track_id[idet,itick,itrue] == -1 or light_sample_inc_true_track_id[idet,itick,itrue] == segment_track_id[itrk]:
-                                        light_sample_inc_true_track_id[idet,itick,itrue] = segment_track_id[itrk]
-                                        light_sample_inc_true_photons[idet,itick,itrue] += photons
-                                        break
+                        profile_tick_float = (profile_time - start_time) / light.LIGHT_TICK_SIZE
+                        profile_tick = int(profile_tick_float)
+
+                        if profile_tick >= 0 and profile_tick < light_sample_inc.shape[1]:
+                            # Base photon rate
+                            base_photons = light_inc['n_photons_det'][itrk,op_channel[idet]] / light.LIGHT_TICK_SIZE
+
+                            # Apply PDE correction (data/MC efficiency ratio) to LUT visibility
+                            base_photons *= pde_correction[idet]
+
+                            if base_photons > 0:
+                                # Apply scintillation convolution
+                                conv_offset = itick - profile_tick
+                                if conv_offset >= 0 and conv_offset < scint_model_heavy.shape[0]:
+                                    if use_heavy_model:
+                                        tick_weight = scint_model_heavy[conv_offset]
+                                    else:
+                                        tick_weight = scint_model_light[conv_offset]
+
+                                    photons = base_photons * tick_weight
+
+                                    if photons > 0:
+                                        light_sample_inc[idet,itick] += photons
+
+                                        if photons > sim.MC_TRUTH_THRESHOLD:
+                                            # get truth information for time tick
+                                            for itrue in range(light_sample_inc_true_track_id.shape[-1]):
+                                                if light_sample_inc_true_track_id[idet,itick,itrue] == -1 or light_sample_inc_true_track_id[idet,itick,itrue] == segment_track_id[itrk]:
+                                                    light_sample_inc_true_track_id[idet,itick,itrue] = segment_track_id[itrk]
+                                                    light_sample_inc_true_photons[idet,itick,itrue] += photons
+                                                    break
 
 @nb.njit
 def scintillation_model(time_tick):
     """
-    Calculates the fraction of scintillation photons emitted 
+    Calculates the fraction of scintillation photons emitted
     during time interval `time_tick` to `time_tick + 1`
-    
+    Supports both double and triple exponential models.
+
     Args:
         time_tick(int): time tick relative to t0
-    
+
     Returns:
         float: fraction of scintillation photons
     """
-    p1 = light.SINGLET_FRACTION * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_S) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_S))
-    p3 = (1 - light.SINGLET_FRACTION) * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_T) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_T))
-    return (p1 + p3) * (time_tick >= 0)
+    if light.USE_TRIPLE_EXPONENTIAL:
+        # Triple exponential model with fast, intermediate, and slow components
+        p1 = light.FAST_FRACTION * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_FAST) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_FAST))
+        p2 = light.INTERMEDIATE_FRACTION * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_INTERMEDIATE) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_INTERMEDIATE))
+        slow_fraction = 1.0 - light.FAST_FRACTION - light.INTERMEDIATE_FRACTION
+        p3 = slow_fraction * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_SLOW) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_SLOW))
+        return (p1 + p2 + p3) * (time_tick >= 0)
+    else:
+        # Standard double exponential model (singlet and triplet)
+        p1 = light.SINGLET_FRACTION * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_S) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_S))
+        p3 = (1 - light.SINGLET_FRACTION) * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_T) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_T))
+        return (p1 + p3) * (time_tick >= 0)
 
 @nb.njit
-def scintillation_array(scint_model):
+def scintillation_array(scint_model, singlet_fraction=None, fast_fraction=None, intermediate_fraction=None):
     """
-    Calculates the fraction of scintillation photons emitted 
+    Calculates the fraction of scintillation photons emitted
     during time interval `time_tick` to `time_tick + 1` for
-    the entire input array.
-    
+    the entire input array. Supports both double and triple exponential models.
+
     Args:
         scint_model: array to store result
-    """
-    for time_tick in range(scint_model.shape[0]):
-        p1 = light.SINGLET_FRACTION * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_S) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_S))
-        p3 = (1 - light.SINGLET_FRACTION) * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_T) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_T))
-        scint_model[time_tick] = p1 + p3
+        singlet_fraction: optional singlet fraction for double exponential (defaults to light.SINGLET_FRACTION)
+        fast_fraction: optional fast fraction for triple exponential (defaults to light.FAST_FRACTION)
+        intermediate_fraction: optional intermediate fraction for triple exponential (defaults to light.INTERMEDIATE_FRACTION)
+    """
+    if light.USE_TRIPLE_EXPONENTIAL:
+        # Triple exponential model with fast, intermediate, and slow components
+        if fast_fraction is None:
+            fast_fraction = light.FAST_FRACTION
+        if intermediate_fraction is None:
+            intermediate_fraction = light.INTERMEDIATE_FRACTION
+        slow_fraction = 1.0 - fast_fraction - intermediate_fraction
+
+        for time_tick in range(scint_model.shape[0]):
+            p1 = fast_fraction * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_FAST) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_FAST))
+            p2 = intermediate_fraction * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_INTERMEDIATE) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_INTERMEDIATE))
+            p3 = slow_fraction * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_SLOW) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_SLOW))
+            scint_model[time_tick] = p1 + p2 + p3
+    else:
+        # Standard double exponential model (singlet and triplet)
+        if singlet_fraction is None:
+            singlet_fraction = light.SINGLET_FRACTION
+        for time_tick in range(scint_model.shape[0]):
+            p1 = singlet_fraction * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_S) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_S))
+            p3 = (1 - singlet_fraction) * exp(-time_tick * light.LIGHT_TICK_SIZE / light.TAU_T) * (1 - exp(-light.LIGHT_TICK_SIZE / light.TAU_T))
+            scint_model[time_tick] = p1 + p3
 
 @cuda.jit
 def calc_scintillation_effect(light_sample_inc, light_sample_inc_true_track_id, light_sample_inc_true_photons, light_sample_inc_scint, light_sample_inc_scint_true_track_id, light_sample_inc_scint_true_photons, scint_model):