diff --git a/straxen/plugins/events/__init__.py b/straxen/plugins/events/__init__.py index b71827041..9098ed043 100644 --- a/straxen/plugins/events/__init__.py +++ b/straxen/plugins/events/__init__.py @@ -13,6 +13,9 @@ from . import event_area_per_channel from .event_area_per_channel import * +from . import event_top_bottom_params +from .event_top_bottom_params import * + from . import event_basics from .event_basics import * diff --git a/straxen/plugins/events/event_top_bottom_params.py b/straxen/plugins/events/event_top_bottom_params.py new file mode 100644 index 000000000..ad96d9e08 --- /dev/null +++ b/straxen/plugins/events/event_top_bottom_params.py @@ -0,0 +1,94 @@ +import numpy as np +import strax +import straxen +from warnings import warn +export, __all__ = strax.exporter() + +@export +class EventTopBottomParams(strax.Plugin): + """ + Pluging that computes timing characteristics of top and bottom waveforms + based on waveforms stored at event level for main/alt S1/S2 + """ + depends_on = ('event_info', 'event_area_per_channel') + provides = 'event_top_bottom_params' + __version__ = '0.0.0' + def infer_dtype(self): + ## Populating data type information + infoline = {'s1': 'main S1', + 's2': 'main S2', + 'alt_s1': 'alternative S1', + 'alt_s2': 'alternative S2', + } + ev_info_fields = self.deps['event_info'].dtype.fields + dtype = [] + # populating APC and waveform samples + self.ptypes = ['s1', 's2', 'alt_s1', 'alt_s2'] + self.arrs = ['top', 'bot'] + for type_ in self.ptypes: + for arr_ in self.arrs: + dtype+=[((f'Central time for {infoline[type_]} for {arr_} PMTs [ ns ]', + f'{type_}_center_time_{arr_}'), + ev_info_fields[f'{type_}_center_time'][0])] + dtype+=[((f'Time between 10% and 50% area quantiles for {infoline[type_]} for {arr_} PMTs [ns]', + f'{type_}_rise_time_{arr_}'), + ev_info_fields[f'{type_}_rise_time'][0])] + dtype+=[((f'Width (in ns) of the central 50% area of the peak for {arr_} PMTs of {infoline[type_]}', + f'{type_}_range_50p_area_{arr_}'), + ev_info_fields[f'{type_}_range_50p_area'][0])] + dtype+=[((f'Width (in ns) of the central 90% area of the peak for {arr_} PMTs of {infoline[type_]}', + f'{type_}_range_90p_area_{arr_}'), + ev_info_fields[f'{type_}_range_90p_area'][0])] + dtype+=[((f'Difference between center times of top and bottom arrays for {infoline[type_]} [ ns ]', + f'{type_}_center_time_diff_top_bot'), + ev_info_fields[f'{type_}_center_time'][0])] + dtype += strax.time_fields + return dtype + + def compute(self, events): + result = np.zeros(events.shape, dtype=self.dtype) + result['time'], result['endtime'] = events['time'], strax.endtime(events) + peak_dtype = strax.peak_dtype(n_channels=straxen.n_tpc_pmts, digitize_top=False) + for type_ in self.ptypes: + for arr_ in self.arrs: + # in order to reuse the same definitions as in other parts, we create "fake peaks" + # based only on data from corresponding array + fpeaks_ = np.zeros(events.shape[0], dtype=peak_dtype) + if arr_ == 'top': + fpeaks_['data']=events[f'{type_}_data_top'] + fpeaks_['area']=events[f'{type_}_area']*events[f'{type_}_area_fraction_top'] + elif arr_ == 'bot': + fpeaks_['data']=(events[f'{type_}_data']-events[f'{type_}_data_top']) + fpeaks_['area']=events[f'{type_}_area']*(1.-events[f'{type_}_area_fraction_top']) + elif arr_ == 'tot': + # This one is ony + fpeaks_['data']=events[f'{type_}_data'] + fpeaks_['area']=events[f'{type_}_area'] + else: + raise RuntimeError(f'Received unknown array type : '+ arr_) + fpeaks_['length']=events[f'{type_}_length'] + fpeaks_['dt']=events[f'{type_}_dt'] + # computing central times + # note that here we ignore 1/2 sample length to be consistent with other definitions + with np.errstate(divide='ignore', invalid='ignore'): + recalc_ctime = np.sum(fpeaks_['data']*(np.arange(0, fpeaks_['data'].shape[1])), axis=1 ) + recalc_ctime/=fpeaks_['area'] + recalc_ctime*=fpeaks_['dt'] + recalc_ctime[~(fpeaks_['area']>0)]=0.0 + # setting central times in the same way as inside peak processing + mask=(fpeaks_['area']>0) + result[f'{type_}_center_time_{arr_}']=events[f'{type_}_time'] + result[f'{type_}_center_time_{arr_}'][mask]+=recalc_ctime[mask].astype(int) + # computing widths ## + # zero or undefined area peaks should have nans + strax.compute_widths(fpeaks_) + result[f'{type_}_rise_time_{arr_}'][:]=np.nan + result[f'{type_}_rise_time_{arr_}'][mask]= -fpeaks_['area_decile_from_midpoint'][mask][:, 1] + result[f'{type_}_range_50p_area_{arr_}'][:]=np.nan + result[f'{type_}_range_50p_area_{arr_}'][mask] = fpeaks_['width'][mask][:, 5] + result[f'{type_}_range_90p_area_{arr_}'][:]=np.nan + result[f'{type_}_range_90p_area_{arr_}'][mask] = fpeaks_['width'][mask][:, 9] + # Difference between center times of top and bottom arrays + result[f'{type_}_center_time_diff_top_bot'] = (result[f'{type_}_center_time_top'] - + result[f'{type_}_center_time_bot']) + return result diff --git a/straxen/plugins/peaks/__init__.py b/straxen/plugins/peaks/__init__.py index 4e87b15e4..42e712438 100644 --- a/straxen/plugins/peaks/__init__.py +++ b/straxen/plugins/peaks/__init__.py @@ -1,3 +1,6 @@ +from . import peak_top_bottom_params +from .peak_top_bottom_params import * + from . import peak_ambience from .peak_ambience import * diff --git a/straxen/plugins/peaks/peak_top_bottom_params.py b/straxen/plugins/peaks/peak_top_bottom_params.py new file mode 100644 index 000000000..3682c1b92 --- /dev/null +++ b/straxen/plugins/peaks/peak_top_bottom_params.py @@ -0,0 +1,78 @@ +import numpy as np +import strax +import straxen +from warnings import warn +export, __all__ = strax.exporter() + +@export +class PeakTopBottomParams(strax.Plugin): + """ + Pluging that computes timing characteristics of top and bottom waveforms + based on waveforms stored at peak level + """ + depends_on = ('peaks', 'peak_basics') + provides = 'peak_top_bottom_params' + __version__ = '0.0.0' + + def infer_dtype(self): + dtype = [] + peak_basics_fields = self.deps['peak_basics'].dtype.fields + self.arrs = ['top', 'bot'] + for arr_ in self.arrs: + dtype+=[((f'Central time for {arr_} PMTs [ ns ]', + f'center_time_{arr_}'), + peak_basics_fields['center_time'][0])] + dtype+=[((f'Time between 10% and 50% area quantiles for {arr_} PMTs [ns]', + f'rise_time_{arr_}'), + peak_basics_fields['rise_time'][0])] + dtype+=[((f'Width (in ns) of the central 50% area of the peak for {arr_} PMTs', + f'range_50p_area_{arr_}'), + peak_basics_fields['range_50p_area'][0])] + dtype+=[((f'Width (in ns) of the central 90% area of the peak for {arr_} PMTs', + f'range_90p_area_{arr_}'), + peak_basics_fields['range_90p_area'][0])] + dtype+=[((f'Difference between center times of top and bottom arrays [ ns ]', + f'center_time_diff_top_bot'), + peak_basics_fields[f'center_time'][0])] + dtype += strax.time_fields + return dtype + + def compute(self, peaks): + result = np.zeros(peaks.shape, dtype=self.dtype) + peak_dtype = self.deps['peaks'].dtype + for arr_ in self.arrs: + fpeaks_ = np.zeros(peaks.shape[0], dtype=peak_dtype) + if arr_ == 'top': + fpeaks_['data']=peaks['data_top'] + fpeaks_['area']=peaks['area']*peaks[f'area_fraction_top'] + elif arr_ == 'bot': + fpeaks_['data']=(peaks['data']-peaks['data_top']) + fpeaks_['area']=peaks['area']*(1.-peaks['area_fraction_top']) + elif arr_ == 'tot': + # This one is ony + fpeaks_['data']=peaks[f'{type_}_data'] + fpeaks_['area']=peaks[f'{type_}_area'] + else: + raise RuntimeError(f'Received unknown array type : '+ arr_) + fpeaks_['length']=peaks[f'length'] + fpeaks_['dt']=peaks[f'dt'] + mask=(fpeaks_['area']>0) + # computing center times + with np.errstate(divide='ignore', invalid='ignore'): + recalc_ctime = np.sum(fpeaks_['data']*(np.arange(0, fpeaks_['data'].shape[1])), axis=1 ) + recalc_ctime/=fpeaks_['area'] + recalc_ctime*=fpeaks_['dt'] + recalc_ctime[~mask]=0.0 + result[f'center_time_{arr_}']=peaks['time'] + result[f'center_time_{arr_}'][mask]+=recalc_ctime[mask].astype(int) + # computing widths times + strax.compute_widths(fpeaks_) + result[f'rise_time_{arr_}'][:]=np.nan + result[f'rise_time_{arr_}'][mask]= -fpeaks_['area_decile_from_midpoint'][mask][:, 1] + result[f'range_50p_area_{arr_}'][:]=np.nan + result[f'range_50p_area_{arr_}'][mask] = fpeaks_['width'][mask][:, 5] + result[f'range_90p_area_{arr_}'][:]=np.nan + result[f'range_90p_area_{arr_}'][mask] = fpeaks_['width'][mask][:, 9] + result['center_time_diff_top_bot'] = (result[f'center_time_top'] - result[f'center_time_bot']) + result['time'], result['endtime'] = peaks['time'], strax.endtime(peaks) + return result