Source code for tesfdmtools.methods.filterect

#!/usr/bin/env python3
#
'''
.. module:: filterect
   :synopsis: filter eventlist to exclude positive electric crosstalk pulses
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''

import numpy

from tesfdmtools.utils.hdfview.convolve import convolve

def filterect(hdf,channel,freq,indx,blines1,blines2):
   '''
   Filter eventlist to exclude records with positive electric crosstalk pulses

   Args:
      *     `hdf` = HDF input file object
      * `channel` = channel number being processed
      *    `freq` = frequency (pixel) being processed
      *    `indx` = index of valid event records
      * `blines1` = baseline levels at start of records
      * `blines2` = baseline levels at end of records 

   Returns:
      * `selindx` = new index of valid records without postive electric crosstalk pulses 
      *   `posxt` = list of filterparameters which sense the positive pulse
      *  `sellim` = exclusion threshold in parameter values 
   '''

   rec0=hdf.channel[channel].freq[freq][indx[0]][:,0]
   sample_rate=float(hdf.channel[channel].freq[freq][indx[0]].attrs['sample_rate'])
   rlen=rec0.size
   xx=numpy.arange(rlen,dtype=float)
   
   posxt=numpy.zeros(indx.size)
   for j,i in enumerate(indx):
      idata=hdf.channel[channel].freq[freq][i][:,0].astype(float)
      baseline=blines1[i]+xx/float(rlen)*(blines2[i]-blines1[i])
      diff=idata-baseline
      diffc=convolve(diff,500.0,sample_rate)			# convolve
      pindx,=numpy.where(diffc > 0)				# look at positive points only
      sdev=numpy.sqrt(numpy.sum(diffc[pindx]**2)/pindx.size)	# compute standard deviation
      pindx,=numpy.where(diffc > 5.0*sdev)			# select only points larger than 5x standard deviation
      if pindx.size == 0:
         posxt[j]=0.0
      else:
         posxt[j]=numpy.sum(diffc[pindx]**2)			# filterparameter
      if ( j % 1000 ) == 0:
         print("process electric crosstalk filter record: ",i)

   histxt,hbinb=numpy.histogram(posxt,bins=50)			# make histogram of filterparameter
   hbins=(hbinb[0:-1]+hbinb[1:])/2.0

   hmax=numpy.argmax(histxt)					# find maximum in histogram
   for i in numpy.arange(hmax,(histxt.size-1)):			# find local minimum, starting at maximum
      if ( histxt[i] < histxt[i+1] ) or (histxt[i] == 0):
         break
   sellim=hbins[i]
   selindx,=numpy.where(posxt < sellim)			# records to keep
   print("Filterect freq %d: %d number of event records discarded" % (freq,(indx.size-selindx.size)) )
   print("                   %d number of event records remain" % (selindx.size) )

   return (selindx,posxt,sellim)