diff --git a/hit_finder.py b/hit_finder.py index 869e3a4..6e98bca 100644 --- a/hit_finder.py +++ b/hit_finder.py @@ -3,28 +3,46 @@ import lar_param as lar import numpy as np +import numba as nb def hit_search(data, view, daq_chan, start, dt_min, thr1, thr2, thr3): - if(cf.view_type[view] == "Collection"): - return hit_search_collection(data, view, daq_chan, start, dt_min, thr1, thr2) - elif(cf.view_type[view] == "Induction"): - return hit_search_induction(data, view, daq_chan, start, dt_min, thr3) - else : + ll = [] + if(cf.view_type[view] != "Collection" and cf.view_type[view] != "Induction"): print(cf.view_type[view], " is not recognized") sys.exit() + elif(cf.view_type[view] == "Collection"): + h_num, h_start, h_stop, h_charge_int, h_max_t, h_max_adc, h_min_t, h_min_adc = hit_search_collection_nb(data, view, daq_chan, start, dt_min, thr1, thr2) + elif(cf.view_type[view] == "Induction"): + h_num, h_start, h_stop, h_charge_int, h_max_t, h_max_adc, h_min_t, h_min_adc = hit_search_induction_nb(data, view, daq_chan, start, dt_min, thr3) + if(h_num == 0): return ll + for start, stop, charge_int, max_t, max_adc, min_t, min_adc in zip(h_start, h_stop, h_charge_int, h_max_t, h_max_adc, h_min_t, h_min_adc): + h = dc.hits(view, daq_chan, start, stop, charge_int, max_t, max_adc, min_t, min_adc) + ll.append(h) + if(len(ll) == h_num): break + return ll -def hit_search_induction(data, view, daq_chan, start, dt_min, thr): +@nb.njit('Tuple((int64,int16[:],int16[:],float64[:],int16[:],float64[:],int16[:],float64[:]))(float64[:],int64,int64,int64,int64,float64)') +def hit_search_induction_nb(data, view, daq_chan, start, dt_min, thr): """ very basic induction-like hit finder """ """ WARNING : CANNOT FIND OVERLAPPING HITS """ + + npts = len(data) - ll = [] + h_num = 0 # store number of found hits + # list of hits parameters to be returned by this numba function (see dc.hits) + h_start = np.zeros(npts,dtype=np.int16) + h_stop = np.zeros(npts,dtype=np.int16) + h_charge_int= np.zeros(npts) + h_max_t = np.zeros(npts,dtype=np.int16) + h_max_adc = np.zeros(npts) + h_min_t = np.zeros(npts,dtype=np.int16) + h_min_adc = np.zeros(npts) - npts = len(data) hitPosFlag = False hitNegFlag = False @@ -33,8 +51,7 @@ def hit_search_induction(data, view, daq_chan, start, dt_min, thr): posSamp = 0 negSamp = 0 - h = dc.hits(view, daq_chan, start, 0, 0., 0., 0., 0., 0.) - + h_start[h_num] = start while(i= thr and hitNegFlag == False): + val = data[i] + while(i < npts and val >= thr and hitNegFlag == False): val = data[i] it = i+start posSamp += 1 @@ -51,19 +69,19 @@ def hit_search_induction(data, view, daq_chan, start, dt_min, thr): if(hitPosFlag == False): hitPosFlag = True - h.start = it + h_start[h_num] = it - h.max_t = it - h.max_adc = val + h_max_t[h_num] = it + h_max_adc[h_num] = val """ update the maximum case """ - if(val > h.max_adc): - h.max_t = it - h.max_adc = val + if(val > h_max_adc[h_num]): + h_max_t[h_num] = it + h_max_adc[h_num] = val if(hitPosFlag): - h.charge_int += val + h_charge_int[h_num] += val i+=1 @@ -71,13 +89,15 @@ def hit_search_induction(data, view, daq_chan, start, dt_min, thr): hitPosFlag = False posSamp = 0 - while(i < npts and hitPosFlag and hitNegFlag == False and np.fabs(data[i]) < thr): - h.charge_int += np.fabs(val) + val = data[i] + while(i < npts and hitPosFlag and hitNegFlag == False and np.fabs(val) < thr): + h_charge_int[h_num] += np.fabs(val) i += 1 - + val = data[i] """ now the negative part """ - while(i < npts and hitPosFlag and data[i] < -1.*thr): + val = data[i] + while(i < npts and hitPosFlag and val < -1.*thr): val = data[i] it = i+start negSamp += 1 @@ -86,19 +106,19 @@ def hit_search_induction(data, view, daq_chan, start, dt_min, thr): if(hitNegFlag == False): hitNegFlag = True - h.min_t = it - h.min_adc = val + h_min_t[h_num] = it + h_min_adc[h_num] = val """ update the minimum case """ - if(val < h.min_adc): - h.min_t = it - h.min_adc = val + if(val < h_min_adc[h_num]): + h_min_t[h_num] = it + h_min_adc[h_num] = val if(hitNegFlag): - h.charge_int += -1.*val + h_charge_int[h_num] += -1.*val - h.stop = it + h_stop[h_num] = it i+=1 if(negSamp < dt_min): @@ -106,20 +126,28 @@ def hit_search_induction(data, view, daq_chan, start, dt_min, thr): negSamp = 0 if(hitPosFlag and hitNegFlag): - ll.append(h) + h_num += 1 break - - - return ll + return h_num, h_start, h_stop, h_charge_int, h_max_t, h_max_adc, h_min_t, h_min_adc -def hit_search_collection(data, view, daq_chan, start, dt_min, thr1, thr2): +@nb.njit('Tuple((int64,int16[:],int16[:],float64[:],int16[:],float64[:],int16[:],float64[:]))(float64[:],int64,int64,int64,int64,float64,float64)') +def hit_search_collection_nb(data, view, daq_chan, start, dt_min, thr1, thr2): """search hit-shape in a list of points""" """algorithm from qscan""" - - ll = [] npts = len(data) + + h_num = 0 # store number of found hits + # list of hits parameters to be returned by this numba function (see dc.hits) + h_start = np.zeros(npts,dtype=np.int16) + h_stop = np.zeros(npts,dtype=np.int16) + h_charge_int= np.zeros(npts) + h_max_t = np.zeros(npts,dtype=np.int16) + h_max_adc = np.zeros(npts) + h_min_t = np.zeros(npts,dtype=np.int16) + h_min_adc = np.zeros(npts) + hitFlag = False i=0 @@ -127,8 +155,6 @@ def hit_search_collection(data, view, daq_chan, start, dt_min, thr1, thr2): minSamp = -1 singleHit = True - - while(i= thr1): val = data[i] @@ -138,42 +164,59 @@ def hit_search_collection(data, view, daq_chan, start, dt_min, thr1, thr2): hitFlag = True singleHit = True - h = dc.hits(view,daq_chan,it,0,0.,it,val, 0., 0.) + h_start[h_num] = it + h_stop[h_num] = 0 + h_charge_int[h_num]= 0. + h_max_t[h_num] = it + h_max_adc[h_num] = val + h_min_t[h_num] = 0. + h_min_adc[h_num] = 0. + + #h = dc.hits(view,daq_chan,it,0,0.,it,val, 0., 0.) minSamp = -1 - if(it > h.max_t and val < h.max_adc - thr2 and (minSamp==-1 or minimum >= val)): + if(it > h_max_t[h_num] and val < h_max_adc[h_num] - thr2 and (minSamp==-1 or minimum >= val)): minSamp = it minimum = val - if(minSamp >= 0 and it > minSamp and val > minimum + thr2 and (it-h.start) >= dt_min): - h.stop = minSamp-1 - ll.append(h) + if(minSamp >= 0 and it > minSamp and val > minimum + thr2 and (it-h_start[h_num]) >= dt_min): + h_stop[h_num] = minSamp-1 + h_num += 1 hitFlag = True singleHit = False - h = dc.hits(view,daq_chan,minSamp,0,0,it,val, 0., 0.) + + h_start[h_num] = minSamp + h_stop[h_num] = 0 + h_charge_int[h_num]= 0. + h_max_t[h_num] = it + h_max_adc[h_num] = val + h_min_t[h_num] = 0. + h_min_adc[h_num] = 0. + minSamp = -1 - if(h.stop == 0 and val > h.max_adc): - h.max_t = it - h.max_adc = val + if(h_stop[h_num] == 0 and val > h_max_adc[h_num]): + h_max_t[h_num] = it + h_max_adc[h_num] = val if(minSamp >= 0): minSamp = -1 minimum = val - h.charge_int += val + h_charge_int[h_num] += val i+=1 if(hitFlag == True): hitFlag = False - h.stop = it-1 + h_stop[h_num] = it-1 - #if((singleHit and (h.stop-h.start >= dt_min)) or not singleHit): - if(h.stop-h.start >= dt_min): - ll.append(h) + #if((singleHit and (h_stop[h_num]-h_start[h_num] >= dt_min)) or not singleHit): + if(h_stop[h_num]-h_start[h_num] >= dt_min): + h_num += 1 i+=1 - return ll + return h_num, h_start, h_stop, h_charge_int, h_max_t, h_max_adc, h_min_t, h_min_adc + def recompute_hit_charge(hit): daq_ch, start, stop = hit.daq_channel, hit.start, hit.stop @@ -184,9 +227,6 @@ def recompute_hit_charge(hit): val += np.fabs(dc.data_daq[daq_ch, t] - mean) hit.charge_int = val - - - def find_hits(pad_left, pad_right, dt_min, n_sig_coll_1, n_sig_coll_2, n_sig_ind): """ get boolean roi based on mask and alive channels """ diff --git a/lardon.py b/lardon.py index fdfb9d8..6726847 100644 --- a/lardon.py +++ b/lardon.py @@ -134,9 +134,11 @@ #ped.set_mask_wf_rms_all() tp = time.time() + for i in range(2): ped.compute_pedestal(noise_type='filt') ped.update_mask(pars.ped_amp_sig_fst) + #plot.plot_noise_daqch(noise_type='filt',option='fft', vrange=pars.plt_noise_zrange) #plot.plot_noise_vch(noise_type='filt', vrange=pars.plt_noise_zrange,option='fft')#,to_be_shown=True) @@ -158,6 +160,7 @@ #ped.set_mask_wf_rms_all() tpm = time.time() + for i in range(2): ped.compute_pedestal(noise_type='filt') ped.update_mask(pars.ped_amp_sig_oth) diff --git a/pedestals.py b/pedestals.py index 651057a..2264f5b 100644 --- a/pedestals.py +++ b/pedestals.py @@ -59,6 +59,171 @@ def update_mask(thresh): dc.mask_daq = ne.evaluate( "where((abs(data) > thresh*rms) | ~alive_chan, 0, 1)", global_dict={'data':dc.data_daq, 'alive_chan':dc.alive_chan, 'rms':dc.evt_list[-1].noise_filt.ped_rms[:,None]}).astype(bool) +def refine_mask(pars,debug=False): + if(debug == True): import matplotlib.pyplot as plt + + mask = dc.mask_daq + for ch in range(len(mask)): + data = dc.data_daq[ch] + rms = dc.evt_list[-1].noise_filt.ped_rms[ch] + view = dc.chmap[ch].view + + if(view == 2): mask_collection_signal(ch,mask,data,pars.ped_rise_thr[view],pars.ped_ramp_thr[view],rms*pars.ped_amp_thr[view],pars.ped_dt_thr) + else: mask_induction_signal(ch,mask,data,pars.ped_rise_thr[view],pars.ped_ramp_thr[view],rms*pars.ped_amp_thr[view],pars.ped_dt_thr,pars.ped_zero_cross_thr) + if(debug==True and ch > 1 and (ch % 50) == 0): + nrows, ncols = 10, 5 + fig = plt.figure() + for it,channel in enumerate(range(ch-50,ch)): + title = 'v'+str(view)+'-ch'+str(channel) + ax = fig.add_subplot(nrows, ncols, it+1) + ax.set_title(title) + ax.plot(dc.data_daq[channel],'o',markersize=0.2) + ax.axhline(y=dc.evt_list[-1].noise_filt.ped_rms[channel]*pars.ped_amp_thr[dc.chmap[channel].view], color='gray') + ax.axhline(y=0, color='lightgray') + ax.axhline(y=-dc.evt_list[-1].noise_filt.ped_rms[channel]*pars.ped_amp_thr[dc.chmap[channel].view], color='gray') + ax.plot(np.max(dc.data_daq[channel, :])*dc.mask_daq[channel,:],linestyle='-',linewidth=1,color='r') + ax.set(xlabel=None,ylabel=None) + ax.axis('off') + plt.subplots_adjust(left=0.01, bottom=0.01, right=0.99, top=0.95, wspace=0.05, hspace=1) + plt.show() + + + dc.mask_daq = mask + +@nb.jit('(int64,boolean[:,:],float64[:],int64,int64,float64,int64)',nopython = True) +def mask_collection_signal(ch,mask,data,rise_thr,ramp_thr,amp_thr,dt_thr): + + tmp_start = 0 # store candidate signal start to be registered when stop is found + trig_pos = 0 # store position at which sample exceed amp_thrs*ped_rms - 0 otherwise + find_pos = 0 # counter for consecutive positive sample + ramp = 0 # counter 3 consecutive rising samples + ongoing = 0 # flag triggered when signal start found + oldval1 = 0 # store previous sample + oldval2 = 0 # store pre-previous sample + zero_cross_check = 0 + + for x,val in np.ndenumerate(data): + + # set positive signal trigger + if(val >= amp_thr): trig_pos = x[0] + # set back positive trigger to zero if a certain period is over + elif(x[0] - trig_pos > dt_thr): trig_pos = 0 + + # if rising edge found - register begining of signal + if(trig_pos > 0 and find_pos >= rise_thr and ramp > ramp_thr and ongoing == 0): + tmp_start = x[0] - find_pos + ongoing = 1 + + # signal found - look when amp. reaches noise floor (neg. value) + if(find_pos >= rise_thr and ongoing == 1 and val < 0): + start = tmp_start + stop = x[0] + for it in range(start,stop): mask[ch][it] = 0 + + # reset counters + zero_cross_check = 0 + find_pos = 0 + ongoing = 0 + trig_pos = 0 + + if(val > 0): + # start counter for positive samples + find_pos += 1 + # identify ramping + if(val > oldval1 and val > oldval2): ramp += 1 + else: ramp -= 1 + else: + zero_cross_check += 1 + + # reset counters when negative sample is found (deep in the noise) for col. + if(zero_cross_check >= 1): + zero_cross_check = 0 + find_pos = 0 + ramp = 0 + ongoing = 0 + + oldval2 = oldval1 + oldval1 = val + +@nb.jit('(int64,boolean[:,:],float64[:],int64,int64,float64,int64,int64)',nopython = True) +def mask_induction_signal(ch,mask,data,rise_thr,ramp_thr,amp_thr,dt_thr,zero_cross_thr): + + tmp_start = 0 # store candidate signal start to be registered when stop is found + trig_pos = 0 # store position at which sample exceed amp_thrs*ped_rms - 0 otherwise + trig_neg = 0 # store position at which sample exceed -amp_thrs*ped_rms - 0 otherwise + find_pos = 0 # counter for consecutive positive sample + find_neg = 0 # counter for consecutive negative sample + ramp = 0 # counter 3 consecutive rising samples + ongoing = 0 # flag triggered when signal start found + oldval1 = 0 # store previous sample + oldval2 = 0 # store pre-previous sample + inv = 0 # invert flag for induction signal + zero_cross_check = 0 + + for x,val in np.ndenumerate(data): + + # set positive signal trigger + if(val >= amp_thr): trig_pos = x[0] + # set back positive trigger to zero if a certain period is over + elif(x[0] - trig_pos > dt_thr): trig_pos = 0 + + # check for sample above threshold for induction signal start + if(trig_pos > 0 and find_pos > rise_thr and val > 0 and ongoing == 0): + tmp_start = x[0] - find_pos + ongoing = 1 + + if(trig_pos > 0 and find_pos >= rise_thr and ongoing == 1): + # set negative signal trigger + if(val <= -amp_thr): trig_neg = x[0] + elif(x[0] - trig_neg > dt_thr): trig_neg = 0 + + # set invert flag when induction crosses zero + if(val < 0 and inv == 0): + inv = 1 + x_zero_cross = x[0] + + # set induction signal stop once a positive sample is found - 5 samples after zero crossing required + if(val > 0 and inv == 1 and x[0] > x_zero_cross+zero_cross_thr): + + # require a minimal signal length/negative amplitude/duration to remove spikes + if(x[0] - tmp_start > dt_thr/4 and trig_pos - tmp_start > 0 and trig_neg >0 and find_neg >= rise_thr): + start= tmp_start + stop = x[0] + for it in range(start,stop): mask[ch][it] = 0 + + # reset counters + zero_cross_check= 0 + find_pos = 0 + find_neg = 0 + ongoing = 0 + trig_pos = 0 + trig_neg = 0 + inv = 0 + + if(val > 0 and inv == 0): + # start counter for positive samples + find_pos += 1 + # identify ramping + if(val > oldval1 and val > oldval2): ramp += 1 + else: ramp -= 1 + elif(val < 0 and inv == 1): + # start counter for negative samples + find_neg += 1 + else: + zero_cross_check += 1 + + # reset counters when negative/positive sample is found (deep in the noise) for col. and ind. resp. + if(zero_cross_check >= 5): + zero_cross_check = 0 + find_pos = 0 + find_neg = 0 + ramp = 0 + ongoing = 0 + inv = 0 + + oldval2 = oldval1 + oldval1 = val + def set_mask_wf_rms_all(): # set signal mask bins from threshold - simplest method