Save to file after each event

cli/simulate_pixels.py

@@ -61,6 +61,7 @@ def swap_coordinates(tracks):
     return tracks
 
 def maybe_create_rng_states(n, seed=0, rng_states=None):
+    """Create or extend random states for CUDA kernel"""
     if rng_states is None:
         return create_xoroshiro128p_states(n, seed=seed)
     elif n > len(rng_states):
@@ -153,7 +154,8 @@ def run_simulation(input_filename,
 
     # Makes an empty array to store data from lightlut
     if light.LIGHT_SIMULATED:
-        light_sim_dat = np.zeros([len(tracks), light.N_OP_CHANNEL*2], dtype=[('n_photons_det','f4'),('t0_det','f4')])
+        light_sim_dat = np.zeros([len(tracks), light.N_OP_CHANNEL*2],
+                                 dtype=[('n_photons_det','f4'),('t0_det','f4')])
 
     if tracks.size == 0:
         print("Empty input dataset, exiting")
@@ -180,35 +182,35 @@ def run_simulation(input_filename,
     print("******************\nRUNNING SIMULATION\n******************")
     # We calculate the number of electrons after recombination (quenching module)
     # and the position and number of electrons after drifting (drifting module)
-    print("Quenching electrons...",end='')
+    print("Quenching electrons..." , end="")
     start_quenching = time()
     quenching.quench[BPG,TPB](tracks, physics.BIRKS)
     end_quenching = time()
     print(f" {end_quenching-start_quenching:.2f} s")
 
-    print("Drifting electrons...",end='')
+    print("Drifting electrons...", end="")
     start_drifting = time()
     drifting.drift[BPG,TPB](tracks)
     end_drifting = time()
     print(f" {end_drifting-start_drifting:.2f} s")
 
     if light.LIGHT_SIMULATED:
-        print("Calculating optical responses...",end='')
-        start_lightLUT = time()
+        print("Calculating optical responses...", end="")
+        start_light_time = time()
         lut = cp.load(light_lut_filename)
         TPB = 256
         BPG = ceil(tracks.shape[0] / TPB)
         lightLUT.calculate_light_incidence[BPG,TPB](tracks, lut, light_sim_dat)
-        end_lightLUT = time()
-        print(f" {end_lightLUT-start_lightLUT:.2f} s")
+        print(f" {time()-start_light_time:.2f} s")
 
-    # initialize lists to collect results from GPU
-    event_id_list = []
-    adc_tot_list = []
-    adc_tot_ticks_list = []
-    track_pixel_map_tot = []
-    unique_pix_tot = []
-    current_fractions_tot = []
+    with h5py.File(output_filename, 'a') as output_file:
+        output_file.create_dataset("tracks", data=tracks)
+        if light.LIGHT_SIMULATED:
+            output_file.create_dataset('light_dat', data=light_sim_dat)
+        if input_has_trajectories:
+            output_file.create_dataset("trajectories", data=trajectories)
+        if input_has_vertices:
+            output_file.create_dataset("vertices", data=vertices)
 
     # create a lookup table that maps between unique event ids and the segments in the file
     tot_evids = np.unique(tracks['eventID'])
@@ -217,14 +219,20 @@ def run_simulation(input_filename,
     end_idx = len(tracks['eventID']) - 1 - rev_idx
 
     # We divide the sample in portions that can be processed by the GPU
-    tracks_batch_runtimes = []
     step = 1
-    tot_events = 0
 
     # pre-allocate some random number states
     rng_states = maybe_create_rng_states(1024*256, seed=0)
+    t0 = 0
     for ievd in tqdm(range(0, tot_evids.shape[0], step), desc='Simulating events...', ncols=80, smoothing=0):
-        start_tracks_batch = time()
+
+        event_id_list = []
+        adc_tot_list = []
+        adc_tot_ticks_list = []
+        track_pixel_map_tot = []
+        unique_pix_tot = []
+        current_fractions_tot = []
+
         first_event = tot_evids[ievd]
         last_event = tot_evids[min(ievd+step, tot_evids.shape[0]-1)]
 
@@ -328,7 +336,9 @@ def run_simulation(input_filename,
             BPG_Z = ceil(signals.shape[2] / TPB[2])
             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]))
+            pixels_tracks_signals = cp.zeros((len(unique_pix),
+                                              len(detector.TIME_TICKS),
+                                              track_pixel_map.shape[1]))
             detsim.sum_pixel_signals[BPG,TPB](pixels_signals,
                                               signals,
                                               track_starts,
@@ -366,46 +376,34 @@ def run_simulation(input_filename,
             RangePop()
 
             event_id_list.append(adc_event_ids)
-            adc_tot_list.append(cp.asnumpy(adc_list))
-            adc_tot_ticks_list.append(cp.asnumpy(adc_ticks_list))
-            unique_pix_tot.append(cp.asnumpy(unique_pix))
-            current_fractions_tot.append(cp.asnumpy(current_fractions))
+            adc_tot_list.append(adc_list)
+            adc_tot_ticks_list.append(adc_ticks_list)
+            unique_pix_tot.append(unique_pix)
+            current_fractions_tot.append(current_fractions)
             track_pixel_map[track_pixel_map != -1] += first_trk_id + itrk
-            track_pixel_map_tot.append(cp.asnumpy(track_pixel_map))
-
-        tot_events += step
-
-        end_tracks_batch = time()
-        tracks_batch_runtimes.append(end_tracks_batch - start_tracks_batch)
-
-    print("*************\nSAVING RESULT\n*************")
-    RangePush("Exporting to HDF5")
-    # Here we export the result in a HDF5 file.
-    event_id_list = np.concatenate(event_id_list, axis=0)
-    adc_tot_list = np.concatenate(adc_tot_list, axis=0)
-    adc_tot_ticks_list = np.concatenate(adc_tot_ticks_list, axis=0)
-    unique_pix_tot = np.concatenate(unique_pix_tot, axis=0)
-    current_fractions_tot = np.concatenate(current_fractions_tot, axis=0)
-    track_pixel_map_tot = np.concatenate(track_pixel_map_tot, axis=0)
-    fee.export_to_hdf5(event_id_list,
-                       adc_tot_list,
-                       adc_tot_ticks_list,
-                       unique_pix_tot,
-                       current_fractions_tot,
-                       track_pixel_map_tot,
-                       output_filename,
-                       bad_channels=bad_channels)
-    RangePop()
+            track_pixel_map_tot.append(track_pixel_map)
+
+        if event_id_list and adc_tot_list:
+            event_id_list_batch = np.concatenate(event_id_list, axis=0)
+            adc_tot_list_batch = np.concatenate(adc_tot_list, axis=0)
+            adc_tot_ticks_list_batch = np.concatenate(adc_tot_ticks_list, axis=0)
+            unique_pix_tot_batch = np.concatenate(unique_pix_tot, axis=0)
+            current_fractions_tot_batch = np.concatenate(current_fractions_tot, axis=0)
+            track_pixel_map_tot_batch = np.concatenate(track_pixel_map_tot, axis=0)
+            _, _, last_time = fee.export_to_hdf5(event_id_list_batch,
+                                                adc_tot_list_batch,
+                                                adc_tot_ticks_list_batch,
+                                                cp.asnumpy(unique_pix_tot_batch),
+                                                cp.asnumpy(current_fractions_tot_batch),
+                                                cp.asnumpy(track_pixel_map_tot_batch),
+                                                output_filename,
+                                                t0=t0,
+                                                bad_channels=bad_channels)
+            t0 = last_time
 
     with h5py.File(output_filename, 'a') as output_file:
-        output_file.create_dataset("tracks", data=tracks)
-        if light.LIGHT_SIMULATED:
-            output_file.create_dataset('light_dat', data=light_sim_dat)
-        if input_has_trajectories:
-            output_file.create_dataset("trajectories", data=trajectories)
-        if input_has_vertices:
-            output_file.create_dataset("vertices", data=vertices)
-        output_file['configs'].attrs['pixel_layout'] = pixel_layout
+        if 'configs' in output_file.keys():
+            output_file['configs'].attrs['pixel_layout'] = pixel_layout
 
     print("Output saved in:", output_filename)
 

larndsim/detector_properties/module0.yaml

@@ -1,5 +1,3 @@
-tpc_offsets: # cm
-  - [0, -21.8236, 0]
 temperature: 87.17 # K
 e_field: 0.50 # kV/cm
 lifetime: 2.2e+3 # us
@@ -8,12 +6,16 @@ long_diff: 4.0e-6 # cm * cm / us
 tran_diff: 8.8e-6 # cm * cm / us
 drift_length: 30.27225 # cm
 response_sampling: 0.1 # us
-response_bin_size: 0.04434 # cm
-tile_map: 
+reponse_bin_size: 0.04434 # cm
+time_padding: 190 # us
+time_window: 189.1 # us
+tpc_offsets: # cm
+  - [0, -21.8236, 0]
+tile_map:
   - [[7,5,3,1],[8,6,4,2]]
   - [[16,14,12,10],[15,13,11,9]]
-n_op_channel: 48
-op_channel_efficiency: [1, 1, 1, 1, 1, 1, 14.9, 14.9, 14.9, 14.9, 14.9, 14.9, 1, 1, 1, 1, 1, 1, 14.9, 14.9, 14.9, 14.9, 14.9, 14.9, 1, 1, 1, 1, 1, 1, 14.9, 14.9, 14.9, 14.9, 14.9, 14.9, 1, 1, 1, 1, 1, 1, 14.9, 14.9, 14.9, 14.9, 14.9, 14.9]
-lut_vox_div: [14, 26, 8]
 module_to_io_groups:
   1: [1,2,3,4]
+n_op_channel: 48
+op_channel_efficiency: [1, 1, 1, 1, 1, 1, 14.9, 14.9, 14.9, 14.9, 14.9, 14.9, 1, 1, 1, 1, 1, 1, 14.9, 14.9, 14.9, 14.9, 14.9, 14.9, 1, 1, 1, 1, 1, 1, 14.9, 14.9, 14.9, 14.9, 14.9, 14.9, 1, 1, 1, 1, 1, 1, 14.9, 14.9, 14.9, 14.9, 14.9, 14.9]
+lut_vox_div: [14, 26, 8]

larndsim/fee.py

@@ -14,11 +14,12 @@
 from larpix.packet import Packet_v2, TimestampPacket, TriggerPacket, SyncPacket
 from larpix.packet import PacketCollection
 from larpix.format import hdf5format
-from tqdm import tqdm
 from .consts import detector, physics
 
 from .pixels_from_track import id2pixel
 
+np.random.seed(1)
+
 #: Maximum number of ADC values stored per pixel
 MAX_ADC_VALUES = 10
 #: Discrimination threshold
@@ -90,6 +91,7 @@ def export_to_hdf5(event_id_list,
                    current_fractions,
                    track_ids,
                    filename,
+                   t0=0,
                    bad_channels=None):
     """
     Saves the ADC counts in the LArPix HDF5 format.
@@ -115,13 +117,14 @@ def export_to_hdf5(event_id_list,
     packets_mc = []
     packets_frac = []
 
-    packets.append(TimestampPacket())
-    packets_mc.append([-1] * track_ids.shape[1])
-    packets_frac.append([0] * current_fractions.shape[2])
-    for io_group in io_groups:
-        packets.append(SyncPacket(sync_type=b'S', timestamp=0, io_group=io_group))
+    if t0 == 0:
+        packets.append(TimestampPacket())
         packets_mc.append([-1] * track_ids.shape[1])
         packets_frac.append([0] * current_fractions.shape[2])
+        for io_group in io_groups:
+            packets.append(SyncPacket(sync_type=b'S', timestamp=0, io_group=io_group))
+            packets_mc.append([-1] * track_ids.shape[1])
+            packets_frac.append([0] * current_fractions.shape[2])
 
     packets_mc_ds = []
     last_event = -1
@@ -133,14 +136,21 @@ def export_to_hdf5(event_id_list,
     unique_events, unique_events_inv = np.unique(event_id_list[...,0], return_inverse=True)
     event_start_time = np.random.exponential(scale=EVENT_RATE, size=unique_events.shape).astype(int)
     event_start_time = np.cumsum(event_start_time)
+    event_start_time += t0
     event_start_time_list = event_start_time[unique_events_inv]
+
     rollover_count = 0
-    for itick, adcs in enumerate(tqdm(adc_list, desc="Writing to HDF5...", ncols=80)):
+    for itick, adcs in enumerate(adc_list):
         ts = adc_ticks_list[itick]
         pixel_id = unique_pix[itick]
 
         pix_x, pix_y, plane_id = id2pixel(pixel_id)
         module_id = plane_id//2+1
+
+        if module_id not in detector.MODULE_TO_IO_GROUPS.keys():
+            logger.warning("Pixel ID not valid %i" % module_id)
+            continue
+
         tile_x = int(pix_x//detector.N_PIXELS_PER_TILE[0])
         tile_y = int(pix_y//detector.N_PIXELS_PER_TILE[1])
         anode_id = 0 if plane_id % 2 == 0 else 1
@@ -154,7 +164,6 @@ def export_to_hdf5(event_id_list,
                     event = event_id_list[itick,iadc]
                     event_t0 = event_start_time_list[itick]
                     time_tick = int(np.floor(t / CLOCK_CYCLE + event_t0))
-
                     if event_t0 > 2**31 - 1 or time_tick > 2**31 - 1:
                         # 31-bit rollover
                         rollover_count += 1
@@ -187,7 +196,7 @@ def export_to_hdf5(event_id_list,
                                                                                 pix_y % detector.N_PIXELS_PER_TILE[1]),
                                                                    tile_id)]
                 except KeyError:
-                    logger.warning("Pixel ID not valid", pixel_id)
+                    logger.warning("Pixel ID not valid %i" % pixel_id)
                     continue
 
                 p.dataword = int(adc)
@@ -222,27 +231,27 @@ def export_to_hdf5(event_id_list,
             else:
                 break
 
-    packet_list = PacketCollection(packets, read_id=0, message='')
-
-    hdf5format.to_file(filename, packet_list)
-
     if packets:
+        packet_list = PacketCollection(packets, read_id=0, message='')
+        hdf5format.to_file(filename, packet_list)
         packets_mc_ds = np.empty(len(packets), dtype=dtype)
         packets_mc_ds['track_ids'] = packets_mc
         packets_mc_ds['fraction'] = packets_frac
 
-    with h5py.File(filename, 'a') as f:
-        if "mc_packets_assn" in f.keys():
-            del f['mc_packets_assn']
-        f.create_dataset("mc_packets_assn", data=packets_mc_ds)
+        with h5py.File(filename, 'a') as f:
+            if t0 == 0:
+                f.create_dataset("mc_packets_assn", data=packets_mc_ds, maxshape=(None,))
+            else:
+                f['mc_packets_assn'].resize((f['mc_packets_assn'].shape[0] + packets_mc_ds.shape[0]), axis=0)
+                f['mc_packets_assn'][-packets_mc_ds.shape[0]:] = packets_mc_ds
 
-        f['configs'].attrs['vdrift'] = detector.V_DRIFT
-        f['configs'].attrs['long_diff'] = detector.LONG_DIFF
-        f['configs'].attrs['tran_diff'] = detector.TRAN_DIFF
-        f['configs'].attrs['lifetime'] = detector.ELECTRON_LIFETIME
-        f['configs'].attrs['drift_length'] = detector.DRIFT_LENGTH
+            f['configs'].attrs['vdrift'] = detector.V_DRIFT
+            f['configs'].attrs['long_diff'] = detector.LONG_DIFF
+            f['configs'].attrs['tran_diff'] = detector.TRAN_DIFF
+            f['configs'].attrs['lifetime'] = detector.ELECTRON_LIFETIME
+            f['configs'].attrs['drift_length'] = detector.DRIFT_LENGTH
 
-    return packets, packets_mc_ds
+    return packets, packets_mc_ds, event_start_time_list[-1]
 
 
 def digitize(integral_list):
@@ -354,7 +363,7 @@ def get_adc_values(pixels_signals,
                         q += curre[ic] * detector.TIME_SAMPLING
                         for itrk in range(current_fractions.shape[2]):
                             current_fractions[ip][iadc][itrk] += pixels_signals_tracks[ip][ic][itrk] * detector.TIME_SAMPLING
-                            
+
                     q_sum += q
                     ic+=1
 
@@ -388,7 +397,7 @@ def get_adc_values(pixels_signals,
                 ic += round(RESET_CYCLES * CLOCK_CYCLE / detector.TIME_SAMPLING)
                 last_reset = ic
                 adc_busy = round(ADC_BUSY_DELAY * CLOCK_CYCLE / detector.TIME_SAMPLING)
-                
+
                 q_sum = xoroshiro128p_normal_float32(rng_states, ip) * RESET_NOISE_CHARGE * physics.E_CHARGE
 
                 iadc += 1

larndsim/pixels_from_track.py

@@ -84,6 +84,7 @@ def get_pixels(tracks, active_pixels, neighboring_pixels, n_pixels_list, radius)
     itrk = cuda.grid(1)
     if itrk < tracks.shape[0]:
         t = tracks[itrk]
+
         this_border = TPC_BORDERS[int(t["pixel_plane"])]
         start_pixel = (
             int((t["x_start"] - this_border[0][0]) // PIXEL_PITCH),
@@ -127,7 +128,7 @@ def get_num_active_pixels(x0, y0, x1, y1, plane_id):
     err = dx + dy
 
     n = 0
-    if 0 <= x0 < N_PIXELS[0] and 0 <= y0 < N_PIXELS[1] and plane_id < TPC_BORDERS.shape[0]:
+    if 0 <= x0 < N_PIXELS[0] and 0 <= y0 < N_PIXELS[1] and 0 <= plane_id < TPC_BORDERS.shape[0]:
         n += 1
 
     while x0 != x1 or y0 != y1:
@@ -141,7 +142,7 @@ def get_num_active_pixels(x0, y0, x1, y1, plane_id):
             err += dx
             y0 += sy
 
-        if 0 <= x0 < N_PIXELS[0] and 0 <= y0 < N_PIXELS[1] and plane_id < TPC_BORDERS.shape[0]:
+        if 0 <= x0 < N_PIXELS[0] and 0 <= y0 < N_PIXELS[1] and 0 <= plane_id < TPC_BORDERS.shape[0]:
             n += 1
 
     return n
@@ -171,7 +172,8 @@ def get_active_pixels(x0, y0, x1, y1, plane_id, tot_pixels):
     err = dx + dy
 
     i = 0
-    if 0 <= x0 < N_PIXELS[0] and 0 <= y0 < N_PIXELS[1] and plane_id < TPC_BORDERS.shape[0]:
+
+    if 0 <= x0 < N_PIXELS[0] and 0 <= y0 < N_PIXELS[1] and 0 <= plane_id < TPC_BORDERS.shape[0]:
         tot_pixels[i] = pixel2id(x0, y0, plane_id)
 
     while x0 != x1 or y0 != y1:
@@ -186,7 +188,7 @@ def get_active_pixels(x0, y0, x1, y1, plane_id, tot_pixels):
             err += dx
             y0 += sy
 
-        if 0 <= x0 < N_PIXELS[0] and 0 <= y0 < N_PIXELS[1] and plane_id < TPC_BORDERS.shape[0]:
+        if 0 <= x0 < N_PIXELS[0] and 0 <= y0 < N_PIXELS[1] and 0 <= plane_id < TPC_BORDERS.shape[0]:
             tot_pixels[i] = pixel2id(x0, y0, plane_id)
 
 @cuda.jit(device=True)
@@ -221,7 +223,7 @@ def get_neighboring_pixels(active_pixels, radius, neighboring_pixels):
                 new_x, new_y = active_x + x_r, active_y + y_r
                 is_unique = True
 
-                if 0 <= new_x < N_PIXELS[0] and 0 <= new_y < N_PIXELS[1] and plane_id < TPC_BORDERS.shape[0]:
+                if 0 <= new_x < N_PIXELS[0] and 0 <= new_y < N_PIXELS[1] and 0 <= plane_id < TPC_BORDERS.shape[0]:
                     new_pixel = pixel2id(new_x, new_y, plane_id)
 
                     for ipix in range(neighboring_pixels.shape[0]):