#!/usr/bin/env python3
#
'''
.. module:: HDF_NEPutils
:synopsis: Methods and utility functions for the HDF_NEP script
.. moduleauthor:: Cor de Vries <c.p.de.vries@sron.nl>
'''
import sys,os
import numpy
from tesfdmtools.hdf.HMUX import HMUX
from tesfdmtools.utils.widgets.filechooser import getfile
from tesfdmtools.methods.IQrotate import IQphase,IQphrot
from tesfdmtools.methods.HDF_Eventutils import getrisetime
from tesfdmtools.methods.positivepulse import positivepulse
from tesfdmtools.utils.cu import cu
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
[docs]class ShowPlots:
'''
Make combined plot of the NEP data
Kwargs:
* `logfreq` = not used
* `title` = subtitle for plot
'''
def __init__(self,logfreq=False,title=None):
self.logfreq=logfreq
self.fig = plt.figure()
self.plots={}
if title is not None:
plt.suptitle(title,size=12)
self.plots[1] = plt.subplot2grid((3, 2), (0, 0))
self.plots[2] = plt.subplot2grid((3, 2), (0, 1))
self.plots[3] = plt.subplot2grid((3, 2), (1, 0))
self.plots[4] = plt.subplot2grid((3, 2), (1, 1))
self.plots[5] = plt.subplot2grid((3, 2), (2, 0), colspan=2)
[docs] def addplot(self,num,y,x=None,title="",xlabel="Sample # (/1000)",ylabel="Y"):
'''
Add a plot to the combined plot
Args:
* `num` = index number within the plot (1-5)
* `y` = array of Y-values
Kwargs:
* `x` = array of X-values. If not specified 'x' will be the index number
* `title` = plot title
* `xlabel` = label for the X-axis
* `ylabel` = label for the Y-axis
'''
if x is None:
x=numpy.arange(y.size)/1000.0
self.plots[num].plot(x,y)
self.plots[num].set_title(title)
self.plots[num].set_xlabel(xlabel)
self.plots[num].set_ylabel(ylabel)
if xlabel.find('Frequency') >= 0:
if int(matplotlib.__version__[0]) >= 3:
self.plots[num].set_xscale('symlog',linthresh=0.01)
else:
self.plots[num].set_xscale('symlog',linthreshx=0.01)
if num >= 3:
self.plots[num].set_yscale('log')
[docs] def annotate(self,num,text):
'''
Annotate plot with text
Args:
* `num` = index number within the plot (1-5)
* `text` = text string (can be multiline)
'''
rr=self.plots[num].axis()
xx=rr[0]
yy=rr[3]
ttt='\n '+text
self.plots[num].text(xx,yy,ttt,verticalalignment='top')
[docs] def show(self,psfile=None):
'''
Show the plot, or make pdf file
Kwargs:
* `psfile` = Name of pdf output file. If 'None', output to screen
'''
if psfile is not None:
fpage = plt.gcf()
fpage.set_size_inches(8.27, 10.75)
plt.tight_layout()
plt.subplots_adjust(top=0.94)
plt.savefig(psfile,format='pdf',papertype='a4')
else:
plt.tight_layout()
plt.show()
plt.close('all')
[docs]def ap2iq(record,scale):
'''
Convert amplitude-phase data record to I-Q record
Args:
* `record` = FDM HDF data record
* `scale` = scale factor for phase
'''
if scale is None:
return record
cc=record[:,0]*numpy.exp(1j*record[:,1]*scale)
crecord=numpy.empty_like(record)
crecord[:,0]=cc.real
crecord[:,1]=cc.imag
return crecord
[docs]def NEPfiles(xfile=None,nfile=None,fpath='/stage/ltsboldat5/LC_XFDM',diag=False,names=False):
'''
Get xray and noise files
Kwargs:
* `xfile` = filename of xrays. If not provided will be asked for.
* `nfile` = filename of noise data. If not provided, accompanying noise file with xray file will be used.
* `fpath` = path for files
* `diag` = if True, print some diagnostics info
* `names` = if True, return filenames instead of file objects
Returns:
* `xrays` = HDF file object with x-ray data
* `noise` = HDF file object with noise data
'''
xx=False
if xfile is None:
xfile=getfile(pattern='*.h5',path=fpath)
if xfile is None:
sys.exit("No hdf5 file selected")
else:
xx=True
if not os.path.exists(xfile):
sys.exit("File does not exist: %s" % xfile)
if type(nfile) is bool:
if nfile == True:
nfile=getfile(pattern='*noise*.h5',path=fpath)
if nfile is None:
sys.exit("No hdf5 noise file selected")
if nfile is None:
nfile=xfile.replace('xray','noise')
if not os.path.exists(nfile):
nfile=xfile
if ( not names ) or diag:
xrays=HMUX(filename=xfile)
noise=HMUX(filename=nfile)
if diag:
print("file: ",os.path.basename(xrays.hdffile))
for channel in xrays.channels:
print(" channel: ",channel)
fff=" freqs#: "+len(xrays.channel[channel].freqs)*' %8d'
print(fff % tuple(xrays.channel[channel].freqs))
hz=xrays.channel[channel].freq_conf['freq']
fff=" kHz: "+len(hz)*' %8.2f'
print(fff % tuple(hz))
else:
print("X-ray file: %s" % os.path.basename(xfile))
print("Noise file: %s" % os.path.basename(nfile))
if names:
return xfile,nfile
else:
return xrays,noise
[docs]def NEP(xrays=None,noise=None,fpath='/stage/ltsboldat5/LC_XFDM',\
xexclude=0.1,xlexclude=0.2,ninclude=0.2,cutfreq=30000.0,diag=True,pdf=None,\
channel=None,freq=None,logfreq=False,pulsepos=True,rot=False,ascii=False,\
dtype=False,prlen=None,prclip='end'):
'''
Compute NEP for xray-data file and noise file combination.
Kwargs:
* `xrays` = HDF file object with xrays. If not provided will be asked for.
* `noise` = HDF file object with noise data. If not provided, accompanying noise file with xray file will be used.
* `fpath` = path for files
* `xexclude` = fraction highest intensity of xray records to exclude. Should be tuned to exclude doubel and triple events in record.
* `xlexclude` = fraction lowest intensity of xray records to exclude. Should be tuned to exclude low events in record.
* `ninclude` = fraction of lowest noise records to use.
* `cutfreq` = frequency cutoff (Hz) for NEP computation
* `diag` = print and plot diagnostics information
* `pdf` = if not 'None', make pdf files, otherwise, show plots interactively
* `channel` = channel to use. If not provided, use first available channel
* `freq` = frequency to use. If not provided, use first avialable frequency
* `logfreq` = not used
* `pulsepos` = if True, only use events at proper trigger position, otherwise use all events
* `rot` = rotate record in I/Q plane for minimum Q signal
* `ascii` = if True, write ascii files of averaged pulse and pulse/noise spectra and NEP
* `dtype` = if True, convert data to I-Q when file dtype is "ampl-phase"
* `prlen` = if not 'None', processing record length
* `prclip` = set to clip record at start or end, when `prlen` is set
Returns:
* `FWHM` = the computed NEP (eV)
* `freq` = the configured frequency (kHz) of pixel
* `amplt` = the configured amplitude of the pixel
* `NEP` = the NEP (power) as function of frequency
* `freq` = the frequency axis for the NEP
* `gbwp` = the configured gbwp for the pixel
* `risetime` = the average pulse rise time
* `falltime` = the average pulse fall time
'''
ps=False
if pdf is not None:
ps=True
if ( xrays is None) or ( noise is None ):
xrays,noise=NEPfiles(fpath=fpath)
fnam=os.path.basename(xrays.hdffile)
nnam=os.path.basename(noise.hdffile)
if diag:
print("X-ray file: ",fnam)
print("Noise file: ",nnam)
if channel is None:
xchan=xrays.channels[0]
nchan=noise.channels[0]
else:
xchan=channel
nchan=channel
if freq is None:
xfreq=xrays.channel[xchan].freqs[0]
nfreq=noise.channel[nchan].freqs[0]
else:
xfreq=freq
nfreq=freq
conf_tab=noise.channel[nchan].freq_conf
# find index (ifreq) for given frequency (pixel) number
for ifreq,fff in enumerate(noise.channel[nchan].freqs):
if fff == nfreq:
break
flip=False
khz=conf_tab['freq'][ifreq]
if khz == 0.0:
print("Data for channel %d, freq %d has zero frequency. Data are DC data" % (channel,freq))
flip=positivepulse(xrays.channel[xchan].freq[xfreq],None) # flip data in record for DC
# (frequency=0.0) positive pulse data
# names for ascii output files
if ascii:
bnam=os.path.basename(xrays.hdffile).rsplit('.',1)[0]
fa=('_pix%d' % xfreq).strip()
plsascii='NEPpulse_'+bnam+fa+'.txt'
frqascii='NEPspec_'+bnam+fa+'.txt'
#
xsamples=xrays.channel[xchan].freq[xfreq][0][:,0].size
xrate=int(xrays.channel[xchan].freq[xfreq].attrs["sample_rate"])
xrecords=int(xrays.channel[xchan].attrs["Nevents"])
nsamples=noise.channel[nchan].freq[nfreq][0][:,0].size
nrate=int(noise.channel[nchan].freq[nfreq].attrs["sample_rate"])
nrecords=int(noise.channel[nchan].attrs["Nevents"])
# check processing record length
if prlen is not None:
if ( prlen > xsamples ) or ( prlen > nsamples ):
print("Set processing record length: ",prlen)
print(" Xray record length: ",xsamples)
print(" Noise record length: ",nsamples)
sys.exit("Error, set processing record length is larger than real record lengths")
xsamples=prlen
nsamples=prlen
if diag:
print("Overview for channel %d, frequency %d" % (xchan,xfreq))
print("Number of samples in X-ray record: %d, rate: %d" % (xsamples,xrate))
print("Number of samples in noise record: %d, rate: %d" % (nsamples,nrate))
else:
print("Process channel %d, frequency (pixel) %d" % (xchan,xfreq))
freq=numpy.arange(nsamples//2)/float(nsamples)*nrate # frequency scale
freqbin=freq[2]-freq[1] # frequency bin
timax=numpy.arange(xsamples)/xrate # time axis
if xrate != nrate:
sys.exit("Error, sample rates of files are not equal")
if xsamples != nsamples:
print(" Xray record length: ",xsamples)
print(" Noise record length: ",nsamples)
sys.exit("Error, number of record samples of files are not equal")
tottime=float(nsamples)/float(nrate) # total time of noise record
# find data_type (ampl-phase or I-Q )
dattype=cu(xrays.channel[xchan].freq[xfreq].attrs['data_type'])
ampphase=False
scaleph=None
if ( dattype == 'ampl-phase' ) and ( dtype == True ) :
ampphase=True
scaleph=float(cu(xrays.channel[xchan].freq[xfreq].attrs['scale_phase_rad']))
print('data_type: ',dtype,scaleph)
#
# if rotation needed, find proper phase rotation from first 100 records
#
if rot:
aphase=0.0
for i,rrecord in enumerate(xrays.channel[xchan].freq[xfreq]):
if ampphase:
rrecord=ap2iq(rrecord,scaleph)
aphase=aphase+(IQphase(rrecord[:,0],rrecord[:,1])+2.0*numpy.pi)
if i > 99:
break
aphase=aphase/float(i+1)-2.0*numpy.pi
print("Average phase correction: ",aphase)
aphase=-aphase # phase correction
#
# first get index of records with real X-ray events in them
#
xpos=numpy.zeros(xrecords)
xptp=numpy.zeros(xrecords)
bl=50 # number of bins for baseline level
dc=100 # seperation for multi peaks
for i,record in enumerate(xrays.channel[xchan].freq[xfreq]):
if ampphase:
record=ap2iq(record,scaleph)
if (i % 1000) == 0:
print(" read xray record: %d" % i)
if rot:
irecord,qr=IQphrot(record[:,0],record[:,1],aphase)
else:
irecord=record[:,0]
if flip or ( irecord[0] < 0):
irecord=-irecord
if prlen is not None:
if prclip == 'end':
irecord=irecord[:prlen]
else:
irecord=irecord[-prlen:]
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
indx,=numpy.where(irecord < tt) # record below threshld
if indx.size < 10:
continue
dd=indx[1:]-indx[0:-1] # peak seperation at threshold
if dd.max() < dc: # only select single peaks
xptp[i]=numpy.ptp(irecord) # look at maximum difference within record
xpos[i]=numpy.argmin(irecord) # position of pulse
vindx,=numpy.where(xptp > 300.0) # only valid records
hh,bb=numpy.histogram(xptp[vindx],bins=5000) # make histogram
hc=numpy.cumsum(hh) # cumulative sum
tindx,=numpy.where( hc > (0.75*vindx.size) ) # look for 0.50 percent fraction
cutoff=bb[tindx[0]] # set cutoff threshold
print("cutoff for valid X-ray records: ",cutoff)
# sptp=numpy.sort(xptp)
# ic=int(xrecords/len(xrays.channel[xchan].freq)/10.)+1 # cutoff in event record-ptp (freq=number of pixels)
# cutoff=0.75*sptp[-ic] # cutoff at 0.75 the max ptp
hpos,ppos=numpy.histogram(xpos[vindx],bins=xsamples)
ppos=(ppos[0:-1]+ppos[1:])/2
mpos=ppos[numpy.argmax(hpos)] # position of triggered pulse
if pulsepos:
cindx,=numpy.where( ( xptp > cutoff ) &
( ( xpos > (mpos-5)) & ( xpos < (mpos+5) ) ) ) # index of records with proper X-ray events
else:
cindx,=numpy.where( ( xptp > cutoff ) &
( ( xpos > 100 ) & ( xpos < int(0.8*nsamples ) ) ) ) # index of records with proper X-ray events
if cindx.size == 0:
print("Insufficient X-ray events found")
return (None,conf_tab['freq'][ifreq],conf_tab['ampl'][ifreq],None,None,conf_tab['gbwp'][ifreq])
else:
print("Number of valid X-ray records: ",cindx.size)
if diag:
f,ax=plt.subplots(2,sharex=False)
ax[0].set_xlabel('ptp')
ax[0].set_ylabel('N')
ax[0].set_title(fnam)
maxh=numpy.amax(xptp)
minh=numpy.amin(xptp)
ax[0].hist(xptp,bins=200,range=[minh,maxh],color='b')
ax[0].hist(xptp[cindx],bins=200,range=[minh,maxh],color='r')
ax[1].set_xlabel('pos')
ax[1].set_ylabel('N')
ax[1].set_title(fnam)
maxh=numpy.amax(xpos)
minh=numpy.amin(xpos)
ax[1].hist(xpos,bins=200,range=[minh,maxh],color='b')
ax[1].hist(xpos[cindx],bins=200,range=[minh,maxh],color='r')
print(" Total number of records: ",xrecords)
print("Number of records with X-ray events: ",cindx.size)
print("Cutoff: ",cutoff)
plt.show()
plt.close('all')
#
# get x-ray records with x-fraction of lowest intensities (to prevent double events from being used)
# always exclude lowest 20% fraction
#
print("reread lowest intensity fraction")
xint=numpy.zeros(cindx.size,dtype=float)
for i,indx in enumerate(cindx):
if (i % 1000) == 0:
print(" read record: %d" % cindx[i])
rrecord=xrays.channel[xchan].freq[xfreq][indx]
if ampphase:
rrecord=ap2iq(rrecord,scaleph)
if rot:
record,qr=IQphrot(rrecord[:,0],rrecord[:,1],aphase) # rotated I-signal
else:
record=rrecord[:,0] # only I-signal
if flip or ( record[0] < 0):
record=-record
if prlen is not None:
if prclip == 'end':
record=record[:prlen]
else:
record=record[-prlen:]
# xmax=record.max() # get maximum as measure of baseline
xmax=(record[0:20].sum()+record[-21:-1].sum())/40.0 # get baseline from outermost 20 points
xint[i]=(xmax-record).sum() # integrate baseline with event subtracted
ss=numpy.sort(xint) # sort maxima
nt=int(xint.size*xexclude) # excluded fraction in number of records
maxthreshold=ss[-nt-1]+10.0 # threshold in maxima
mt=int(xint.size*xlexclude) # always exclude lowest `xlexclude` %
if mt > (xint.size-1):
minthreshold=maxthreshold-20.0
else:
minthreshold=ss[mt]-1.0
print("X-ray summed event selection thresholds (min,max): %12.5e %12.5e" % (minthreshold,maxthreshold))
xindx,=numpy.where( ( xint < maxthreshold ) & \
( xint > minthreshold ) ) # list of record numbers to be used
if xindx.size == 0:
print("Insufficient single X-ray events found")
return (None,conf_tab['freq'][ifreq],conf_tab['ampl'][ifreq],None,None,conf_tab['gbwp'][ifreq])
print("Number of valid X-ray events found: ",xindx.size)
#
# retrieve average of selected x-ray records
#
apulse=numpy.zeros(xsamples,dtype=float)
qpulse=numpy.zeros(xsamples,dtype=float)
xsamples2=xsamples//2
pspectrum=numpy.zeros(xsamples2,dtype=float)
for indx in xindx:
rrecord=xrays.channel[xchan].freq[xfreq][cindx[indx]]
if ampphase:
rrecord=ap2iq(rrecord,scaleph)
if rot:
record,qr=IQphrot(rrecord[:,0],rrecord[:,1],aphase) # rotated I-signal
else:
record=rrecord[:,0] # I-signal, used for NEP
if flip:
qr=numpy.zeros(record.size)
else:
qr=rrecord[:,1] # Q-signal
if flip or ( record[0] ) < 0:
record=-record
qr=-qr
if prlen is not None:
if prclip == 'end':
record=record[:prlen]
qr=qr[:prlen]
else:
record=record[-prlen:]
qr=qr[-prlen:]
pp=record.astype(float)
apulse=apulse+pp
qpulse=qpulse+qr.astype(float)
pspectrum=pspectrum+numpy.absolute(numpy.fft.fft(pp))[0:xsamples2] # get half of the spectrum, leave the symmetric part
apulse=apulse/float(xindx.size)
qpulse=qpulse/float(xindx.size)
pspectrum=pspectrum/float(xindx.size)
pspectrum[0]=0.0 # ignore integral value
if diag or ps:
ptitle="%s ch=%d, px=%d" % (fnam.split('.')[0],xchan,xfreq)
plots=ShowPlots(logfreq=logfreq,title=ptitle)
plots.addplot(1,apulse,title="average pulse",ylabel="Pulse I-data")
plots.addplot(3,pspectrum,x=freq/1000.0,title="average pulse spectum",\
ylabel="Pulse power",xlabel="Frequency (kHz)")
if ascii:
plsampl=numpy.sqrt(apulse**2+qpulse**2)
plsphas=numpy.arctan2(qpulse,apulse)
plstxt=open(plsascii,'w')
hh=['Time','Amplitude','Phase']
head='#'
for ih in hh:
head=head+('%15.15s' % ih)
plstxt.write(head+'\n')
for i,time in enumerate(timax):
ll=' %14.5e %14.5e %14.5e\n' % (time,plsampl[i],plsphas[i])
plstxt.write(ll)
plstxt.close()
#
# compute rise and fall time of average pulse
#
nbck=15
bckindx=numpy.concatenate((numpy.arange(nbck),numpy.arange(apulse.size-nbck,apulse.size)))
a,b=numpy.polyfit(bckindx,apulse[bckindx],1)
pulsex=numpy.arange(apulse.size)
spulse=apulse-(a*pulsex+b)
risetime,falltime=getrisetime(pulsex,spulse)
risetime=risetime/xrate
falltime=falltime/xrate
print("Pixel: %d :: risetime: %10.7f (ms), falltime: %10.7f (ms)" % (xfreq,risetime*1000.0,falltime*1000.0))
#
# get noise records, using only lowest included fraction
#
print("read noise records")
next=numpy.zeros(nrecords)
for i,rrecord in enumerate(noise.channel[nchan].freq[nfreq]):
if ampphase:
rrecord=ap2iq(rrecord,scaleph)
if (i % 1000) == 0:
print(" read record: %d" % i)
if rot:
record,qr=IQphrot(rrecord[:,0],rrecord[:,1],aphase) # rotated I-signal
else:
record=rrecord[:,0] # only I-signal
if flip or ( record[0] ) < 0:
record=-record
if prlen is not None:
if prclip == 'end':
record=record[:prlen]
else:
record=record[-prlen:]
next[i]=numpy.ptp(record) # look at maximum difference within record
vindx,=numpy.where(next != 0.0) # discard dummy records which have only zeroes
ss=numpy.sort(next[vindx]) # sort maxima
mt=int(vindx.size*ninclude) # only include lowest 10%
minthreshold=ss[mt]
nindx,=numpy.where( ( next <= minthreshold ) & ( next != 0.0 ) ) # list of record numbers to be used
#
# retrieve selected noise records
#
nsamples2=nsamples//2
nspectrum=numpy.zeros(nsamples2)
anoise=numpy.zeros(nsamples)
print("reread selected noise records")
for i,j in enumerate(nindx):
if ( i % 1000 ) == 0:
print(" read record no: %d" % i)
rrecord=noise.channel[nchan].freq[nfreq][j]
if ampphase:
rrecord=ap2iq(rrecord,scaleph)
if rot:
record,qr=IQphrot(rrecord[:,0],rrecord[:,1],aphase) # rotated I-signal
else:
record=rrecord[:,0] # only I-signal
if flip or ( record[0] ) < 0:
record=-record
if prlen is not None:
if prclip == 'end':
record=record[:prlen]
else:
record=record[-prlen:]
anoise=anoise+record
nspectrum=nspectrum+(numpy.absolute(numpy.fft.fft(record))[0:nsamples2])**2 # get half of the power spectrum, leave the symmetric part
if nindx.size == 0:
print("Insufficient valid noise records found")
sys.exit()
print("Number of valid noise record found: ",nindx.size)
anoise=anoise/float(nindx.size)
nspectrum=nspectrum/float(nindx.size)
nspectrum[0]=0.0
if diag or ps:
print(" Number of noise records used: ",nindx.size)
plots.addplot(2,anoise,title="average noise",ylabel="Noise I-data")
plots.addplot(4,nspectrum,x=freq/1000.0,title="average noise power",\
ylabel="Power",xlabel="Frequency (kHz)")
#
# now the compute NEP using the X-ray pulse spectrum and the noise spectrum:
#
# starting points:
# pspectrum = average spectrum of pulse ( abs(fft(data) )
# nspectrum = average power spectrum of noise ( abs(fft(data))^2 )
# tottime = total time of noise record ( number_of_samples/sampling_rate )
#
pspectrum[0]=1.0 # prevent divide by zero
nspec=numpy.sqrt(nspectrum*2.0*tottime) # noise average power density
# multiply by 2, to get full power of complete spectrum
# (spectrum is only the first half symmetric part)
photonenergy=(5.89875*16.2+5.88765*8.2+6.49045*2.85)/(16.2+8.2+2.85)*1000.0 # average photon energy (alpha+beta lines)
si=pspectrum/photonenergy*tottime # divide by photon energy to get scale
NEPev=nspec/si # NEP
if ascii:
frqtxt=open(frqascii,'w')
hh=['Freq(Hz)','NEP','Pulse','Noise']
head='#'
for ih in hh:
head=head+('%15.15s' % ih)
frqtxt.write(head+'\n')
for i,fff in enumerate(freq):
ll=' %14.5e %14.5e %14.5e %14.5e\n' % (fff,NEPev[i]*(1.6E-19),pspectrum[i],nspectrum[i])
frqtxt.write(ll)
frqtxt.close()
cindx,=numpy.where(freq < cutfreq) # compute NEP only until cutoff frequency
NEPev=NEPev[cindx]
if diag or ps: # plot NEP in SI unit; multiply by electron charge
plots.addplot(5,NEPev*(1.6E-19),title="",ylabel="NEP",x=freq[cindx]/1000.0,xlabel="Frequency (kHz)")
dErms=1.0/numpy.sqrt(numpy.sum(4.0/NEPev[1:]**2*freqbin)) # integrate 4.0/NEP^2
dEfwhm=2.0*numpy.sqrt(2.0*numpy.log(2.0))*dErms # energy FWHM
print(" dE(rms) = %.2f eV" % dErms)
print("dE(fwhm) = %.2f eV" % dEfwhm)
if diag or ps:
atxt=' NEP= %.2f eV\n freq(%d)=%7.2f (kHz)\n Ampl=%7.2f' % \
(dEfwhm,nfreq,conf_tab['freq'][ifreq],conf_tab['ampl'][ifreq])
plots.annotate(5,atxt)
if ps:
# psnam='NEP_'+fnam.split('.')[0]+('_ch%2.2d_fr%3.3d_' % (xchan,xfreq))+'.ps'
plots.show(psfile=pdf)
else:
plots.show()
return (dEfwhm,conf_tab['freq'][ifreq],conf_tab['ampl'][ifreq],(NEPev*(1.6E-19)),freq[cindx],\
conf_tab['gbwp'][ifreq],risetime,falltime)
[docs]def NEPs(xrays=None,noise=None,fpath='/stage/ltsboldat5/LC_XFDM',\
xexclude=0.3,xlexclude=0.2,ninclude=0.2,cutfreq=30000.0,diag=True,pdf=None,\
channel=None,freqs=None,logfreq=False,summary=False,pulsepos=True,sumret=False,\
rot=False,ascii=False,dtype=False,prlen=None,prclip='end'):
'''
Get NEP for all frequencies (pixels)
Kwargs:
* `xrays` = HDF file object with xrays. If not provided will be asked for.
* `noise` = HDF file object with noise data. If not provided, accompanying noise file with xray file will be used.
* `fpath` = path for files
* `xexclude` = fraction highest intensity of xray records to exclude. Should be tuned to exclude doubel and triple events in record.
* `xlexclude` = fraction lowest intensity of xray records to exclude. Should be tuned to exclude low events in record.
* `ninclude` = fraction of lowest noise records to use.
* `cutfreq` = frequency cutoff (Hz) for NEP computation
* `diag` = print and plot diagnostics information
* `pdf` = if not 'None', make pdf files, otherwise, show plots interactively
* `channel` = channel to use. If not provided, use first available channel
* `freqs` = frequencies to use. If not provided, process all frequencies
* `logfreq` = not used
* `summary` = if True, provide summary of results in text file
* `pulsepos` = if True, only use events at proper trigger position, otherwise use all events
* `sumret` = if True, return summary as additional output of routine; do not write to file
(only applicable if freqs contains only one frequency)
* `rot` = if True, rotate I/Q records for minimum Q signal, based on average phase of 100 records
* `ascii` = if True, write ascii files of averaged pulse and pulse/noise spectra and NEP
* `dtype` = if True, convert records to I-Q when the file dtype is "ampl-phase"
* `prlen` = if not 'None', processing record length
* `prclip` = set to clip record at start or end, when `prlen` is set
Returns:
* `NEPdata` = list of the computed NEP (eV) for the given frequencies (pixels)
'''
ps=False
if pdf is not None:
ps=True
# get NEP data for all freqs
if channel is None:
channel=xrays.channels[0]
if freqs is None:
freqs=[None]
nepdata=numpy.zeros(len(freqs),dtype=float)
khz=numpy.zeros(len(freqs),dtype=float)
ampl=numpy.zeros(len(freqs),dtype=float)
nepc={}
nepx={}
gbwp={}
risetimes={}
falltimes={}
if summary:
sname='NEP_'+os.path.basename(xrays.hdffile).split('.')[0]+'.txt'
sfile=open(sname,'w')
sfile.write('#file: %s\n#\n' % (os.path.basename(xrays.hdffile).split('.')[0]))
sfile.write('# i px freq(khz) V_bias gbwp NEP(fwhm) (eV) risetime(ms) falltime(ms)\n#\n')
for j,freq in enumerate(freqs):
nepdata[j],khz[j],ampl[j],nepc[j],nepx[j],gbwp[j],risetimes[j],falltimes[j]=NEP(xrays=xrays,noise=noise,fpath=fpath,\
xexclude=xexclude,xlexclude=xlexclude,ninclude=ninclude,\
cutfreq=cutfreq,diag=diag,pdf=pdf,\
channel=channel,freq=freq,logfreq=logfreq,pulsepos=pulsepos,rot=rot,\
ascii=ascii,dtype=dtype,prlen=prlen,prclip=prclip)
sline='%3d %3d %10.2f %10.2f %10.4f %11.2f %14.6f %12.6f\n' % \
(j,freq,khz[j],ampl[j],gbwp[j],nepdata[j],risetimes[j]*1000.0,falltimes[j]*1000.0)
if summary:
sfile.write(sline)
if summary:
sfile.close()
if (diag or ps) and (len(freqs) > 1):
fnam=os.path.basename(xrays.hdffile)
f,ax=plt.subplots(1)
ax.set_title(fnam)
ax.set_xlabel('Frequency (kHz)')
ax.set_ylabel('NEP')
for j,freq in enumerate(freqs):
if ( nepx[j] is not None ) and ( nepc[j] is not None ):
xindx,=numpy.where(( nepx[j] > 10.0 ) & ( nepx[j] < 10000.0) )
ax.plot((nepx[j][xindx]/1000.),nepc[j][xindx],label=('px: %d%8.1f (kHz)' % (freq,khz[j])))
ax.legend(loc='upper left')
if int(matplotlib.__version__[0]) >= 3:
ax.set_xscale('symlog',linthresh=0.01)
else:
ax.set_xscale('symlog',linthreshx=0.01)
ax.set_yscale('log')
if pdf is not None:
plt.savefig(pdf,format='pdf',papertype='a4')
else:
plt.show()
plt.close('all')
f,ax=plt.subplots(1)
ax.set_xlabel('Frequency (kHz)')
ax.set_ylabel('dE(fwhm) (eV)')
ax.set_title(fnam)
ax.plot(khz,nepdata,'bo')
axx=list(ax.axis())
ax.axis([(min(khz)-100.0),(max(khz)+100.0),axx[2],axx[3]])
if pdf is not None:
plt.savefig(pdf,format='pdf',papertype='a4')
else:
plt.show()
plt.close('all')
if sumret:
return nepdata,sline
else:
return nepdata,"\n"