#!/usr/bin/env python3
#
'''
.. module:: HDF_Eventutils
:synopsis: Methods and utility functions for the HDF_Eventpar script
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import os
import numpy
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from tesfdmtools.methods.holzer import holzer,holzerfit
from tesfdmtools.methods.spectrumfit import spectrumfit
from tesfdmtools.methods.IQrotate import IQphase,IQphrot
from tesfdmtools.methods.xcorplot import xcorplot
from tesfdmtools.hdf.HMUX import HMUX
from tesfdmtools.utils.widgets.filechooser import getfile
from tesfdmtools.utils.cu import cu
[docs]def dellist(delist,nrecs):
'''
make logical array specifying which event records to use from HDF file
Args:
* `delist` = full file name which hold the records which decribe which events to ignore
* `nrecs` = size of array to make
Returns:
* evuse = logical array with True values for events to use
'''
if delist is not None:
if not os.path.exists(delist):
print("Warning, event delete list file '%s' does not exist" % delist)
print("Will use all events")
delist=None
evuse=numpy.ones(nrecs,dtype=bool) # all True array; use all events
if delist is not None:
ff=open(delist,'r')
for j,line in enumerate(ff):
lll=line.strip('\n')
ll="".join(lll.split()) # remove all spaces
nn=ll.split('-')
if len(nn) != 2:
print('Illegal format in event delete list line %d: "%s"' % ((j+1),lll))
print('Use: "from - to" with (from,to) integer values')
else:
try:
d1=min([(nrecs-1),max([0,int(nn[0])])])
d2=min([(nrecs-1),max([0,int(nn[1])])])
dd1=min([d1,d2])
dd2=max([d1,d2])+1
evuse[numpy.arange(dd1,dd2,dtype=int)]=False
except ValueError:
print('Illegal format in event delete list line %d: "%s"' % ((j+1),lll))
print('Use: "from - to" with (from,to) integer values')
ff.close()
return evuse
[docs]def eventpars(hdf,channel=None,freq=None,threshold=2000,bsec=0.05,nrecs=None,\
rotation=0.0,bpt=False,delist=None,flip=False):
'''
get basic event record parameters for further event selection. Only the I-signal is considered
Args:
* `hdf` = HDF5 input file object
Kwargs:
* `channel` = channel number to process. If Nont, take first channel [default: None]
* `freq` = frequency (pixel) number to process. If None, take first frequency [default: None]
* `threshold` = consider only records which have a max-min larger than threshold [default: None]
* `bsec` = part of record (start and end) from which baselines are derived. (This part of the record should not contain a pulse or partial pulse) [default: 0.05]
* `nrecs` = number of records to process. If None, process all records. [default: None]
* `rotation` = if not zero, rotate events with this amount [default: 0.0]
* `bpt` = if True, select pretrigger section of event only for baseline [default: False]
* `flip` = flip record data
Returns:
* `channel` = channel number which was processed
* `freq` = frequency (pixel number) which was processed
* `blines1` = background levels at start of record
* `blines2` = background levels at end of record
* `positions` = position of pulse minimum (pulse runs negative)
* `peaks` = value of pulse minimum (value at maximum pulse signal)
* `surf` = surface of pulse minus background. Bacground is linear fit between blines1 and blines2
* `ptp` = maximum signal difference within record (max - min)
* `evttim` = event times (time of event peak value in seconds)
* `rsize` = length of record (number of samples in record)
'''
if channel is None:
channel=hdf.channels[0]
if freq is None:
freq=hdf.channel[channel].freq_conf['pixel_index'][0]
print("read records from channel %d, frequency number %d" % (channel,freq))
nrecords=cu(hdf.channel[channel].attrs['Nevents'])
if ( type(nrecs) is str ) and ( nrecs == 'all' ):
nrecs=nrecords
else:
nrecs=min(nrecords,int(nrecs))
mrecs=nrecs-1
bl=100 # baseline length for double peak check
# initialize data arrays
blines1=numpy.zeros(nrecs,dtype=float)
blines2=numpy.zeros(nrecs,dtype=float)
positions=numpy.zeros(nrecs,dtype=int)
peaks=numpy.zeros(nrecs,dtype=int)
surf=numpy.zeros(nrecs,dtype=float)
ptp=numpy.zeros(nrecs,dtype=int)
evttim=numpy.zeros(nrecs,dtype=float)
evuse=dellist(delist,nrecs)
# find record size, sample frequency from first record and initialize x-axis
sirecord=hdf.channel[channel].freq[freq][0][:,0]
samplefreq=float(hdf.channel[channel].freq[freq].attrs['sample_rate']) # sample frequency
bl=int(sirecord.size*bsec) # part of record used for background
dc=int(sirecord.size/100.) # length for double peak check
xax=numpy.arange(sirecord.size)
if bpt:
bxax=xax[0:bl] # use only pretrigger section start
else:
bxax=numpy.concatenate((xax[0:bl],xax[-bl:])) # use both start and end of record
if rotation != 0.0:
print("Eventpars phase rotation: ",rotation)
# go through all data records
for i,record in enumerate(hdf.channel[channel].freq[freq]):
irecord=record[:,0] # use I-signal only
if rotation != 0.0:
irecord,qr=IQphrot(record[:,0],record[:,1],rotation)
if flip or ( irecord[0] < 0 ): # when wrong rotation, inverse of irecord
irecord=-irecord
ptp[i]=numpy.ptp(irecord) # retrieve max-min for record
if ptp[i] > threshold and evuse[i]: # if there is a pulse ( max-min > threshold ) and valid record
rectim=float(record.attrs['event_time']) # time of record
bb=(numpy.sum(irecord[0:bl])+numpy.sum(irecord[-bl:]))/(bl*2.0) # baseline level
pp=irecord.min()
tt=bb-0.3*(bb-pp) # threshold 1
indx1,=numpy.where(irecord < tt) # record below threshld 1
tt=bb-0.5*(bb-pp) # threshold 2
indx2,=numpy.where(irecord < tt) # record below threshold 2
dd1=indx1[1:]-indx1[0:-1] # peak seperation threshold 1
dd2=indx2[1:]-indx2[0:-1] # peak seperation threshold 2
position=numpy.argmin(irecord) # position of pulse peak (minimum)
evttim[i]=rectim+position/samplefreq # time of the event position in seconds
if ( dd1.size > 1 ) and ( dd2.size > 1 ) and \
( dd1.max() < dc ) and ( dd2.max() < dc ): # check for double peaks
positions[i]=position # position of the pulse (minimum)
peaks[i]=irecord[positions[i]] # value of the pulse minimum
pp=numpy.polyfit(bxax,irecord[bxax],1) # fit linear background
blines1[i]=pp[0]*xax[0]+pp[1] # background at start of record
blines2[i]=pp[0]*xax[-1]+pp[1] # background at end of record
surf[i]=numpy.sum((irecord-(pp[0]*xax+pp[1]))) # surface of pulse minus background
if (i % 1000) == 0:
print('read record: ',i)
if i >= mrecs:
break
return (channel,freq,blines1,blines2,positions,peaks,surf,ptp,evttim,sirecord.size)
[docs]def plotepar(hdf,channel,freq,bline,positions,peaks,surf,show=None,pdf=None,indx=None,ectpar=None):
'''
Plot basic event record parameters
Args:
* `hdf` = HDF5 input file object
* `channel` = channel number
* `freq` = frequency (pixel) number
* `bline` = baseline level
* `positions` = pulse positions
* `peaks` = pulse minima (pulse peaks)
* `surf` = pulse surface
Kwargs:
* `show` = plot mode [default: None] :
* 'X' : both show on screen and print to pdf file
* 'screen' : show on screen only
* 'pdf' : print to pdf file
* None : do not plot
* `pdf` = pdf object for plotting [default: None]
* `indx` = indexnumbers of records used
* `ectpar` = electric crosstalk parameters (tuple):
* `ectpar[0]` = actual list of selected records
* `ectpar[1]` = actual list of parameters for the indexed records
* `ectpar[2]` = actual threshold applied for the electric crosstalk
Returns:
none
'''
fnam=os.path.basename(hdf.hdffile).split('.')[0]
if show is not None:
if indx is None:
rax=numpy.arange(bline.size)
else:
rax=indx
nplots=4
if ectpar is not None:
if len(ectpar) == 3:
if ( ectpar[0] is not None ) and \
( ectpar[1] is not None ) and \
( ectpar[2] is not None ):
nplots=5
ptitle='%s ch=%d,px=%d' % (fnam,channel,freq)
f, ax = plt.subplots(nplots,sharex=True)
plt.suptitle(ptitle,size=12)
# ax[0].plot(rax,bline,'b.')
xcorplot(ax[0],rax,bline)
ax[0].set_ylabel('Baseline')
# ax[1].plot(rax,positions,'b.')
xcorplot(ax[1],rax,positions)
ax[1].set_ylabel('Event position')
# ax[2].plot(rax,peaks,'b.')
xcorplot(ax[2],rax,peaks)
ax[2].set_ylabel('Event minimum')
# ax[3].plot(rax,surf,'b.')
xcorplot(ax[3],rax,surf)
if nplots == 4:
ax[3].set_xlabel('Event number')
ax[3].set_ylabel('Event surface')
if nplots == 5:
# ax[4].plot(rax,ectpar[1],'b.')
xcorplot(ax[4],rax,ectpar[1])
axrr=ax[4].axis()
ax[4].plot([axrr[0],axrr[1]],[ectpar[2],ectpar[2]],'r-')
ax[4].set_xlabel('Event number')
ax[4].set_ylabel('El. crosst. par.')
if show != 'screen':
# psfile="Events_px%d_%s.ps" % (freq,fnam)
fpage = plt.gcf()
fpage.set_size_inches(8.27, 10.75)
plt.tight_layout()
plt.subplots_adjust(top=0.94)
plt.savefig(pdf,format='pdf',papertype='a4')
if ( show =='X' ) or ( show == 'screen' ):
plt.show()
plt.close('all')
return
[docs]def selecthist(data,tlevel=0.01,debug=False,label=''):
'''
compute data range which selects the main body of the data, leaving
out data values where data become sparse
Args:
* `data` = array of input values
Kwargs:
* `tlevel` = select only data where data density is above this level [default: 0.01]
* `debug` = if True, plot debug info [default: False]
* `label` = plot label for debug plot [default: ""]
Returns:
* `lowb` = low range value
* `higb` = high rang evalue
'''
n=data.size # number of total data
threshold=n*tlevel # data density threshold level
hist,xx=numpy.histogram(data,bins=100) # make histogram of data
indx,=numpy.where(hist > threshold) # select only data higher than given density
dx=xx[1]-xx[0]
lowb=xx[indx[0]]-dx # lower boundary of range
higb=xx[indx[-1]+1]+dx # high boundary of range
if debug:
f, ax = plt.subplots(1)
ax.plot(((xx[0:-1]+xx[1:])/2.0),hist,'b-')
rr=ax.axis()
yr=[rr[2],rr[3]]
ax.plot([lowb,lowb],yr,'r-')
ax.plot([higb,higb],yr,'r-')
ax.set_xlabel(label)
if int(matplotlib.__version__[0]) >= 3:
ax.set_yscale('symlog',linthresh=1.0)
else:
ax.set_yscale('symlog',linthreshy=1.0)
ax.set_ylabel('N')
plt.show()
plt.close('all')
return lowb,higb
[docs]def subsel(indx,par,prange=[None,None]):
'''
select a subset within the given parameter range
Args:
* `indx` = index of the parameters to be considered
* `par` = list of parameter values
* `prange` = selection range of the parameters
Returns:
* `sindx` = new index selecting only the relevant items
'''
if ( prange[0] is None) or ( prange[1] is None):
return indx
ss,=numpy.where( ( par[indx] > prange[0] ) & ( par[indx] < prange[1] ) )
print('Number of records in subselection: ',ss.size)
return indx[ss]
[docs]def selectevs(hdf,channel,freq,bline,positions,peaks,surf,rclen,show=None,debug=False,\
nopos=False,pdf=None,nopulse=False,bsec=0.05):
'''
select the proper events to process from the basic event record parameters
Args:
* `hdf` = HDF5 file which is being used
* `channel` = channel number whch was processed
* `freq` = frequency (pixel) number which was processed
* `bline` = baseline levels of the records
* `positions` = pulse positions for the records
* `peaks` = record minima (pulse peak values)
* `surf` = pulse surfaces
* `rclen` = record length (number of samples in record)
* `nopulse` = do not select on pulse intensity when set **True**
Kwargs:
* `show` = plot mode [default: None]:
* 'X' : both show on screen and print to pdf file
* 'screen' : show on screen only
* 'pdf' : print to pdf file
* None : do not plot
* `debug` = if True, print/plot some debug info [default: False]
* `nopos` = if True, do not select pulses on pulse position [default: False]
* `bsec` = section of record (baseline section) where pulse should not occur when `nopos`
* `pdf` = pdf object for plotting [default: None]
Returns:
* `sindx` = list of selected record numbers to use for further processing
'''
indx,=numpy.where(positions > 0) # select only valid records
print("Number of valid records: ",indx.size)
if ( show is not None ) and ( show != 'pdf' ) :
plotepar(hdf,channel,freq,bline[indx],positions[indx],peaks[indx],surf[indx],show='screen',pdf=pdf)
if nopos:
p1=int(bsec*rclen)
p2=int(0.8*rclen)
pindx,=numpy.where( ( positions > p1 ) & ( positions < p2 ) ) # pulse must not be on record limits
# pindx=indx # take all events, when no position selection is set
else:
# select positions based on maximum in position histogram
pmax=positions.max()
hpos,pedges=numpy.histogram(positions,range=(1,pmax),bins=pmax)
mpos=numpy.argmax(hpos)
if nopulse:
ww=150
else:
ww=max([2,rclen//3000])
# only use records with pulses which have their peak on the proper spot
pindx,=numpy.where(( positions > (pedges[mpos]-ww) ) & ( positions < (pedges[mpos]+ww+1 ) ) )
if debug:
try:
print("Average pulse position: ",pedges[mpos])
print("Selection width: ",ww)
except:
pass
print("Number of records after position selection: ",pindx.size)
plotepar(hdf,channel,freq,bline[pindx],\
positions[pindx],peaks[pindx],surf[pindx],show='screen',pdf=pdf)
blow,bhig=selecthist(bline[pindx],debug=debug,label='baseline') # discard tails of baseline values
if nopulse:
eindx,=numpy.where( ( bline[pindx] > blow ) & ( bline[pindx] < bhig ) ) # only select on baseline
else:
plow,phig=selecthist(peaks[pindx],debug=debug,label='peaks') # discard tails of peak values
# slow,shig=selecthist(surf[pindx],debug=debug,label='surface',tlevel=0.05) # discard tails of surface values
slow,shig=selecthist(surf[pindx],debug=debug,label='surface',tlevel=0.01) # discard tails of surface values
eindx,=numpy.where( ( bline[pindx] > blow ) & ( bline[pindx] < bhig ) & \
( peaks[pindx] > plow ) & ( peaks[pindx] < phig ) & \
( surf[pindx] > slow ) & ( surf[pindx] < shig ) )
if debug:
print("Number of selected pulses: ",eindx.size)
sindx=pindx[eindx] # index of selected records
# if show is not None:
# plotepar(hdf,channel,freq,bline[sindx],\
# positions[sindx],peaks[sindx],surf[sindx],show=show,pdf=pdf)
print("Number of selected records: ",sindx.size)
return sindx
[docs]def noisefile(hdf,channel,freq,nfile=None):
'''
See if there is an associated noise file, or a noise file is given.
If so , open file and retrieve ptp info
Args:
* `hdf` = HDF5 input xray file
* `channel` = channel number to use
* `freq` = frequency number to use
Kwargs:
* `nfile` = noise file to use, or True to ask for noise file [default: None]
Returns:
* `nhdf` = HDF5 file object for noise data
* `ptp` = max-min signal fluctuations for the records
'''
xnam=hdf.hdffile
if type(nfile) is bool:
if nfile == True:
ndir=os.path.dirname(xnam)
nfile=getfile(pattern='*noise*.h5',path=ndir)
else:
if nfile is None:
if os.path.basename(xnam).find('xray') > 0:
nnam=xnam.replace('xray','noise')
if os.path.exists(nnam):
nfile=nnam
if nfile is None:
print("No dedicated noise file found. Get noise from xray file")
return (None,None)
print("noise file: ",os.path.basename(nfile))
nhdf=HMUX(filename=nfile)
nrecords=cu(nhdf.channel[channel].attrs['Nevents'])
ptp=numpy.zeros(nrecords,dtype=int)
for i,record in enumerate(nhdf.channel[channel].freq[freq]):
ptp[i]=numpy.ptp(record[:,0])
return (nhdf,ptp)
[docs]def noisespec(hdf,channel,freq,ptp,fraction=0.1,debug=False,absolute=False,prlen=None,flip=False):
'''
Retrieve noise spectra. Use given fraction of I-signal records with least max-min signal fluctuations.
Args:
* `hdf` = HDF5 input file object
* `channel` = channel number to use
* `freq` = frequency (pixel) number to use
* `ptp` = max-min signal fluctuations for the records
Kwargs:
* `fraction` = fraction of least signal fluctuations to use [default: 0.1]
* `debug` = if True, print/plot debug info [default: False]
* `prlen` = length of record to take [default: None]
* `flip` = flip record data
Returns:
* `noisspec` = average noise spectrum
* `fax` = frequency axis for noise spectrum
'''
n=ptp.size
if hdf.hdffile.find('noise') > 0:
fraction=min([0.9,(4.0*fraction)]) # for real noise files, take larger fraction
hh,pp=numpy.histogram(ptp,bins=5000) # histogram (max-min) values
shh=numpy.cumsum(hh[1:]) # cumulative distribution (excuding possible 0-value first)
nn=numpy.sum(hh[1:])
mindx,=numpy.where(shh <= (int(fraction*nn+0.5)) ) # take given fraction
if len(mindx) == 0:
mindx=[0]
minptp=pp[mindx[-1]+2] # boundary value for this fraction
indx,=numpy.where( ( ptp < minptp ) & ( ptp != 0 ) ) # lowest fraction
sirecord=hdf.channel[channel].freq[freq][0][:,0] # get record size from first record
pr1=0
pr2=sirecord.size
if prlen is not None:
sirecord=sirecord[0:prlen]
pr2=prlen
samplerate=float(cu(hdf.channel[channel].freq[freq].attrs['sample_rate']))
slen=sirecord.size//2
noisspec=numpy.zeros(slen,dtype=float) # initialize noisespectrum
print("Make noise spectrum")
for i,irec in enumerate(indx): # go through all selected records
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]
irecord=irecord[pr1:pr2]
# print("Noisespec: sirecord.size,slen: ",sirecord.size,slen)
noisspec=noisspec+numpy.absolute(numpy.fft.fft(irecord)[0:slen]) # accumulate signal data to noisespectrum
if ( i % 1000 ) == 0:
print("process record ",irec)
print("number of noise records processed: ",i)
noisspec=noisspec/float(indx.size) # take average
noisspec[0]=0.0 # discard offsets
fax=numpy.arange(slen)/float(slen*2)*samplerate/1000.0 # compute frequency axis using sample rate
if debug:
f, ax = plt.subplots(1)
ax.plot(fax[1:],noisspec[1:],'b-')
ax.set_xlabel('frequency (kHz)')
if int(matplotlib.__version__[0]) >= 3:
ax.set_yscale('symlog',linthresh=1.0)
else:
ax.set_yscale('symlog',linthreshy=1.0)
ax.set_ylabel('Power')
ax.set_title('Noise spectrum')
plt.show()
plt.close('all')
return (noisspec,fax)
[docs]def noisplot(pltax,spectrum,nspectrum,fax,title):
'''
plot noise spectrum
Args:
* `pltax` = plot instance for plotting
* `spectrum` = spectrum data
* `fax` = frequency axis
* `title` = title of plot
'''
# bbs=numpy.mean(spectrum[-100:])
# bbn=numpy.mean(nspectrum[-100:])
# scale=bbs/bbn
pltax.plot(fax[1:]*1000.0,nspectrum[1:],'g-')
pltax.plot(fax[1:]*1000.0,spectrum[1:],'b-')
pltax.set_xlabel('frequency (Hz)')
if int(matplotlib.__version__[0]) >= 3:
pltax.set_yscale('symlog',linthresh=1.00)
pltax.set_xscale('symlog',linthresh=0.01)
else:
pltax.set_yscale('symlog',linthreshy=1.00)
pltax.set_xscale('symlog',linthreshx=0.01)
pltax.set_ylabel('Power')
pltax.set_title(title)
[docs]def getphase(hdf,channel,freq,indx,nevents=100,debug=False):
'''
Get average phase of first 'nevents' events
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
Kwargs:
* `nevents` = number of events to process [default: 100]
* `debug` = if True, show some events [default: False]
Returns:
* `phase` = average phase correction for events
'''
nmax=min([indx.size,nevents])
phase=0.0
for i,irec in enumerate(indx):
if i >= nmax:
break
record=hdf.channel[channel].freq[freq][irec]
phase=phase+IQphase(record[:,0],record[:,1])
phase=phase/float(nmax)
print("Average phase: %7.2f (deg)" % (phase/numpy.pi*180.0))
if debug:
nax=numpy.arange(record[:,0].size)
ioplt={}
qoplt={}
inplt={}
qnplt={}
for i in numpy.arange(4):
record=hdf.channel[channel].freq[freq][indx[i]]
ioplt[i] = plt.subplot2grid((4,2),(i,0))
qoplt[i]= ioplt[i].twinx()
ioplt[i].set_ylabel('ADC')
ioplt[i].plot(nax,record[:,0],'b-')
qoplt[i].plot(nax,record[:,1],'r-')
for tl in qoplt[i].get_yticklabels():
tl.set_color('r')
inplt[i] = plt.subplot2grid((4,2),(i,1))
qnplt[i]= inplt[i].twinx()
inew,qnew=IQphrot(record[:,0],record[:,1],phase)
inplt[i].plot(nax,inew,'b-')
qnplt[i].plot(nax,qnew,'r-')
for tl in qnplt[i].get_yticklabels():
tl.set_color('r')
ioplt[3].set_xlabel('sample #')
inplt[3].set_xlabel('sample #')
plt.show()
plt.close('all')
return -phase # return phase correction
from scipy.optimize import curve_fit
def pulseshape(p,x):
tt=x-p[1]
indx,=numpy.where(tt < 0)
shape=p[0]*numpy.exp(-tt/p[2])*(1.0-numpy.exp(-tt/p[3]))
shape[indx]=0.0
return shape
[docs]def getrisetime(xax,pulse,debug=False,plotfit=None):
'''
compute rise time for pulse, based on linear fit of log values
Args:
* `xax` = X-axis of pulse expressed in sample #
* `pulse` = pulse shape, with background subtracted
Kwargs:
* `debug` = if True, plot each pulse and fit [default: False]
* `plotfit` = plotinstance for log plot of fitted pulse [default: None]
Returns:
* `rtime` = rise time of pulse in sample units
* `ftime` = fall time of pulse in sample units
'''
hfunc=lambda x, p1, p2, p3, p4: pulseshape((p1,p2,p3,p4),x)
pp=numpy.argmin(pulse) # position of pulse top (minimum)
top=pulse[pp] # pulse minimum
indxr,=numpy.where( ( pulse[0:pp] < 0.15*top ) & ( pulse[0:pp] > 0.8*top ) ) # rising flank
if len(indxr) > 1:
ll=numpy.log(-pulse[indxr]) # fit on log pulse
lfit=numpy.polyfit(xax[indxr],ll,1)
rtime=1.0/lfit[0]
else:
return (0.0,0.0)
position=xax[indxr[0]]
indxf,=numpy.where( ( pulse[pp:] < 0.15*top ) & ( pulse[pp:] > 0.8 * top ) ) # pulse decay
if len(indxf) > 1:
indxf=indxf+pp
lld=numpy.log(-pulse[indxf]) # fit on log pulse
lfit=numpy.polyfit(xax[indxf],lld,1)
ftime=-1.0/lfit[0]
else:
return (0.0,0.0)
p=numpy.array([top,position,ftime,rtime],dtype=float)
if debug:
print("initial fitpars: ",p)
# sigma=numpy.ones(pulse.size,dtype=float)
findx,=numpy.where( ( pulse < 0.1*top ) & ( pulse > 0.8 * top ) )
sigma=numpy.ones(findx.size,dtype=float)
# indx,=numpy.where( ( pulse[findx] > 0.1*top ) | ( pulse[findx] < 0.8 * top ) )
# sigma[indx]=sigma[indx]*100.0
try:
p1, covar = curve_fit(hfunc, xax[findx], pulse[findx], p.copy(), sigma)
except:
if debug:
print("fit failed")
return (0.0,0.0)
if debug:
print("fitted pars: ",p1)
rtime=p1[3] # rise time
ftime=p1[2] # fall time
if debug:
f, ax = plt.subplots(4,sharex=False)
ax[0].set_xlabel("Sample #")
ax[0].set_ylabel("pulse ADC")
ax[0].set_title("Rise time: %6.3f" % rtime)
ax[1].set_xlabel("Sample #")
ax[1].set_ylabel("Log pulsefit")
ax[1].set_title("")
dp=indxr[-1]-indxr[0]
pp1=indxr[0]-2*dp
pp2=indxr[-1]+3*dp
ax[0].plot(xax[pp1:pp2],pulse[pp1:pp2])
ax[0].plot(xax[indxr],pulseshape(p1,xax[indxr]))
ax[1].plot(xax[indxr],ll)
ax[1].plot(xax[indxr],numpy.log(-pulseshape(p1,xax[indxr])))
ax[2].set_xlabel("Sample #")
ax[2].set_ylabel("pulse ADC")
ax[2].set_title("Fall time: %6.3f" % ftime)
ax[3].set_xlabel("Sample #")
ax[3].set_ylabel("Log pulsefit")
ax[3].set_title("")
pp1=indxr[-1]
pp2=indxf[-1]+3*dp
ax[2].plot(xax[pp1:pp2],pulse[pp1:pp2])
ax[2].plot(xax[indxf],pulseshape(p1,xax[indxf]))
ax[3].plot(xax[indxf],lld)
ax[3].plot(xax[indxf],numpy.log(-pulseshape(p1,xax[indxf])))
plt.show()
plt.close("all")
if plotfit is not None:
plotfit.plot(xax[indxr],-numpy.log10(-pulseshape(p1,xax[indxr])))
plotfit.plot(xax[indxf],-numpy.log10(-pulseshape(p1,xax[indxf])))
return (rtime,ftime)
[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,flip=False,showpspec=None,\
**kwargs):
'''
perform optimal filtering fit of individual pulses
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]
* `flip` = flip record data
* `showpspec` = interactive plot of a given set of pulse spectra
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:
if freqcutoff > samplerate:
print("Warning: given frequency cutoff is larger than samplerate ",samplerate)
print("Frequency cutoff set to samplerate")
slen=sirecord.size//2
else:
slen=int(float(freqcutoff)/float(samplerate)*sirecord.size)//2 # 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
# showpspec=1000
if showpspec is not None: # show set of pulse spectra used for opt. filter
f,ax=plt.subplots(1)
sint=showpspec//10
spax=numpy.arange((slen-1))/(sirecord.size-1)*samplerate
for k in numpy.arange(0,showpspec,sint):
ax.plot(spax,fpulses[k,:])
ax.plot(spax,apulse,'k-')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Intensity')
ax.set_title('Pulses spectra')
ax.set_yscale('symlog',linthresh=1.0)
ax.set_xscale('symlog',linthresh=1.0)
plt.show()
plt.close('all')
nns=(noisspec**2)[1:slen]
# print("size: noisspec: ",noisspec.size)
# print("sizes: apulse,nns,slen:",apulse.size,nns.size,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',color='blue')
ax[0].set_xlabel('Energy (arb. units)')
ax[0].set_ylabel('Counts')
if risetime:
ax[2].hist(rtimes*1000,bins=512,align='mid',color='blue')
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)
[docs]def limdist(cc,ext=0.5,bins=2500):
'''
compute range of main data distribution in array of values
Args:
* `cc` = array of values
Kwargs:
* `bins` = number of bins to use for histogram
* `ext` = range extension [default: 0.5]
Returns:
* `xmin` = minimum of range
* `xmax` = maximum of range
'''
hh,xx=numpy.histogram(cc,bins=bins) # make histogram
imax=numpy.argmax(hh)
i=hh.size-1
for i in numpy.arange((imax+1),hh.size): # find position of first zero bin below histogram maximum
if hh[i] == 0:
break
xmax=xx[i]
for i in ((imax-1)-numpy.arange(imax)): # find position of first zero bin above histogram maximum
if hh[i] == 0:
break
xmin=xx[i+1]
ww=xmax-xmin
xmin=xmin-ext*ww # extend range; add half the distribution width
xmax=xmax+ext*ww
return (xmin,xmax)
[docs]def showcorr(c1,c2,xlabel='X',ylabel='Y',title='',center=None,ytilt=None):
'''
plot correlation between two parameters and the fitted linear correlation relation
Args:
* `c1` = array of first parameter
* `c2` = array of second parameter
Kwargs:
* `xlabel` = label for x-axis [default: 'X']
* `ylabel` = label for y-axis [default: 'Y']
* `title` = title for plot [default: '']
* `center` = fitted center of distribution [default: None]
* `ytilt` = fitted linear slope between parameters [default: None]
Returns:
none
'''
f, ax = plt.subplots(1)
# ax.plot(c1,c2,'b.') # plot data points
xcorplot(ax,c1,c2) # plot data points
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
if ( center is not None ) and ( ytilt is not None ):
rr=ax.axis()
y1=center[1]+(rr[0]-center[0])*ytilt # plot fitted linear correlation
y2=center[1]+(rr[1]-center[0])*ytilt
ax.plot([rr[0],rr[1]],[y1,y2],'r-',label=('%5.3f' % (ytilt*1000.0)))
ax.legend()
plt.show()
plt.close('all')
# sf=open('crossdat.txt','w')
# sf.write('%d\n' % c1.size )
# for i,d in enumerate(c1):
# sf.write('%15.8e %15.8e\n' % (d,c2[i]) )
# sf.close()
return
[docs]def timplt(ax,hdf,channel,freq,c1,c2,frange=None,xlabel='c1',ylabel='c2'):
'''
time evolution plot of pulse fit parameter
Args:
* `ax` = a plotting instance for the plot output
* `hdf` = HDF5 input file object
* `channel` = channel number being processed
* `freq` = frequency number (pixel) being processed
* `c1` = array of time parameter (event record number)
* `c2` = array of pulse fits
Kwargs:
* `frange` = range of fits parameter values to plot [default: None]
* `xlabel` = label for x-axis [default: 'c1']
* `ylabel` = label for y-axis [default: 'c2']
Returns:
none
'''
# freqindx=hdf.channel[channel].freqs.index(freq)
fffx,=numpy.where( hdf.channel[channel].freq_conf['pixel_index'] == freq )
freqindx=fffx[0]
khz=hdf.channel[channel].freq_conf['freq'][freqindx]
title='ch=%d, f=%d: %7.2f (kHz)' % (channel,freq,khz) # make plot title
emin,emax=limdist(c2,ext=0.3)
indx,=numpy.where( ( c2 > emin ) & ( c2 < emax ) )
# ax.plot(c1[indx],c2[indx],'b.') # plot fitted linear correlation
xcorplot(ax,c1[indx],c2[indx]) # plot fitted linear correlation
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
return
[docs]def crosscplt(ax,hdf,channel,freq,c1,c2,xlabel='c1',ylabel='c2',center=None,ytilt=None):
'''
plot correlation between two parameters and the fitted linear correlation relation
Args:
* `ax` = a plotting instance for the plot output
* `hdf` = HDF5 input file object
* `channel` = channel number being processed
* `freq` = frequency number (pixel) being processed
* `c1` = array of first parameter
* `c2` = array of second parameter
Kwargs:
* `xlabel` = label for x-axis [default: 'c1']
* `ylabel` = label for y-axis [default: 'c12']
* `center` = fitted center of distribution [default: None]
* `ytilt` = fitted linear slope between parameters [default: 'None']
Returns:
none
'''
# freqindx=hdf.channel[channel].freqs.index(freq)
fffx,=numpy.where( hdf.channel[channel].freq_conf['pixel_index'] == freq )
freqindx=fffx[0]
khz=hdf.channel[channel].freq_conf['freq'][freqindx]
title='ch=%d, f=%d: %7.2f (kHz)' % (channel,freq,khz) # make plot title
xmin,xmax=limdist(c2) # compute proper x-axis range
hh,yy=numpy.histogram(c1,bins=10000) # histogram for y values, to compute proper range
cc=numpy.cumsum(hh)
efrac=0.05 # fraction of y values to leave out
cndx,=numpy.where( (cc > int(efrac*c1.size)) & (cc < int((1.0-efrac)*c1.size)) )
if cndx.size < 10:
cndx=numpy.arange(cc.size)
xmin=numpy.amin(c2)
xmax=numpy.amax(c2)
ymin=yy[cndx[0]] # proper y-axis range
ymax=yy[cndx[-1]+1]
indx,=numpy.where( ( c2 > xmin ) & ( c2 < xmax ) &
( c1 > ymin ) & ( c1 < ymax ) )
# ax.plot(c1[indx],c2[indx],'b.') # plot fitted linear correlation
xcorplot(ax,c1[indx],c2[indx]) # plot fitted linear correlation
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
if ( center is not None ) and ( ytilt is not None ):
rr=ax.axis()
y1=center[1]+(rr[0]-center[0])*ytilt # plot fitted linear correlation
y2=center[1]+(rr[1]-center[0])*ytilt
pt=ax.plot([rr[0],rr[1]],[y1,y2],'r-',label=('%f5.3' % (ytilt*1000.0)))
ax.legend()
return
[docs]def makespectrum(pulsefits,sbins=5000,debug=False,efrac=0.05):
'''
make spectrum of fitted optimal filtering scale factors, and determine proper
range (leave out the outliers)
Args:
* `pulsefits` = array of pulse fit paramters (scale factors)
Kwargs:
* `sbins` = number of bins for spectrum [default: 5000]
* `debug` = if True, plot/print debug output [default: False]
* `efrac` = fraction of the events at low and hig end defined as outliers [default: 0.05]
Returns:
* `spectrum` = spectrum of the pulsefits
* `bins` = centers of the spectral bins
'''
print("makespectrum")
spectrum,bins=numpy.histogram(pulsefits,bins=10000) # make initial spectrum
bins=(bins[0:-1]+bins[1:])/2.0 # compute bin centers
cc=numpy.cumsum(spectrum) # cumulative sum
indx,=numpy.where((cc > int(efrac*pulsefits.size)) & (cc < int((1.0-efrac)*pulsefits.size)) ) # discard the outliers
if len(indx) == 0:
b1=bins[0]
b2=bins[-1]
else:
b1=bins[indx[0]]
b2=bins[indx[-1]]
bb=b2-b1
b1=b1-0.1*bb # extend range a bit
b2=b2+0.1*bb
spectrum,bins=numpy.histogram(pulsefits,bins=sbins,range=[b1,b2]) # final spectrum
bins=(bins[0:-1]+bins[1:])/2.0 # bin centers
if debug:
f, axarr = plt.subplots(1)
axarr.set_xlabel('Spectral bin')
axarr.set_ylabel('Counts')
axarr.set_title('Binned spectrum')
axarr.plot(bins,spectrum)
plt.show()
plt.close('all')
return (spectrum,bins)
[docs]def fitspectrum(pltdata,iplt,hdf,channel,freq,spectrum,bins,debug=False,afitrange=[5.860,5.920]):
'''
fit the Holzer alpha and beta Mn54 lines to the spectrum
Args:
* `pltdata` = dictionary for the fit results
* `iplt` = number of the plot to process
* `hdf` = HDF5 input file object
* `channel` = channel number being processed
* `freq` = frequency number (pixel) being processed
* `spectrum` = spectrum of the pulse fit parameters
* `bins` = centers of the spectrum bins
Kwargs:
* `debug` = if True, plot debug information [default: False]
* `afitrange` = range (keV) for spectrum fit [default: 5.860,5.920]
Results:
* `pltdata[iplot]` = spectrum fit results
Returns:
* `fwhm` = fitted resolution inFWHM
* `fitpars` = additional fitted parameters
'''
fffx,=numpy.where( hdf.channel[channel].freq_conf['pixel_index'] == freq )
freqindx=fffx[0]
khz=hdf.channel[channel].freq_conf['freq'][freqindx]
fid="ch=%d, px=%d %7.2f (kHz)" % (channel,freq,khz) # plot title
nppp,spectrum,xxx,holy,lratio=spectrumfit(spectrum,bins,debug=debug,afitrange=afitrange)
if ( nppp is None ) or ( nppp[0] == 0.0 and nppp[1] == 0.0):
pltdata[iplt]=None
fitpars={'khz':khz,'gain':0.0,'counts':0.0,'lratio':0.0,'fwhmerr':0.0}
return (0.0,fitpars)
fitpars={'khz':khz,'gain':nppp[1],'counts':numpy.sum(spectrum),'lratio':lratio,\
'fwhmerr':(nppp[9]*1000.0)}
fwhmtxt="%5.5s %12.4g +/- %12.4g (eV)" % ("FWHM: ",(nppp[3]*1000.0),(nppp[9]*1000.0)) # FWHM + 1-sigma error
# txt="%5.5s %12.4g\n" % ("Shift: ",nppp[0]) # text for plot
# txt=txt+( "%5.5s %12.4g\n" % ("Gain: ",nppp[1]) )
txt=' '.join(fwhmtxt.split())+'\n'
# txt=txt+( "%5.5s %12.4g\n" % ("Corr: ",((bpos-avapos)/(bpos-(avapos-avashift/1000.0)))) )
txt=txt+( "%d counts\n" % (fitpars['counts']) )
txt=txt+( "lin.ratio:%7.4f\n" % lratio )
txt=txt+( "C-stat: %5.1f\n" % nppp[4] )
txt=txt+( "degs of frdm: %d\n" % nppp[5] )
print("FWHM: %12.4g +/- %10.4g (eV)" % (nppp[3]*1000.0,nppp[9]*1000.0)) # FWHM with 1-sigma error
yrr=numpy.sqrt(spectrum) # 1-sigma errors
yerr=numpy.where(yrr <= 0.0,1.0,yrr) # for zero counts insert one-count error
pltdata[iplt]={'xax':xxx,'spectrum':spectrum,'yerr':yerr,'holy':holy,\
'ptitle':fid,'ptxt':txt} # spectral fit data
return ((nppp[3]*1000.0),fitpars)
[docs]def makeplot(pdf,ps=False):
'''
finish plot
Args:
* `pdf` = plot output file
Kwargs:
* `ps` = if False, also put output on screen [default: False]
'''
fpage = plt.gcf()
fpage.set_size_inches(8.27, 10.75)
plt.tight_layout()
plt.subplots_adjust(top=0.94)
plt.savefig(pdf,format='pdf',papertype='a4')
if not ps:
plt.show()
plt.close('all')
[docs]def xoutl(values):
'''
find reasonable index which excludes outliers
Args:
* `values` = array of values
Returns:
* `rindx` = index which excludes outliers
'''
rindx,=numpy.where( (values > 0.0) & ( values < values.max() ) )
if rindx.size == 0:
return numpy.arange(values.size)
hh,bb=numpy.histogram(values[rindx],bins=500)
bb=(bb[0:-1]+bb[1:])/2.0
hindx,=numpy.where(hh > 1)
rindx,=numpy.where( (values > 0.0) & ( values <= bb[hindx[-1]] ) )
if rindx.size == 0:
return numpy.arange(values.size)
hh,bb=numpy.histogram(values[rindx],bins=500)
bb=(bb[0:-1]+bb[1:])/2.0
hindx,=numpy.where(hh > 2)
if hindx.size < 1:
return numpy.arange(values.size)
h1=max([hindx[0]-10,0])
h2=min([hindx[-1]+10,(bb.size-1)])
rindx,=numpy.where( ( values > bb[h1] ) & ( values < bb[h2] ) )
if rindx.size < 2:
rindx=numpy.arange(values.size)
return rindx
[docs]def plotptimes(rtimes,ftimes,nfit,avpulse,ptitle,pdf,ps=False,outtxt=False,erange=[None,None],\
avbline=None):
'''
Make plot of pulse times (rise and fall) and correlations with energy + average pulse profile
Also write ASCII data of average pulse profile when outtxt=True
Args:
* `rtimes` = pulse rise times for the events
* `ftimes` = pulse fall times for the events
* `nfit` = optimum filter fits for the events (scales with energy)
* `avpulse` = the average pulse profile
* `ptitle` = title of the plot page
* `pdf` = pdf plot instance
Kwargs:
* `ps` = if False, also show plot interactively [default: False]
* `outtxt` = if True, write ASCII data of average pulse profile [default: False]
* `erange` = range of the nfit parameter to plot in the energy/(rise/fall)time correlation [default: None,None]
* `avbline` = average baseline level for the average pulse. To be included in the ASCII pulse profile [default: None]
'''
plt.suptitle(ptitle,size=12)
rplot={}
for i in [0,1]:
for j in [0,1]:
rplot[i*2+j]=plt.subplot2grid((3,2),(j,i))
rplot[4]=plt.subplot2grid((3,1),(2,0))
if ( erange[0] is not None ) and ( erange[1] is not None ):
eindx,=numpy.where( ( nfit > erange[0] ) & ( nfit < erange[1] ) )
tindx=xoutl(rtimes[eindx])
rindx=eindx[tindx]
tindx=xoutl(ftimes[eindx])
findx=eindx[tindx]
else:
rindx=xoutl(rtimes)
findx=xoutl(ftimes)
rplot[0].hist(rtimes[rindx]*1E6,bins=500,color='blue')
rplot[0].set_xlabel(r'Rise time ($\mu$s)')
rplot[0].set_ylabel("N")
rplot[0].set_title("Histogram of rise times")
# rplot[1].plot(rtimes[rindx]*1E6,nfit[rindx],'b.')
xcorplot(rplot[1],rtimes[rindx]*1E6,nfit[rindx])
rplot[1].set_xlabel(r'Rise time ($\mu$s)')
rplot[1].set_ylabel("Energy (arb. units)")
rplot[1].set_title("Correlation")
rplot[2].set_xlabel(r'Fall time ($\mu$s)')
rplot[2].set_ylabel("N")
rplot[2].set_title("Histogram of fall times")
rplot[2].hist(ftimes[findx]*1E6,bins=500,color='blue')
# rplot[3].plot(ftimes[findx]*1E6,nfit[findx],'b.')
xcorplot(rplot[3],ftimes[findx]*1E6,nfit[findx])
rplot[3].set_xlabel(r'Fall time ($\mu$s)')
rplot[3].set_ylabel("Energy (arb. units)")
rplot[3].set_title("Correlation")
pax=numpy.arange(avpulse.size)
rplot[4].set_title("Average pulse profile")
rplot[4].set_xlabel("Record sample #")
rplot[4].set_ylabel("-Log(-ADC) units")
indx,=numpy.where(avpulse < -0.01)
pmin=numpy.argmin(avpulse)
# p1=max([0,(indx[0]-20)])
# p2=min([(avpulse.size-1),(indx[-1]+20)])
dw=int(500./16384.*avpulse.size)
p1=max([0,pmin-dw])
p2=min([pmin+10*dw,avpulse.size])
logpulse=numpy.empty(avpulse.size)
logpulse[:]=numpy.nan
logpulse[indx]=-numpy.log10(-avpulse[indx])
rplot[4].plot(pax[p1:p2],logpulse[p1:p2],'b-')
prt,pft=getrisetime(pax,avpulse,plotfit=rplot[4])
makeplot(pdf,ps=ps)
if not outtxt:
return
pp=ptitle.split(',')
pix=pp[-1].split('=')[-1].strip(' ')
if pix.isdigit():
pixel=int(pix)
px='_px%3.3d' % pixel
else:
px='_px___'
pixel=-1
chl=pp[0].split()[-1].split('=')[-1].strip(' ')
if chl.isdigit():
channel=int(chl)
else:
channel=-1
txtnam='Avpulse_'+ptitle.strip().split()[0]+px+'.txt'
ofile=open(txtnam,'w')
ofile.write('# File_id: %s\n' % ptitle.strip().split()[0])
ofile.write('# channel: %d\n' % channel)
ofile.write('# pixel: %d\n' % pixel)
if avbline is not None:
ofile.write('# average baseline level: %10.2f\n' % avbline)
ofile.write('# x ADC\n')
# for i,x in enumerate(pax[p1:p2]):
# ofile.write('%7d %14.5f\n' % (x,avpulse[p1+i]) )
for x in pax:
ofile.write('%7d %14.5f\n' % (x,avpulse[x]) )
ofile.close()
return
[docs]def eventlist_out(wfile,freq,recordno,optfit,baseline,peaks,hdffilename,channel=None,
rtimes=None,ftimes=None,evttim=None,eefit=None,surf=None):
'''
write eventlist to file
Args:
* `wfile` = write eventlist or not
* `freq` = frequency number in file
* `recordno` = list of record numbers
* `optfit` = list of optimal filtering results for each event
* `baseline` = list of baseline levels for each event
* `peaks` = list of pulse peak values for each event
* `hdffilename` = full filename (including path) of processed HDF file
Kwargs:
* `channel` = channel number used
* `rtimes` = pulses rise times
* `ftimes` = pulses fall times
* `evttim` = time of the events
* `eefit` = list of (calibrated) energies for the events
* `surf` = list of pulse surfaces for each event
'''
if not wfile:
return
ws=False
if surf is not None:
ws=True
wt=False
if ( rtimes is not None ) and ( ftimes is not None ):
wt=True
we=False
if evttim is not None:
we=True
ifnam=os.path.basename(hdffilename)
ifbas=ifnam.rsplit('.',1)[0]
fnam='eventlist_'+ifbas+('_px%3.3d' % freq)+'.lst'
ofile=open(fnam,'w')
hdfaname=os.path.abspath(hdffilename)
if hdfaname.find('/stage/') >= 0:
hdfnam=hdfaname
else:
hdfnam=hdfaname.replace('/xmmdat10/','/stage/xmmdat10/')
ofile.write("#inputfile: %s\n" % hdfnam)
if channel is not None:
ofile.write("# channel: %d\n" % channel)
ofile.write("# pixel: %d\n" % freq)
ofile.write("#record-no opt.fit baseline peak-value")
if eefit is not None:
ofile.write(" energy(eV)")
if we:
ofile.write(" event-time")
if wt:
ofile.write(' risetime falltime')
if ws:
ofile.write(' event-surface')
ofile.write('\n')
if eefit is not None:
if wt:
if we:
if ws:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %13.2f %14.4f %14.1f %14.1f %14.4e\n" \
% (recno,optfit[i],baseline[i],peaks[i],eefit[i],\
evttim[i],rtimes[i]*1E6,ftimes[i]*1E6,surf[i]) )
else:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %13.2f %14.4f %14.1f %14.1f\n" \
% (recno,optfit[i],baseline[i],peaks[i],eefit[i],\
evttim[i],rtimes[i]*1E6,ftimes[i]*1E6) )
else:
if ws:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %13.2f %14.1f %14.1f %14.4e\n" \
% (recno,optfit[i],baseline[i],peaks[i],eefit[i],\
rtimes[i]*1E6,ftimes[i]*1E6,surf[i]) )
else:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %13.2f %14.1f %14.1f\n" \
% (recno,optfit[i],baseline[i],peaks[i],eefit[i],\
rtimes[i]*1E6,ftimes[i]*1E6) )
else:
if we:
if ws:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %13.2f %14.4f %14.4e\n" \
% (recno,optfit[i],baseline[i],peaks[i],eefit[i],\
evttim[i],surf[i]) )
else:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %13.2f %14.4f\n" \
% (recno,optfit[i],baseline[i],peaks[i],eefit[i],\
evttim[i]) )
else:
if ws:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %13.2f %14.4e\n" \
% (recno,optfit[i],baseline[i],peaks[i],eefit[i],surf[i]) )
else:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %13.2f\n" \
% (recno,optfit[i],baseline[i],peaks[i],eefit[i]) )
else:
if wt:
if we:
if ws:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %14.4f %14.1f %14.1f %14.4e\n" \
% (recno,optfit[i],baseline[i],peaks[i],evttim[i],rtimes[i]*1E6,\
ftimes[i]*1E6,surf[i]) )
else:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %14.4f %14.1f %14.1f\n" \
% (recno,optfit[i],baseline[i],peaks[i],evttim[i],rtimes[i]*1E6,\
ftimes[i]*1E6) )
else:
if ws:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %14.1f %14.1f %14.4e\n" \
% (recno,optfit[i],baseline[i],peaks[i],rtimes[i]*1E6,\
ftimes[i]*1E6,surf[i]) )
else:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %14.1f %14.1f\n" \
% (recno,optfit[i],baseline[i],peaks[i],rtimes[i]*1E6,ftimes[i]*1E6) )
else:
if we:
if ws:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %14.4f %14.4e\n" \
% (recno,optfit[i],baseline[i],peaks[i],evttim[i],surf[i]) )
else:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %14.4f\n" \
% (recno,optfit[i],baseline[i],peaks[i],evttim[i]) )
else:
if ws:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f %14.4e\n" \
% (recno,optfit[i],baseline[i],peaks[i],surf[i]) )
else:
for i,recno in enumerate(recordno):
ofile.write("%10d %14.8f %14.4f %14.f\n" \
% (recno,optfit[i],baseline[i],peaks[i]) )
ofile.close()
return
[docs]def selectdtim(indx,evttim,etprev=[None,None],etnext=[None,None]):
'''
Select events based on the delta times to the previous and next events
Args:
* `indx` = index of valid events
* `evttim` = time tages of the events
Kwargs:
* `etprev` = selection window of delta time to the previous event [min,max]
* `etnext` = selection window of delta time to the next event [min,max]
'''
if ( etprev[0] is None ) and ( etnext[0] is None ):
return indx
vindx,=numpy.where( evttim != 0 )
vvttim=evttim[vindx]
xevttim=numpy.concatenate((numpy.array([vvttim[0]]),vvttim,numpy.array([vvttim[-1]])))
pprevtim=vvttim-xevttim[:vvttim.size] # delta time to previous event
nnexttim=xevttim[2:]-vvttim # delta time to next event
prevtim=numpy.zeros(evttim.size) # expand to comply with evttim size
nexttim=numpy.zeros(evttim.size)
prevtim[vindx]=pprevtim
nexttim[vindx]=nnexttim
if etprev[0] is not None:
pindx,=numpy.where( ( prevtim[indx] > etprev[0] ) & ( prevtim[indx] < etprev[1] ) )
indx=indx[pindx]
if etnext[0] is not None:
nindx,=numpy.where( ( nexttim[indx] > etnext[0] ) & ( nexttim[indx] < etnext[1] ) )
indx=indx[nindx]
return indx