Skip to content

Commit

Permalink
Merge pull request #91 from soleti/cleanup
Browse files Browse the repository at this point in the history 
Cleaning up some code
  • Loading branch information
@soleti
soleti authored Aug 25, 2022
2 parents b5c971b + 0d75a6a commit a9477f8
Show file tree
Hide file tree
Showing 6 changed files with 104 additions and 259 deletions.
112 changes: 0 additions & 112 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,115 +61,3 @@ simulate_pixels.py \
The `response_38.npy` is a file containing an array of induced currents for several $(x,y)$ positions on a pixel with a 38 mm pitch. It is calculated externally to `larnd-sim`. Two versions, one with 44 mm pitch and one with 38 mm pitch, are available in the `larndsim/bin` directory.

The output file will contain the datasets described in the [LArPix HDF5 documentation](https://larpix-control.readthedocs.io/en/stable/api/format/hdf5format.html), plus a dataset `tracks` containing the _true_ energy depositions in the detector, and a dataset `mc_packets_assn`, which has a list of indeces corresponding to the true energy deposition associated to each packet.

## Step-by-step simulation

Here we will describe, step-by-step, how we perform the simulation. A full example for the ND-LAr detector is available in the `examples/NDLAr example.ipynb` notebook.

### Quenching and drifting stage

The particles that interact in the TPC ionize the argon atoms. Some of the resulting electrons will immediately recombine with the atoms. This effect is simulated in the `quenching` module.

The remaining electrons travel towards the anode and their spatial distribution is affected by longitudinal and transverse diffusion. The presence of impurities reduces the amount of electrons that reach the anode. These effects are simulated in the `drifting` module.

These two modules modify in place the input array, here called `input_tracks`.

```python
from larndsim import quenching, drifting

TPB = 256
BPG = ceil(tracks.shape[0] / TPB)
quenching.quench[BPG,TPB](input_tracks, consts.birks)
drifting.drift[BPG,TPB](input_tracks)
```

### Pixel simulation stage

Once we have calculated the number and the position of the electrons reaching the anode, we can calculate the current induced on each pixel.
First, we find the pixels interesected by the projection of each track segment on the anode plane using the [Bresenham's line algorithm](https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm) implented in the `pixels_from_track` module. Due to diffusion, we consider also the neighboring pixels.

```python
from larndsim import pixels_from_track
...
pixels_from_track.get_pixels[BPG,TPB](input_tracks,
active_pixels,
neighboring_pixels,
n_pixels_list,
max_radius)
```

Finally, we calculate the current induced on each pixel using the `tracks_current` function in the `detsim` module. The induced current is stored in the `signals` array, which is a three-dimensional array, where the dimensions correspond to the track segment, the pixel, and the time tick, respectively: `signals[0][1][2]` will contain the current induced by the track `0`, for the pixel `1`, at the time tick `2`.

```python
from larndsim import detsim
...
detsim.tracks_current[BPG,TPB](signals,
neighboring_pixels,
input_tracks,
response,
rng_states)
```

Here, `response` is a Numpy array containing a look-up table with a pre-calculated field response. The file valid for Module0 and SingleCube LArPix tiles is availabe at `larndsim/response-44.npy`. For ND-LAr the file is `larndsim/response_38.npy`.

### Accessing the signals

The three-dimensional array can contain more than one signal for each pixel at different times. If we want to plot the full induced signal on the pixel, we need to join the signals corresponding to the same pixel. First, we find the start time of each signal with `time_intervals`:

```python
from larndsim import detsim
...
detsim.time_intervals[blockspergrid,threadsperblock](track_starts,
max_length,
input_tracks)
```

Then, we join them using `sum_pixel_signals`:

```python
from larndsim import detsim
...
detsim.sum_pixel_signals[blockspergrid,threadsperblock](pixels_signals,
signals,
track_starts,
pixel_index_map,
track_pixel_map,
pixels_tracks_signals)
```

### Electronics simulation

Once we have the induced current for each active pixel in the detector we can apply our electronics simulation, which will calculate the ADC values for each pixel:

```python
from larndsim import fee
from numba.cuda.random import create_xoroshiro128p_states
...
fee.get_adc_values[BPG,TPB](pixels_signals,
pixels_tracks_signals,
time_ticks,
integral_list,
adc_ticks_list,
0,
rng_states,
current_fractions,
pixel_thresholds)
```

where the random states `rng_states` are needed for the noise simulation.

### Export

The final output is then exported to the [LArPix HDF5 format](https://larpix-control.readthedocs.io/en/stable/api/format/hdf5format.html) after each event in the input file:

```python
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)
```
15 changes: 10 additions & 5 deletions cli/simulate_pixels.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,13 +68,16 @@ def swap_coordinates(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):

if n > len(rng_states):
new_states = device_array(n, dtype=rng_states.dtype)
new_states[:len(rng_states)] = rng_states
new_states[len(rng_states):] = create_xoroshiro128p_states(n - len(rng_states), seed=seed)
return new_states

return rng_states

def run_simulation(input_filename,
Expand Down Expand Up @@ -273,7 +276,8 @@ def run_simulation(input_filename,
# pre-allocate some random number states
rng_states = maybe_create_rng_states(1024*256, seed=0)
last_time = 0
for ievd in tqdm(range(0, tot_evids.shape[0], step), desc='Simulating events...', ncols=80, smoothing=0):
for ievd in tqdm(range(0, tot_evids.shape[0], step),
desc='Simulating events...', ncols=80, smoothing=0):

first_event = tot_evids[ievd]
last_event = tot_evids[min(ievd+step, tot_evids.shape[0]-1)]
Expand All @@ -288,7 +292,8 @@ def run_simulation(input_filename,
evt_tracks = track_subset[event_mask]
first_trk_id = np.where(track_subset['eventID'] == evt_tracks['eventID'][0])[0][0] + min(start_idx[ievd:ievd + step])

for itrk in tqdm(range(0, evt_tracks.shape[0], BATCH_SIZE), desc=' Simulating event %i batches...' % ievd, leave=False, ncols=80):
for itrk in tqdm(range(0, evt_tracks.shape[0], BATCH_SIZE),
desc=' Simulating event %i batches...' % ievd, leave=False, ncols=80):
selected_tracks = evt_tracks[itrk:itrk+BATCH_SIZE]
RangePush("event_id_map")
event_ids = selected_tracks['eventID']
Expand Down Expand Up @@ -441,7 +446,7 @@ def run_simulation(input_filename,

TPB = (1,64)
BPG = (ceil(light_sample_inc.shape[0] / TPB[0]),
ceil(light_sample_inc.shape[1] / TPB[1]))
ceil(light_sample_inc.shape[1] / TPB[1]))
light_sim.sum_light_signals[BPG, TPB](
selected_tracks, track_light_voxel[track_slice][event_mask][itrk:itrk+BATCH_SIZE], selected_track_id,
light_inc, lut, light_t_start, light_sample_inc, light_sample_inc_true_track_id,
Expand Down Expand Up @@ -512,7 +517,7 @@ def run_simulation(input_filename,
light_waveforms_true_track_id_list_batch = np.concatenate(light_waveforms_true_track_id_list, axis=0)
light_waveforms_true_photons_list_batch = np.concatenate(light_waveforms_true_photons_list, axis=0)
light_trigger_modules = cp.array([detector.TPC_TO_MODULE[tpc] for tpc in light.OP_CHANNEL_TO_TPC[light_op_channel_idx_list_batch.get()][:,0]])

fee.export_to_hdf5(event_id_list_batch,
adc_tot_list_batch,
adc_tot_ticks_list_batch,
Expand Down
20 changes: 11 additions & 9 deletions larndsim/consts/light.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,11 +80,11 @@ def set_light_properties(detprop_file):
global ENABLE_LUT_SMEARING
global LIGHT_TICK_SIZE
global LIGHT_WINDOW

global SINGLET_FRACTION
global TAU_S
global TAU_T

global LIGHT_GAIN
global SIPM_RESPONSE_MODEL
global LIGHT_RESPONSE_TIME
Expand All @@ -106,34 +106,35 @@ def set_light_properties(detprop_file):
LUT_VOX_DIV = np.array(detprop['lut_vox_div'])
N_OP_CHANNEL = detprop['n_op_channel']
OP_CHANNEL_EFFICIENCY = np.array(detprop['op_channel_efficiency'])

tpc_to_op_channel = detprop['tpc_to_op_channel']
OP_CHANNEL_TO_TPC = np.zeros((N_OP_CHANNEL,), int)
TPC_TO_OP_CHANNEL = np.zeros((len(tpc_to_op_channel), len(tpc_to_op_channel[0])), int)
for itpc in range(len(tpc_to_op_channel)):
TPC_TO_OP_CHANNEL[itpc] = np.array(tpc_to_op_channel[itpc])
for idet in tpc_to_op_channel[itpc]:
for itpc, tpc_to_op in enumerate(tpc_to_op_channel):
TPC_TO_OP_CHANNEL[itpc] = np.array(tpc_to_op)
for idet in tpc_to_op:
OP_CHANNEL_TO_TPC[idet] = itpc

ENABLE_LUT_SMEARING = bool(detprop.get('enable_lut_smearing', ENABLE_LUT_SMEARING))
LIGHT_TICK_SIZE = float(detprop.get('light_tick_size', LIGHT_TICK_SIZE))
LIGHT_WINDOW = tuple(detprop.get('light_window', LIGHT_WINDOW))
assert len(LIGHT_WINDOW) == 2

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))

LIGHT_GAIN = np.array(detprop.get('light_gain', np.full(OP_CHANNEL_EFFICIENCY.shape, LIGHT_GAIN)))

if LIGHT_GAIN.size == 1:
LIGHT_GAIN = np.full(OP_CHANNEL_EFFICIENCY.shape, LIGHT_GAIN)

assert LIGHT_GAIN.shape == OP_CHANNEL_EFFICIENCY.shape
SIPM_RESPONSE_MODEL = int(detprop.get('sipm_response_model', SIPM_RESPONSE_MODEL))
assert SIPM_RESPONSE_MODEL in (0,1)
LIGHT_DET_NOISE_SAMPLE_SPACING = float(detprop.get('light_det_noise_sample_spacing', LIGHT_DET_NOISE_SAMPLE_SPACING))
LIGHT_RESPONSE_TIME = float(detprop.get('light_response_time', LIGHT_RESPONSE_TIME))
LIGHT_OSCILLATION_PERIOD = float(detprop.get('light_oscillation_period', LIGHT_OSCILLATION_PERIOD))
impulse_model_filename = str(detprop.get('impulse_model', ''))

if impulse_model_filename and SIPM_RESPONSE_MODEL == 1:
print('Light impulse model:', impulse_model_filename)
try:
Expand All @@ -145,6 +146,7 @@ def set_light_properties(detprop_file):
IMPULSE_MODEL = np.load(os.path.join(os.path.dirname(__file__), '../../') + impulse_model_filename)
except FileNotFoundError:
print("Impulse model file not found:", impulse_model_filename)

IMPULSE_TICK_SIZE = float(detprop.get('impulse_tick_size', IMPULSE_TICK_SIZE))

OP_CHANNEL_PER_TRIG = int(detprop.get('op_channel_per_det', OP_CHANNEL_PER_TRIG))
Expand All @@ -154,7 +156,7 @@ def set_light_properties(detprop_file):
LIGHT_DIGIT_SAMPLE_SPACING = float(detprop.get('light_digit_sample_spacing', LIGHT_DIGIT_SAMPLE_SPACING))
LIGHT_NBIT = int(detprop.get('light_nbit', LIGHT_NBIT))



except KeyError:
LIGHT_SIMULATED = False
Loading

0 comments on commit a9477f8

Please sign in to comment.