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:

    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
    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
    SetDataFolder root:
    Variable i
    Variable j
        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)
            String currentPeak = "Peak "+num2str(j)+" Coefs"
            Wave peakInfo1 = $currentPeak
            tempEnergy[j] = peakInfo1[0]
            tempWidth[j] = peakInfo1[1]
            print currentPeak,peakInfo1[0],peakInfo1[1]
    SetDataFolder root:
    SetDataFolder root:OscillatorStrengths:
    Variable l
        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
    SetDataFolder root:
    Variable m
        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
    SetDataFolder root:


Sample output for binEnergiesPrep


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

        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
            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))     
        print eWaveName,OSnameWave,peakEnergyName, peakWidthName,currentPeakEnergy,currentPeakWidth
    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
    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
            String currentPeakName = "Peak"+num2str(m)
            NewDataFolder/O $currentPeakName
    Make/O/N=(numPeaks[k]) $currentPeakName
    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

        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     

Output of BinWaves

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
        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
            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
                    binEnergy[k] = currentEnergyWave[k+lowPeakEpoint]
                    binOS[k] = currentOSwave[k+lowPeakEpoint]
                    //print k,totalPoints              
        print currentFolder
        print currentEnergyWaveName,currentPeakEnergy,currentPeakWidth,peakEnergyPoint
        print energyVal,lowEnergy,highEnergy       
        //Make sure peak widths are positive!!
        SetDataFolder root:

Output for BinEnergies

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.



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".