diff --git a/data/proto_nd_flow/filter_model_pars.pkl b/data/proto_nd_flow/filter_model_pars.pkl deleted file mode 100644 index 03850577..00000000 Binary files a/data/proto_nd_flow/filter_model_pars.pkl and /dev/null differ diff --git a/src/proto_nd_flow/reco/charge/calib_noise_filter.py b/src/proto_nd_flow/reco/charge/calib_noise_filter.py deleted file mode 100644 index fac0ad51..00000000 --- a/src/proto_nd_flow/reco/charge/calib_noise_filter.py +++ /dev/null @@ -1,292 +0,0 @@ -import numpy as np -import numpy.ma as ma -import logging -from collections import defaultdict -import sys -from h5flow.core import H5FlowStage, resources -from sklearn.neighbors import KernelDensity -from pickle import dump, load - -from proto_nd_flow.reco.charge.calib_prompt_hits import CalibHitBuilder - -## Some useful functions for the filter classes -## -def unique_to_io_group(unique): - return ((unique // (100*1000*1000)) % 1000) - -def unique_channel_id(d): - return ((d['io_group'].astype(int)*1000+d['io_channel'].astype(int))*1000 \ - + d['chip_id'].astype(int))*100 + d['channel_id'].astype(int) - -def unique_to_channel_id(unique): - return (unique % 100) - -## - -## Filter classes -## - -class low_current_filter: - ''' - Filters out hits with low indunced current. Specify threshold to remove hitsnwith Q1 and hits.shape[0]>1: - return self.filter_chunked(hits) - lrs=np.zeros(hits.shape) - un = unique_channel_id(hits) - un_chip = (un//100)*100 - - for chip in set(un_chip[~un_chip.mask]): - - mask = un_chip==chip - - n_chip_hits = np.sum(mask) - if np.sum(mask)<1: continue - - min_ts = np.min(hits[mask]['ts_pps']) - - - for chan in set(hits[mask]['channel_id']): - - cc = unique_to_channel_id(chan) - - if not cc in [6,7,24]: continue - - m = np.logical_and(mask, hits['channel_id']==chan) - - chan_nhit = np.sum(m) - - m = np.logical_and(m, hits['Q']self.RANGE_Q[0]) - - ts = hits['ts_pps'][m].astype(int)-min_ts - qs = hits['Q'][m] - sumqs = np.array([np.sum(hits['Q'][m])]*ts.shape[0]) - if np.sum(m)<1: continue - - lrs[m] = self.get_lr(n_chip_hits, np.array([ts, qs, sumqs]).transpose(), cc, chan_nhit) - - return lrs > self.lrs_cut - - def get_lr(self, n_chip_hits, X, chan, chan_nhit): - - #scale data - X[:, 0] = X[:, 0]/self.SCALE_T - X[:, 1] = X[:, 1]/self.SCALE_Q - if X.shape[1]==3: X[:, 2] = X[:, 2]/self.SCALE_Q - - nhit_lr=1 - chip_nhit=n_chip_hits-chan_nhit - - nhit_bins = self.nhit_bins - - if chip_nhit > 80: chip_nhit=80 - - kdes = self.kdes - kdes_cpth=self.kdes_cpth6 - if chan==6: - nhit_lr = self.vals6[chan_nhit]/self.vals0[chan_nhit] - nhit_lr *= self.chvals6[chan_nhit]/self.chvals0[chan_nhit] - kdes_cpth=self.kdes_cpth6 - elif chan==7: - nhit_lr = self.vals7[chan_nhit]/self.vals0[chan_nhit] - nhit_lr *= self.chvals7[chan_nhit]/self.chvals0[chan_nhit] - kdes_cpth=self.kdes_cpth7 - - elif chan==24: - nhit_lr = self.vals24[chan_nhit]/self.vals0[chan_nhit] - nhit_lr *= self.chvals24[chan_nhit]/self.chvals0[chan_nhit] - kdes_cpth=self.kdes_cpth24 - - if n_chip_hits > nhit_bins[0][0] and n_chip_hits <= nhit_bins[0][1] : - y = np.exp(kdes[0].score_samples(X))+1e-3 - y_cpt = np.exp(kdes_cpth[0].score_samples(X))+1e-5 - return nhit_lr*y_cpt/y - - if n_chip_hits > nhit_bins[1][0] and n_chip_hits <= nhit_bins[1][1] : - y = np.exp(kdes[1].score_samples(X))+1e-3 - y_cpt = np.exp(kdes_cpth[1].score_samples(X))+1e-5 - return nhit_lr*y_cpt/y - - if n_chip_hits > nhit_bins[2][0] and n_chip_hits <= nhit_bins[2][1] : - y = np.exp(kdes[2].score_samples(X))+1e-3 - y_cpt = np.exp(kdes_cpth[2].score_samples(X))+1e-5 - return nhit_lr*y_cpt/y - - else: - y = np.exp(kdes[3].score_samples(X))+1e-3 - y_cpt = np.exp(kdes_cpth[3].score_samples(X))+1e-5 - return nhit_lr*y_cpt/y - - -class hot_pixel_filter: - - def __init__(self, max_n_hits=35): - self.max_n_hits = int(max_n_hits) - - def filter(self, hits): - - un = unique_channel_id(hits) - - chans, counts = np.unique( un, return_counts=True ) - - hot_pixels = chans[counts>self.max_n_hits] - - mask = np.zeros( hits.shape ).astype(bool) - - for pix in hot_pixels: - chan_mask = un==pix - event_mask = np.sum( chan_mask, axis=-1 ) > self.max_n_hits - - mask = mask | np.dot(event_mask, chan_mask) - - - - return mask -#Main class defintion -## - -class CalibNoiseFilter(H5FlowStage): - ''' - Noise Filter for charge readout. - ''' - class_version = '0.0.0' - defaults = dict( - events_dset_name = 'charge/events', - hits_name = 'charge/calib_prompt_hits', - hit_charge_name = 'charge/calib_prompt_hits', - merged_name = 'charge/hits/calib_merged_hits', - low_current_filter__threshold=6.0, - hot_pixel_filter__max_n_hits=35, - filter_function_names = ['hot_pixel_filter'] - ) - valid_filter_functions = ['low_current_filter', 'correlated_post_trigger_filter', 'hot_pixel_filter'] - - merged_dtype = CalibHitBuilder.calib_hits_dtype - - def __init__(self, **params): - super(CalibNoiseFilter, self).__init__(**params) - for key in self.defaults: - setattr(self, key, params.get(key, self.defaults[key])) - for f in self.filter_function_names: - assert f in self.valid_filter_functions, f'invalid filter function name: {f}' - - def init(self, source_name): - super(CalibNoiseFilter, self).init(source_name) - self.data_manager.create_dset(self.merged_name, dtype=self.merged_dtype) - self.data_manager.create_ref(self.hits_name, self.merged_name) - self.data_manager.create_ref(source_name, self.merged_name) - self.data_manager.create_ref(self.events_dset_name, self.merged_name) - - self.init_filter_functions() - - def init_filter_functions(self): - self.filter_functions=[] - for filter_name in self.filter_function_names: - params={} - for p in self.__dict__.keys(): - if filter_name in p: params[p.split('__')[-1]]=getattr(self, p) - self.filter_functions.append( getattr(sys.modules[__name__], filter_name )( **params ) ) - - def default_filter_function(self, hits): - return hits['Q']>np.inf - - #@staticmethod - def filter_hits(self, hits, seg_fracs): - ''' - - :param hits: original hits array, shape: (N,M) - - :param fracs: fractional contributions of true segments per packet - - :returns: new hit array, shape: (N,m), new hit charge array, shape: (N,m), and an index array with shape (L,2), [:,0] being the index into the original hit array and [:,1] being the flattened index into the compressed new array - - ''' - - mask = hits.mask['id'].copy() - new_hits = hits.data.copy() - old_ids = hits.data['id'].copy()[...,np.newaxis] - old_id_mask = hits.mask['id'].copy()[...,np.newaxis] - filter_mask = self.default_filter_function(new_hits) - for f in self.filter_functions: - filter_mask = filter_mask | f.filter(hits) - - mask = filter_mask | mask - - new_ids = old_ids[~mask] - back_track = None - return ( - ma.array(new_hits, mask=mask), - np.c_[new_ids, np.array(range(new_ids.shape[0]))], - ma.array(back_track, mask=mask) if back_track is not None else None - ) - - def run(self, source_name, source_slice, cache): - super(CalibNoiseFilter, self).run(source_name, source_slice, cache) - - event_id = np.r_[source_slice] - packet_frac_bt = cache['packet_frac_backtrack'] - hits = cache[self.hits_name] - - merged, ref, back_track = self.filter_hits(hits, seg_fracs=packet_frac_bt) - - merged_mask = merged.mask['id'] - - # first write the new merged hits to the file - new_nhit = int((~merged_mask).sum()) - #print('hits after filter:', merged.shape, new_nhit) - - merge_slice = self.data_manager.reserve_data(self.merged_name, new_nhit) - merge_idx = np.r_[merge_slice].astype(merged.dtype['id']) - if new_nhit > 0: - ref[:,1] += merge_idx[0] # offset references based on reserved region in output file - np.place(merged['id'], ~merged_mask, merge_idx) - - self.data_manager.write_data(self.merged_name, merge_idx, merged[~merged_mask]) - - # sort based on the ID of the prompt hit, to make analysis more convenient - ref = ref[np.argsort(ref[:, 0])] - self.data_manager.write_ref(self.hits_name, self.merged_name, ref) - - ev_ref = np.c_[(np.indices(merged_mask.shape)[0] + source_slice.start)[~merged_mask], merge_idx] - self.data_manager.write_ref(source_name, self.merged_name, ev_ref) - self.data_manager.write_ref(self.events_dset_name, self.merged_name, ev_ref) diff --git a/src/proto_nd_flow/reco/charge/calib_prompt_hits.py b/src/proto_nd_flow/reco/charge/calib_prompt_hits.py index 350cdc87..01f3689e 100644 --- a/src/proto_nd_flow/reco/charge/calib_prompt_hits.py +++ b/src/proto_nd_flow/reco/charge/calib_prompt_hits.py @@ -56,12 +56,10 @@ class CalibHitBuilder(H5FlowStage): x f8, pixel x location [cm] y f8, pixel y location [cm] z f8, pixel z location [cm] - t_drift u8, drift time [ticks???] - ts_pps f8, PPS packet timestamp [ns] + t_drift u8, drift time [tick = 100ns] + ts_pps f8, PPS packet timestamp [tick = 100ns] io_group u8, io group ID (PACMAN number) io_channel u8, io channel ID (related to PACMAN number & PACMAN UART Number) - chip_id u8, chip_id on tile - channel_id u8, channel_id on single chip (0-63) Q f8, hit charge [ke-] E f8, hit energy [MeV] @@ -88,8 +86,6 @@ class CalibHitBuilder(H5FlowStage): ('ts_pps', 'u8'), ('io_group', 'u8'), ('io_channel', 'u8'), - ('chip_id', 'u8'), - ('channel_id', 'u8'), ('Q', 'f8'), ('E', 'f8') ]) @@ -106,6 +102,11 @@ def __init__(self, **params): self.t0_dset_name = params.get('t0_dset_name') self.pedestal_file = params.get('pedestal_file', '') self.configuration_file = params.get('configuration_file', '') + self.pedestal_mv = params.get('pedestal_mv', 580.0) + self.vref_mv = params.get('vref_mv', 1568.0) + self.vcm_mv = params.get('vcm_mv', 478.1) + self.adc_counts = params.get('adc_counts', 256) + self.gain = params.get('gain', 4.522) def init(self, source_name): super(CalibHitBuilder, self).init(source_name) @@ -120,6 +121,7 @@ def run(self, source_name, source_slice, cache): if resources['RunData'].is_mc: packet_frac_bt = cache['packet_frac_backtrack'] packet_seg_bt = cache['packet_seg_backtrack'] + t0_data = cache[self.t0_dset_name] raw_hits = cache[self.raw_hits_dset_name] @@ -215,12 +217,19 @@ def run(self, source_name, source_slice, cache): + packets_arr['chip_id'].astype(int)*100 + packets_arr['channel_id'].astype(int)) hit_uniqueid_str = hit_uniqueid.astype(str) - vref = np.array( - [self.configuration[unique_id]['vref_mv'] for unique_id in hit_uniqueid_str]) - vcm = np.array([self.configuration[unique_id]['vcm_mv'] - for unique_id in hit_uniqueid_str]) - ped = np.array([self.pedestal[unique_id]['pedestal_mv'] - for unique_id in hit_uniqueid_str]) + if self.configuration_file != '': + vref = np.array( + [self.configuration[unique_id]['vref_mv'] for unique_id in hit_uniqueid_str]) + vcm = np.array([self.configuration[unique_id]['vcm_mv'] + for unique_id in hit_uniqueid_str]) + else: + vref = np.full(len(hit_uniqueid_str), self.vref_mv) + vcm = np.full(len(hit_uniqueid_str), self.vcm_mv) + if self.pedestal_file != '': + ped = np.array([self.pedestal[unique_id]['pedestal_mv'] + for unique_id in hit_uniqueid_str]) + else: + ped = np.full(len(hit_uniqueid_str), self.pedestal_mv) calib_hits_arr['id'] = calib_hits_slice.start + np.arange(n, dtype=int) calib_hits_arr['x'] = x if has_mc_truth: @@ -231,14 +240,12 @@ def run(self, source_name, source_slice, cache): calib_hits_arr['t_drift'] = drift_t calib_hits_arr['io_group'] = packets_arr['io_group'] calib_hits_arr['io_channel'] = packets_arr['io_channel'] - calib_hits_arr['chip_id'] = packets_arr['chip_id'] - calib_hits_arr['channel_id'] = packets_arr['channel_id'] - hits_charge = self.charge_from_dataword(packets_arr['dataword'],vref,vcm,ped) # ke- + hits_charge = self.charge_from_dataword(packets_arr['dataword'], vref, vcm, ped, self.adc_counts, self.gain) # ke- calib_hits_arr['Q'] = hits_charge # ke- #FIXME supply more realistic dEdx in the recombination; also apply measured electron lifetime calib_hits_arr['E'] = hits_charge * (1000 * units.e) / resources['LArData'].ionization_recombination(mode=2,dEdx=2) * (resources['LArData'].ionization_w / units.MeV) # MeV if has_mc_truth: - true_recomb = resources['LArData'].ionization_recombination(mode=1,dEdx=packet_seg_bt_arr['dEdx']) + true_recomb = resources['LArData'].ionization_recombination(mode=2,dEdx=packet_seg_bt_arr['dEdx']) calib_hits_arr['E_true_recomb_elife'] = np.divide(hits_charge.reshape((hits_charge.shape[0],1)) * (1000 * units.e), true_recomb, out=np.zeros_like(true_recomb), where=true_recomb!=0) / resources['LArData'].charge_reduction_lifetime(t_drift=drift_t_true) * (resources['LArData'].ionization_w / units.MeV) # MeV # if back tracking information was available, write the merged back tracking @@ -272,8 +279,8 @@ def run(self, source_name, source_slice, cache): @staticmethod - def charge_from_dataword(dw, vref, vcm, ped): - return (dw / 256. * (vref - vcm) + vcm - ped) / 4. # hardcoding 4 mV/ke- conv. + def charge_from_dataword(dw, vref, vcm, ped, adc_counts, gain): + return (dw / adc_counts * (vref - vcm) + vcm - ped) / gain def load_pedestals(self): if self.pedestal_file != '' and not resources['RunData'].is_mc: diff --git a/yamls/proto_nd_flow/reco/charge/CalibNoiseFilter.yaml b/yamls/proto_nd_flow/reco/charge/CalibNoiseFilter.yaml deleted file mode 100644 index c92b023a..00000000 --- a/yamls/proto_nd_flow/reco/charge/CalibNoiseFilter.yaml +++ /dev/null @@ -1,22 +0,0 @@ -classname: CalibNoiseFilter # reco/charge/calib_noise_filter -path: proto_nd_flow.reco.charge.calib_noise_filter -requires: - - 'charge/events' - - 'charge/calib_prompt_hits' - - name: 'packet_frac_backtrack' - path: ['charge/calib_prompt_hits','mc_truth/calib_prompt_hit_backtrack'] - -params: - # inputs - events_dset_name: 'charge/events' - hits_name: 'charge/calib_prompt_hits' - filter_function_names: ['hot_pixel_filter','correlated_post_trigger_filter', 'low_current_filter'] #which filter functions to apply (this is an ordered list, but all filters will be "OR'd") - - - #params specific to each cut - low_current_filter__threshold: 6.0 - hot_pixel_filter__max_n_hits: 35 - # outputs - merged_name: 'charge/calib_final_hits' - mc_hit_frac_dset_name: 'mc_truth/calib_final_hit_backtrack' - diff --git a/yamls/proto_nd_flow/workflows/charge/final_calibration_data.yaml b/yamls/proto_nd_flow/workflows/charge/final_calibration_data.yaml index b6b8a7e5..41fbdee8 100644 --- a/yamls/proto_nd_flow/workflows/charge/final_calibration_data.yaml +++ b/yamls/proto_nd_flow/workflows/charge/final_calibration_data.yaml @@ -3,7 +3,9 @@ flow: source: raw_events - stages: [calib_noise_filter] + #source: calib_prompt_hits + #stages: [temp_hit_builder, calib_hit_merger] + stages: [calib_hit_merger] drop: [] @@ -19,6 +21,6 @@ raw_events: params: chunk_size: 32 -calib_noise_filter: - !include yamls/proto_nd_flow/reco/charge/CalibNoiseFilter.yaml +calib_hit_merger: + !include yamls/proto_nd_flow/reco/charge/CalibHitMerger.yaml diff --git a/yamls/proto_nd_flow/workflows/charge/final_calibration_mc.yaml b/yamls/proto_nd_flow/workflows/charge/final_calibration_mc.yaml index e151a143..d7323527 100644 --- a/yamls/proto_nd_flow/workflows/charge/final_calibration_mc.yaml +++ b/yamls/proto_nd_flow/workflows/charge/final_calibration_mc.yaml @@ -3,7 +3,9 @@ flow: source: raw_events - stages: [calib_noise_filter] + #source: calib_prompt_hits + #stages: [temp_hit_builder, calib_hit_merger] + stages: [calib_hit_merger] drop: [] @@ -19,6 +21,6 @@ raw_events: params: chunk_size: 32 -calib_noise_filter: - !include yamls/proto_nd_flow/reco/charge/CalibNoiseFilter.yaml +calib_hit_merger: + !include yamls/proto_nd_flow/reco/charge/CalibHitMerger.yaml