From ce4f1edf39cc6c9728e66d4a4a55253ac2d99db4 Mon Sep 17 00:00:00 2001 From: liviocali Date: Wed, 1 Nov 2023 16:25:58 -0500 Subject: [PATCH 1/5] Update Light MC event generator for new trigger logic. Check configs in yaml to change back --- .../reco/light/mc_event_generator.py | 16 ++++++++-------- .../reco/light/LightEventGeneratorMC.yaml | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/proto_nd_flow/reco/light/mc_event_generator.py b/src/proto_nd_flow/reco/light/mc_event_generator.py index 83c5235f..a27b14ab 100644 --- a/src/proto_nd_flow/reco/light/mc_event_generator.py +++ b/src/proto_nd_flow/reco/light/mc_event_generator.py @@ -22,8 +22,8 @@ class LightEventGeneratorMC(H5FlowGenerator): - ``wvfm_dset_name`` : ``str``, required, path to dataset to store raw waveforms - ``n_adcs`` : ``int``, number of ADC serial numbers - ``n_channels`` : ``int``, number of channels per ADC - - ``n_sipms`` : ``int``, total number of SiPM channels in simulation data - - ``n_modules`` : ``int``, number of proto ND modules + - ``n_sipms_per_module`` : ``int``, total number of SiPM channels per module + - ``n_modules_per_trig`` : ``int``, number of proto ND modules per trigger - ``adc_sn`` : ``list`` of ``int``, serial number of each ADC - ``channel_map``: ``list`` of ``list`` of ``int``, mapping from simulation optical detector to adc, channel, ``-1`` values indicate channel is not connected - ``busy_channel``: ``list`` of ``int``, channel used for busy signal on each ADC (if relevant) @@ -79,8 +79,8 @@ class LightEventGeneratorMC(H5FlowGenerator): mc_truth_dset_name='mc_truth/light', n_adcs=2, n_channels=64, - n_sipms=96, - n_modules=4, + n_sipms_per_module= 96, + n_modules_per_trig= 4, busy_delay=123, busy_ampl=20e3, chunk_size=32 @@ -167,8 +167,8 @@ def init(self): class_version=self.class_version, n_adcs=self.n_adcs, n_channels=self.n_channels, - n_sipms=self.n_sipms, - n_modules=self.n_modules, + n_sipms_per_module=self.n_sipms_per_module, + n_modules_per_trig=self.n_modules_per_trig, n_samples=self.n_samples, chunk_size=self.chunk_size, busy_delay=self.busy_delay, @@ -210,9 +210,9 @@ def next(self): next_startch = [tr_ev[0][0] for tr_ev in next_trig] # convert channel map - tmp_wvfms = np.zeros((next_wvfms.shape[0], self.n_sipms, next_wvfms.shape[-1]),dtype=next_wvfms.dtype) + tmp_wvfms = np.zeros((next_wvfms.shape[0], next_wvfms.shape[1], next_wvfms.shape[-1]),dtype=next_wvfms.dtype) for ev in range(next_wvfms.shape[0]): - tmp_wvfms[ev,next_startch[ev]:next_startch[ev]+(self.n_sipms//self.n_modules),:] = next_wvfms[ev] + tmp_wvfms[ev,next_startch[ev]:next_startch[ev]+(self.n_sipms_per_module*self.n_modules_per_trig),:] = next_wvfms[ev] next_wvfms=tmp_wvfms del tmp_wvfms remapped_wvfms = self._remap_array(self.channel_map, -next_wvfms, axis=-2) diff --git a/yamls/proto_nd_flow/reco/light/LightEventGeneratorMC.yaml b/yamls/proto_nd_flow/reco/light/LightEventGeneratorMC.yaml index dd6e334f..3ab90a56 100644 --- a/yamls/proto_nd_flow/reco/light/LightEventGeneratorMC.yaml +++ b/yamls/proto_nd_flow/reco/light/LightEventGeneratorMC.yaml @@ -8,8 +8,8 @@ params: # configuration parameters n_adcs: 8 n_channels: 64 - n_sipms: 384 - n_modules: 4 + n_sipms_per_module: 96 + n_modules_per_trig: 4 adc_sn: - 0 - 1 From d86d389bac90f24ff9574c297da7ab14555118ad Mon Sep 17 00:00:00 2001 From: Karolina Wresilo Date: Thu, 2 Nov 2023 14:09:16 -0700 Subject: [PATCH 2/5] add yamls for proto_nd c+l matching --- src/proto_nd_flow/reco/charge/charge2light.py | 222 ++++++++++++++++++ .../reco/charge/Charge2LightAssociation.yaml | 15 ++ .../workflows/charge/charge_light_assoc.yaml | 14 ++ 3 files changed, 251 insertions(+) create mode 100644 src/proto_nd_flow/reco/charge/charge2light.py create mode 100644 yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml create mode 100644 yamls/proto_nd_flow/workflows/charge/charge_light_assoc.yaml diff --git a/src/proto_nd_flow/reco/charge/charge2light.py b/src/proto_nd_flow/reco/charge/charge2light.py new file mode 100644 index 00000000..06861c74 --- /dev/null +++ b/src/proto_nd_flow/reco/charge/charge2light.py @@ -0,0 +1,222 @@ +import numpy as np +import numpy.ma as ma +import numpy.lib.recfunctions as rfn +import logging + +from h5flow.core import H5FlowStage, resources +from h5flow import H5FLOW_MPI +import module0_flow.util.units as units + + +class Charge2LightAssociation(H5FlowStage): + ''' + Generate references between charge events and light events. In general, + matches a given light event to a given charge event if:: + + |light_unix_ts_second - charge_unix_ts_second| <= unix_ts_window + AND + |light_ts_10MHz - charge_ts_10MHz| <= ts_window + + where ``*_unix_ts_second`` is the unix timestamp of the event in seconds and + ``*_ts_10MHz`` is the timestamp in 10MHz ticks since SYNC / PPS. Creates + references from both external triggers to light events as well as references + from charge events to light events. + + Requires the ``ext_trigs_dset`` in the data cache as well as its indices + (stored under the name ``ext_trigs_dset + '_idcs'``). + + Also requires RunData resource in workflow. + + Example config:: + + charge_light_associator: + classname: Charge2LightAssociation + requires: + - 'charge/ext_trigs' + - name: 'charge/ext_trigs_idcs' + path: 'charge/ext_trigs' + index_only: True + params: + light_event_dset_name: 'light/events' + ext_trigs_dset_name: 'charge/ext_trigs' + unix_ts_window: 3 + ts_window: 10 + + ''' + class_version = '0.0.1' + + default_unix_ts_window = 1 # how big of a symmetric window to use with unix timestamps (0=exact match, 1=±1 second, ...) [s] + default_ts_window = 1000 # how big of a symmetric window to use with PPS timestamps (0=exact match, 10=±10 ticks, ...) [ticks] + + def __init__(self, **params): + super(Charge2LightAssociation, self).__init__(**params) + + self.light_event_dset_name = params.get('light_event_dset_name') + self.ext_trigs_dset_name = params.get('ext_trigs_dset_name') + self.events_dset_name = None # put off until init stage + + self.unix_ts_window = params.get('unix_ts_window', self.default_unix_ts_window) + self.ts_window = params.get('ts_window', self.default_ts_window) + + self.total_charge_events = 0 + self.total_charge_triggers = 0 + self.total_light_events = 0 + self.total_matched_triggers = 0 + self.total_matched_events = 0 + self.matched_light = np.zeros((0,), dtype=bool) + self.total_matched_light = 0 + + def init(self, source_name): + super(Charge2LightAssociation, self).init(source_name) + + # save all config info + self.events_dset_name = source_name + self.data_manager.set_attrs(self.events_dset_name, + charge_to_light_assoc_classname=self.classname, + charge_to_light_assoc_class_version=self.class_version, + light_event_dset=self.light_event_dset_name, + charge_to_light_assoc_unix_ts_window=self.unix_ts_window, + charge_to_light_assoc_ts_window=self.ts_window + ) + + # then set up new datasets + self.data_manager.create_ref(self.events_dset_name, self.light_event_dset_name) + self.data_manager.create_ref(self.ext_trigs_dset_name, self.light_event_dset_name) + + # load in light system timestamps (use max to get non-null timestamp entries) + self.light_event_id = self.data_manager.get_dset(self.light_event_dset_name)['id'][:] + self.light_event_mask = self.data_manager.get_dset(self.light_event_dset_name)['wvfm_valid'][:].astype(bool) + self.light_unix_ts = self.data_manager.get_dset(self.light_event_dset_name)['utime_ms'][:] + # reshape unix ts array to use with mask + self.light_unix_ts = self.light_unix_ts[:, :, np.newaxis] + self.light_unix_ts = np.where(self.light_event_mask, self.light_unix_ts, 0) + self.light_unix_ts = ma.array(self.light_unix_ts, mask=~self.light_event_mask).mean(axis=-1).mean(axis=-1) + self.light_unix_ts = self.light_unix_ts * (units.ms / units.s) # convert ms -> s + self.light_ts = self.data_manager.get_dset(self.light_event_dset_name)['tai_ns'][:] + # reshape tai_ns array as above + self.light_ts = self.light_ts[:, :, np.newaxis] + self.light_ts = np.where(self.light_event_mask, self.light_ts, 0) + self.light_ts = ma.array(self.light_ts, mask=~self.light_event_mask).mean(axis=-1).mean(axis=-1) + if not resources['RunData'].is_mc: + self.light_ts = self.light_ts % int(1e9) + self.light_ts = self.light_ts * (units.ns / resources['RunData'].crs_ticks) # convert ns -> larpix clock ticks + + self.light_unix_ts_start = self.light_unix_ts.min() + self.light_unix_ts_end = self.light_unix_ts.max() + self.total_light_events = len(self.light_unix_ts) + self.matched_light = np.zeros((self.total_light_events,), dtype=bool) + + def finish(self, source_name): + super(Charge2LightAssociation, self).finish(source_name) + + if H5FLOW_MPI: + self.total_charge_events = self.comm.reduce(self.total_charge_events, root=0) + self.total_charge_triggers = self.comm.reduce(self.total_charge_triggers, root=0) + self.total_matched_triggers = self.comm.reduce(self.total_matched_triggers, root=0) + self.total_matched_events = self.comm.reduce(self.total_matched_events, root=0) + self.matched_light = self.comm.reduce(self.matched_light, root=0) + + if self.rank == 0: + self.total_matched_light = self.matched_light.clip(0,1).sum() + trigger_eff = self.total_matched_triggers/max(self.total_charge_triggers, 1) + event_eff = self.total_matched_events/max(self.total_charge_events, 1) + light_eff = self.total_matched_light/max(self.total_light_events, 1) + print(f'Total charge trigger matching: {self.total_matched_triggers}/{self.total_charge_triggers} ({trigger_eff:0.04f})') + print(f'Total charge event matching: {self.total_matched_events}/{self.total_charge_events} ({event_eff:0.04f})') + print(f'Total light event matching: {self.total_matched_light}/{self.total_light_events} ({light_eff:0.04f})') + + + def match_on_timestamp(self, charge_unix_ts, charge_pps_ts): + unix_ts_start = charge_unix_ts.min() + unix_ts_end = charge_unix_ts.max() + + if self.light_unix_ts_start >= unix_ts_end + self.unix_ts_window or \ + self.light_unix_ts_end <= unix_ts_start - self.unix_ts_window: + # no overlap, short circuit + return np.empty((0, 2), dtype=int) + + # subselect only portion of light events that overlaps with unix timestamps + i_min = np.argmax((self.light_unix_ts >= unix_ts_start - self.unix_ts_window)) + i_max = len(self.light_unix_ts) - 1 - np.argmax((self.light_unix_ts <= unix_ts_end + self.unix_ts_window)[::-1]) + sl = slice(i_min, i_max) + + assoc_mat = (np.abs(self.light_unix_ts[sl].reshape(1, -1) - charge_unix_ts.reshape(-1, 1)) <= self.unix_ts_window) \ + & (np.abs(self.light_ts[sl].reshape(1, -1) - charge_pps_ts.reshape(-1, 1)) <= self.ts_window) + idcs = np.argwhere(assoc_mat) + if len(idcs): + idcs[:, 1] = self.light_event_id[sl][idcs[:, 1]] # idcs now contains ext trigger index <-> global light event id + else: + idcs = np.empty((0,2), dtype=int) + + return idcs + + + def match_on_timestamp(self, charge_unix_ts, charge_pps_ts): + unix_ts_start = charge_unix_ts.min() + unix_ts_end = charge_unix_ts.max() + + if self.light_unix_ts_start >= unix_ts_end + self.unix_ts_window or \ + self.light_unix_ts_end <= unix_ts_start - self.unix_ts_window: + # no overlap, short circuit + return np.empty((0, 2), dtype=int) + + # subselect only portion of light events that overlaps with unix timestamps + i_min = np.argmax((self.light_unix_ts >= unix_ts_start - self.unix_ts_window)) + i_max = len(self.light_unix_ts) - 1 - np.argmax((self.light_unix_ts <= unix_ts_end + self.unix_ts_window)[::-1]) + sl = slice(i_min, i_max) + + assoc_mat = (np.abs(self.light_unix_ts[sl].reshape(1, -1) - charge_unix_ts.reshape(-1, 1)) <= self.unix_ts_window) \ + & (np.abs(self.light_ts[sl].reshape(1, -1) - charge_pps_ts.reshape(-1, 1)) <= self.ts_window) + idcs = np.argwhere(assoc_mat) + if len(idcs): + idcs[:, 1] = self.light_event_id[sl][idcs[:, 1]] # idcs now contains ext trigger index <-> global light event id + else: + idcs = np.empty((0,2), dtype=int) + + return idcs + + def run(self, source_name, source_slice, cache): + super(Charge2LightAssociation, self).run(source_name, source_slice, cache) + + event_data = cache[self.events_dset_name] + ext_trigs_data = cache[self.ext_trigs_dset_name] + ext_trigs_idcs = cache[self.ext_trigs_dset_name + '_idcs'] + ext_trigs_mask = ~rfn.structured_to_unstructured(ext_trigs_data.mask).any(axis=-1) + + nevents = len(event_data) + + ev_id = np.arange(source_slice.start, source_slice.stop, dtype=int) + ext_trig_ref = np.empty((0, 2), dtype=int) + ev_ref = np.empty((0, 2), dtype=int) + + # check match on external triggers + if nevents: + ext_trigs_mask = ~rfn.structured_to_unstructured(ext_trigs_data.mask).any(axis=-1) + if np.any(ext_trigs_mask): + ext_trigs_all = ext_trigs_data.data[ext_trigs_mask] + ext_trigs_idcs = ext_trigs_idcs.data[ext_trigs_mask] + ext_trigs_unix_ts = np.broadcast_to(event_data['unix_ts'].reshape(-1, 1), ext_trigs_data.shape)[ext_trigs_mask] + ext_trigs_ts = ext_trigs_all['ts'] + + idcs = self.match_on_timestamp(ext_trigs_unix_ts, ext_trigs_ts) + + if len(idcs): + ext_trig_ref = np.append(ext_trig_ref, np.c_[ext_trigs_idcs[idcs[:, 0]], idcs[:, 1]], axis=0) + ev_id_bcast = np.broadcast_to(ev_id[:,np.newaxis], ext_trigs_mask.shape) + ev_ref = np.unique(np.append(ev_ref, np.c_[ev_id_bcast[ext_trigs_mask][idcs[:, 0]], idcs[:, 1]], axis=0), axis=0) + + logging.info(f'found charge/light match on {len(ext_trig_ref)}/{ext_trigs_mask.sum()} triggers') + logging.info(f'found charge/light match on {len(ev_ref)}/{len(event_data)} events') + self.total_charge_triggers += ext_trigs_mask.sum() + self.total_matched_triggers += len(np.unique(ext_trig_ref[:,0])) + self.total_matched_events += len(np.unique(ev_ref[:,0])) + self.matched_light[np.unique(ext_trig_ref[:,1])] = True + + # write references + # ext trig -> light event + self.data_manager.write_ref(self.ext_trigs_dset_name, self.light_event_dset_name, ext_trig_ref) + + # charge event -> light event + self.data_manager.write_ref(self.events_dset_name, self.light_event_dset_name, ev_ref) + + self.total_charge_events += len(event_data) diff --git a/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml b/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml new file mode 100644 index 00000000..8c7e7184 --- /dev/null +++ b/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml @@ -0,0 +1,15 @@ +classname: Charge2LightAssociation # reco/charge/charge2light.py +path: proto_nd_flow.reco.charge.charge2light +requires: + - 'charge/ext_trigs' + - name: 'charge/ext_trigs_idcs' + path: 'charge/ext_trigs' + index_only: True +params: + # inputs + light_event_dset_name: 'light/events' + ext_trigs_dset_name: 'charge/ext_trigs' + + # configuration parameters + unix_ts_window: 1 + ts_window: 20 diff --git a/yamls/proto_nd_flow/workflows/charge/charge_light_assoc.yaml b/yamls/proto_nd_flow/workflows/charge/charge_light_assoc.yaml new file mode 100644 index 00000000..02d1c8ff --- /dev/null +++ b/yamls/proto_nd_flow/workflows/charge/charge_light_assoc.yaml @@ -0,0 +1,14 @@ +# Generates associations between charge events and light events + +flow: + source: 'charge/events' + stages: [charge_light_associator] + drop: [] + + +resources: + - !include yamls/proto_nd_flow/resources/RunData.yaml + + +charge_light_associator: + !include yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml From ac753fce937e3ad86504c75857341434c5661f9e Mon Sep 17 00:00:00 2001 From: Karolina Wresilo Date: Thu, 2 Nov 2023 14:13:33 -0700 Subject: [PATCH 3/5] update charge2light --- src/proto_nd_flow/reco/charge/charge2light.py | 25 ------------------- 1 file changed, 25 deletions(-) diff --git a/src/proto_nd_flow/reco/charge/charge2light.py b/src/proto_nd_flow/reco/charge/charge2light.py index 06861c74..2930ef54 100644 --- a/src/proto_nd_flow/reco/charge/charge2light.py +++ b/src/proto_nd_flow/reco/charge/charge2light.py @@ -150,31 +150,6 @@ def match_on_timestamp(self, charge_unix_ts, charge_pps_ts): return idcs - - def match_on_timestamp(self, charge_unix_ts, charge_pps_ts): - unix_ts_start = charge_unix_ts.min() - unix_ts_end = charge_unix_ts.max() - - if self.light_unix_ts_start >= unix_ts_end + self.unix_ts_window or \ - self.light_unix_ts_end <= unix_ts_start - self.unix_ts_window: - # no overlap, short circuit - return np.empty((0, 2), dtype=int) - - # subselect only portion of light events that overlaps with unix timestamps - i_min = np.argmax((self.light_unix_ts >= unix_ts_start - self.unix_ts_window)) - i_max = len(self.light_unix_ts) - 1 - np.argmax((self.light_unix_ts <= unix_ts_end + self.unix_ts_window)[::-1]) - sl = slice(i_min, i_max) - - assoc_mat = (np.abs(self.light_unix_ts[sl].reshape(1, -1) - charge_unix_ts.reshape(-1, 1)) <= self.unix_ts_window) \ - & (np.abs(self.light_ts[sl].reshape(1, -1) - charge_pps_ts.reshape(-1, 1)) <= self.ts_window) - idcs = np.argwhere(assoc_mat) - if len(idcs): - idcs[:, 1] = self.light_event_id[sl][idcs[:, 1]] # idcs now contains ext trigger index <-> global light event id - else: - idcs = np.empty((0,2), dtype=int) - - return idcs - def run(self, source_name, source_slice, cache): super(Charge2LightAssociation, self).run(source_name, source_slice, cache) From d35071f41d9899d0fe3e45a96fd374aab6d14345 Mon Sep 17 00:00:00 2001 From: Karolina Wresilo Date: Tue, 7 Nov 2023 04:33:08 -0800 Subject: [PATCH 4/5] working version of cl matching + matching window change --- src/proto_nd_flow/reco/charge/charge2light.py | 23 ++++++++++--------- .../reco/charge/Charge2LightAssociation.yaml | 2 +- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/src/proto_nd_flow/reco/charge/charge2light.py b/src/proto_nd_flow/reco/charge/charge2light.py index 2930ef54..72f02106 100644 --- a/src/proto_nd_flow/reco/charge/charge2light.py +++ b/src/proto_nd_flow/reco/charge/charge2light.py @@ -5,7 +5,7 @@ from h5flow.core import H5FlowStage, resources from h5flow import H5FLOW_MPI -import module0_flow.util.units as units +import proto_nd_flow.util.units as units class Charge2LightAssociation(H5FlowStage): @@ -87,16 +87,18 @@ def init(self, source_name): self.light_event_id = self.data_manager.get_dset(self.light_event_dset_name)['id'][:] self.light_event_mask = self.data_manager.get_dset(self.light_event_dset_name)['wvfm_valid'][:].astype(bool) self.light_unix_ts = self.data_manager.get_dset(self.light_event_dset_name)['utime_ms'][:] + self.light_unix_ts = self.light_unix_ts.mean(axis=-1) # reshape unix ts array to use with mask - self.light_unix_ts = self.light_unix_ts[:, :, np.newaxis] - self.light_unix_ts = np.where(self.light_event_mask, self.light_unix_ts, 0) - self.light_unix_ts = ma.array(self.light_unix_ts, mask=~self.light_event_mask).mean(axis=-1).mean(axis=-1) + # self.light_unix_ts = self.light_unix_ts[:, :, np.newaxis] + # self.light_unix_ts = np.where(self.light_event_mask, self.light_unix_ts, 0) + # self.light_unix_ts = ma.array(self.light_unix_ts, mask=~self.light_event_mask).mean(axis=-1).mean(axis=-1) self.light_unix_ts = self.light_unix_ts * (units.ms / units.s) # convert ms -> s self.light_ts = self.data_manager.get_dset(self.light_event_dset_name)['tai_ns'][:] + self.light_ts = self.light_ts.mean(axis=-1) # reshape tai_ns array as above - self.light_ts = self.light_ts[:, :, np.newaxis] - self.light_ts = np.where(self.light_event_mask, self.light_ts, 0) - self.light_ts = ma.array(self.light_ts, mask=~self.light_event_mask).mean(axis=-1).mean(axis=-1) + # self.light_ts = self.light_ts[:, :, np.newaxis] + # self.light_ts = np.where(self.light_event_mask, self.light_ts, 0) + # self.light_ts = ma.array(self.light_ts, mask=~self.light_event_mask).mean(axis=-1).mean(axis=-1) if not resources['RunData'].is_mc: self.light_ts = self.light_ts % int(1e9) self.light_ts = self.light_ts * (units.ns / resources['RunData'].crs_ticks) # convert ns -> larpix clock ticks @@ -123,8 +125,7 @@ def finish(self, source_name): light_eff = self.total_matched_light/max(self.total_light_events, 1) print(f'Total charge trigger matching: {self.total_matched_triggers}/{self.total_charge_triggers} ({trigger_eff:0.04f})') print(f'Total charge event matching: {self.total_matched_events}/{self.total_charge_events} ({event_eff:0.04f})') - print(f'Total light event matching: {self.total_matched_light}/{self.total_light_events} ({light_eff:0.04f})') - + print(f'Total light event matching: {self.total_matched_light}/{self.total_light_events} ({light_eff:0.04f})') def match_on_timestamp(self, charge_unix_ts, charge_pps_ts): unix_ts_start = charge_unix_ts.min() @@ -143,6 +144,7 @@ def match_on_timestamp(self, charge_unix_ts, charge_pps_ts): assoc_mat = (np.abs(self.light_unix_ts[sl].reshape(1, -1) - charge_unix_ts.reshape(-1, 1)) <= self.unix_ts_window) \ & (np.abs(self.light_ts[sl].reshape(1, -1) - charge_pps_ts.reshape(-1, 1)) <= self.ts_window) idcs = np.argwhere(assoc_mat) + #idcs = 1000 if len(idcs): idcs[:, 1] = self.light_event_id[sl][idcs[:, 1]] # idcs now contains ext trigger index <-> global light event id else: @@ -159,7 +161,7 @@ def run(self, source_name, source_slice, cache): ext_trigs_mask = ~rfn.structured_to_unstructured(ext_trigs_data.mask).any(axis=-1) nevents = len(event_data) - + print('nevents') ev_id = np.arange(source_slice.start, source_slice.stop, dtype=int) ext_trig_ref = np.empty((0, 2), dtype=int) ev_ref = np.empty((0, 2), dtype=int) @@ -172,7 +174,6 @@ def run(self, source_name, source_slice, cache): ext_trigs_idcs = ext_trigs_idcs.data[ext_trigs_mask] ext_trigs_unix_ts = np.broadcast_to(event_data['unix_ts'].reshape(-1, 1), ext_trigs_data.shape)[ext_trigs_mask] ext_trigs_ts = ext_trigs_all['ts'] - idcs = self.match_on_timestamp(ext_trigs_unix_ts, ext_trigs_ts) if len(idcs): diff --git a/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml b/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml index 8c7e7184..a880506b 100644 --- a/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml +++ b/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml @@ -12,4 +12,4 @@ params: # configuration parameters unix_ts_window: 1 - ts_window: 20 + ts_window: 1000 From 7c455467f15ddb70d9de872dd3d4ee9ce4f9f82e Mon Sep 17 00:00:00 2001 From: Karolina Wresilo Date: Fri, 15 Dec 2023 09:50:27 -0800 Subject: [PATCH 5/5] 2x2 CL matching --- scripts/proto_nd_scripts/run_CL_matching.sh | 41 +++++++++++++++++++ src/proto_nd_flow/reco/charge/charge2light.py | 9 ++-- .../reco/charge/Charge2LightAssociation.yaml | 2 +- 3 files changed, 45 insertions(+), 7 deletions(-) create mode 100644 scripts/proto_nd_scripts/run_CL_matching.sh diff --git a/scripts/proto_nd_scripts/run_CL_matching.sh b/scripts/proto_nd_scripts/run_CL_matching.sh new file mode 100644 index 00000000..ea87c1b1 --- /dev/null +++ b/scripts/proto_nd_scripts/run_CL_matching.sh @@ -0,0 +1,41 @@ +#!/bin/bash +# Runs proto_nd_flow on an example file. +# Before using this script, use +# >> source get_proto_nd_input.sh +# to download all the necessary inputs into the correct directories +# + +INPUT_FILE=$1 + +OUTPUT_DIR='/global/cfs/cdirs/dune/users/kwresilo/data/MiniRun4' +OUTPUT_NAME=(${INPUT_FILE//"/"/ }) +OUTPUT_NAME=${OUTPUT_NAME[-1]} +OUTPUT_FILE="${OUTPUT_DIR}/${OUTPUT_NAME}" +OUTPUT_FILE=${OUTPUT_FILE//.h5/.matched_taco.h5} +echo ${OUTPUT_FILE} + +# for running on a login node +H5FLOW_CMD='h5flow' +# for running on a single compute node with 32 cores +#H5FLOW_CMD='srun -n32 h5flow' + +# run all stages +WORKFLOW1='yamls/proto_nd_flow/workflows/charge/charge_light_assoc.yaml' +HERE=`pwd` +#cd ndlar_flow +# assumes this is being run from ndlar_flow/scripts/proto_nd_flow: +cd ../../ + +# avoid being asked if we want to overwrite the file if it exists. +# this is us answering "yes". +if [ -e $OUTPUT_FILE ]; then + rm $OUTPUT_FILE +fi + +$H5FLOW_CMD -c $WORKFLOW1 -i $INPUT_FILE -o $OUTPUT_FILE + +echo "Done!" +echo "Output can be found at $OUTPUT_FILE" + +cd ${HERE} + diff --git a/src/proto_nd_flow/reco/charge/charge2light.py b/src/proto_nd_flow/reco/charge/charge2light.py index 72f02106..da842ce1 100644 --- a/src/proto_nd_flow/reco/charge/charge2light.py +++ b/src/proto_nd_flow/reco/charge/charge2light.py @@ -129,8 +129,8 @@ def finish(self, source_name): def match_on_timestamp(self, charge_unix_ts, charge_pps_ts): unix_ts_start = charge_unix_ts.min() - unix_ts_end = charge_unix_ts.max() - + unix_ts_end = charge_unix_ts.max() + if self.light_unix_ts_start >= unix_ts_end + self.unix_ts_window or \ self.light_unix_ts_end <= unix_ts_start - self.unix_ts_window: # no overlap, short circuit @@ -138,13 +138,11 @@ def match_on_timestamp(self, charge_unix_ts, charge_pps_ts): # subselect only portion of light events that overlaps with unix timestamps i_min = np.argmax((self.light_unix_ts >= unix_ts_start - self.unix_ts_window)) - i_max = len(self.light_unix_ts) - 1 - np.argmax((self.light_unix_ts <= unix_ts_end + self.unix_ts_window)[::-1]) + i_max = len(self.light_unix_ts) - np.argmax((self.light_unix_ts <= unix_ts_end + self.unix_ts_window)[::-1]) sl = slice(i_min, i_max) - assoc_mat = (np.abs(self.light_unix_ts[sl].reshape(1, -1) - charge_unix_ts.reshape(-1, 1)) <= self.unix_ts_window) \ & (np.abs(self.light_ts[sl].reshape(1, -1) - charge_pps_ts.reshape(-1, 1)) <= self.ts_window) idcs = np.argwhere(assoc_mat) - #idcs = 1000 if len(idcs): idcs[:, 1] = self.light_event_id[sl][idcs[:, 1]] # idcs now contains ext trigger index <-> global light event id else: @@ -161,7 +159,6 @@ def run(self, source_name, source_slice, cache): ext_trigs_mask = ~rfn.structured_to_unstructured(ext_trigs_data.mask).any(axis=-1) nevents = len(event_data) - print('nevents') ev_id = np.arange(source_slice.start, source_slice.stop, dtype=int) ext_trig_ref = np.empty((0, 2), dtype=int) ev_ref = np.empty((0, 2), dtype=int) diff --git a/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml b/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml index a880506b..7bb67661 100644 --- a/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml +++ b/yamls/proto_nd_flow/reco/charge/Charge2LightAssociation.yaml @@ -12,4 +12,4 @@ params: # configuration parameters unix_ts_window: 1 - ts_window: 1000 + ts_window: 1000000