Binning Results from MultiPeak Fitting
Hi there,
I have simulated x-ray absorbance data that comes from DFT calculations for a given molecule on which I perform MultiPeak Fitting Analysis. For example, if I have a molecule that has 6 atoms, my DFT calculation would give me 6 individual spectra at the end and I would perform MultiPeak Fitting(MPF) on each spectra to obtain the peak energies(i.e. peak position) and peak widths which are the parameters I intend to use for the "binning" process that follows.
For the "binning " process I want to use the peak energy obtained from MPF and match it to the peak energy obtained from my DFT calculations. Next, I want to use the peak width from MPF as my "bin width" for the DFT calculations in order to match the energies that are within the peak width.
In order to start this process I first run the function below that I made called "binEnergiesPrep" that essentially cycles through the various MPF folders and counts how many Peaks each folder contains. This is what the first loop does. The next part of the function involves placing the peak energies and peak widths for all peaks for every MPF run into corresponding waves that I iteratively named as peakEnergy and peakWidth. This is what the second loop does. The last two loops can be ignored at the moment since they are not relevant to the problem I'm having at the moment. This function prints out the peak Coefficients for every peak, for every MPF run. I've shown an example output below.
Function binEnergiesPrep(fnum) //This function basically prepares folders and waves from MP fitting for binning
Variable fnum //Number of MP fitting runs donw
SetDataFolder root:
Variable a
for(a=1;a<=fnum;a+=1) //This loop will create the waves where the peak energies and peak widths for each MP fitting run are stored
String peakEnergy = "pEnergy"+num2str(a)+"_c"+num2str(a) //Names the peak energy wave for MP run "a"
String peakWidth = "pWidth"+num2str(a)+"_c"+num2str(a) //Names the peak width wave for run "a"
Wave pEnergy = $peakEnergy
Wave pWidth = $peakWidth
Make/O/N=(fnum) numPeaks,$peakEnergy,$peakWidth //. The numPeaks wave tells me how many peaks are present in that particular MP run.Makes the waves according to convention above
endfor
print "Listing of Peaks"
Variable k
for(k=1;k<=fnum;k+=1) //This loop will go through each MPfitting run folder and create a list of the Peak Coefs waves.
String currentMPfolder = "MPF_SetFolder_"+ num2str(k)
SetDataFolder root:Packages:MultiPeakFit2:$currentMPfolder
String numPeaksinFit = WaveList("Peak*Coefs",";","")
Variable nPeaks = ItemsInList(numPeaksinFit,";") //How many peaks are there in this MPfitting run?
numPeaks[k-1] = nPeaks
print numPeaksinFit,nPeaks //Sanity check
endfor
SetDataFolder root:
Variable i
Variable j
for(i=1;i<=fnum;i+=1)
SetDataFolder root:
currentMPfolder = "MPF_SetFolder_"+ num2str(i) //Name of folder containing current atom MP fitting results
peakEnergy = "pEnergy"+num2str(i)+"_c"+num2str(i) //
peakWidth = "pWidth"+num2str(i)+"_c"+num2str(i)
Wave tempEnergy = $peakEnergy
Wave tempWidth = $peakWidth
SetDataFolder root:Packages:MultiPeakFit2:$currentMPfolder
print "c"+num2str(i)
for(j=0;j<numPeaks[i-1]-1;j+=1)
String currentPeak = "Peak "+num2str(j)+" Coefs"
Wave peakInfo1 = $currentPeak
tempEnergy[j] = peakInfo1[0]
tempWidth[j] = peakInfo1[1]
print currentPeak,peakInfo1[0],peakInfo1[1]
endfor
endfor
SetDataFolder root:
SetDataFolder root:OscillatorStrengths:
Variable l
for(l=1;l<=fnum;l+=1)
String currentOSfolder = "c"+num2str(l)
SetDataFolder root:OscillatorStrengths:$currentOSfolder
String currentOSwave = "int_OS_c"+num2str(l)
String currenteVwave = "n_eV_c"+num2str(l)
Wave os_w = $currentOSwave
Wave eV_w = $currenteVwave
Duplicate/O $currentOSwave, root:$currentOSwave
Duplicate/O $currenteVwave, root:$currenteVwave
endfor
SetDataFolder root:
Variable m
for(m=1;m<=fnum;m+=1)
peakEnergy = "pEnergy"+num2str(m)+"_c"+num2str(m)
peakWidth = "pWidth"+num2str(m)+"_c"+num2str(m)
currentOSwave = "OS_c"+num2str(m)
Wave E_current = $peakEnergy
Wave width_current = $peakWidth
Wave os_current = $currentOSwave
endfor
SetDataFolder root:
End
The problem I'm facing comes when I try to organize the data and match the results from MPF with my DFT results using the function below "BinWaves". The purpose of the first loop is to use peak Energy obtained from MPF and match it to the DFT data. I do this by first identifying the point in the DFT wave that matches the energy from the peak energy obtained from MPF using BinarySearchInterp. Then I use the peak width from MPF to get a lower energy bound and an upper energy bound in order to essentially bin the DFT data according to the results from MPF. This is done for every peak for every atom. This part of the procedure kind of works but for some reason the results are not getting saved properly. As a bit of a sanity check, I have it print out the waves on which its operating on and it seems to be doing the right thing(picture below)
Function BinWaves(eWave,osc_strwave,peakEn,peakWid,numPeaks,fnum) //This function sorts data from DFT and MP fitting
String eWave,osc_strwave,peakEn,peakWid
Wave numPeaks
Variable fnum
SetDataFolder root:
Variable i,j,k,l,m,n
for(i=1;i<=fnum;i+=1)
String eWaveName = "n_eV_" + eWave + num2str(i)
String OSnameWave = "OS_" + osc_strwave + num2str(i)
String peakEnergyName ="pEnergy"+num2str(i)+"_" + peakEn + num2str(i)
String peakWidthName ="pWidth"+num2str(i)+"_" + peakWid +num2str(i)
Wave current_eVwave = $eWaveName
Wave current_OSwave = $OSnameWave
Wave peakEnergyWave = $peakEnergyName
Wave peakWidthWave = $peakWidthName
for(j=0;j<=numPeaks[i-1];j+=1)
Variable currentPeakEnergy = peakEnergyWave[j]
Variable currentPeakWidth = peakWidthWave[j]
Variable peakEnergyPoint = round(BinarySearchInterp($eWaveName,currentPeakEnergy))
Variable energyVal = current_eVwave[peakEnergyPoint]
Variable lowEnergy = energyVal - currentPeakWidth
Variable highEnergy = energyVal + currentPeakWidth
Variable lowPeakEpoint = round(BinarySearchInterp(current_eVwave,lowEnergy))
Variable highPeakEpoint = round(BinarySearchInterp(current_eVwave,highEnergy))
endfor
print eWaveName,OSnameWave,peakEnergyName, peakWidthName,currentPeakEnergy,currentPeakWidth
endfor
NewDataFolder/O BinnedWaves
SetDataFolder root:BinnedWaves
String binnedFolder = GetDataFolder(1)
for(l=1;l<=fnum;l+=1) //Make folders to contain binned waves for each atom
String currentAtom = "c"+num2str(l)
NewDataFolder/O $currentAtom
SetDataFolder root:BinnedWaves:$currentAtom
//String currentPeakName = "Peak"+num2str(k)
//NewDataFolder $currentPeakName
SetDataFolder binnedFolder
endfor
for(k=0;k<=fnum;k+=1) //Make folders to contain peak parameters for binned waves
currentAtom = "c"+num2str(k+1)
SetDataFolder root:BinnedWaves:$currentAtom
for(m=0;m<=numPeaks[k]-1;m+=1)
String currentPeakName = "Peak"+num2str(m)
NewDataFolder/O $currentPeakName
endfor
Make/O/N=(numPeaks[k]) $currentPeakName
endfor
SetDataFolder root:
//Need to make wave with set of energies within peak width energy range for each peak
//Need to make wave containing set of oscillator strengths corresponding to energy set for each peak
for(n=1;n<=fnum;n+=1)
currentAtom = "c"+num2str(n)
eWaveName = "eV_" + eWave + num2str(n)
OSnameWave = "OS_" + osc_strwave + num2str(n)
peakEnergyName ="pEnergy"+num2str(n)+"_" + peakEn + num2str(n)
peakWidthName ="pWidth"+num2str(n)+"_" + peakWid +num2str(n)
Wave current_eVwave = $eWaveName
Wave current_OSwave = $OSnameWave
Wave peakEnergyWave = $peakEnergyName
Wave peakWidthWave = $peakWidthName
Duplicate/O $eWaveName, root:BinnedWaves:$(currentAtom):$eWaveName
Duplicate/O $OSnameWave, root:BinnedWaves:$(currentAtom):$OSnameWave
Duplicate/O $peakEnergyName, root:BinnedWaves:$(currentAtom):$peakEnergyName
Duplicate/O $peakWidthName, root:BinnedWaves:$(currentAtom):$peakWidthName
endfor
End
So the procedure is operating on the correct waves but it's not returning the correct parameters for the peak energy and peak width.
The last part of the process involves the function BinEnergies shown below which I want to use to place the actual DFT energies and corresponding intensities( named OS in the procedure) into waves for each peak and for each atom. The output of the function is shown below.
Function binEnergies(fnum)
Variable fnum
Variable i,j,k,m
SetDataFolder root:BinnedWaves
for(i=1;i<=fnum;i+=1)
SetDataFolder root:BinnedWaves:
String currentFolder = "c"+num2str(i)
SetDataFolder root:BinnedWaves:$currentFolder
String currentOSwaveName = "OS_c"+num2str(i)
String currentEnergyWaveName = "eV_c"+num2str(i)
String peakEnergyWaveName = "pEnergy"+num2str(i)+"_c"+num2str(i)
String peakWidthWaveName = "pWidth"+num2str(i)+"_c"+num2str(i)
Wave currentOSwave = $currentOSwaveName
Wave peakEnergyWave = $peakEnergyWaveName
Wave peakWidthWave = $peakWidthWaveName
Wave currentEnergyWave = $currentEnergyWaveName
Wave numPeaks = root:numPeaks
for(j=0;j<=(numPeaks[i-1]-1);j+=1)
String currentPeakFolder = "root:BinnedWaves:"+currentFolder+":Peak"+num2str(j)
SetDataFolder $currentPeakFolder
Variable currentPeakEnergy = peakEnergyWave[j]
Variable currentPeakWidth = abs(peakWidthWave[j])
Variable peakEnergyPoint = round(BinarySearchInterp(currentEnergyWave,currentPeakEnergy))
Variable energyVal = currentEnergyWave[peakEnergyPoint]//Energy at peakEnergyPoint
Variable lowEnergy = energyVal - currentPeakWidth //Lower energy range
Variable highEnergy = energyVal + currentPeakWidth //Upper energy range
Variable lowPeakEpoint = round(BinarySearchInterp(currentEnergyWave,lowEnergy))
Variable highPeakEpoint = round(BinarySearchInterp(currentEnergyWave,highEnergy))
Variable totalPoints = highPeakEPoint-lowPeakEPoint
String binPointsWave ="Peak"+num2str(j)+"_points"
String binEnergyWave ="Peak"+num2str(j)+"_Erange"
String binOSWave = "Peak"+num2str(j)+"_OS"
Make/O/N=(totalPoints) $binPointsWave, $binEnergyWave, $binOSWave
Wave binPoints = $binPointsWave
Wave binEnergy = $binEnergyWave
Wave binOS = $binOSWave
k =0
do
binEnergy[k] = currentEnergyWave[k+lowPeakEpoint]
binOS[k] = currentOSwave[k+lowPeakEpoint]
k+=1
//print k,totalPoints
while(k<totalPoints)
print currentFolder
print currentEnergyWaveName,currentPeakEnergy,currentPeakWidth,peakEnergyPoint
print energyVal,lowEnergy,highEnergy
endfor
//Make sure peak widths are positive!!
SetDataFolder root:
endfor
End
As you can see, the peak energies and peak widths are being correctly identified and returned by the procedure but the waves are not being placed into the appropriate folders.
I think I'm just running in circles now and I'm kind of out of ideas. If you have any suggestions or ideas I'd be more than happy to try them out.
I'm attaching an IGOR file with the data I've been working with in addition to the Procedure containing the functions above. I made a wrapper function that runs all 3 functions but since one of the functions is not quite working it's of limited use.
Thanks!!
That's a lot of information and a lot of code- I haven't followed the whole business. But one recommendation is to use Igor's debugger to examine what's going on. It will allow you to step through all that code and examine what's what while the code runs. If you don't know about the debugger yet, you need to! DisplayHelpTopic "The Debugger".
July 16, 2018 at 04:41 pm - Permalink