Functions | Variables
simulate_pixels Namespace Reference

Functions

def swap_coordinates (tracks)
 
def maybe_create_rng_states (n, seed=0, rng_states=None)
 
def run_simulation (input_filename, pixel_layout, detector_properties, output_filename, response_file='../larndsim/bin/response_44.npy', light_lut_filename='../larndsim/bin/lightLUT.npy', bad_channels=None, n_tracks=None, pixel_thresholds_file=None)
 

Variables

int BATCH_SIZE = 4000
 
string LOGO
 

Detailed Description

Command-line interface to larnd-sim module.

Function Documentation

def simulate_pixels.maybe_create_rng_states (   n,
  seed = 0,
  rng_states = None 
)

Definition at line 63 of file simulate_pixels.py.

63 def maybe_create_rng_states(n, seed=0, rng_states=None):
64  if rng_states is None:
65  return create_xoroshiro128p_states(n, seed=seed)
66  elif n > len(rng_states):
67  new_states = device_array(n, dtype=rng_states.dtype)
68  new_states[:len(rng_states)] = rng_states
69  new_states[len(rng_states):] = create_xoroshiro128p_states(n - len(rng_states), seed=seed)
70  return new_states
71  return rng_states
72 
def maybe_create_rng_states(n, seed=0, rng_states=None)
def simulate_pixels.run_simulation (   input_filename,
  pixel_layout,
  detector_properties,
  output_filename,
  response_file = '../larndsim/bin/response_44.npy',
  light_lut_filename = '../larndsim/bin/lightLUT.npy',
  bad_channels = None,
  n_tracks = None,
  pixel_thresholds_file = None 
)
Command-line interface to run the simulation of a pixelated LArTPC

Args:
    input_filename (str): path of the edep-sim input file
    pixel_layout (str): path of the YAML file containing the pixel
        layout and connection details.
    detector_properties (str): path of the YAML file containing
        the detector properties
    output_filename (str): path of the HDF5 output file. If not specified
        the output is added to the input file.
    response_file (str, optional): path of the Numpy array containing the pre-calculated
        field responses. Defaults to ../larndsim/bin/response_44.npy.
    light_lut_file (str, optional): path of the Numpy array containing the light
        look-up table. Defaults to ../larndsim/bin/lightLUT.npy.
    bad_channels (str, optional): path of the YAML file containing the channels to be
        disabled. Defaults to None
    n_tracks (int, optional): number of tracks to be simulated. Defaults to None
        (all tracks).
    pixel_thresholds_file (str): path to npz file containing pixel thresholds. Defaults
        to None.

Definition at line 81 of file simulate_pixels.py.

