#%TITLE% PSEARCH.MAC
#%NAME%
# Peak searching and fitting with energy calibration
#%DESCRIPTION%
# Psearch.mac is the spec front end to the psearch external program to
# search and fit peaks. The basic algorithm consists of removing first
# all peaks (above a certain statistical threshold), fitting a back ground
# function, isolating groups of peaks and trying to fit multipletts in
# an iterative procedure.
# The macros in this macro set are not directly used by the user, but
# provide the peak fitting to other macro sets.
# The macro set also provides the necessary functions to calibrate
# spectra in energy.
# The initial width estimation and the Energy-channel relation are: %BR%
# E=ea + eb ch + ec ch^2 ; FWHM = sqrt( wa^2 + wb E + (wc E)^2 )
#
#%EXAMPLE%
# These example are not really meant to be typed in by the user. They
# should give you an example how to use the psearch function in your
# macros.
#%PRE%
# %B%ps_calib%B%
# FWHM = sqr[a^2 + b*ch]
# Term (a) in channels (5)?
# Term (b) in channels (0.04)?
# 54 Peaks found
# No Position Intensity FWHM Flag
# 1 51.284000 599.0 3.780 2
# 2 53.000000 368.0 3.645 2
# ......
# 8 232.835007 7014.0 8.290 3
# 9 248.266006 26738.0 6.822 3
# ....
# No of 1. peak (9)?
# Energy at peak 9, channel 248.266006 (Int:26738.0) (8.0441)?
# No of 2. peak (54)?
# Energy at peak 54, channel 3894.892090 (Int:1178503.0) (122.06)?
# No of 3. peak (0 to end) (0)?
# Linear FIT E= a + b*ch
# CHI2: 0.000 ; a: 0.2817817 b: 0.03126613
# Do you accept the fit (YES)?
#
# The energy calibration is done. It is stored in the global
# variables PS_EA , PS_EB, PS_EC. You can use the conversion
# macros :
# ps_calcch <Energy> <Channel>
# ps_calcE <Channel> <Energy>
# ps_Etoch group1 elem1 group2 elem2
# ps_Etoch group1 elem1 group2 elem2
#
# %B%pssearch%B%
# Asks for initial width parameters and displays the peaks in
# Energy.
#%PRE%
#%SETUP%
# The program psearch must be accessible. (in the path)
#%IU% (program , arguments , array_in, array_out)
#%MDESC% like data_pipe but for arrays. This
# This needs the new version of the binary program psearch!!!!!
# Please update this program.
def old_array_pipe(program , arguments , array_in, array_out,nodata) '{
local array temp_in[nodata][1]
local array temp_out[200][5]
local retval
temp_in[][0] = array_in
retval = array_pipe(program, arguments, temp_in, temp_out)
array_out = temp_out
return retval;
}'
def _obsolete_array_pipe(program , arguments , array_in, array_out,nodata) '{
# Ca va sauter - we have to use groups to interface to psearch
local PS_GRPSEARCH PS_GRPDATA PS_GRPFIT
local array temp[nodata][2]
local retval
PS_GRPIN = 125 ; PS_GRPOUT = 126
data_grp(PS_GRPIN,nodata,2) # data : x y E
data_grp(PS_GRPOUT,200,5) # peaks : no position intensity fwhm flag
temp [][1] = array_in
data_read(temp,PS_GRPIN,0,0)
retval=data_pipe(program,arguments,PS_GRPIN,PS_GRPOUT)
# No Channel Counts FWHM Flag Energy FWHM-E
for (i=0;i<=data_info(PS_GRPOUT,"last");i++)
for (j=0;j<data_info(PS_GRPOUT,"elem");j++)
array_out[i][j] = data_get(PS_GRPOUT,i,j)
return retval;
}'
#%IU% (data, data_col, minch, no, a, b, c, wa, wb, wc)
#%MDESC% Does the peak search for the array given in data. The data
# holds only counts. The corresponding channel numbers go from minch
# to minch + no - 1. The a, b and c parameters provide means to input
# an enery calibration (E = a + b * ch + c * ch * ch).
# The routine needs an estimate for the peak width. The parameters for
# this estimate are given by W(E) = sqrt (wa**2 + wb*E + (wc*E)**2).
# The results will be in the global array PS_PEAKS.
# The individual columns are: %BR%
# No Channel Counts FWHM Flag Energy FWHM-E. %BR%
# The calculation will always be done in
# channels, only the width parameters will be ajusted for the given
# calibration.
def ps_psearch (data, minch, no, a, b, c, wa, wb, wc, eunit) '{
local acorr bcorr offset nopts minch
local lwa lwb lwc
array PS_PEAKS[200][7]
if (eunit) {
lwa = sqrt(fabs(wa*wa + wb*a + wc * wc * c *c))/b
lwb = wb/b + 2*wc *a /b
} else {
lwa = wa
lwb = wb
}
PS_NOPEAKS = old_array_pipe("psearch",sprintf("%g %g %g %g %g %g %d",\
minch,1,0,lwa,lwb,wc,no),data,PS_PEAKS,no)
ps_chtoE (PS_PEAKS[][1],PS_PEAKS[][5])
PS_PEAKS[][6]= b * PS_PEAKS[][3]
return PS_NOPEAKS
} '
#%UU% (data arr, first channel, no data, wa, wb, wc, usee, minint, maxerror)
#%MDESC%
# Searches the current data array for peaks.
# The Energy calibration must have been done before.
# The result is put into PS_PEAKS.
# Can be used interactively by the user or from other macros
# with given width parameters. Only if called with no parameter
# it will produce output.
# The number of found peaks is stored in PS_NOPEAKS.
def pssearch (data, firstch, nodata, wa, wb, wc, usee ,minint, maxerror) '{
userinput = length(wa)
if (!userinput) {
printf ("Current Calibration: E = %10g+%10g*ch+%10g*ch^2\n",\
PS_EA,PS_EB,PS_EC)
if (PS_USEE = yesno ("Estimate width as function of Energy",PS_USEE)) {
printf("FWHM = sqrt[a^2 + b * E + (c * E)^2]\n")
} else {
printf("FWHM = sqrt[a^2 + b * channel + (c * channel)^2]\n")
}
PS_WA=getval("Term (a) ",PS_WA)
PS_WB=getval("Term (b) ",PS_WB)
PS_WC=getval("Term (c) ",PS_WC)
PS_MININT =getval ("Minimum Intensity (counts)",PS_MININT)
PS_MAXERROR = getval ("Max Error in FWHM in %",PS_MAXERROR)
} else {
PS_WA = wa ; PS_WB = wb ;PS_WC = wc ; PS_USEE = usee
PS_MININT= minint ; PS_MAXERROR= maxerror
}
PS_NOPEAKS = ps_psearch (data, firstch, nodata,\
PS_EA,PS_EB,PS_EC,PS_WA,PS_WB,PS_WC,PS_USEE )
ps_reject (PS_MININT,PS_WA,PS_WB,PS_WC,PS_MAXERROR,0,0)
if (!userinput) {
ps_dump
}
return PS_NOPEAKS;
}'
#%IU% ( minint, wa, wb, wc, %error, minE, maxE )
#%MDESC%
# FWHM = sqr[wa^2 + wb*E + (wc*E)^2] %BR%
# Deletes all peaks from the PS_PEAKS which do not
# fullfill the given conditions
#%UL%
# %LI% minint < I
# %LI% minE < E < maxE
# %LI% FWHM is in %error of specified FWHM
#%XUL%
def ps_reject (minint, wa, wb, wc, perror, minE, maxE) '{
local fwhm
if (0 == PS_NOPEAKS) return 0
local array reject[PS_NOPEAKS]
for (ii=0;ii<PS_NOPEAKS;ii++) {
if (minint > 0)
reject[ii] = (PS_PEAKS[ii][2] < minint)
else
reject[ii] = 0
if (minE > 0)
reject[ii] += (PS_PEAKS[ii][5] < minE)
if (maxE > 0)
reject[ii] += (PS_PEAKS[ii][5] > maxE)
if (perror > 0) {
fwhm = sqrt(pow(wa,2) + wb*PS_PEAKS[ii][1] + pow(wc * PS_PEAKS[ii][1],2))
reject[ii] += fabs(PS_PEAKS[ii][3] - fwhm) > fwhm / 100 * perror
}
}
for (ii=0, jj=0; ii< PS_NOPEAKS; ii++)
if (!reject[ii])
PS_PEAKS[jj++][] = PS_PEAKS[ii][]
PS_NOPEAKS=jj
array_op("fill",PS_PEAKS[0:PS_NOPEAKS][0],1)
PS_PEAKS[0:PS_NOPEAKS][0] += 1
return PS_NOPEAKS
}'
#%UU%
#%MDESC%
# Displays a long list of peaks found
def ps_dump '
printf ("%d Peaks found\n",PS_NOPEAKS)
printf ("%5s %15s %10s %10s %10s %10s %3s\n",\
"No","Channel","Int Counts","FWHM (ch)","E","FWHM (E)","Flag")
for (ii=0;ii<PS_NOPEAKS;ii++) {
printf("%5d %15f %10.1f %10.3f %10.3f %10.4f %3d\n",\
PS_PEAKS[ii][0],PS_PEAKS[ii][1],PS_PEAKS[ii][2],PS_PEAKS[ii][3],\
PS_PEAKS[ii][5],PS_PEAKS[ii][6],PS_PEAKS[ii][4])
}
'
def ps_testdata '
data_read("ps_testdata.dat",0,0,0)
array PS_DATA[8192][2]
for (ii=0 ; ii< 2047 ; ii ++) {
PS_DATA[ii][0]=data_get(0,ii,0)
PS_DATA[ii][1]=data_get(0,ii,1)
}
data_plot(PS_DATA)
PS_NODATA = 2048
'
#%UU%
#%MDESC%
# Displays a short list of peaks found
def ps_sdump '
for (ii=0;ii<PS_NOPEAKS;ii++) {
printf ("%9.4f %9.1f ",PS_PEAKS[ii][1],PS_PEAKS[ii][2])
if (int((ii+1)/4) == (ii+1)/4) print
}
print
'
#%UU%
#%MDESC%
# Displays the calibration parameters, number of peaks and a short
# list of the last found peaks.
def ps_show '
printf ("Current Calibration: E = %10g+%10g*ch+%10g*ch^2\n",\
PS_EA,PS_EB,PS_EC)
printf ("Peaks found : (Energy Intensity) \n")
ps_sdump
'
#%UU% (data, firstch, nodata)
#%MDESC%
# Does a peaksearch on array data. Then
# asks the user for calibration energy of some peaks and fits a
# linear or quadratic function.
def ps_calib (data, firstch, nodata) '{
pssearch (data, firstch, nodata)
if (PS_NOPEAKS > 0) {
ps_calibration
}
}'
#%IU%
#%MDESC% see ps_calib
def ps_calibration '{
global PEAKN PEAKC PEAKE
local fitno
for (ii=0 ; ; ii++) {
PEAKN[ii] = getval(sprintf("No of %d. peak %s",ii+1,\
(ii>1)?"(0 to end)":""),PEAKN[ii]+1)-1
if (PEAKN[ii] == -1 ) {
fitno = ii
break;
}
PEAKC[ii] = PS_PEAKS[PEAKN[ii]][1]
PEAKE[ii] = getval(sprintf("Energy at peak %d, channel %f (Int:%.1f)",\
PEAKN[ii]+1,PEAKC[ii],PS_PEAKS[PEAKN[ii]][2]),PEAKE[ii])
}
if (fitno < 2) {
print "ERROR : I need at least two points to calibrate"
exit
}
ps_ofit (PEAKE,PEAKC,fitno, (fitno>2) ? 3 : 2 ,1)
if (yesno("Do you accept the fit",1)) {
PS_EA=PS_PAR[0]; PS_EB=PS_PAR[1] ; PS_EC=PS_PAR[2]
}
}'
#%UU%
#%MDESC%
# Shows the calibration parameters and lets the user change them
def ps_parinput '
{
printf ("Current Calibration: E = %10g+%10g*ch+%10g*ch^2\n",\
PS_EA,PS_EB,PS_EC)
PS_EA = getval("Parameter A :",PS_EA)
PS_EB = getval("Parameter B (*ch):",PS_EB)
PS_EC = getval("Parameter C (*ch^2):",PS_EC)
}
'
#%IU%
#%MDESC%
# Same functionality as ps_fit but for old type arrays
def ps_ofit (ya,xa,no,order,verbose) '{
local array x[no] , y[no]
for (ii = 0 ; ii < no ; ii++ ) {
x[ii] = xa[ii]
y[ii] = ya[ii]
}
ps_fit (y,x,no,order,verbose)
}'
#%IU% (y, x, no of data points, order, verbose flag)
#%MDESC%
# Do a linear or quadratic fit y = a + b * x + c * x**2. The
# (a,b,c) parameters are returned in the global array PS_PAR. The return
# value of the function is the chi2 of the fit. If the verbose flag is
# set, the fitvalues are printed to the screen.
def ps_fit (y,x,no,order,verbose) '{
global PS_PAR
local chi2
local array fitarray[4][no]
fitarray[0][] = y ; fitarray[1][] = 1 ; fitarray[2][] = x ; fitarray[3][] = x*x
p fitarray
array_op("col_wise",fitarray,1)
if (order == 2) {
chi2 = array_fit (PS_PAR,fitarray[0:2][])
PS_PAR[2]=0
if (verbose) {
printf ("Linear FIT E= a + b*ch \n")
printf ("CHI2: %10.3f ; a: %15.7g b: %15.7g\n",chi2,PS_PAR[0],PS_PAR[1])
}
} else {
chi2 = array_fit (PS_PAR,fitarray)
if (verbose) {
printf ("Quadr. FIT E= a + b*ch + c*ch^2 \n")
printf ("CHI2: %10.3f ; a: %15.7g b: %15.7g c: %15.7g \n",\
chi2,PS_PAR[0],PS_PAR[1],PS_PAR[2])
}
}
return chi2
}'
# Transformation function which can be used by other programs
#%UU% <Channel> <Energy>
#%MDESC%
# Converts from Energy to channel (channel to Energy). The first
# parameter is the value to convert and the second parameter a variable
# name where to put the converted value.
def ps_calcE (channel) '{
return PS_EA + PS_EB * channel + PS_EC * channel * channel
}'
#%UU% (Energy)
#%MDESC% return the corresponding channel. Suppose that there is only
# a small quadratic correction
def ps_calcch (energy) '{
if (PS_EC == 0) {
return ((energy - PS_EA) / PS_EB)
} else if (PS_EB > 0 ) {
return ((-PS_EB+sqrt(pow(PS_EB,2)-4*(PS_EA-energy)*PS_EC))/(2*PS_EC))
} else {
return ((-PS_EB-sqrt(pow(PS_EB,2)-4*(PS_EA-energy)*PS_EC))/(2*PS_EC))
}
}'
#%UU% (channel array, Energy array)
#%MDESC%
# Converts all the channels in source array[][column1] to energies and puts
# these values into destination array[][column2]. The calibration function
# E=a+b ch + c ch^2 is used.
def ps_chtoE (charr, Earr) '{
Earr = PS_EA + charr * (PS_EB + PS_EC * charr)
}'
#%UU% (Energy array, channel array)
#%MDESC%
# Converts all the energies in sorce array[][column1] to channels and puts
# these values into destination array[][column2]. The taylor expansion for
# the parabolic calibration function E=a+b ch + c ch^2 is used.
def ps_Etoch (Earr, charr) '{
if (PS_EC == 0) {
charr = ((Earr - PS_EA) / PS_EB)
} else if (PS_EB > 0 ) {
charr = ((-PS_EB+sqrt(pow(PS_EB,2)-4*(PS_EA-Earr)*PS_EC))/(2*PS_EC))
} else {
charr = ((-PS_EB-sqrt(pow(PS_EB,2)-4*(PS_EA-Earr)*PS_EC))/(2*PS_EC))
}
}'
# %END%
# Thats just A(y-a) + B(y-a)**2 from Taylor , trans. to save group
# aa = -PS_EA/PS_EB(PS_EA*PS_EC/PS_EB/PS_EB + 1)
# bb = (1+2*PS_EC*PS_EA/PS_EB/PS_EB)/PS_EB
# cc = -PS_EC/PS_EB/PS_EB/PS_EB
# charr = aa + Earr * (bb + cc * charr)
#%UU%
#%MDESC%
# defines the global variables and looks for psearch executable
def pssetup '
global PS_EA PS_EB PS_EC PS_WA PS_WB PS_WC
global PS_MININT PS_MAXERROR PS_USEE PS_NOPEAKS
if (unix(sprintf("test -x %s/data_pipe/psearch",SPECD))) {
printf("Install psearch in %s/data_pipe first\n",SPECD)
exit
}
'
#%MACROS%
#%IMACROS%
#%ATTENTION%
# make sure psearch is in the path
#
#%DEPENDENCIES%
# The file psearch.mac has to be read in !done by: startup script
#%AUTHOR%
# PSEARCH.MAC - JK 2.94
## DATE Fri Apr 12 11:19:16 METDST 1996 REV 1.4;
#%TOC%
|
