From 088142dfeff2ed60a50273909c218e90e94d6009 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Mon, 17 Apr 2023 07:17:50 +0000 Subject: [PATCH 01/15] added plugins --- straxen/plugins/events/event_bipo_basics.py | 156 ++++++++++++++++++ straxen/plugins/events/event_bipo_matching.py | 100 +++++++++++ 2 files changed, 256 insertions(+) create mode 100644 straxen/plugins/events/event_bipo_basics.py create mode 100644 straxen/plugins/events/event_bipo_matching.py diff --git a/straxen/plugins/events/event_bipo_basics.py b/straxen/plugins/events/event_bipo_basics.py new file mode 100644 index 000000000..f763dd1f1 --- /dev/null +++ b/straxen/plugins/events/event_bipo_basics.py @@ -0,0 +1,156 @@ +import strax +import numpy as np +import numba +import straxen + + +export, __all__ = strax.exporter() + +@export +class BiPoVariables(strax.LoopPlugin): + """ + Compute: + - peak properties + - peak positions + of the first three main (in area) S1 and ten S2. + + The standard PosRec algorithm and the three different PosRec algorithms (mlp, gcn, cnn) + are given for the five S2. + """ + + __version__ = '3.0.0' + + depends_on = ('events', + 'peak_basics', + 'peak_positions', + 'peak_proximity') + + # TODO change name + provides = 'bi_po_variables' + data_kind = 'events' + loop_over = 'events' + + max_n_s1 = straxen.URLConfig(default=3, infer_type=False, + help='Number of S1s to consider') + + max_n_s2 = straxen.URLConfig(default=10, infer_type=False, + help='Number of S2s to consider') + + def setup(self): + + self.peak_properties = ( + # name dtype comment + ('time', np.int64, 'start time since unix epoch [ns]'), + ('center_time', np.int64, 'weighted center time since unix epoch [ns]'), + ('endtime', np.int64, 'end time since unix epoch [ns]'), + ('area', np.float32, 'area, uncorrected [PE]'), + ('n_channels', np.int32, 'count of contributing PMTs'), + ('n_competing', np.float32, 'number of competing PMTs'), + ('max_pmt', np.int16, 'PMT number which contributes the most PE'), + ('max_pmt_area', np.float32, 'area in the largest-contributing PMT (PE)'), + ('range_50p_area', np.float32, 'width, 50% area [ns]'), + ('range_90p_area', np.float32, 'width, 90% area [ns]'), + ('rise_time', np.float32, 'time between 10% and 50% area quantiles [ns]'), + ('area_fraction_top', np.float32, 'fraction of area seen by the top PMT array') + ) + self.to_store = [name for name, _, _ in peak_properties] + + self.pos_rec_labels = ['cnn', 'gcn', 'mlp'] # sorted alphabetically + self.posrec_save = [(xy + algo, xy + algo) for xy in ['x_', 'y_'] for algo in pos_rec_labels] # ???? + + def infer_dtype(self): + + # Basic event properties + basics_dtype = [] + basics_dtype += strax.time_fields + basics_dtype += [('n_peaks', np.int32, 'Number of peaks in the event'), + ('n_incl_peaks_s1', np.int32, 'Number of included S1 peaks in the event'), + ('n_incl_peaks_s2', np.int32, 'Number of included S2 peaks in the event')] + + # For S1s and S2s + for p_type in [1, 2]: + if p_type == 1: + max_n = self.max_n_s1 + if p_type == 2: + max_n = self.max_n_s2 + for n in range(max_n): + # Peak properties + for name, dt, comment in self.peak_properties: + basics_dtype += [(f's{p_type}_{name}_{n}', dt, f'S{p_type}_{n} {comment}'), ] + + if p_type == 2: + # S2 Peak positions + for algo in self.pos_rec_labels: + basics_dtype += [(f's2_x_{algo}_{n}', + np.float32, f'S2_{n} {algo}-reconstructed X position, uncorrected [cm]'), + (f's2_y_{algo}_{n}', + np.float32, f'S2_{n} {algo}-reconstructed Y position, uncorrected [cm]')] + + return basics_dtype + + def compute_loop(self, event, peaks): + + result = dict(time=event['time'], + endtime=strax.endtime(event)) + result['n_peaks'] = len(peaks) + + if not len(peaks): + return result + + ######## + # S1 # + ######## + + mask_s1s = (peaks['type']==1) + mask_s1s &= (peaks['area']>100) + + if not len(peaks[mask_s1s]): + return result + + ## Save the biggest peaks + max_s1s = min(self.max_n_s1, len(peaks[mask_s1s])) + + # Need to initialize them to be able to use them in S2 mask without errors + result['s1_time_0'], result['s1_time_1'] = float('nan'), float('nan') + for i, p in enumerate(reversed(np.sort(peaks[mask_s1s], order='area'))): + + for prop in self.to_store: + result[f's1_{prop}_{i}'] = p[prop] + if i == self.max_n_s1 - 1: + break + + result['n_incl_peaks_s1'] = max_s1s + + ######## + ## S2 # + ######## + + # TODO + # all this mask thingis should me moved to the next plugin + # you can have a minimal one but should be very basic + # the complicated stuff should be in the matching plugin + # and this one can be used generally (not only for bipos) + # same for the S1s + + mask_s2s = peaks['type']==2 + mask_s2s &= peaks['area'] > 1500 # low area limit + mask_s2s &= peaks['area_fraction_top'] > 0.5 # to remove S1 afterpulses + mask_s2s &= np.abs(peaks['time'] - result['s1_time_0']) > 1000 # again to remove afterpulses + mask_s2s &= np.abs(peaks['time'] - result['s1_time_1']) > 1000 # and again to remove afterpulses + mask_s2s &= peaks['time']-result['s1_time_0'] < 5000000 # 5000mus, S2 is too far in time, not related to Po + + if not len(peaks[mask_s2s]): + return result + + max_s2s = min(self.max_n_s2, len(peaks[mask_s2s])) + for i, p in enumerate(reversed(np.sort(peaks[mask_s2s], order='area'))): + for prop in self.to_store: + result[f's2_{prop}_{i}'] = p[prop] + for name_alg in self.posrec_save: + result[f's2_{name_alg[0]}_{i}'] = p[name_alg[1]] + if i == self.max_n_s2 - 1: + break + + result['n_incl_peaks_s2'] = max_s2s + + return result diff --git a/straxen/plugins/events/event_bipo_matching.py b/straxen/plugins/events/event_bipo_matching.py new file mode 100644 index 000000000..a23b9c6a2 --- /dev/null +++ b/straxen/plugins/events/event_bipo_matching.py @@ -0,0 +1,100 @@ +import strax +import numpy as np +import numba +import straxen + + +export, __all__ = strax.exporter() + +@export +class BiPo214Matching(strax.Plugin): + + depends_on=('bi_po_variables', ) + provides = 'bi_po_214_matching' + + __version__ = "2.1.4" + + rechunk_on_save=False + + def infer_dtype(self): + + dtype = strax.time_fields + [ + (f's2_bi_match', np.int, + f'Index of the S2 reconstructed as a Bismuth peak'), + (f's2_po_match', np.int, + f'Index of the S2 reconstructed as a Polonium peak')] + + return dtype + + def setup(self): + + self.tol = 3000 # 2 mus tolerance window + + def compute(self,events): + + result = np.ones(len(events), dtype=self.dtype) + result['time'] = events['time'] + result['endtime'] = events['endtime'] + + # Calculate the delta_t between the two S1s + dt = events['s1_center_time_0'] - events['s1_center_time_1'] + + L = [] + for n in range(10): + L.append(str(n)) + + combs = list(itertools.combinations(L,2)) + + ds1_u = (dt + self.tol) + ds1_l = (dt - self.tol) + for i in np.where(ds1_l < self.tol): + ds1_l[i] = 0 + + flag = np.repeat('-1_-1', len(dt)) + overall_cross_mask = np.full(len(dt), False) + + # try every possible combination + for comb in combs: + + # compute dt for this combination + dtcomb = events['s2_center_time_'+comb[1]] - events['s2_center_time_'+comb[0]] + dtcomb = np.abs(dtcomb) + + # mask true if match + mask = (ds1_l < dtcomb) & (ds1_u > dtcomb) + + #mask &= (events['s2_area_fraction_top_'+comb[0]] > .55) + #mask &= (events['s2_area_fraction_top_'+comb[1]] > .55) + + # check if we found a match for something that was already matched + # if two combinations are possible discard the event + cross_mask = mask & (flag != '-1_-1') + overall_cross_mask |= cross_mask + + flag = np.where(mask , comb[0]+'_'+comb[1], flag) + + flag = np.where(overall_cross_mask , '-2_-2', flag) + + # This part is very badly coded but ok.. + + flag_bi = [] + flag_po = [] + + for i, f in enumerate(flag): + c = list(map(int, f.split('_'))) + if c[0] >= 0: + if events['s2_center_time_%i'%c[1]][i] - events['s2_center_time_%i'%c[0]][i] < 0: + flag_bi.append(c[1]) + flag_po.append(c[0]) + else: + flag_bi.append(c[0]) + flag_po.append(c[1]) + else: + flag_bi.append(c[0]) + flag_po.append(c[1]) + + + result['s2_bi_match'] = flag_bi + result['s2_po_match'] = flag_po + + return result \ No newline at end of file From edd2fcd51d6cbc7fef5ad5f125fedab415c3fc0a Mon Sep 17 00:00:00 2001 From: cfuselli Date: Mon, 1 May 2023 08:36:59 +0000 Subject: [PATCH 02/15] still testing --- straxen/plugins/events/__init__.py | 6 + straxen/plugins/events/event_bipo_basics.py | 147 ++-------------- .../plugins/events/event_bipo_variables.py | 159 ++++++++++++++++++ 3 files changed, 181 insertions(+), 131 deletions(-) create mode 100644 straxen/plugins/events/event_bipo_variables.py diff --git a/straxen/plugins/events/__init__.py b/straxen/plugins/events/__init__.py index 0079976a5..ea1e15960 100644 --- a/straxen/plugins/events/__init__.py +++ b/straxen/plugins/events/__init__.py @@ -60,3 +60,9 @@ from . import multi_scatter from .multi_scatter import * + +from . import event_bipo_variables +from .event_bipo_variables import * + +from . import event_bipo_matching +from .event_bipo_matching import * diff --git a/straxen/plugins/events/event_bipo_basics.py b/straxen/plugins/events/event_bipo_basics.py index f763dd1f1..73cad9613 100644 --- a/straxen/plugins/events/event_bipo_basics.py +++ b/straxen/plugins/events/event_bipo_basics.py @@ -7,150 +7,35 @@ export, __all__ = strax.exporter() @export -class BiPoVariables(strax.LoopPlugin): +class EventBiPoBasics(straxen.EventBasics): """ - Compute: - - peak properties - - peak positions - of the first three main (in area) S1 and ten S2. - - The standard PosRec algorithm and the three different PosRec algorithms (mlp, gcn, cnn) - are given for the five S2. + Carlo explain please """ - __version__ = '3.0.0' + __version__ = '0.0.1' depends_on = ('events', + 'bi_po_variables', + 'bi_po_214_matching', 'peak_basics', 'peak_positions', 'peak_proximity') # TODO change name - provides = 'bi_po_variables' + provides = 'event_basics' data_kind = 'events' loop_over = 'events' - max_n_s1 = straxen.URLConfig(default=3, infer_type=False, - help='Number of S1s to consider') - - max_n_s2 = straxen.URLConfig(default=10, infer_type=False, - help='Number of S2s to consider') - - def setup(self): - - self.peak_properties = ( - # name dtype comment - ('time', np.int64, 'start time since unix epoch [ns]'), - ('center_time', np.int64, 'weighted center time since unix epoch [ns]'), - ('endtime', np.int64, 'end time since unix epoch [ns]'), - ('area', np.float32, 'area, uncorrected [PE]'), - ('n_channels', np.int32, 'count of contributing PMTs'), - ('n_competing', np.float32, 'number of competing PMTs'), - ('max_pmt', np.int16, 'PMT number which contributes the most PE'), - ('max_pmt_area', np.float32, 'area in the largest-contributing PMT (PE)'), - ('range_50p_area', np.float32, 'width, 50% area [ns]'), - ('range_90p_area', np.float32, 'width, 90% area [ns]'), - ('rise_time', np.float32, 'time between 10% and 50% area quantiles [ns]'), - ('area_fraction_top', np.float32, 'fraction of area seen by the top PMT array') - ) - self.to_store = [name for name, _, _ in peak_properties] - - self.pos_rec_labels = ['cnn', 'gcn', 'mlp'] # sorted alphabetically - self.posrec_save = [(xy + algo, xy + algo) for xy in ['x_', 'y_'] for algo in pos_rec_labels] # ???? - def infer_dtype(self): - - # Basic event properties - basics_dtype = [] - basics_dtype += strax.time_fields - basics_dtype += [('n_peaks', np.int32, 'Number of peaks in the event'), - ('n_incl_peaks_s1', np.int32, 'Number of included S1 peaks in the event'), - ('n_incl_peaks_s2', np.int32, 'Number of included S2 peaks in the event')] + def fill_events(self, result_buffer, events, split_peaks): + """Loop over the events and peaks within that event""" + for event_i, _ in enumerate(events): + peaks_in_event_i = split_peaks[event_i] + n_peaks = len(peaks_in_event_i) + result_buffer[event_i]['n_peaks'] = n_peaks - # For S1s and S2s - for p_type in [1, 2]: - if p_type == 1: - max_n = self.max_n_s1 - if p_type == 2: - max_n = self.max_n_s2 - for n in range(max_n): - # Peak properties - for name, dt, comment in self.peak_properties: - basics_dtype += [(f's{p_type}_{name}_{n}', dt, f'S{p_type}_{n} {comment}'), ] - - if p_type == 2: - # S2 Peak positions - for algo in self.pos_rec_labels: - basics_dtype += [(f's2_x_{algo}_{n}', - np.float32, f'S2_{n} {algo}-reconstructed X position, uncorrected [cm]'), - (f's2_y_{algo}_{n}', - np.float32, f'S2_{n} {algo}-reconstructed Y position, uncorrected [cm]')] - - return basics_dtype - - def compute_loop(self, event, peaks): - - result = dict(time=event['time'], - endtime=strax.endtime(event)) - result['n_peaks'] = len(peaks) - - if not len(peaks): - return result - - ######## - # S1 # - ######## - - mask_s1s = (peaks['type']==1) - mask_s1s &= (peaks['area']>100) - - if not len(peaks[mask_s1s]): - return result - - ## Save the biggest peaks - max_s1s = min(self.max_n_s1, len(peaks[mask_s1s])) - - # Need to initialize them to be able to use them in S2 mask without errors - result['s1_time_0'], result['s1_time_1'] = float('nan'), float('nan') - for i, p in enumerate(reversed(np.sort(peaks[mask_s1s], order='area'))): - - for prop in self.to_store: - result[f's1_{prop}_{i}'] = p[prop] - if i == self.max_n_s1 - 1: - break - - result['n_incl_peaks_s1'] = max_s1s - - ######## - ## S2 # - ######## - - # TODO - # all this mask thingis should me moved to the next plugin - # you can have a minimal one but should be very basic - # the complicated stuff should be in the matching plugin - # and this one can be used generally (not only for bipos) - # same for the S1s - - mask_s2s = peaks['type']==2 - mask_s2s &= peaks['area'] > 1500 # low area limit - mask_s2s &= peaks['area_fraction_top'] > 0.5 # to remove S1 afterpulses - mask_s2s &= np.abs(peaks['time'] - result['s1_time_0']) > 1000 # again to remove afterpulses - mask_s2s &= np.abs(peaks['time'] - result['s1_time_1']) > 1000 # and again to remove afterpulses - mask_s2s &= peaks['time']-result['s1_time_0'] < 5000000 # 5000mus, S2 is too far in time, not related to Po - - if not len(peaks[mask_s2s]): - return result - - max_s2s = min(self.max_n_s2, len(peaks[mask_s2s])) - for i, p in enumerate(reversed(np.sort(peaks[mask_s2s], order='area'))): - for prop in self.to_store: - result[f's2_{prop}_{i}'] = p[prop] - for name_alg in self.posrec_save: - result[f's2_{name_alg[0]}_{i}'] = p[name_alg[1]] - if i == self.max_n_s2 - 1: - break - - result['n_incl_peaks_s2'] = max_s2s + if not n_peaks: + raise ValueError(f'No peaks within event?\n{events[event_i]}') - return result + if event_i['s2_bi_match']>0 and event_i['s2_po_match']>0: + self.fill_result_i(result_buffer[event_i], peaks_in_event_i) \ No newline at end of file diff --git a/straxen/plugins/events/event_bipo_variables.py b/straxen/plugins/events/event_bipo_variables.py new file mode 100644 index 000000000..98ee5fba4 --- /dev/null +++ b/straxen/plugins/events/event_bipo_variables.py @@ -0,0 +1,159 @@ +import strax +import numpy as np +import numba +import straxen + + +export, __all__ = strax.exporter() + +@export +class BiPoVariables(strax.LoopPlugin): + """ + Compute: + - peak properties + - peak positions + of the first three main (in area) S1 and ten S2. + + The standard PosRec algorithm and the three different PosRec algorithms (mlp, gcn, cnn) + are given for the five S2. + """ + + __version__ = '3.0.1' + + depends_on = ('events', + 'peak_basics', + 'peak_positions', + 'peak_proximity') + + # TODO change name + provides = 'bi_po_variables' + data_kind = 'events' + loop_over = 'events' + + max_n_s1 = straxen.URLConfig(default=3, infer_type=False, + help='Number of S1s to consider') + + max_n_s2 = straxen.URLConfig(default=10, infer_type=False, + help='Number of S2s to consider') + + + peak_properties = ( + # name dtype comment + ('time', np.int64, 'start time since unix epoch [ns]'), + ('center_time', np.int64, 'weighted center time since unix epoch [ns]'), + ('endtime', np.int64, 'end time since unix epoch [ns]'), + ('area', np.float32, 'area, uncorrected [PE]'), + ('n_channels', np.int32, 'count of contributing PMTs'), + ('n_competing', np.float32, 'number of competing PMTs'), + ('max_pmt', np.int16, 'PMT number which contributes the most PE'), + ('max_pmt_area', np.float32, 'area in the largest-contributing PMT (PE)'), + ('range_50p_area', np.float32, 'width, 50% area [ns]'), + ('range_90p_area', np.float32, 'width, 90% area [ns]'), + ('rise_time', np.float32, 'time between 10% and 50% area quantiles [ns]'), + ('area_fraction_top', np.float32, 'fraction of area seen by the top PMT array') + ) + + pos_rec_labels = ['cnn', 'gcn', 'mlp'] + + def setup(self): + + self.posrec_save = [(xy + algo, xy + algo) for xy in ['x_', 'y_'] for algo in self.pos_rec_labels] # ???? + self.to_store = [name for name, _, _ in self.peak_properties] + + + def infer_dtype(self): + + # Basic event properties + basics_dtype = [] + basics_dtype += strax.time_fields + basics_dtype += [('n_peaks', np.int32, 'Number of peaks in the event'), + ('n_incl_peaks_s1', np.int32, 'Number of included S1 peaks in the event'), + ('n_incl_peaks_s2', np.int32, 'Number of included S2 peaks in the event')] + + # For S1s and S2s + for p_type in [1, 2]: + if p_type == 1: + max_n = self.max_n_s1 + if p_type == 2: + max_n = self.max_n_s2 + for n in range(max_n): + # Peak properties + for name, dt, comment in self.peak_properties: + basics_dtype += [(f's{p_type}_{name}_{n}', dt, f'S{p_type}_{n} {comment}'), ] + + if p_type == 2: + # S2 Peak positions + for algo in self.pos_rec_labels: + basics_dtype += [(f's2_x_{algo}_{n}', + np.float32, f'S2_{n} {algo}-reconstructed X position, uncorrected [cm]'), + (f's2_y_{algo}_{n}', + np.float32, f'S2_{n} {algo}-reconstructed Y position, uncorrected [cm]')] + + return basics_dtype + + def compute_loop(self, event, peaks): + + result = dict(time=event['time'], + endtime=strax.endtime(event)) + result['n_peaks'] = len(peaks) + + if not len(peaks): + return result + + ######## + # S1 # + ######## + + mask_s1s = (peaks['type']==1) + mask_s1s &= (peaks['area']>100) + + if not len(peaks[mask_s1s]): + return result + + ## Save the biggest peaks + max_s1s = min(self.max_n_s1, len(peaks[mask_s1s])) + + # Need to initialize them to be able to use them in S2 mask without errors + result['s1_time_0'], result['s1_time_1'] = float('nan'), float('nan') + for i, p in enumerate(reversed(np.sort(peaks[mask_s1s], order='area'))): + + for prop in self.to_store: + result[f's1_{prop}_{i}'] = p[prop] + if i == self.max_n_s1 - 1: + break + + result['n_incl_peaks_s1'] = max_s1s + + ######## + ## S2 # + ######## + + # TODO + # all this mask thingis should me moved to the next plugin + # you can have a minimal one but should be very basic + # the complicated stuff should be in the matching plugin + # and this one can be used generally (not only for bipos) + # same for the S1s + + mask_s2s = peaks['type']==2 + mask_s2s &= peaks['area'] > 1500 # low area limit + mask_s2s &= peaks['area_fraction_top'] > 0.5 # to remove S1 afterpulses + mask_s2s &= np.abs(peaks['time'] - result['s1_time_0']) > 1000 # again to remove afterpulses + mask_s2s &= np.abs(peaks['time'] - result['s1_time_1']) > 1000 # and again to remove afterpulses + mask_s2s &= peaks['time']-result['s1_time_0'] < 5000000 # 5000mus, S2 is too far in time, not related to Po + + if not len(peaks[mask_s2s]): + return result + + max_s2s = min(self.max_n_s2, len(peaks[mask_s2s])) + for i, p in enumerate(reversed(np.sort(peaks[mask_s2s], order='area'))): + for prop in self.to_store: + result[f's2_{prop}_{i}'] = p[prop] + for name_alg in self.posrec_save: + result[f's2_{name_alg[0]}_{i}'] = p[name_alg[1]] + if i == self.max_n_s2 - 1: + break + + result['n_incl_peaks_s2'] = max_s2s + + return result From 9accd88d121d5a3d4c90e142793e4117648c0a0b Mon Sep 17 00:00:00 2001 From: cfuselli Date: Mon, 1 May 2023 15:18:35 +0000 Subject: [PATCH 03/15] amazing refactoring --- straxen/plugins/events/event_bipo_basics.py | 45 ++++- straxen/plugins/events/event_bipo_matching.py | 161 +++++++++++------- .../plugins/events/event_bipo_variables.py | 115 ++++++------- 3 files changed, 188 insertions(+), 133 deletions(-) diff --git a/straxen/plugins/events/event_bipo_basics.py b/straxen/plugins/events/event_bipo_basics.py index 73cad9613..25e52b813 100644 --- a/straxen/plugins/events/event_bipo_basics.py +++ b/straxen/plugins/events/event_bipo_basics.py @@ -12,7 +12,7 @@ class EventBiPoBasics(straxen.EventBasics): Carlo explain please """ - __version__ = '0.0.1' + __version__ = '0.0.3' depends_on = ('events', 'bi_po_variables', @@ -37,5 +37,44 @@ def fill_events(self, result_buffer, events, split_peaks): if not n_peaks: raise ValueError(f'No peaks within event?\n{events[event_i]}') - if event_i['s2_bi_match']>0 and event_i['s2_po_match']>0: - self.fill_result_i(result_buffer[event_i], peaks_in_event_i) \ No newline at end of file + self.fill_result_i(result_buffer[event_i], peaks_in_event_i, _['s2_bi_match'], _['s2_po_match']) + + + def fill_result_i(self, event, peaks, bi_i, po_i): + """For a single event with the result_buffer""" + + if (bi_i>=0) and (po_i>=0): + + largest_s2s, s2_idx = self.get_largest_sx_peaks(peaks, s_i=2, number_of_peaks=0) + + bipo_mask = [bi_i, po_i] + s2_idx = s2_idx[bipo_mask] + largest_s2s = largest_s2s[bipo_mask] + + largest_s1s, s1_idx = self.get_largest_sx_peaks( + peaks, + s_i=1, + number_of_peaks=2) + + largest_s1s = largest_s1s[::-1] + s1_idx = s1_idx[::-1] + + self.set_sx_index(event, s1_idx, s2_idx) + self.set_event_properties(event, largest_s1s, largest_s2s, peaks) + + # Loop over S1s and S2s and over main / alt. + for s_i, largest_s_i in enumerate([largest_s1s, largest_s2s], 1): + # Largest index 0 -> main sx, 1 -> alt sx + for largest_index, main_or_alt in enumerate(['s', 'alt_s']): + peak_properties_to_save = [name for name, _, _ in self.peak_properties] + if s_i == 2: + peak_properties_to_save += ['x', 'y'] + peak_properties_to_save += self.posrec_save + field_names = [f'{main_or_alt}{s_i}_{name}' for name in peak_properties_to_save] + self.copy_largest_peaks_into_event(event, + largest_s_i, + largest_index, + field_names, + peak_properties_to_save) + + diff --git a/straxen/plugins/events/event_bipo_matching.py b/straxen/plugins/events/event_bipo_matching.py index a23b9c6a2..73f9d8f6b 100644 --- a/straxen/plugins/events/event_bipo_matching.py +++ b/straxen/plugins/events/event_bipo_matching.py @@ -2,27 +2,42 @@ import numpy as np import numba import straxen - +import itertools export, __all__ = strax.exporter() @export class BiPo214Matching(strax.Plugin): + """Plugin for matching S2 signals reconstructed as bismuth or polonium peaks in a dataset + containing BiPo214 events. + + Provides: + -------- + bi_po_214_matching : numpy.array + Array containing the indices of the S2 signals reconstructed as bismuth or polonium peaks. + The index is -9 if the S2 signal is not found, -2 if multiple S2 signals are found, + or the index of the S2 signal if only one is found. + + Configuration: + ----------- + tol : int + Time tolerance window in ns to match S2 signals to BiPo214 events. Default is 3000 ns. + """ depends_on=('bi_po_variables', ) provides = 'bi_po_214_matching' - __version__ = "2.1.4" - - rechunk_on_save=False - + __version__ = "2.1.7" + def infer_dtype(self): dtype = strax.time_fields + [ (f's2_bi_match', np.int, f'Index of the S2 reconstructed as a Bismuth peak'), (f's2_po_match', np.int, - f'Index of the S2 reconstructed as a Polonium peak')] + f'Index of the S2 reconstructed as a Polonium peak'), + (f'n_incl_peaks_s2', np.int, + f'S2s considered in the combinations, that passed the mask requirements')] return dtype @@ -30,71 +45,85 @@ def setup(self): self.tol = 3000 # 2 mus tolerance window - def compute(self,events): - - result = np.ones(len(events), dtype=self.dtype) + def compute(self, events): + result = np.zeros(len(events), dtype=self.dtype) result['time'] = events['time'] result['endtime'] = events['endtime'] + result['s2_bi_match'] = -9 + result['s2_po_match'] = -9 + + + mask = self.box_mask(events) + + s2_bi_match, s2_po_match, n_s2s = self.find_match(events[mask], self.tol) + + result['s2_bi_match'][mask] = s2_bi_match + result['s2_po_match'][mask] = s2_po_match + result['n_incl_peaks_s2'][mask] = n_s2s + + return result + + + def find_match(self, events, tol): + - # Calculate the delta_t between the two S1s - dt = events['s1_center_time_0'] - events['s1_center_time_1'] - - L = [] - for n in range(10): - L.append(str(n)) - - combs = list(itertools.combinations(L,2)) - - ds1_u = (dt + self.tol) - ds1_l = (dt - self.tol) - for i in np.where(ds1_l < self.tol): - ds1_l[i] = 0 - - flag = np.repeat('-1_-1', len(dt)) - overall_cross_mask = np.full(len(dt), False) + # Compute the time difference between the two S1s + dt_s1 = events['s1_center_time_0'] - events['s1_center_time_1'] + + dt_s1_lower = dt_s1 - tol + dt_s1_upper = dt_s1 + tol + dt_s1_lower[dt_s1_lower < 0] = 0 + + # Prepare arrays to store matched S2s + s2_bi_match = np.full(len(events), -1) + s2_po_match = np.full(len(events), -1) + n_s2s = np.full(len(events), 0) + + # Find matching S2s + for i, event in enumerate(events): + # Create a list of possible S2 pairs to match + s2s_idx = self.consider_s2s(event) + n_s2s[i] = len(s2s_idx) + possible_pairs = list(itertools.combinations(s2s_idx,2)) - # try every possible combination - for comb in combs: - - # compute dt for this combination - dtcomb = events['s2_center_time_'+comb[1]] - events['s2_center_time_'+comb[0]] - dtcomb = np.abs(dtcomb) - - # mask true if match - mask = (ds1_l < dtcomb) & (ds1_u > dtcomb) - - #mask &= (events['s2_area_fraction_top_'+comb[0]] > .55) - #mask &= (events['s2_area_fraction_top_'+comb[1]] > .55) + for pair in possible_pairs: + + p0, p1 = str(pair[0]), str(pair[1]) + t0, t1 = event['s2_center_time_' + p0], event['s2_center_time_' + p1] + + dt_s2 = abs(t0 - t1) + if dt_s1_lower[i] <= dt_s2 <= dt_s1_upper[i]: + if s2_bi_match[i] == -1: + s2_bi_match[i] = pair[np.argmin([t0, t1])] + s2_po_match[i] = pair[np.argmax([t0, t1])] + else: + s2_bi_match[i] = -2 + s2_po_match[i] = -2 + break + + return s2_bi_match, s2_po_match, n_s2s - # check if we found a match for something that was already matched - # if two combinations are possible discard the event - cross_mask = mask & (flag != '-1_-1') - overall_cross_mask |= cross_mask - - flag = np.where(mask , comb[0]+'_'+comb[1], flag) + @staticmethod + def consider_s2s(event): + res = [] + for ip in range(10): + p = '_'+str(ip) + consider = True + consider &= event['s2_area'+p] > 1500 # low area limit + consider &= event['s2_area_fraction_top'+p] > 0.5 # to remove S1 afterpulses + consider &= np.abs(event['s2_time'+p] - event['s1_time_0']) > 1000 # again to remove afterpulses + consider &= np.abs(event['s2_time'+p] - event['s1_time_1']) > 1000 # and again to remove afterpulses + consider &= event['s2_time'+p]-event['s1_time_0'] < 5000000 # 5000mus, S2 is too far in time, not related to Po + if consider: + res.append(ip) + return res - flag = np.where(overall_cross_mask , '-2_-2', flag) - - # This part is very badly coded but ok.. - - flag_bi = [] - flag_po = [] - - for i, f in enumerate(flag): - c = list(map(int, f.split('_'))) - if c[0] >= 0: - if events['s2_center_time_%i'%c[1]][i] - events['s2_center_time_%i'%c[0]][i] < 0: - flag_bi.append(c[1]) - flag_po.append(c[0]) - else: - flag_bi.append(c[0]) - flag_po.append(c[1]) - else: - flag_bi.append(c[0]) - flag_po.append(c[1]) + @staticmethod + def box_mask(events): + # s1_0 == alpha (Po) + mask = (events['s1_area_0'] > 40000) # this is implicit from band cut + mask &= (events['s1_area_0'] < 120000) + mask &= (events['s1_area_fraction_top_0'] < 0.6) - result['s2_bi_match'] = flag_bi - result['s2_po_match'] = flag_po - - return result \ No newline at end of file + return mask \ No newline at end of file diff --git a/straxen/plugins/events/event_bipo_variables.py b/straxen/plugins/events/event_bipo_variables.py index 98ee5fba4..c6334aff3 100644 --- a/straxen/plugins/events/event_bipo_variables.py +++ b/straxen/plugins/events/event_bipo_variables.py @@ -7,7 +7,7 @@ export, __all__ = strax.exporter() @export -class BiPoVariables(strax.LoopPlugin): +class BiPoVariables(strax.Plugin): """ Compute: - peak properties @@ -18,7 +18,7 @@ class BiPoVariables(strax.LoopPlugin): are given for the five S2. """ - __version__ = '3.0.1' + __version__ = '4.0.0' depends_on = ('events', 'peak_basics', @@ -60,7 +60,6 @@ def setup(self): self.posrec_save = [(xy + algo, xy + algo) for xy in ['x_', 'y_'] for algo in self.pos_rec_labels] # ???? self.to_store = [name for name, _, _ in self.peak_properties] - def infer_dtype(self): # Basic event properties @@ -90,70 +89,58 @@ def infer_dtype(self): np.float32, f'S2_{n} {algo}-reconstructed Y position, uncorrected [cm]')] return basics_dtype + + @staticmethod + def set_nan_defaults(buffer): + """ + When constructing the dtype, take extra care to set values to + np.Nan / -1 (for ints) as 0 might have a meaning + """ + for field in buffer.dtype.names: + if np.issubdtype(buffer.dtype[field], np.integer): + buffer[field][:] = -1 + else: + buffer[field][:] = np.nan + + @staticmethod + def get_largest_sx_peaks(peaks, + s_i, + number_of_peaks=2): + """Get the largest S1/S2. For S1s allow a min coincidence and max time""" + # Find all peaks of this type (S1 or S2) + + s_mask = peaks['type'] == s_i + selected_peaks = peaks[s_mask] + s_index = np.arange(len(peaks))[s_mask] + largest_peaks = np.argsort(selected_peaks['area'])[-number_of_peaks:][::-1] + return selected_peaks[largest_peaks], s_index[largest_peaks] - def compute_loop(self, event, peaks): - - result = dict(time=event['time'], - endtime=strax.endtime(event)) - result['n_peaks'] = len(peaks) - if not len(peaks): - return result - - ######## - # S1 # - ######## - - mask_s1s = (peaks['type']==1) - mask_s1s &= (peaks['area']>100) + def compute(self, events, peaks): - if not len(peaks[mask_s1s]): - return result - - ## Save the biggest peaks - max_s1s = min(self.max_n_s1, len(peaks[mask_s1s])) - - # Need to initialize them to be able to use them in S2 mask without errors - result['s1_time_0'], result['s1_time_1'] = float('nan'), float('nan') - for i, p in enumerate(reversed(np.sort(peaks[mask_s1s], order='area'))): - - for prop in self.to_store: - result[f's1_{prop}_{i}'] = p[prop] - if i == self.max_n_s1 - 1: - break + result = np.zeros(len(events), dtype=self.dtype) + self.set_nan_defaults(result) - result['n_incl_peaks_s1'] = max_s1s + split_peaks = strax.split_by_containment(peaks, events) - ######## - ## S2 # - ######## - - # TODO - # all this mask thingis should me moved to the next plugin - # you can have a minimal one but should be very basic - # the complicated stuff should be in the matching plugin - # and this one can be used generally (not only for bipos) - # same for the S1s - - mask_s2s = peaks['type']==2 - mask_s2s &= peaks['area'] > 1500 # low area limit - mask_s2s &= peaks['area_fraction_top'] > 0.5 # to remove S1 afterpulses - mask_s2s &= np.abs(peaks['time'] - result['s1_time_0']) > 1000 # again to remove afterpulses - mask_s2s &= np.abs(peaks['time'] - result['s1_time_1']) > 1000 # and again to remove afterpulses - mask_s2s &= peaks['time']-result['s1_time_0'] < 5000000 # 5000mus, S2 is too far in time, not related to Po - - if not len(peaks[mask_s2s]): - return result - - max_s2s = min(self.max_n_s2, len(peaks[mask_s2s])) - for i, p in enumerate(reversed(np.sort(peaks[mask_s2s], order='area'))): - for prop in self.to_store: - result[f's2_{prop}_{i}'] = p[prop] - for name_alg in self.posrec_save: - result[f's2_{name_alg[0]}_{i}'] = p[name_alg[1]] - if i == self.max_n_s2 - 1: - break - - result['n_incl_peaks_s2'] = max_s2s + result['time'] = events['time'] + result['endtime'] = events['endtime'] + + for event_i, _ in enumerate(events): + + peaks_in_event_i = split_peaks[event_i] + + largest_s1s, s1_idx = self.get_largest_sx_peaks(peaks_in_event_i, s_i=1, number_of_peaks=self.max_n_s1) + largest_s2s, s2_idx = self.get_largest_sx_peaks(peaks_in_event_i, s_i=2, number_of_peaks=self.max_n_s2) + + for i, p in enumerate(largest_s1s): + for prop in self.to_store: + result[event_i][f's1_{prop}_{i}'] = p[prop] + + for i, p in enumerate(largest_s2s): + for prop in self.to_store: + result[event_i][f's2_{prop}_{i}'] = p[prop] + for name_alg in self.posrec_save: + result[event_i][f's2_{name_alg[0]}_{i}'] = p[name_alg[1]] - return result + return result \ No newline at end of file From 923d74a9efcce3e82a070cd16fb73cbbbfeb7b3d Mon Sep 17 00:00:00 2001 From: cfuselli Date: Tue, 2 May 2023 07:51:27 +0000 Subject: [PATCH 04/15] rename plugin --- straxen/plugins/events/event_basics.py | 2 -- .../{event_bipo_variables.py => event_basics_multi.py} | 6 +++--- straxen/plugins/events/event_bipo_basics.py | 2 +- straxen/plugins/events/event_bipo_matching.py | 2 +- 4 files changed, 5 insertions(+), 7 deletions(-) rename straxen/plugins/events/{event_bipo_variables.py => event_basics_multi.py} (98%) diff --git a/straxen/plugins/events/event_basics.py b/straxen/plugins/events/event_basics.py index 179d104c1..f3983664e 100644 --- a/straxen/plugins/events/event_basics.py +++ b/straxen/plugins/events/event_basics.py @@ -3,10 +3,8 @@ import numba import straxen - export, __all__ = strax.exporter() - @export class EventBasics(strax.Plugin): """ diff --git a/straxen/plugins/events/event_bipo_variables.py b/straxen/plugins/events/event_basics_multi.py similarity index 98% rename from straxen/plugins/events/event_bipo_variables.py rename to straxen/plugins/events/event_basics_multi.py index c6334aff3..0bf04a118 100644 --- a/straxen/plugins/events/event_bipo_variables.py +++ b/straxen/plugins/events/event_basics_multi.py @@ -7,7 +7,7 @@ export, __all__ = strax.exporter() @export -class BiPoVariables(strax.Plugin): +class EventBasicsMulti(strax.Plugin): """ Compute: - peak properties @@ -15,7 +15,7 @@ class BiPoVariables(strax.Plugin): of the first three main (in area) S1 and ten S2. The standard PosRec algorithm and the three different PosRec algorithms (mlp, gcn, cnn) - are given for the five S2. + are given for the S2s. """ __version__ = '4.0.0' @@ -26,7 +26,7 @@ class BiPoVariables(strax.Plugin): 'peak_proximity') # TODO change name - provides = 'bi_po_variables' + provides = 'event_basics_multi' data_kind = 'events' loop_over = 'events' diff --git a/straxen/plugins/events/event_bipo_basics.py b/straxen/plugins/events/event_bipo_basics.py index 25e52b813..e4c0609e2 100644 --- a/straxen/plugins/events/event_bipo_basics.py +++ b/straxen/plugins/events/event_bipo_basics.py @@ -15,7 +15,7 @@ class EventBiPoBasics(straxen.EventBasics): __version__ = '0.0.3' depends_on = ('events', - 'bi_po_variables', + 'event_basics_multi', 'bi_po_214_matching', 'peak_basics', 'peak_positions', diff --git a/straxen/plugins/events/event_bipo_matching.py b/straxen/plugins/events/event_bipo_matching.py index 73f9d8f6b..4160e7cdd 100644 --- a/straxen/plugins/events/event_bipo_matching.py +++ b/straxen/plugins/events/event_bipo_matching.py @@ -24,7 +24,7 @@ class BiPo214Matching(strax.Plugin): Time tolerance window in ns to match S2 signals to BiPo214 events. Default is 3000 ns. """ - depends_on=('bi_po_variables', ) + depends_on=('event_basics_multi', ) provides = 'bi_po_214_matching' __version__ = "2.1.7" From 304c8782ac088a92f4fc9d062490953b62705758 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Tue, 2 May 2023 07:52:51 +0000 Subject: [PATCH 05/15] rename plugin --- straxen/plugins/events/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/straxen/plugins/events/__init__.py b/straxen/plugins/events/__init__.py index ea1e15960..8cba45d84 100644 --- a/straxen/plugins/events/__init__.py +++ b/straxen/plugins/events/__init__.py @@ -61,8 +61,8 @@ from . import multi_scatter from .multi_scatter import * -from . import event_bipo_variables -from .event_bipo_variables import * +from . import event_basics_multi +from .event_basics_multi import * from . import event_bipo_matching from .event_bipo_matching import * From 1365d99ee2c1d30541e2f22e4672f85cac4c65b9 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Tue, 2 May 2023 15:33:47 +0000 Subject: [PATCH 06/15] shape fit --- straxen/plugins/events/__init__.py | 3 + .../events/event_basics_multi_s2_shape_fit.py | 133 ++++++++++++++++++ 2 files changed, 136 insertions(+) create mode 100644 straxen/plugins/events/event_basics_multi_s2_shape_fit.py diff --git a/straxen/plugins/events/__init__.py b/straxen/plugins/events/__init__.py index 8cba45d84..7f132ed4c 100644 --- a/straxen/plugins/events/__init__.py +++ b/straxen/plugins/events/__init__.py @@ -64,5 +64,8 @@ from . import event_basics_multi from .event_basics_multi import * +from . import event_basics_multi_s2_shape_fit +from .event_basics_multi_s2_shape_fit import * + from . import event_bipo_matching from .event_bipo_matching import * diff --git a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py new file mode 100644 index 000000000..a61523cd4 --- /dev/null +++ b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py @@ -0,0 +1,133 @@ +import strax +import numpy as np +import numba +import straxen +from scipy.optimize import curve_fit + +export, __all__ = strax.exporter() + +@export +class EventBasicsMultiS2ShapeFit(strax.Plugin): + """ + Compute: + + """ + + __version__ = '4.0.0' + + depends_on = ('events', + 'peaks') + + # TODO change name + provides = 'event_basics_multi_s2_shape_fit' + data_kind = 'events' + + max_n_s1 = straxen.URLConfig(default=3, infer_type=False, + help='Number of S1s to consider') + + max_n_s2 = straxen.URLConfig(default=10, infer_type=False, + help='Number of S2s to consider') + + peak_properties = ( + # name dtype comment + ('chi2', np.float32, 'chi2'), + ('mean', np.float32, 'mean'), + ('sigma', np.float32, 'sigma'), + ('area', np.float32, 'area of fitted gaussian') + ) + + + def setup(self): + + self.posrec_save = [(xy + algo, xy + algo) for xy in ['x_', 'y_'] for algo in self.pos_rec_labels] # ???? + self.to_store = [name for name, _, _ in self.peak_properties] + + def infer_dtype(self): + + # Basic event properties + basics_dtype = [] + basics_dtype += strax.time_fields + + + + # For S2s + p_type == 2: + max_n = self.max_n_s2 + for n in range(max_n): + # Peak properties + for name, dt, comment in self.peak_properties: + basics_dtype += [(f's{p_type}_fit_{name}_{n}', dt, f'S{p_type}_{n} fit {comment}'), ] + + return basics_dtype + + @staticmethod + def set_nan_defaults(buffer): + """ + When constructing the dtype, take extra care to set values to + np.Nan / -1 (for ints) as 0 might have a meaning + """ + for field in buffer.dtype.names: + if np.issubdtype(buffer.dtype[field], np.integer): + buffer[field][:] = -1 + else: + buffer[field][:] = np.nan + + @staticmethod + def get_largest_sx_peaks(peaks, + s_i, + number_of_peaks=2): + """Get the largest S1/S2. For S1s allow a min coincidence and max time""" + # Find all peaks of this type (S1 or S2) + + s_mask = peaks['type'] == s_i + selected_peaks = peaks[s_mask] + s_index = np.arange(len(peaks))[s_mask] + largest_peaks = np.argsort(selected_peaks['area'])[-number_of_peaks:][::-1] + return selected_peaks[largest_peaks], s_index[largest_peaks] + + @staticmethod + def gaussian(x, a, x0, sigma): + return a * np.exp(-(x - x0)**2 / (2 * sigma**2)) + + def compute(self, events, peaks): + + result = np.zeros(len(events), dtype=self.dtype) + self.set_nan_defaults(result) + split_peaks = strax.split_by_containment(peaks, events) + + result['time'] = events['time'] + result['endtime'] = events['endtime'] + + for event_i, _ in enumerate(events): + + peaks_in_event_i = split_peaks[event_i] + + largest_s2s, s2_idx = self.get_largest_sx_peaks(peaks_in_event_i, s_i=2, number_of_peaks=self.max_n_s2) + + for i, p in enumerate(largest_s2s): + + # Define the data to be fit + x = np.arange(200)[:p['length']] + y = p['data'][:p['length']]/max(p['data']) # assuming you want to fit the first peak + + # Fit the data with one Gaussian + try: + popt, pcov = curve_fit(self.gaussian, x, y, p0=[.1, 80, 10]) + except Exception as e: + popt = [0,0,0] + + # Calculate the chi-square goodness of fit statistic + residuals = y - self.gaussian(x, *popt) + chi2 = np.sum(residuals**2) + chi2_red = chi2 / (len(x) - len(popt)) + + + result[event_i][f's2_fit_chi2_{i}'] = chi2_red + result[event_i][f's2_fit_mean_{i}'] = popt[1] + result[event_i][f's2_fit_sigma_{i}'] = popt[2] + result[event_i][f's2_fit_area_{i}'] = max(p['data'])*popt[0]*sigma/(1/np.sqrt(2*np.pi)) + + return result + + + From 6a95eeae7d3ef18524e4c364d17441db83188911 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Tue, 2 May 2023 15:35:04 +0000 Subject: [PATCH 07/15] shape fit --- straxen/plugins/events/event_basics_multi_s2_shape_fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py index a61523cd4..3b7073256 100644 --- a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py +++ b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py @@ -51,7 +51,7 @@ def infer_dtype(self): # For S2s - p_type == 2: + p_type == 2 max_n = self.max_n_s2 for n in range(max_n): # Peak properties From 9384ea38c08d1ed002ef6c1c1e67a3858d5a6d0d Mon Sep 17 00:00:00 2001 From: cfuselli Date: Tue, 2 May 2023 15:37:24 +0000 Subject: [PATCH 08/15] shape fit --- straxen/plugins/events/event_basics_multi_s2_shape_fit.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py index 3b7073256..3e397edee 100644 --- a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py +++ b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py @@ -48,10 +48,8 @@ def infer_dtype(self): basics_dtype = [] basics_dtype += strax.time_fields - - # For S2s - p_type == 2 + p_type = 2 max_n = self.max_n_s2 for n in range(max_n): # Peak properties From f631780894bf35dbc9d335b80c2ae6e8cf6561ef Mon Sep 17 00:00:00 2001 From: cfuselli Date: Tue, 2 May 2023 15:38:50 +0000 Subject: [PATCH 09/15] shape fit --- straxen/plugins/events/event_basics_multi_s2_shape_fit.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py index 3e397edee..510f356b4 100644 --- a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py +++ b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py @@ -37,11 +37,6 @@ class EventBasicsMultiS2ShapeFit(strax.Plugin): ) - def setup(self): - - self.posrec_save = [(xy + algo, xy + algo) for xy in ['x_', 'y_'] for algo in self.pos_rec_labels] # ???? - self.to_store = [name for name, _, _ in self.peak_properties] - def infer_dtype(self): # Basic event properties From e6aa7bfbbf341e71ff9b13f09fc8e3126bc4f9fc Mon Sep 17 00:00:00 2001 From: cfuselli Date: Tue, 2 May 2023 15:40:51 +0000 Subject: [PATCH 10/15] shape fit --- straxen/plugins/events/event_basics_multi_s2_shape_fit.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py index 510f356b4..dbd9ef701 100644 --- a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py +++ b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py @@ -31,6 +31,7 @@ class EventBasicsMultiS2ShapeFit(strax.Plugin): peak_properties = ( # name dtype comment ('chi2', np.float32, 'chi2'), + ('ampl', np.float32, 'ampl'), ('mean', np.float32, 'mean'), ('sigma', np.float32, 'sigma'), ('area', np.float32, 'area of fitted gaussian') @@ -116,9 +117,10 @@ def compute(self, events, peaks): result[event_i][f's2_fit_chi2_{i}'] = chi2_red + result[event_i][f's2_fit_ampl_{i}'] = popt[0]*max(p['data']) result[event_i][f's2_fit_mean_{i}'] = popt[1] result[event_i][f's2_fit_sigma_{i}'] = popt[2] - result[event_i][f's2_fit_area_{i}'] = max(p['data'])*popt[0]*sigma/(1/np.sqrt(2*np.pi)) + result[event_i][f's2_fit_area_{i}'] = max(p['data'])*popt[0]*popt[2]/(1/np.sqrt(2*np.pi)) return result From 08ef1b9b279dc149f9f819aef5abce75613c0fc0 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Tue, 2 May 2023 15:53:49 +0000 Subject: [PATCH 11/15] shape fit --- straxen/plugins/events/event_basics_multi_s2_shape_fit.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py index dbd9ef701..9487b5c59 100644 --- a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py +++ b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py @@ -13,7 +13,7 @@ class EventBasicsMultiS2ShapeFit(strax.Plugin): """ - __version__ = '4.0.0' + __version__ = '5.0.0' depends_on = ('events', 'peaks') @@ -31,6 +31,7 @@ class EventBasicsMultiS2ShapeFit(strax.Plugin): peak_properties = ( # name dtype comment ('chi2', np.float32, 'chi2'), + ('dt',np.int,'dt of peak'), ('ampl', np.float32, 'ampl'), ('mean', np.float32, 'mean'), ('sigma', np.float32, 'sigma'), @@ -120,6 +121,7 @@ def compute(self, events, peaks): result[event_i][f's2_fit_ampl_{i}'] = popt[0]*max(p['data']) result[event_i][f's2_fit_mean_{i}'] = popt[1] result[event_i][f's2_fit_sigma_{i}'] = popt[2] + result[event_i][f's2_fit_dt_{i}'] = p['dt'] result[event_i][f's2_fit_area_{i}'] = max(p['data'])*popt[0]*popt[2]/(1/np.sqrt(2*np.pi)) return result From 8c1f60c11bc36283dfc47843b151d19bb7a19822 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Wed, 3 May 2023 07:52:26 +0000 Subject: [PATCH 12/15] improve fit --- straxen/plugins/events/event_basics_multi.py | 6 ++---- .../plugins/events/event_basics_multi_s2_shape_fit.py | 11 ++++++++--- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/straxen/plugins/events/event_basics_multi.py b/straxen/plugins/events/event_basics_multi.py index 0bf04a118..866f8e45c 100644 --- a/straxen/plugins/events/event_basics_multi.py +++ b/straxen/plugins/events/event_basics_multi.py @@ -18,7 +18,7 @@ class EventBasicsMulti(strax.Plugin): are given for the S2s. """ - __version__ = '4.0.0' + __version__ = '5.0.0' depends_on = ('events', 'peak_basics', @@ -65,9 +65,7 @@ def infer_dtype(self): # Basic event properties basics_dtype = [] basics_dtype += strax.time_fields - basics_dtype += [('n_peaks', np.int32, 'Number of peaks in the event'), - ('n_incl_peaks_s1', np.int32, 'Number of included S1 peaks in the event'), - ('n_incl_peaks_s2', np.int32, 'Number of included S2 peaks in the event')] + basics_dtype += [('n_peaks', np.int32, 'Number of peaks in the event'),] # For S1s and S2s for p_type in [1, 2]: diff --git a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py index 9487b5c59..14b02a161 100644 --- a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py +++ b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py @@ -13,7 +13,7 @@ class EventBasicsMultiS2ShapeFit(strax.Plugin): """ - __version__ = '5.0.0' + __version__ = '6.0.0' depends_on = ('events', 'peaks') @@ -105,18 +105,23 @@ def compute(self, events, peaks): x = np.arange(200)[:p['length']] y = p['data'][:p['length']]/max(p['data']) # assuming you want to fit the first peak + _p = x[y > np.exp(-0.5)*y.max()] + guess_sigma = 0.5*(_p.max() - _p.min()) + + p0 = [1, np.argmax(y), guess_sigma] + # Fit the data with one Gaussian try: - popt, pcov = curve_fit(self.gaussian, x, y, p0=[.1, 80, 10]) + popt, pcov = curve_fit(self.gaussian, x, y, p0=p0) except Exception as e: popt = [0,0,0] # Calculate the chi-square goodness of fit statistic + popt[2] = np.abs(popt[2]) residuals = y - self.gaussian(x, *popt) chi2 = np.sum(residuals**2) chi2_red = chi2 / (len(x) - len(popt)) - result[event_i][f's2_fit_chi2_{i}'] = chi2_red result[event_i][f's2_fit_ampl_{i}'] = popt[0]*max(p['data']) result[event_i][f's2_fit_mean_{i}'] = popt[1] From c6e97622c571717e601bc61943c684138d2f42f8 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Wed, 3 May 2023 12:55:56 +0000 Subject: [PATCH 13/15] shape only bipo --- .../events/event_basics_multi_s2_shape_fit.py | 61 ++++++++++--------- 1 file changed, 32 insertions(+), 29 deletions(-) diff --git a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py index 14b02a161..ca50e73b3 100644 --- a/straxen/plugins/events/event_basics_multi_s2_shape_fit.py +++ b/straxen/plugins/events/event_basics_multi_s2_shape_fit.py @@ -13,9 +13,10 @@ class EventBasicsMultiS2ShapeFit(strax.Plugin): """ - __version__ = '6.0.0' + __version__ = '7.0.0' depends_on = ('events', + 'bi_po_214_matching', 'peaks') # TODO change name @@ -93,41 +94,43 @@ def compute(self, events, peaks): result['time'] = events['time'] result['endtime'] = events['endtime'] - for event_i, _ in enumerate(events): + for event_i, e in enumerate(events): - peaks_in_event_i = split_peaks[event_i] + if e['s2_bi_match']>-5: - largest_s2s, s2_idx = self.get_largest_sx_peaks(peaks_in_event_i, s_i=2, number_of_peaks=self.max_n_s2) + peaks_in_event_i = split_peaks[event_i] - for i, p in enumerate(largest_s2s): + largest_s2s, s2_idx = self.get_largest_sx_peaks(peaks_in_event_i, s_i=2, number_of_peaks=self.max_n_s2) - # Define the data to be fit - x = np.arange(200)[:p['length']] - y = p['data'][:p['length']]/max(p['data']) # assuming you want to fit the first peak + for i, p in enumerate(largest_s2s): - _p = x[y > np.exp(-0.5)*y.max()] - guess_sigma = 0.5*(_p.max() - _p.min()) + # Define the data to be fit + x = np.arange(200)[:p['length']] + y = p['data'][:p['length']]/max(p['data']) # assuming you want to fit the first peak - p0 = [1, np.argmax(y), guess_sigma] + _p = x[y > np.exp(-0.5)*y.max()] + guess_sigma = 0.5*(_p.max() - _p.min()) - # Fit the data with one Gaussian - try: - popt, pcov = curve_fit(self.gaussian, x, y, p0=p0) - except Exception as e: - popt = [0,0,0] - - # Calculate the chi-square goodness of fit statistic - popt[2] = np.abs(popt[2]) - residuals = y - self.gaussian(x, *popt) - chi2 = np.sum(residuals**2) - chi2_red = chi2 / (len(x) - len(popt)) - - result[event_i][f's2_fit_chi2_{i}'] = chi2_red - result[event_i][f's2_fit_ampl_{i}'] = popt[0]*max(p['data']) - result[event_i][f's2_fit_mean_{i}'] = popt[1] - result[event_i][f's2_fit_sigma_{i}'] = popt[2] - result[event_i][f's2_fit_dt_{i}'] = p['dt'] - result[event_i][f's2_fit_area_{i}'] = max(p['data'])*popt[0]*popt[2]/(1/np.sqrt(2*np.pi)) + p0 = [1, np.argmax(y), guess_sigma] + + # Fit the data with one Gaussian + try: + popt, pcov = curve_fit(self.gaussian, x, y, p0=p0) + except Exception as e: + popt = [0,0,0] + + # Calculate the chi-square goodness of fit statistic + popt[2] = np.abs(popt[2]) + residuals = y - self.gaussian(x, *popt) + chi2 = np.sum(residuals**2) + chi2_red = chi2 / (len(x) - len(popt)) + + result[event_i][f's2_fit_chi2_{i}'] = chi2_red + result[event_i][f's2_fit_ampl_{i}'] = popt[0]*max(p['data']) + result[event_i][f's2_fit_mean_{i}'] = popt[1] + result[event_i][f's2_fit_sigma_{i}'] = popt[2] + result[event_i][f's2_fit_dt_{i}'] = p['dt'] + result[event_i][f's2_fit_area_{i}'] = max(p['data'])*popt[0]*popt[2]/(1/np.sqrt(2*np.pi)) return result From caf260d5511069e8623ac4abb518769058c411cb Mon Sep 17 00:00:00 2001 From: cfuselli Date: Fri, 19 May 2023 02:48:14 -0500 Subject: [PATCH 14/15] fix area correction Po --- straxen/plugins/events/event_bipo_basics.py | 30 +++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/straxen/plugins/events/event_bipo_basics.py b/straxen/plugins/events/event_bipo_basics.py index e4c0609e2..e2fd654bf 100644 --- a/straxen/plugins/events/event_bipo_basics.py +++ b/straxen/plugins/events/event_bipo_basics.py @@ -77,4 +77,34 @@ def fill_result_i(self, event, peaks, bi_i, po_i): field_names, peak_properties_to_save) + @staticmethod + @numba.njit + def set_event_properties(result, largest_s1s, largest_s2s, peaks): + """Get properties like drift time and area before main S2""" + # Compute drift times only if we have a valid S1-S2 pair + if len(largest_s1s) > 0 and len(largest_s2s) > 0: + result['drift_time'] = largest_s2s[0]['center_time'] - largest_s1s[0]['center_time'] + + # Correcting alt S1 and S2 based on BiPo + + if len(largest_s1s) > 1: + result['alt_s1_interaction_drift_time'] = largest_s2s[1]['center_time'] - largest_s1s[1]['center_time'] + result['alt_s1_delay'] = largest_s1s[1]['center_time'] - largest_s1s[0]['center_time'] + if len(largest_s2s) > 1: + result['alt_s2_interaction_drift_time'] = largest_s2s[1]['center_time'] - largest_s1s[1]['center_time'] + result['alt_s2_delay'] = largest_s2s[1]['center_time'] - largest_s2s[0]['center_time'] + + # areas before main S2 + if len(largest_s2s): + peaks_before_ms2 = peaks[peaks['time'] < largest_s2s[0]['time']] + result['area_before_main_s2'] = np.sum(peaks_before_ms2['area']) + + s2peaks_before_ms2 = peaks_before_ms2[peaks_before_ms2['type'] == 2] + if len(s2peaks_before_ms2) == 0: + result['large_s2_before_main_s2'] = 0 + else: + result['large_s2_before_main_s2'] = np.max(s2peaks_before_ms2['area']) + return result + + From f2ce2722ac6bbbfedbb0ef9aba40886667d885e1 Mon Sep 17 00:00:00 2001 From: cfuselli Date: Sat, 20 May 2023 02:03:32 -0500 Subject: [PATCH 15/15] bump ev basics --- straxen/plugins/events/event_bipo_basics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/straxen/plugins/events/event_bipo_basics.py b/straxen/plugins/events/event_bipo_basics.py index e2fd654bf..115bfebe0 100644 --- a/straxen/plugins/events/event_bipo_basics.py +++ b/straxen/plugins/events/event_bipo_basics.py @@ -12,7 +12,7 @@ class EventBiPoBasics(straxen.EventBasics): Carlo explain please """ - __version__ = '0.0.3' + __version__ = '1.0.0' depends_on = ('events', 'event_basics_multi',