diff --git a/cli/simulate_pixels.py b/cli/simulate_pixels.py index a71ccaae..d94a0ad6 100755 --- a/cli/simulate_pixels.py +++ b/cli/simulate_pixels.py @@ -187,6 +187,8 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o results['unique_pix'], results['current_fractions'], results['track_pixel_map'], + results['traj_pixel_map'], + #results['segidx2trackid'], output_filename, # defined earlier in script uniq_event_times, is_first_batch=is_first_batch, @@ -454,17 +456,22 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o # First of all we load the edep-sim output with h5py.File(input_filename, 'r') as f: tracks = np.array(f['segments']) - if 'segment_id' in tracks.dtype.names: - segment_ids = tracks['segment_id'] - else: - dtype = tracks.dtype.descr - dtype = [('segment_id','u4')] + dtype - new_tracks = np.empty(tracks.shape, dtype=np.dtype(dtype, align=True)) - new_tracks['segment_id'] = np.arange(tracks.shape[0], dtype='u4') - for field in dtype[1:]: - new_tracks[field[0]] = tracks[field[0]] - tracks = new_tracks - segment_ids = tracks['segment_id'] + #if 'segment_id' in tracks.dtype.names: + # segment_ids = tracks['segment_id'] + # trajectory_ids = tracks['file_traj_id'] + #else: + # dtype = tracks.dtype.descr + # dtype = [('segment_id','u4')] + dtype + # new_tracks = np.empty(tracks.shape, dtype=np.dtype(dtype, align=True)) + # new_tracks['segment_id'] = np.arange(tracks.shape[0], dtype='u4') + # for field in dtype[1:]: + # new_tracks[field[0]] = tracks[field[0]] + # tracks = new_tracks + # segment_ids = tracks['segment_id'] + + segment_ids = tracks['segment_id'] + trajectory_ids = tracks['file_traj_id'] + try: trajectories = np.array(f['trajectories']) input_has_trajectories = True @@ -506,6 +513,7 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o print(f'Selecting only the first {n_events} events for simulation.') max_eventID = np.unique(tracks[sim.EVENT_SEPARATOR])[n_events-1] segment_ids = segment_ids[tracks[sim.EVENT_SEPARATOR] <= max_eventID] + trajectory_ids = trajectory_ids[tracks[sim.EVENT_SEPARATOR] <= max_eventID] tracks = tracks[tracks[sim.EVENT_SEPARATOR] <= max_eventID] if input_has_trajectories: @@ -604,6 +612,7 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o # First copy all tracks and segment_ids all_mod_tracks = tracks all_mod_segment_ids = segment_ids + all_mod_trajectory_ids = trajectory_ids if mod2mod_variation == None or mod2mod_variation == False: mod_ids = [-1] # Sub-select only segments in active volumes @@ -613,6 +622,7 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o active_tracks_mask = active_volume.select_active_volume(all_mod_tracks, detector.TPC_BORDERS) tracks = all_mod_tracks[active_tracks_mask] segment_ids = all_mod_segment_ids[active_tracks_mask] + trajectory_ids = all_mod_trajectory_ids[active_tracks_mask] end_mask = time() print(f" {end_mask-start_mask:.2f} s") else: @@ -622,6 +632,7 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o # Convention module counting start from 1 for i_mod in mod_ids: + print(f'Simulating module {i_mod}') if mod2mod_variation: consts.detector.set_detector_properties(detector_properties, pixel_layout, i_mod) # Currently shouln't be necessary to reload light props, but if @@ -646,6 +657,7 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o module_tracks_mask = active_volume.select_active_volume(all_mod_tracks, module_borders) tracks = all_mod_tracks[module_tracks_mask] segment_ids = all_mod_segment_ids[module_tracks_mask] + trajectory_ids = all_mod_trajectory_ids[module_tracks_mask] RangePop() # find the module that triggers @@ -763,6 +775,8 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o track_ids = cp.asarray(np.arange(segment_ids.shape[0], dtype=int)) segment_ids_arr = cp.asarray(segment_ids) + trajectory_ids_arr = cp.asarray(trajectory_ids) + # We divide the sample in portions that can be processed by the GPU is_first_batch = True logger.start() @@ -771,8 +785,16 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o i_trig = 0 sync_start = event_times[0] // (fee.CLOCK_RESET_PERIOD * fee.CLOCK_CYCLE) * (fee.CLOCK_RESET_PERIOD * fee.CLOCK_CYCLE) + (fee.CLOCK_RESET_PERIOD * fee.CLOCK_CYCLE) det_borders = module_borders if mod2mod_variation else detector.TPC_BORDERS + #print('Detector boarders',det_borders) + i_batch=0 + #for batch_mask in tqdm(batching.TPCBatcher(all_mod_tracks, tracks, sim.EVENT_SEPARATOR, tpc_batch_size=sim.EVENT_BATCH_SIZE, tpc_borders=det_borders), + # desc='Simulating batches...', ncols=80, smoothing=0): + # i_batch = i_batch+1 + # print(i_batch,np.unique(all_mod_tracks[sim.EVENT_SEPARATOR])) + i_batch = 0 for batch_mask in tqdm(batching.TPCBatcher(all_mod_tracks, tracks, sim.EVENT_SEPARATOR, tpc_batch_size=sim.EVENT_BATCH_SIZE, tpc_borders=det_borders), desc='Simulating batches...', ncols=80, smoothing=0): + #print('Batch index',i_batch) i_batch = i_batch+1 # grab only tracks from current batch track_subset = tracks[batch_mask] @@ -804,13 +826,19 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o del null_light_results_acc['light_event_id'] # Nothing to simulate for charge readout? continue - + #print('Looping over',evt_tracks.shape[0],'tracks') for itrk in tqdm(range(0, evt_tracks.shape[0], sim.BATCH_SIZE), delay=1, desc=' Simulating event %i batches...' % ievd, leave=False, ncols=80): if itrk > 0: warnings.warn(f"Entered sub-batch loop, results may not be accurate! Consider increasing batch_size (currently {sim.BATCH_SIZE}) in the simulation_properties file.") selected_tracks = evt_tracks[itrk:itrk+sim.BATCH_SIZE] + #segidx2trackid = [] + #print('Tracks in batch',selected_tracks.dtype.fields.keys()) + #for i,seg in enumerate(selected_tracks): + # #print(i,seg['file_traj_id'],seg['segment_id']) + # segidx2trackid.append(seg['file_traj_id']) + #segidx2trackid=np.array(segidx2trackid) RangePush("event_id_map") event_ids = selected_tracks[sim.EVENT_SEPARATOR] @@ -835,13 +863,14 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o active_pixels = cp.full((selected_tracks.shape[0], max_pixels[0]), -1, dtype=np.int32) neighboring_pixels = cp.full((selected_tracks.shape[0], max_neighboring_pixels), -1, dtype=np.int32) + neighboring_radius = cp.full((selected_tracks.shape[0], max_neighboring_pixels), -1, dtype=np.int32) n_pixels_list = cp.zeros(shape=(selected_tracks.shape[0])) if not active_pixels.shape[1] or not neighboring_pixels.shape[1]: if light.LIGHT_SIMULATED and (light.LIGHT_TRIG_MODE == 0 or light.LIGHT_TRIG_MODE == 1): null_light_results_acc['light_event_id'].append(cp.full(1, ievd)) # one event save_results(event_times, is_first_batch, null_light_results_acc, i_trig, i_mod, light_only=True) - i_trig += 1 # add to the trigger counter + i_trig += 1 # add to the trigger countern_max_pixels del null_light_results_acc['light_event_id'] continue @@ -849,6 +878,7 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o pixels_from_track.get_pixels[BPG,TPB](selected_tracks, active_pixels, neighboring_pixels, + neighboring_radius, n_pixels_list, max_radius) RangePop() @@ -860,6 +890,34 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o unique_pix = unique_pix[(unique_pix != -1)] RangePop() + def invert_array_map(in_map,pix_set): + ''' + Invert the map of unique segment id => a set of unique pixel IDs to a map of unique + pixel index => a set of segment IDs. + + Args: + in_map (:obj:`numpy.ndarray`): 2D array where segment index => list of pixel IDs + pix_set (:obj:`numpy.ndarray`): 1D array containing all unique pixel IDs + Returns: + ndarray: 2D array where pixel index => list of segment index + ''' + pixids,counts=cp.unique(in_map[in_map>=0].flatten(),return_counts=True) + + pix_id2idx = {val.item():i for i,val in enumerate(pix_set)} + + mymap=cp.full(shape=(pix_set.shape[0],counts.max().item()),fill_value=-1,dtype=int) + curr_idx=cp.zeros(shape=(len(pix_id2idx),),dtype=int) + for seg_idx in range(in_map.shape[0]): + ass = in_map[seg_idx] + for pixid in ass: + if pixid<0: break + pix_idx = pix_id2idx[pixid.item()] + mymap[pix_idx][curr_idx[pix_idx]]=seg_idx + curr_idx[pix_idx] += 1 + return mymap + + assmap_pix2seg = invert_array_map(active_pixels,unique_pix) + if not unique_pix.shape[0]: if light.LIGHT_SIMULATED and (light.LIGHT_TRIG_MODE == 0 or light.LIGHT_TRIG_MODE == 1): null_light_results_acc['light_event_id'].append(cp.full(1, ievd)) # one event @@ -900,10 +958,18 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o RangePush("track_pixel_map") # Mapping between unique pixel array and track array index - track_pixel_map = cp.full((unique_pix.shape[0], detsim.MAX_TRACKS_PER_PIXEL), -1) + max_segments_to_trace = max(assmap_pix2seg.shape[1],detsim.MAX_TRACKS_PER_PIXEL) + track_pixel_map = cp.full((unique_pix.shape[0], max_segments_to_trace), -1) TPB = 32 + #BPG = max(ceil(unique_pix.shape[0] / TPB),1) BPG = max(ceil(unique_pix.shape[0] / TPB),1) - detsim.get_track_pixel_map[BPG, TPB](track_pixel_map, unique_pix, neighboring_pixels) + detsim.get_track_pixel_map2[BPG, TPB](track_pixel_map, + unique_pix, + #active_pixels, + neighboring_pixels, + neighboring_radius, + neighboring_radius.max().item()+1, + ) RangePop() RangePush("sum_pixels_signals") @@ -915,16 +981,32 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o BPG = (BPG_X, BPG_Y, BPG_Z) pixels_signals = cp.zeros((len(unique_pix), len(detector.TIME_TICKS))) pixels_tracks_signals = cp.zeros((len(unique_pix), - len(detector.TIME_TICKS), - track_pixel_map.shape[1])) + len(detector.TIME_TICKS), + track_pixel_map.shape[1],) + ) detsim.sum_pixel_signals[BPG,TPB](pixels_signals, signals, track_starts, pixel_index_map, track_pixel_map, - pixels_tracks_signals) + pixels_tracks_signals, + ) RangePop() + #np.savez('wf.npz', + # actives=cp.asnumpy(active_pixels), + # signals=cp.asnumpy(signals), + # pix2idx=cp.asnumpy(pixel_index_map), + # trk2pix=cp.asnumpy(track_pixel_map), + # neighbors=cp.asnumpy(neighboring_pixels), + # distances=cp.asnumpy(neighboring_radius), + # uniquepix=cp.asnumpy(unique_pix), + # track_starts=cp.asnumpy(track_starts/detector.TIME_SAMPLING), + # pixels_tracks_signals=cp.asnumpy(pixels_tracks_signals), + # assmap_pix2seg=cp.asnumpy(assmap_pix2seg), + # ) + + RangePush("get_adc_values") # Here we simulate the electronics response (the self-triggering cycle) and the signal digitization time_ticks = cp.linspace(0, len(unique_eventIDs) * detector.TIME_INTERVAL[1], pixels_signals.shape[1]+1) @@ -948,7 +1030,6 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o rng_states, current_fractions, pixel_thresholds) - # get list of adc values if pixel_gains_file is not None: pixel_gains = cp.array(pixel_gains_lut[unique_pix.ravel()]) @@ -965,10 +1046,15 @@ def save_results(event_times, is_first_batch, results, i_trig, i_mod=-1, light_o results_acc['adc_tot_ticks'].append(adc_ticks_list) results_acc['unique_pix'].append(unique_pix) results_acc['current_fractions'].append(current_fractions) - #track_pixel_map[track_pixel_map != -1] += first_trk_id + itrk + #results_acc['segidx2trackid'].append(segidx2trackid) + traj_pixel_map = cp.full(track_pixel_map.shape,-1) + traj_pixel_map[:] = track_pixel_map + traj_pixel_map[traj_pixel_map != -1] = trajectory_ids_arr[batch_mask][traj_pixel_map[traj_pixel_map != -1] + itrk] track_pixel_map[track_pixel_map != -1] = segment_ids_arr[batch_mask][track_pixel_map[track_pixel_map != -1] + itrk] + results_acc['traj_pixel_map'].append(traj_pixel_map) results_acc['track_pixel_map'].append(track_pixel_map) + # ~~~ Light detector response simulation ~~~ if light.LIGHT_SIMULATED: RangePush("sum_light_signals") diff --git a/larndsim/detsim.py b/larndsim/detsim.py index 6d11739c..2dfd4e9c 100644 --- a/larndsim/detsim.py +++ b/larndsim/detsim.py @@ -4,7 +4,7 @@ """ from math import pi, ceil, sqrt, erf, exp, log - +import cupy as cp import numba as nb from numba import cuda @@ -14,7 +14,7 @@ from .consts import detector from .pixels_from_track import id2pixel -MAX_TRACKS_PER_PIXEL = 20 +MAX_TRACKS_PER_PIXEL = 40 MIN_STEP_SIZE = 0.001 # cm MC_SAMPLE_MULTIPLIER = 1 @@ -469,7 +469,13 @@ def sign(x): return 1 if x >= 0 else -1 @cuda.jit -def sum_pixel_signals(pixels_signals, signals, track_starts, pixel_index_map, track_pixel_map, pixels_tracks_signals): +def sum_pixel_signals(pixels_signals, + signals, + track_starts, + pixel_index_map, + track_pixel_map, + pixels_tracks_signals, + ): """ This function sums the induced current signals on the same pixel. @@ -489,20 +495,27 @@ def sum_pixel_signals(pixels_signals, signals, track_starts, pixel_index_map, tr pixels_tracks_signals (:obj:`numpy.ndarray`): 3D array that will contain the waveforms for each pixel and each track that induced current on the pixel. """ - itrk, ipix, itick = cuda.grid(3) + # itrk = segment index in signals collection, same as the index in all-segments (axis 0) + # ipix = pixel index within signals collection (axis 1) + # track_idx = segment index in track_pixel_map collection (axis 1) + # counter = segment index in pixels_tracks_signals collection (axis 2), same as track_idx ordering + + itrk, ipix, itick = cuda.grid(3) if itrk < signals.shape[0] and ipix < signals.shape[1]: pixel_index = pixel_index_map[itrk][ipix] start_tick = round(track_starts[itrk] / detector.TIME_SAMPLING) if pixel_index >= 0: - counter = 0 + counter = -1 for track_idx in range(track_pixel_map[pixel_index].shape[0]): if itrk == -1: continue if itrk == int(track_pixel_map[pixel_index][track_idx]): counter = track_idx + #if counter == 0: + # print(pixel_index,track_idx,itrk) break if itick < signals.shape[2]: @@ -511,12 +524,16 @@ def sum_pixel_signals(pixels_signals, signals, track_starts, pixel_index_map, tr cuda.atomic.add(pixels_signals, (pixel_index, itime), signals[itrk][ipix][itick]) - cuda.atomic.add(pixels_tracks_signals, - (pixel_index, itime, counter), - signals[itrk][ipix][itick]) + #if (counter==0) and (pixel_index==64): + # print(itrk,ipix,itick,signals[itrk][ipix][itick],pixels_tracks_signals[pixel_index][itime][counter]) + + if (counter >= 0): + cuda.atomic.add(pixels_tracks_signals, + (pixel_index, itime, counter), + signals[itrk][ipix][itick]) @cuda.jit -def get_track_pixel_map(track_pixel_map, unique_pix, pixels): +def get_track_pixel_map(track_pixel_map, unique_pix, pixels, distances): """ This kernel fills a 2D array which contains, for each unique pixel, an array with the track indeces associated to that pixel. @@ -546,3 +563,48 @@ def get_track_pixel_map(track_pixel_map, unique_pix, pixels): if imap < track_pixel_map.shape[1]: track_pixel_map[index][imap] = itrk + +@cuda.jit +def get_track_pixel_map2(track_pixel_map, unique_pix, pixels, distances, max_distance): + """ + This kernel fills a 2D array which contains, for each unique pixel, + an array with the track indeces associated to that pixel. + + Args: + track_pixel_map_col (:obj:`numpy.ndarray`): 2D array that will contain the + association between the unique pixels array and the track indeces + unique_pix (:obj:`numpy.ndarray`): 1D array containing the unique pixels + pixels (:obj:`numpy.ndarray`): 2D array containing the pixels for each + track. + """ + # index of unique_pix array + index = cuda.grid(1) + + upix = unique_pix[index] + + for target_dist in range(max_distance): + + for itrk in range(pixels.shape[0]): + + for ipix in range(pixels.shape[1]): + pID = pixels[itrk][ipix] + dist = distances[itrk][ipix] + + if (upix == pID): + if (dist == target_dist): + + imap = 0 + #while imap < track_pixel_map.shape[1] and track_pixel_map[index][imap] != -1 and track_pixel_map[index][imap] != itrk: + while imap < track_pixel_map.shape[1]: + if track_pixel_map[index][imap] == itrk: + imap = -1 + break + if track_pixel_map[index][imap] == -1: + break + else: + imap += 1 + + if (imap >= 0) and (imap < track_pixel_map.shape[1]): + track_pixel_map[index][imap] = itrk + + break \ No newline at end of file diff --git a/larndsim/fee.py b/larndsim/fee.py index e1facfa8..581e3a51 100644 --- a/larndsim/fee.py +++ b/larndsim/fee.py @@ -23,7 +23,7 @@ from .consts import units #: Number of back-tracked segments to be recorded -ASSOCIATION_COUNT_TO_STORE = 10 +ASSOCIATION_COUNT_TO_STORE = 20 #: Maximum number of ADC values stored per pixel MAX_ADC_VALUES = 20 #: Discrimination threshold in e- @@ -134,6 +134,8 @@ def export_to_hdf5(event_id_list, unique_pix, current_fractions, track_ids, + traj_ids, + #segidx2trackid, filename, event_start_times, is_first_batch, @@ -171,6 +173,7 @@ def export_to_hdf5(event_id_list, packets = [] packets_mc_evt = [] packets_mc_trk = [] + packets_mc_trj = [] packets_frac = [] #if is_first_batch: @@ -252,11 +255,13 @@ def export_to_hdf5(event_id_list, packets[-1].chip_key = Key(io_group,0,0) packets_mc_evt.append([-1]) packets_mc_trk.append([-1] * track_ids.shape[1]) + packets_mc_trj.append([-1] * traj_ids.shape[1]) packets_frac.append([0] * current_fractions.shape[2]) packets.append(SyncPacket(sync_type=b'S', timestamp=time_tick , io_group=io_group)) packets_mc_evt.append([-1]) packets_mc_trk.append([-1] * track_ids.shape[1]) + packets_mc_trj.append([-1] * traj_ids.shape[1]) packets_frac.append([0] * current_fractions.shape[2]) trig_mask = light_trigger_event_id == event @@ -268,6 +273,7 @@ def export_to_hdf5(event_id_list, packets.append(TriggerPacket(io_group=io_group, trigger_type=b'\x02', timestamp=t_trig)) packets_mc_evt.append([-1]) packets_mc_trk.append([-1] * track_ids.shape[1]) + packets_mc_trj.append([-1] * traj_ids.shape[1]) packets_frac.append([0] * current_fractions.shape[2]) # redundant here elif light.LIGHT_TRIG_MODE == 1: @@ -276,6 +282,7 @@ def export_to_hdf5(event_id_list, packets.append(TriggerPacket(io_group=io_group, trigger_type=b'\x02', timestamp=t_trig)) packets_mc_evt.append([-1]) packets_mc_trk.append([-1] * track_ids.shape[1]) + packets_mc_trj.append([-1] * traj_ids.shape[1]) packets_frac.append([0] * current_fractions.shape[2]) last_event = event @@ -307,7 +314,7 @@ def export_to_hdf5(event_id_list, if channel in bad_channels_list[chip_key]: logger.info(f"Channel {channel} on chip {chip_key} disabled") continue - + p.pixel_id = pixel_id p.chip_key = chip_key p.channel_id = channel p.receipt_timestamp = time_tick @@ -324,9 +331,11 @@ def export_to_hdf5(event_id_list, packets[-1].chip_key = Key(io_group,0,0) packets_mc_evt.append([-1]) packets_mc_trk.append([-1] * (ASSOCIATION_COUNT_TO_STORE * 2)) + packets_mc_trj.append([-1] * (ASSOCIATION_COUNT_TO_STORE * 2)) packets_frac.append([0] * (ASSOCIATION_COUNT_TO_STORE*2)) packets_mc_evt.append([event]) packets_mc_trk.append(track_ids[itick]) + packets_mc_trj.append(traj_ids[itick]) packets_frac.append(current_fractions[itick][iadc]) packets.append(p) @@ -339,25 +348,30 @@ def export_to_hdf5(event_id_list, hdf5format.to_file(filename, packet_list, workers=1) dtype = np.dtype([('event_ids',f'(1,)i8'), ('segment_ids',f'({ASSOCIATION_COUNT_TO_STORE},)i8'), - ('fraction', f'({ASSOCIATION_COUNT_TO_STORE},)f8')]) + ('fraction', f'({ASSOCIATION_COUNT_TO_STORE},)f8'), + ('file_traj_ids',f'({ASSOCIATION_COUNT_TO_STORE},)i8'), + ('fraction_traj',f'({ASSOCIATION_COUNT_TO_STORE},)f8'),] + ) packets_mc_ds = np.empty(len(packets), dtype=dtype) # First, sort the back-tracking information by the magnitude of the fraction packets_frac = np.array(packets_frac) - packets_mc_trk = np.array(packets_mc_trk) - packets_mc_evt = np.array(packets_mc_evt) + packets_mc_trk = np.array(packets_mc_trk) + packets_mc_trj = np.array(packets_mc_trj) + packets_mc_evt = np.array(packets_mc_evt) - frac_order = np.flip(np.argsort(np.abs(packets_frac),axis=1),axis=1) - ass_track_ids = np.take_along_axis(packets_mc_trk, frac_order, axis=1) + frac_order = np.flip(np.argsort(packets_frac,axis=1),axis=1) + ass_segment_ids = np.take_along_axis(packets_mc_trk, frac_order, axis=1) + ass_trajectory_ids = np.take_along_axis(packets_mc_trj, frac_order, axis=1) ass_fractions = np.take_along_axis(packets_frac, frac_order, axis=1) # Second, only store the relevant portion. - if ass_track_ids.shape[1] >= ASSOCIATION_COUNT_TO_STORE: - packets_mc_ds['segment_ids'] = ass_track_ids[:,:ASSOCIATION_COUNT_TO_STORE] + if ass_segment_ids.shape[1] >= ASSOCIATION_COUNT_TO_STORE: + packets_mc_ds['segment_ids'] = ass_segment_ids[:,:ASSOCIATION_COUNT_TO_STORE] packets_mc_ds['fraction' ] = ass_fractions[:,:ASSOCIATION_COUNT_TO_STORE] else: - num_to_pad = ASSOCIATION_COUNT_TO_STORE - ass_track_ids.shape[1] - packets_mc_ds['segment_ids'] = np.pad(ass_track_ids, + num_to_pad = ASSOCIATION_COUNT_TO_STORE - ass_segment_ids.shape[1] + packets_mc_ds['segment_ids'] = np.pad(ass_segment_ids, pad_width=((0,0),(0,num_to_pad)), mode='constant', constant_values=-1) @@ -365,6 +379,32 @@ def export_to_hdf5(event_id_list, pad_width=((0,0),(0,num_to_pad)), mode='constant', constant_values=0.) + + + #print(segidx2trackid.shape,segidx2trackid.min(),segidx2trackid.max()) + # fractions for file_traj_ids + ass_track_ids = np.full(ass_trajectory_ids.shape,fill_value=-1,dtype=np.int32) + ass_fractions_track = np.full(ass_fractions.shape,fill_value=0.,dtype=np.float32) + for pidx, tids in enumerate(ass_trajectory_ids): + mask = tids > -1 + for tidx, unique_tid in enumerate(np.unique(tids[mask])): + ass_track_ids[pidx][tidx] = unique_tid + ass_fractions_track[pidx][tidx] = np.sum(ass_fractions[pidx][mask][tids[mask]==unique_tid]) + + if ass_segment_ids.shape[1] >= ASSOCIATION_COUNT_TO_STORE: + packets_mc_ds['file_traj_ids'] = ass_track_ids[:,:ASSOCIATION_COUNT_TO_STORE] + packets_mc_ds['fraction_traj'] = ass_fractions_track[:,:ASSOCIATION_COUNT_TO_STORE] + else: + num_to_pad = ASSOCIATION_COUNT_TO_STORE - ass_track_ids.shape[1] + packets_mc_ds['file_traj_ids'] = np.pad(ass_track_ids, + pad_width=((0,0),(0,num_to_pad)), + mode='constant', + constant_values=-1) + packets_mc_ds['fraction_traj' ] = np.pad(ass_fractions_track, + pad_width=((0,0),(0,num_to_pad)), + mode='constant', + constant_values=0.) + packets_mc_ds['event_ids'] = packets_mc_evt with h5py.File(filename, 'a') as f: diff --git a/larndsim/pixels_from_track.py b/larndsim/pixels_from_track.py index cdff25c0..e5546067 100644 --- a/larndsim/pixels_from_track.py +++ b/larndsim/pixels_from_track.py @@ -8,6 +8,8 @@ #from .consts.detector import PIXEL_PITCH, N_PIXELS, TPC_BORDERS from .consts import detector +MAX_NEIGHBOR_BACKTRACK_DISTANCE=4 + @nb.njit def pixel2id(pixel_x, pixel_y, pixel_plane): """ @@ -63,7 +65,7 @@ def max_pixels(tracks, n_max_pixels): cuda.atomic.max(n_max_pixels, 0, n_active_pixels) @cuda.jit -def get_pixels(tracks, active_pixels, neighboring_pixels, n_pixels_list, radius): +def get_pixels(tracks, active_pixels, neighboring_pixels, neighboring_radius, n_pixels_list, radius): """ For all tracks, takes the xy start and end position and calculates all impacted pixels by the track segment @@ -77,6 +79,9 @@ def get_pixels(tracks, active_pixels, neighboring_pixels, n_pixels_list, radius) neighboring_pixels (:obj:`numpy.ndarray`): array where we store the IDs of the pixels directly below the projection of the segments and the ones next to them + neighboring_radius (:obj:`numpy.ndarray`): array where we store + the distance to the pixels directly below the projection of + the segments n_pixels_list (:obj:`numpy.ndarray`): number of total involved pixels radius (int): number of pixels around the active pixels that @@ -100,7 +105,8 @@ def get_pixels(tracks, active_pixels, neighboring_pixels, n_pixels_list, radius) t["pixel_plane"], active_pixels[itrk]) n_pixels_list[itrk] = get_neighboring_pixels(active_pixels[itrk], radius, - neighboring_pixels[itrk]) + neighboring_pixels[itrk], + neighboring_radius[itrk]) @nb.njit def get_num_active_pixels(x0, y0, x1, y1, plane_id): @@ -193,7 +199,7 @@ def get_active_pixels(x0, y0, x1, y1, plane_id, tot_pixels): tot_pixels[i] = pixel2id(x0, y0, plane_id) @cuda.jit(device=True) -def get_neighboring_pixels(active_pixels, radius, neighboring_pixels): +def get_neighboring_pixels(active_pixels, radius, neighboring_pixels, neighboring_radius): """ For each active_pixel, it includes all neighboring pixels within a specified radius @@ -207,6 +213,9 @@ def get_neighboring_pixels(active_pixels, radius, neighboring_pixels): neighboring_pixels (:obj:`numpy.ndarray`): array where we store the IDs of the pixels directly below the projection of the segments and the ones next to them + neighboring_radius (:obj:`numpy.ndarray`): array where we store + the distance to the pixels directly below the projection of + the segments Returns: int: number of total involved pixels @@ -234,6 +243,31 @@ def get_neighboring_pixels(active_pixels, radius, neighboring_pixels): if is_unique: neighboring_pixels[count] = new_pixel + dist=pow(x_r**2+y_r**2,0.5) + + dx,dy = abs(x_r),abs(y_r) + dmax,dmin = max(dx,dy),min(dx,dy) + dsum = dmax+dmin + if dsum > MAX_NEIGHBOR_BACKTRACK_DISTANCE: + dist = -1 + elif dsum <=1: + dist = dsum + elif dsum == 2: + dist = 2 if dmax==1 else 3 + elif dsum == 3: + dist = 4 if dmax==2 else 5 + elif dsum == 4: + if dmax == 2: + dist = 6 + elif dmax == 3: + dist = 7 + elif dmax == 4: + dist = 8 + else: + print('Unsupported dsum',dsum) + dist=-1 + #raise ValueError(f'Found neighbor distance larger than set threshold {MAX_NEIGHBOR_BACKGRACK_DISTANCE}') + neighboring_radius[count] = dist count += 1 return count