#!/usr/bin/env python3
#
'''
.. module:: Moptfilter
:synopsis: Eventprofiles optimal filtering using multiple templates
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import os
import numpy
import matplotlib.pyplot as plt
from tesfdmtools.methods.IQrotate import IQphase,IQphrot
from tesfdmtools.utils.cu import cu
# ============================================================================================
def write_optfilter(templates,tnorms,hdf,channel,freq):
'''
write a file with the optimal filter data
Args:
* `templates` = numpy array with template filter data
* `tnorms` = average amplitudes for the filters
* `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='MOptFilters_%s_ch%2.2d_px%2.2d.txt' % (bnam,channel,freq)
olst=open(fnam,'w')
olst.write("#Optimal multiple filter data\n")
olst.write("#Inputfile: %s\n" % hdf.hdffile)
olst.write("#Channel: %d\n" % channel)
olst.write("#Freq: %d\n#\n" % freq)
fmt='%17.8e'
for i in numpy.arange(1,tnorms.size):
fmt=fmt+' %17.8e'
fmt=fmt+'\n'
olst.write("#Filter scaling identifiers: \n")
olst.write(fmt % tuple(tnorms))
olst.write("#Templates: \n")
for j in numpy.arange(templates[0,:].size):
olst.write(fmt % tuple(templates[:,j]))
olst.close()
print("Templates filter file written: ",fnam)
return
# ============================================================================================
[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,ntemplates=4,flip=False,**kwargs):
'''
perform optimal filtering fit of individual pulses using multiple templates
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 templates to file, for further analysis
* `ntemplates` = number of templates to use
* `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)
print(" number of templates to use: %d" % ntemplates)
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
amplitudes=numpy.empty(indx.size,dtype=float) # amplitudes of pulses
for i,irec in enumerate(indx): # go through all records 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
# amplitudes[i]=fpulses[i,:].max()
amplitudes[i]=fpulses[i,0] # power spectrum amplitude from first frequency
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
nns=(noisspec**2)[1:slen] # noise spectrum normalisation
# make multiple templates
ampindx=numpy.argsort(amplitudes) # sort amplitudes
fntemp=amplitudes.size/float(ntemplates) # number of amplitudes for one template
templates=numpy.empty((ntemplates,(slen-1)),dtype=float) # to keep templates
tempx=numpy.arange(ntemplates,dtype=float) # index coordinate for templates
norms=numpy.empty(ntemplates)
tnorms=numpy.empty(ntemplates) # templates normalization
itx=0.0
ntref=ntemplates//2
for i in numpy.arange(ntemplates):
it1=int(itx)
itx=itx+fntemp
it2=min(((int(itx)+1),indx.size))
tnorms[i]=numpy.mean(amplitudes[ampindx[it1:it2]]) # template scaling identifier
templates[i,:]=numpy.mean(fpulses[ampindx[it1:it2],:],axis=0) # template
norms[i]=numpy.sum((templates[i,:]**2)/nns)
tfits=numpy.empty(ntemplates)
tfits[0]=numpy.sum(templates[0,:]*templates[0,:]/nns)/norms[0]
for i in numpy.arange(1,ntemplates): # factor of template with respect to previous
tfits[i]=numpy.sum(templates[i,:]*templates[(i-1),:]/nns)/norms[(i-1)]
for i in numpy.arange(2,ntemplates): # normalize factors with respect to first template
tfits[i]=tfits[i]*tfits[i-1]
tfits=tfits/tfits[ntref] # normalize with respect to reference
for i in numpy.arange(ntemplates): # mutual normalisation of the templates
templates[i,:]=templates[i,:]/tfits[i]
norms[i]=numpy.sum((templates[i,:]**2)/nns)
#
if wrtfilter: # write the templates to file
write_optfilter(templates,tnorms,hdf,channel,freq)
#
if debug:
f, ax = plt.subplots(1)
ax.set_yscale('log')
ppx=numpy.arange(slen-1)
for i in numpy.arange(ntemplates):
ax.plot(ppx,templates[i,:])
ax.set_title('Templates')
plt.show()
plt.close('all')
#
ifit=numpy.zeros(indx.size,dtype=float) # initialize array to store optimal filter fit parameters
nlast=ntemplates-1
for i in numpy.arange(indx.size):
ccx=numpy.interp(amplitudes[i],tnorms,tempx) # interpolated coordinate for template
icx=int(ccx)
if ccx < 0.0: # coordinate before first template
apulse=templates[0,:] # use first template
norm=norms[0] # and first normalization
elif ccx >= float(nlast): # coordinate beyond last template
apulse=templates[nlast,:] # use last template
norm=norms[nlast] # and last normalization
else:
apulse=templates[icx,:]+(ccx-icx)*(templates[(icx+1),:]-templates[icx,:]) # interpolated template
norm=norms[icx]+(ccx-icx)*(norms[icx+1]-norms[icx]) # interpolated normalization
# if ccx < 0.0: # coordinate before first template
# icx=0 # extrapolate from first template
# elif ccx >= float(nlast): # coordinate beyond last template
# icx=nlast-2 # extrapolate from last template
# apulse=templates[icx,:]+(ccx-icx)*(templates[(icx+1),:]-templates[icx,:]) # interpolated template
# norm=norms[icx]+(ccx-icx)*(norms[icx+1]-norms[icx]) # interpolated normalization
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)