From 9f26884f3bf2dde3fdee155d84cf90878f6554a8 Mon Sep 17 00:00:00 2001 From: James Vincent Mead Date: Tue, 16 Dec 2025 06:22:37 -0800 Subject: [PATCH 1/3] Implement per-segment LAr scintillation with particle-dependent parameters MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Modified light simulation to apply LAr pulse shape convolution per segment before summing, allowing for particle-by-particle variation in pulse shape - Added variable singlet fraction based on PDG code: * Neutrons and nuclei (including alphas): singlet_fraction = 0.7 * All other particles: singlet_fraction = 0.3 - Implemented optional triple exponential scintillation model: * New constants: USE_TRIPLE_EXPONENTIAL, FAST_FRACTION, INTERMEDIATE_FRACTION, TAU_FAST, TAU_INTERMEDIATE, TAU_SLOW * Default remains double exponential for compatibility * Can be enabled via detector properties YAML - Order of operations: LAr convolution (per-segment) → Sum → Poisson → SPE - Updated sum_light_signals to sum_light_signals_with_scintillation kernel - Truth backtracking logic verified and maintained correctly 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.5 --- cli/simulate_pixels.py | 40 +++++---- larndsim/consts/light.py | 33 ++++++++ larndsim/light_sim.py | 175 +++++++++++++++++++++++++++++---------- 3 files changed, 189 insertions(+), 59 deletions(-) diff --git a/cli/simulate_pixels.py b/cli/simulate_pixels.py index 54a9d60a..9d2b0101 100644 --- a/cli/simulate_pixels.py +++ b/cli/simulate_pixels.py @@ -1353,33 +1353,45 @@ 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) + 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) 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) @@ -1387,7 +1399,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..244785ec 100644 --- a/larndsim/consts/light.py +++ b/larndsim/consts/light.py @@ -27,6 +27,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 +38,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,)) @@ -83,9 +100,15 @@ def set_light_properties(detprop_file): 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 @@ -141,10 +164,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 222475bb..16eae9a5 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): """ - 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,13 @@ 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 """ idet,itick = cuda.grid(2) @@ -89,75 +94,155 @@ 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 + + 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 + + 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): From 2cf8fe83a85485d0c46348848ec0d073ef3c0dec Mon Sep 17 00:00:00 2001 From: James Vincent Mead Date: Tue, 16 Dec 2025 06:53:37 -0800 Subject: [PATCH 2/3] Add PDE correction feature for in-situ detector efficiency corrections MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Implemented optional PDE correction to apply data/MC efficiency ratios to LUT visibility before convolutions and Poisson fluctuations - Added constants ENABLE_PDE_CORRECTION and OP_CHANNEL_PDE_CORRECTION - Supports loading corrections from: * Numpy file (pde_correction_file in YAML) * Direct YAML specification (single value or array) * Default of 1.0 (no correction) when disabled - Modified sum_light_signals_with_scintillation kernel to apply corrections immediately after LUT lookup, before scintillation convolution - Corrections apply uniformly across all channels in each detector unit (2 channels for LCMs, 6 channels for ACLs) - Order: LUT → PDE correction → LAr convolution → Sum → Poisson → SPE - Validated: syntax checks pass, shape validation included 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.5 --- cli/simulate_pixels.py | 9 ++++++++- larndsim/consts/light.py | 36 ++++++++++++++++++++++++++++++++++++ larndsim/light_sim.py | 10 +++++++++- 3 files changed, 53 insertions(+), 2 deletions(-) diff --git a/cli/simulate_pixels.py b/cli/simulate_pixels.py index 8e78d0ee..4d051ed6 100644 --- a/cli/simulate_pixels.py +++ b/cli/simulate_pixels.py @@ -1372,6 +1372,13 @@ def save_results(event_times, results, i_trig, i_mod=-1, light_only=False): 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)) @@ -1381,7 +1388,7 @@ def save_results(event_times, results, i_trig, i_mod=-1, light_only=False): 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, scint_model_heavy, scint_model_light) + 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)") diff --git a/larndsim/consts/light.py b/larndsim/consts/light.py index 244785ec..4bbb86c5 100644 --- a/larndsim/consts/light.py +++ b/larndsim/consts/light.py @@ -17,6 +17,14 @@ 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_op_channel,) - one correction per optical channel +#: Applied uniformly across all channels in the same detector unit +#: Default value of 1.0 means no correction +OP_CHANNEL_PDE_CORRECTION = np.ones(0) +#: 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` @@ -95,6 +103,8 @@ def set_light_properties(detprop_file): global OP_CHANNEL_EFFICIENCY global OP_CHANNEL_TO_TPC global TPC_TO_OP_CHANNEL + global OP_CHANNEL_PDE_CORRECTION + global ENABLE_PDE_CORRECTION global ENABLE_LUT_SMEARING global LIGHT_TICK_SIZE @@ -142,6 +152,32 @@ 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) + ENABLE_PDE_CORRECTION = bool(detprop.get('enable_pde_correction', ENABLE_PDE_CORRECTION)) + pde_correction_file = str(detprop.get('pde_correction_file', '')) + + if ENABLE_PDE_CORRECTION and pde_correction_file: + # Load from file (numpy array) + try: + # First try to load from current directory + OP_CHANNEL_PDE_CORRECTION = np.load(pde_correction_file) + except FileNotFoundError: + # Then try from larnd-sim base directory + try: + OP_CHANNEL_PDE_CORRECTION = 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") + OP_CHANNEL_PDE_CORRECTION = np.ones(N_OP_CHANNEL) + else: + # Load from YAML or use default + OP_CHANNEL_PDE_CORRECTION = np.array(detprop.get('op_channel_pde_correction', np.ones(N_OP_CHANNEL))) + + # Ensure correct shape + if OP_CHANNEL_PDE_CORRECTION.size == 1: + OP_CHANNEL_PDE_CORRECTION = np.full(N_OP_CHANNEL, OP_CHANNEL_PDE_CORRECTION) + elif OP_CHANNEL_PDE_CORRECTION.size != N_OP_CHANNEL: + raise ValueError(f"PDE correction array size ({OP_CHANNEL_PDE_CORRECTION.size}) must match N_OP_CHANNEL ({N_OP_CHANNEL})") + try: tpc_to_op_channel = detprop['tpc_to_op_channel'] OP_CHANNEL_TO_TPC = np.zeros((N_OP_CHANNEL,), int) diff --git a/larndsim/light_sim.py b/larndsim/light_sim.py index 41d74568..275deffe 100644 --- a/larndsim/light_sim.py +++ b/larndsim/light_sim.py @@ -56,7 +56,7 @@ def get_active_op_channel(light_incidence): return cp.empty((0,), dtype='i4') @cuda.jit -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): +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, applying LAr scintillation convolution per segment with particle-dependent singlet fraction @@ -77,6 +77,7 @@ def sum_light_signals_with_scintillation(segments, segment_voxel, segment_track_ 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) @@ -123,6 +124,10 @@ def sum_light_signals_with_scintillation(segments, segment_voxel, segment_track_ # 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 @@ -161,6 +166,9 @@ def sum_light_signals_with_scintillation(segments, segment_voxel, segment_track_ # 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 From 95b375d41ab6b92abd7f04546a68f0de7eb47214 Mon Sep 17 00:00:00 2001 From: James Vincent Mead Date: Tue, 16 Dec 2025 07:38:14 -0800 Subject: [PATCH 3/3] Refactor PDE correction to support mixed detector types (2-ch and 6-ch) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Changed from uniform channel-per-detector assumption to flexible mapping - Added DETECTOR_CHANNEL_MAP to explicitly track (TPC, detector, channels) - Supports explicit mapping via YAML config for mixed detector types: * 2-channel LCMs * 6-channel ACLs - Falls back to inferred uniform mapping if not specified - PDE_CORRECTION_2D remains (n_tpc, n_detectors_per_tpc) structure - Expansion to per-channel array now uses detector mapping for correctness - Enables realistic in-situ measurements with varying detector geometries Example config: detector_channel_map: - {tpc: 0, detector: 0, channels: [0,1,2,3,4,5]} # ACL - {tpc: 0, detector: 1, channels: [6,7]} # LCM 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Sonnet 4.5 --- larndsim/consts/light.py | 75 ++++++++++++++++++++++++++++++++-------- 1 file changed, 61 insertions(+), 14 deletions(-) diff --git a/larndsim/consts/light.py b/larndsim/consts/light.py index 4bbb86c5..be925540 100644 --- a/larndsim/consts/light.py +++ b/larndsim/consts/light.py @@ -18,10 +18,17 @@ TPC_TO_OP_CHANNEL = np.zeros((0,0)) #: PDE correction factors (data/MC efficiency ratios) applied to LUT visibility -#: Shape: (n_op_channel,) - one correction per optical channel -#: Applied uniformly across all channels in the same detector unit +#: 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 -OP_CHANNEL_PDE_CORRECTION = np.ones(0) +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 @@ -103,7 +110,9 @@ 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 @@ -153,30 +162,68 @@ def set_light_properties(detprop_file): 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) + # Load from file (numpy array with shape (n_tpc, n_detectors_per_tpc)) try: # First try to load from current directory - OP_CHANNEL_PDE_CORRECTION = np.load(pde_correction_file) + PDE_CORRECTION_2D = np.load(pde_correction_file) except FileNotFoundError: # Then try from larnd-sim base directory try: - OP_CHANNEL_PDE_CORRECTION = np.load(os.path.join(os.path.dirname(__file__), '../../') + pde_correction_file) + 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") - OP_CHANNEL_PDE_CORRECTION = np.ones(N_OP_CHANNEL) + PDE_CORRECTION_2D = np.ones((n_tpc, n_detectors_per_tpc)) else: # Load from YAML or use default - OP_CHANNEL_PDE_CORRECTION = np.array(detprop.get('op_channel_pde_correction', np.ones(N_OP_CHANNEL))) - - # Ensure correct shape - if OP_CHANNEL_PDE_CORRECTION.size == 1: - OP_CHANNEL_PDE_CORRECTION = np.full(N_OP_CHANNEL, OP_CHANNEL_PDE_CORRECTION) - elif OP_CHANNEL_PDE_CORRECTION.size != N_OP_CHANNEL: - raise ValueError(f"PDE correction array size ({OP_CHANNEL_PDE_CORRECTION.size}) must match N_OP_CHANNEL ({N_OP_CHANNEL})") + 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']