#!/usr/bin/env python3
'''
.. module:: TMoptfilter
:synopsis: Optimum filter in the frequency domain, using a predefined optimal filter
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import os
import numpy
import matplotlib
import matplotlib.pyplot as plt
from tesfdmtools.methods.HDF_Eventutils import getrisetime,getphase
from tesfdmtools.methods.IQrotate import IQphase,IQphrot
from tesfdmtools.utils.cu import cu
def write_optfilter(ofilter,hdf,channel,freq):
'''
write a file with the optimal filter data
Args:
* `ofilter` = numpy array with optimal filter data
* `hdf` = the HDF5 object with the event data
* `channel` = channel used inthe HDF5 data
* `freq` = frequency (pixel) data used in the HDF5 data
'''
bnam=os.path.basename(hdf.hdffile).rsplit('.',1)[0]
fnam='OptFilter_%s_ch%2.2d_px%2.2d.txt' % (bnam,channel,freq)
olst=open(fnam,'w')
olst.write("#Optimal filter data\n")
olst.write("#Inputfile: %s\n" % hdf.hdffile)
olst.write("#Channel: %d\n" % channel)
olst.write("#Freq: %d\n" % freq)
for dd in ofilter:
olst.write("%22.8e\n" % dd)
olst.close()
print("Optimal filter file written: ",fnam)
return
def read_optfilter(length,optname,channel,freq):
'''
read an optimal filter file
Args:
* `length` = length of the optimal filter data
* `optname` = filename of the optimal filter data
* `channel` = opt filter file for this channel
* `freq` = opt filter file for this frequency
Returns:
* `ofilter` = optimal filter data
'''
ofilter=numpy.zeros(length,dtype=float)
if os.path.exists(optname):
optfname=optname
else:
if optname.find('.txt') < 0:
optfname="%s_ch%2.2d_px%2.2d.txt" % (optname,channel,freq)
else:
kk=optname.find('_ch')
onam=optname[0:kk]
optfname="%s_ch%2.2d_px%2.2d.txt" % (onam,channel,freq)
print("Open optimal filter file: ",optfname)
olst=open(optfname,'r')
i=0
opt=False
for ll in olst:
line=ll.strip(' \n')
if line[0] == '#':
if line.find("Optimal") >= 0:
opt=True
else:
if not opt:
sys.exit("Given optimal filterfile is not an optimal filter")
if i >= length:
sys.exit("Given optimal filter file is too large")
ofilter[i]=float(line)
i=i+1
if i != length:
sys.exit("Given optimal filter file is too small")
return ofilter
[docs]def optfilter(hdf,channel,freq,indx,base1,base2,noisspec,debug=False,freqcutoff=None,\
rotate=False,risetime=False,bsec=0.05,absolute=False,prlen=None,apulsepos=None,\
wrtfilter=False,usefilter=None,flip=False,**kwargs):
'''
perform optimal filtering fit of individual pulses in the frequency domain
There is the option to write the optimal filter to file
and the option to use an external optimal filter
Args:
* `hdf` = HDF5 input file object
* `channel` = channel number being processed
* `freq` = frequency number (pixel) being processed
* `indx` = index of selected events, to be processed
* `base1` = baseline level at start of record
* `base2` = baseline level at end of record
* `noisspec` = noise spectrum
Kwargs:
* `debug` = if True, plot various fit parameters [default: False]
* `freqcutoff` = cutoff frequency (Hz) for optimal filtering [default: None]
* `rotate` = rotation angle (radians) for I/Q pulse record
* `risetime` = it True, compute rise time of pulses [default: None]
* `bsec` = section of record to take for background (only used for 'rotate' or 'absolute') [default: 0.05]
* `absolute` = use sqrt(I^2+Q^2) signal [default: False]
* `prlen` = length of record to process (None = entire record) [default: None]
* `apulsepos` = average position of pulse maximum; used for limited `prlen`. [default: None]
* `wrtfilter` = write optimal filter to file
* `usefilter` = use this external optimal filter
* `flip` = flip data in record
Returns:
* `ifit` = fitted pulse heigth parameters
* `rtimes` = computed rise times of fitted pulses
* `ftimes` = computed fall times of fitted pulses
* `avpulse` = average pulse profile with baseline subtracted
* `avbline` = average baseline for the average pulse prfile
'''
print("Optimal filter, channel=",channel,' freq=',freq)
aphase=0.0
if type(rotate) is bool:
if rotate:
aphase=getphase(hdf,channel,freq,indx,debug=debug)
print("Optimal fit phase rotation: ",aphase)
else:
if type(rotate) is float:
aphase=rotate
print("Optimal fit phase rotation: ",aphase)
sirecord=hdf.channel[channel].freq[freq][0][:,0] # read first record to get record parameters
recp1=0 # use entire record
recp2=sirecord.size
if ( prlen is not None ) and ( apulsepos is not None ): # limit record length ?
p1=int(apulsepos-100-bsec*sirecord.size)
if p1 > 0:
recp1=p1
p2=recp1+prlen
if p2 < sirecord.size:
recp2=p2
print("Optimal fit record length limited to: ",recp1,recp2)
sirecord=sirecord[recp1:recp2]
bl=int(bsec*sirecord.size) # section for background level
samplerate=float(cu(hdf.channel[channel].freq[freq].attrs['sample_rate'])) # get sample rate
if freqcutoff is None:
slen=sirecord.size//2 # max samplingfrequency
else:
slen=int(float(freqcutoff)/float(samplerate)*sirecord.size) # frequency cutoff
xax=numpy.arange(sirecord.size,dtype=float) # make x-axsis for record
iax=numpy.arange(sirecord.size,dtype=int)
bax=numpy.concatenate((iax[0:bl],iax[-bl:])) # baseline section
fpulses=numpy.zeros((indx.size,(slen-1)),dtype=float) # initialize array to store pulse power spectra
avpulse=numpy.zeros(sirecord.size,dtype=float) # average pulse profile
if risetime:
rtimes=numpy.empty(indx.size,dtype=float) # array of risetimes for events
ftimes=numpy.empty(indx.size,dtype=float) # array of fall times for events
else:
rtimes=None
ftimes=None
for i,irec in enumerate(indx): # go through all record in selection list
bas1=base1[irec]
bas2=base2[irec]
record=hdf.channel[channel].freq[freq][irec]
if absolute:
irecord=numpy.sqrt(record[:,0].astype(float)**2+record[:,1].astype(float)**2)
else:
if flip or ( record[0,0] < 0 ): # wrong rotation, rotate by pi
irecord=-record[:,0]
else:
irecord=record[:,0]
if aphase != 0.0:
irecord,qrecord=IQphrot(record[:,0],record[:,1],aphase)
irecord=irecord[recp1:recp2] # limit record length
if absolute or ( aphase != 0.0 ):
pp=numpy.polyfit(bax,irecord[bax],1)
cirecord=irecord-(pp[0]*xax+pp[1])
else:
cirecord=irecord-(bas1+xax/float(sirecord.size-1)*(bas2-bas1)) # retrieve I-signal and subtract background
avpulse=avpulse+irecord # accumulate average pulse
if risetime:
rt,ft=getrisetime(xax,cirecord,debug=debug) # rise time in seconds
rtimes[i]=rt/samplerate
ftimes[i]=ft/samplerate
fpulses[i,:]=numpy.absolute(numpy.fft.fft(cirecord)[1:slen]) # compute power spectrum
if ( i % 1000 ) == 0:
print("process record: ",irec)
print("number of pulses processed: ",indx.size)
avpulse=avpulse/float(indx.size) # average pulse
bindx=numpy.concatenate((numpy.arange(100,dtype=int),numpy.arange((avpulse.size-100),avpulse.size,dtype=int)))
abfit=numpy.polyfit(xax[bindx],avpulse[bindx],1)
avpulse=avpulse-(abfit[0]*xax+abfit[1])
minav=numpy.argmin(avpulse)
avbline=abfit[0]*minav+abfit[1] # average baseline
apulse=fpulses.sum(axis=0)/float(indx.size) # compute average power spectrum for all selected pulses
if wrtfilter: # write optimal filter data to file?
write_optfilter(apulse,hdf,channel,freq)
else:
if usefilter is not None: # use external optimal filter file?
apulse=read_optfilter(apulse.size,usefilter,channel,freq)
nns=(noisspec**2)[1:slen]
norm=numpy.sum(apulse**2/nns) # average pulse normalization
ifit=numpy.zeros(indx.size,dtype=float) # initialize array to store optimal filter fit parameters
for i in numpy.arange(indx.size):
ifit[i]=numpy.sum(fpulses[i]*apulse/nns)/norm # compute optimal filter fit
if debug:
if risetime:
np=3
else:
np=2
f, ax = plt.subplots(np)
nax=numpy.arange(indx.size)
ax[1].plot(nax,ifit,'b.')
ax[1].set_xlabel('event number')
ax[1].set_ylabel('Pulse intensity')
ax[0].set_title('Optimal filter fits')
ax[0].hist(ifit,bins=500,histtype='stepfilled')
ax[0].set_xlabel('Energy (arb. units)')
ax[0].set_ylabel('Counts')
if risetime:
ax[2].hist(rtimes*1000,bins=512,align='mid')
ax[2].set_xlabel("Rise time (ms)")
ax[2].set_ylabel("N")
ax[2].set_title("Pulse rise times")
plt.show()
plt.close('all')
return (ifit,rtimes,ftimes,avpulse,avbline)