81  pixel_thresholds_file=None):
82  """
83  Command-line interface to run the simulation of a pixelated LArTPC
84 
85  Args:
86  input_filename (str): path of the edep-sim input file
87  pixel_layout (str): path of the YAML file containing the pixel
88  layout and connection details.
89  detector_properties (str): path of the YAML file containing
90  the detector properties
91  output_filename (str): path of the HDF5 output file. If not specified
92  the output is added to the input file.
93  response_file (str, optional): path of the Numpy array containing the pre-calculated
94  field responses. Defaults to ../larndsim/bin/response_44.npy.
95  light_lut_file (str, optional): path of the Numpy array containing the light
96  look-up table. Defaults to ../larndsim/bin/lightLUT.npy.
97  bad_channels (str, optional): path of the YAML file containing the channels to be
98  disabled. Defaults to None
99  n_tracks (int, optional): number of tracks to be simulated. Defaults to None
100  (all tracks).
101  pixel_thresholds_file (str): path to npz file containing pixel thresholds. Defaults
102  to None.
103  """
104  start_simulation = time()
105 
106  RangePush("run_simulation")
107 
108  print(LOGO)
109  print("**************************\nLOADING SETTINGS AND INPUT\n**************************")
110  print("Pixel layout file:", pixel_layout)
111  print("Detector propeties file:", detector_properties)
112  print("edep-sim input file:", input_filename)
113  print("Response file:", response_file)
114  if bad_channels:
115  print("Disabled channel list: ", bad_channels)
116  RangePush("load_detector_properties")
117  consts.load_properties(detector_properties, pixel_layout)
118  from larndsim.consts import light, detector, physics
119  RangePop()
120 
121  RangePush("load_larndsim_modules")
122  # Here we load the modules after loading the detector properties
123  # maybe can be implemented in a better way?
124  from larndsim import quenching, drifting, detsim, pixels_from_track, fee, lightLUT
125  RangePop()
126 
127  RangePush("load_pixel_thresholds")
128  if pixel_thresholds_file is not None:
129  print("Pixel thresholds file:", pixel_thresholds_file)
130  pixel_thresholds_lut = CudaDict.load(pixel_thresholds_file, 256)
131  else:
132  pixel_thresholds_lut = CudaDict(cp.array([fee.DISCRIMINATION_THRESHOLD]), 1, 1)
133  RangePop()
134 
135  RangePush("load_hd5_file")
136  # First of all we load the edep-sim output
137  with h5py.File(input_filename, 'r') as f:
138  tracks = np.array(f['segments'])
139  try:
140  trajectories = np.array(f['trajectories'])
141  input_has_trajectories = True
142  except KeyError:
143  input_has_trajectories = False
144 
145  try:
146  vertices = np.array(f['eventTruth'])
147  input_has_vertices = True
148  except KeyError:
149  print("Input file does not have true vertices info")
150  input_has_vertices = False
151 
152  RangePop()
153 
154  # Makes an empty array to store data from lightlut
155  if light.LIGHT_SIMULATED:
156  light_sim_dat = np.zeros([len(tracks), light.N_OP_CHANNEL*2], dtype=[('n_photons_det','f4'),('t0_det','f4')])
157 
158  if tracks.size == 0:
159  print("Empty input dataset, exiting")
160  return
161 
162  if n_tracks:
163  tracks = tracks[:n_tracks]
164  if light.LIGHT_SIMULATED:
165  light_sim_dat = light_sim_dat[:n_tracks]
166 
167  if 'n_photons' not in tracks.dtype.names:
168  n_photons = np.zeros(tracks.shape[0], dtype=[('n_photons', 'f4')])
169  tracks = rfn.merge_arrays((tracks, n_photons), flatten=True)
170 
171  # Here we swap the x and z coordinates of the tracks
172  # because of the different convention in larnd-sim wrt edep-sim
173  tracks = swap_coordinates(tracks)
174 
175  response = cp.load(response_file)
176 
177  TPB = 256
178  BPG = ceil(tracks.shape[0] / TPB)
179 
180  print("******************\nRUNNING SIMULATION\n******************")
181  # We calculate the number of electrons after recombination (quenching module)
182  # and the position and number of electrons after drifting (drifting module)
183  print("Quenching electrons...",end='')
184  start_quenching = time()
185  quenching.quench[BPG,TPB](tracks, physics.BIRKS)
186  end_quenching = time()
187  print(f" {end_quenching-start_quenching:.2f} s")
188 
189  print("Drifting electrons...",end='')
190  start_drifting = time()
191  drifting.drift[BPG,TPB](tracks)
192  end_drifting = time()
193  print(f" {end_drifting-start_drifting:.2f} s")
194 
195  if light.LIGHT_SIMULATED:
196  print("Calculating optical responses...",end='')
197  start_lightLUT = time()
198  lut = cp.load(light_lut_filename)
199  TPB = 256
200  BPG = ceil(tracks.shape[0] / TPB)
201  lightLUT.calculate_light_incidence[BPG,TPB](tracks, lut, light_sim_dat)
202  end_lightLUT = time()
203  print(f" {end_lightLUT-start_lightLUT:.2f} s")
204 
205  # initialize lists to collect results from GPU
206  event_id_list = []
207  adc_tot_list = []
208  adc_tot_ticks_list = []
209  track_pixel_map_tot = []
210  unique_pix_tot = []
211  current_fractions_tot = []
212 
213  # create a lookup table that maps between unique event ids and the segments in the file
214  tot_evids = np.unique(tracks['eventID'])
215  _, _, start_idx = np.intersect1d(tot_evids, tracks['eventID'], return_indices=True)
216  _, _, rev_idx = np.intersect1d(tot_evids, tracks['eventID'][::-1], return_indices=True)
217  end_idx = len(tracks['eventID']) - 1 - rev_idx
218 
219  # We divide the sample in portions that can be processed by the GPU
220  tracks_batch_runtimes = []
221  step = 1
222  tot_events = 0
223 
224  # pre-allocate some random number states
225  rng_states = maybe_create_rng_states(1024*256, seed=0)
226  for ievd in tqdm(range(0, tot_evids.shape[0], step), desc='Simulating events...', ncols=80, smoothing=0):
227  start_tracks_batch = time()
228  first_event = tot_evids[ievd]
229  last_event = tot_evids[min(ievd+step, tot_evids.shape[0]-1)]
230 
231  if first_event == last_event:
232  last_event += 1
233 
234  # load a subset of segments from the file and process those that are from the current event
235  track_subset = tracks[min(start_idx[ievd:ievd + step]):max(end_idx[ievd:ievd + step])+1]
236  evt_tracks = track_subset[(track_subset['eventID'] >= first_event) & (track_subset['eventID'] < last_event)]
237  first_trk_id = np.where(track_subset['eventID'] == evt_tracks['eventID'][0])[0][0] + min(start_idx[ievd:ievd + step])
238 
239  for itrk in tqdm(range(0, evt_tracks.shape[0], BATCH_SIZE), desc=' Simulating event %i batches...' % ievd, leave=False, ncols=80):
240  selected_tracks = evt_tracks[itrk:itrk+BATCH_SIZE]
241  RangePush("event_id_map")
242  event_ids = selected_tracks['eventID']
243  unique_eventIDs = np.unique(event_ids)
244  RangePop()
245 
246  # We find the pixels intersected by the projection of the tracks on
247  # the anode plane using the Bresenham's algorithm. We also take into
248  # account the neighboring pixels, due to the transverse diffusion of the charges.
249  RangePush("pixels_from_track")
250  max_radius = ceil(max(selected_tracks["tran_diff"])*5/detector.PIXEL_PITCH)
251 
252  TPB = 128
253  BPG = ceil(selected_tracks.shape[0] / TPB)
254  max_pixels = np.array([0])
255  pixels_from_track.max_pixels[BPG,TPB](selected_tracks, max_pixels)
256 
257  # This formula tries to estimate the maximum number of pixels which can have
258  # a current induced on them.
259  max_neighboring_pixels = (2*max_radius+1)*max_pixels[0]+(1+2*max_radius)*max_radius*2
260 
261  active_pixels = cp.full((selected_tracks.shape[0], max_pixels[0]), -1, dtype=np.int32)
262  neighboring_pixels = cp.full((selected_tracks.shape[0], max_neighboring_pixels), -1, dtype=np.int32)
263  n_pixels_list = cp.zeros(shape=(selected_tracks.shape[0]))
264 
265  if not active_pixels.shape[1] or not neighboring_pixels.shape[1]:
266  continue
267 
268  pixels_from_track.get_pixels[BPG,TPB](selected_tracks,
269  active_pixels,
270  neighboring_pixels,
271  n_pixels_list,
272  max_radius)
273  RangePop()
274 
275  RangePush("unique_pix")
276  shapes = neighboring_pixels.shape
277  joined = neighboring_pixels.reshape(shapes[0] * shapes[1])
278  unique_pix = cp.unique(joined)
279  unique_pix = unique_pix[(unique_pix != -1)]
280  RangePop()
281 
282  if not unique_pix.shape[0]:
283  continue
284 
285  RangePush("time_intervals")
286  # Here we find the longest signal in time and we store an array with the start in time of each track
287  max_length = cp.array([0])
288  track_starts = cp.empty(selected_tracks.shape[0])
289  detsim.time_intervals[BPG,TPB](track_starts, max_length, selected_tracks)
290  RangePop()
291 
292  RangePush("tracks_current")
293  # Here we calculate the induced current on each pixel
294  signals = cp.zeros((selected_tracks.shape[0],
295  neighboring_pixels.shape[1],
296  cp.asnumpy(max_length)[0]), dtype=np.float32)
297  TPB = (1,1,64)
298  BPG_X = ceil(signals.shape[0] / TPB[0])
299  BPG_Y = ceil(signals.shape[1] / TPB[1])
300  BPG_Z = ceil(signals.shape[2] / TPB[2])
301  BPG = (BPG_X, BPG_Y, BPG_Z)
302  rng_states = maybe_create_rng_states(int(np.prod(TPB[:2]) * np.prod(BPG[:2])), seed=ievd, rng_states=rng_states)
303  detsim.tracks_current_mc[BPG,TPB](signals, neighboring_pixels, selected_tracks, response, rng_states)
304  RangePop()
305 
306  RangePush("pixel_index_map")
307  # Here we create a map between tracks and index in the unique pixel array
308  pixel_index_map = cp.full((selected_tracks.shape[0], neighboring_pixels.shape[1]), -1)
309  for i_ in range(selected_tracks.shape[0]):
310  compare = neighboring_pixels[i_, ..., cp.newaxis] == unique_pix
311  indices = cp.where(compare)
312  pixel_index_map[i_, indices[0]] = indices[1]
313  RangePop()
314 
315  RangePush("track_pixel_map")
316  # Mapping between unique pixel array and track array index
317  track_pixel_map = cp.full((unique_pix.shape[0], detsim.MAX_TRACKS_PER_PIXEL), -1)
318  TPB = 32
319  BPG = ceil(unique_pix.shape[0] / TPB)
320  detsim.get_track_pixel_map[BPG, TPB](track_pixel_map, unique_pix, neighboring_pixels)
321  RangePop()
322 
323  RangePush("sum_pixels_signals")
324  # Here we combine the induced current on the same pixels by different tracks
325  TPB = (8,8,8)
326  BPG_X = ceil(signals.shape[0] / TPB[0])
327  BPG_Y = ceil(signals.shape[1] / TPB[1])
328  BPG_Z = ceil(signals.shape[2] / TPB[2])
329  BPG = (BPG_X, BPG_Y, BPG_Z)
330  pixels_signals = cp.zeros((len(unique_pix), len(detector.TIME_TICKS)))
331  pixels_tracks_signals = cp.zeros((len(unique_pix),len(detector.TIME_TICKS),track_pixel_map.shape[1]))
332  detsim.sum_pixel_signals[BPG,TPB](pixels_signals,
333  signals,
334  track_starts,
335  pixel_index_map,
336  track_pixel_map,
337  pixels_tracks_signals)
338  RangePop()
339 
340  RangePush("get_adc_values")
341  # Here we simulate the electronics response (the self-triggering cycle) and the signal digitization
342  time_ticks = cp.linspace(0, len(unique_eventIDs) * detector.TIME_INTERVAL[1], pixels_signals.shape[1]+1)
343  integral_list = cp.zeros((pixels_signals.shape[0], fee.MAX_ADC_VALUES))
344  adc_ticks_list = cp.zeros((pixels_signals.shape[0], fee.MAX_ADC_VALUES))
345  current_fractions = cp.zeros((pixels_signals.shape[0], fee.MAX_ADC_VALUES, track_pixel_map.shape[1]))
346 
347  TPB = 128
348  BPG = ceil(pixels_signals.shape[0] / TPB)
349  rng_states = maybe_create_rng_states(int(TPB * BPG), seed=ievd, rng_states=rng_states)
350  pixel_thresholds_lut.tpb = TPB
351  pixel_thresholds_lut.bpg = BPG
352  pixel_thresholds = pixel_thresholds_lut[unique_pix.ravel()].reshape(unique_pix.shape)
353 
354  fee.get_adc_values[BPG, TPB](pixels_signals,
355  pixels_tracks_signals,
356  time_ticks,
357  integral_list,
358  adc_ticks_list,
359  0,
360  rng_states,
361  current_fractions,
362  pixel_thresholds)
363 
364  adc_list = fee.digitize(integral_list)
365  adc_event_ids = np.full(adc_list.shape, unique_eventIDs[0]) # FIXME: only works if looping on a single event
366  RangePop()
367 
368  event_id_list.append(adc_event_ids)
369  adc_tot_list.append(cp.asnumpy(adc_list))
370  adc_tot_ticks_list.append(cp.asnumpy(adc_ticks_list))
371  unique_pix_tot.append(cp.asnumpy(unique_pix))
372  current_fractions_tot.append(cp.asnumpy(current_fractions))
373  track_pixel_map[track_pixel_map != -1] += first_trk_id + itrk
374  track_pixel_map_tot.append(cp.asnumpy(track_pixel_map))
375 
376  tot_events += step
377 
378  end_tracks_batch = time()
379  tracks_batch_runtimes.append(end_tracks_batch - start_tracks_batch)
380 
381  print("*************\nSAVING RESULT\n*************")
382  RangePush("Exporting to HDF5")
383  # Here we export the result in a HDF5 file.
384  event_id_list = np.concatenate(event_id_list, axis=0)
385  adc_tot_list = np.concatenate(adc_tot_list, axis=0)
386  adc_tot_ticks_list = np.concatenate(adc_tot_ticks_list, axis=0)
387  unique_pix_tot = np.concatenate(unique_pix_tot, axis=0)
388  current_fractions_tot = np.concatenate(current_fractions_tot, axis=0)
389  track_pixel_map_tot = np.concatenate(track_pixel_map_tot, axis=0)
390  fee.export_to_hdf5(event_id_list,
391  adc_tot_list,
392  adc_tot_ticks_list,
393  unique_pix_tot,
394  current_fractions_tot,
395  track_pixel_map_tot,
396  output_filename,
397  bad_channels=bad_channels)
398  RangePop()
399 
400  with h5py.File(output_filename, 'a') as output_file:
401  output_file.create_dataset("tracks", data=tracks)
402  if light.LIGHT_SIMULATED:
403  output_file.create_dataset('light_dat', data=light_sim_dat)
404  if input_has_trajectories:
405  output_file.create_dataset("trajectories", data=trajectories)
406  if input_has_vertices:
407  output_file.create_dataset("eventTruth", data=vertices)
408  output_file['configs'].attrs['pixel_layout'] = pixel_layout
409 
410  print("Output saved in:", output_filename)
411 
412  RangePop()
413  end_simulation = time()
414  print(f"Elapsed time: {end_simulation-start_simulation:.2f} s")
415 
def maybe_create_rng_states(n, seed=0, rng_states=None)
static int max(int a, int b)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
def swap_coordinates(tracks)
def simulate_pixels.swap_coordinates (   tracks)
Swap x and z coordinates in tracks.
This is because the convention in larnd-sim is different
from the convention in edep-sim. FIXME.

Args:
    tracks (:obj:`numpy.ndarray`): tracks array.
    
Returns:
    :obj:`numpy.ndarray`: tracks with swapped axes. 

Definition at line 37 of file simulate_pixels.py.

37 def swap_coordinates(tracks):
38  """
39  Swap x and z coordinates in tracks.
40  This is because the convention in larnd-sim is different
41  from the convention in edep-sim. FIXME.
42 
43  Args:
44  tracks (:obj:`numpy.ndarray`): tracks array.
45 
46  Returns:
47  :obj:`numpy.ndarray`: tracks with swapped axes.
48  """
49  x_start = np.copy(tracks['x_start'] )
50  x_end = np.copy(tracks['x_end'])
51  x = np.copy(tracks['x'])
52 
53  tracks['x_start'] = np.copy(tracks['z_start'])
54  tracks['x_end'] = np.copy(tracks['z_end'])
55  tracks['x'] = np.copy(tracks['z'])
56 
57  tracks['z_start'] = x_start
58  tracks['z_end'] = x_end
59  tracks['z'] = x
60 
61  return tracks
62 
def swap_coordinates(tracks)

Variable Documentation

int simulate_pixels.BATCH_SIZE = 4000

Definition at line 25 of file simulate_pixels.py.

string simulate_pixels.LOGO
Initial value:
1 = """
2  _ _ _
3  | | | | (_)
4  | | __ _ _ __ _ __ __| |______ ___ _ _ __ ___
5  | |/ _` | '__| '_ \ / _` |______/ __| | '_ ` _ \\
6  | | (_| | | | | | | (_| | \__ \ | | | | | |
7  |_|\__,_|_| |_| |_|\__,_| |___/_|_| |_| |_|
8 
9 """

Definition at line 27 of file simulate_pixels.py.