bartels-rausch
091cf00e77
* Transitioned nappNorm, nappDetermineShift, and nappCopyShift to use modern MATLAB arguments blocks, removing their legacy functionSignatures.json entries.
* Refactored nappNorm to rely on prior peak fitting results (e.g., from nappFitGaussXps) for normalization rather than performing the fitting internally.
* Decoupled energy calibration: nappDetermineShift now solely calculates internal binding energy shifts, and nappCopyShift distributes them using a unified matching engine ('earlier', 'after', 'spectrumNumber').
* Implemented a Universal Logging Block to track data provenance, automatically appending execution timestamps and caller script names to processed rows.
* Updated nappCombineSpectra to output a logDataSource cell array and cleanly name combined spectra.
* Created nappSaveDataTable for automatic, script-named .mat file saving, complete with a Live Editor interceptor.
* Added a publish-mode silencer to nappFitGaussXps to prevent optimizer spam in generated PDFs.
* Updated readme.md to document the new arguments structure, the 3-step calibration workflow, the matching engine, and the new table variables.
399 lines
19 KiB
Matlab
399 lines
19 KiB
Matlab
function DataTable = nappBackground(DataTable,varargin)
|
|
%NAPPBACKGROUND Removes background from photoelectron or XAS spectra.
|
|
% This function determines and removes the background from PE or XAS
|
|
% spectra using mean subtraction or curve fitting (linear or exponential).
|
|
% Background correction can be performed on kinetic energy, binding energy,
|
|
% or X-ray photon energy scales. Results are optionally normalized by
|
|
% dividing instead of subtracting.
|
|
%
|
|
% SYNTAX
|
|
% ------
|
|
% DataTable = NAPPBACKGROUND(DataTable, energyStart, energyEnd)
|
|
%
|
|
% DataTable = NAPPBACKGROUND(DataTable, energyStart, energyEnd, 'Name', Value)
|
|
%
|
|
% EXAMPLES
|
|
% --------
|
|
% % Simple mean subtraction with kinetic energy scale
|
|
% DataTable = nappBackground(DataTable, 368, 370);
|
|
%
|
|
% % Linear fit on selected rows with binding energy scale
|
|
% DataTable = nappBackground(DataTable, 368, 370, ...
|
|
% 'tableSelection', tmpDataSelection, ...
|
|
% 'energyScale', 'binding', ...
|
|
% 'mode', 'linfit', ...
|
|
% 'showFigures', 'off');
|
|
%
|
|
% INPUT ARGUMENTS
|
|
% ---------------
|
|
% DataTable table
|
|
% Data table containing spectra and energy scales. Must contain
|
|
% one of the following variables: 'spectrum_counts',
|
|
% 'spectrum_countsPerSecond', or 'spectrum_countsNorm'.
|
|
%
|
|
% energyStart — Start energies for excluded ranges
|
|
% Numeric scalar or vector.
|
|
% Defines the lower bound(s) of the data range(s) excluded from
|
|
% background fitting.
|
|
%
|
|
% energyEnd — End energies for excluded ranges
|
|
% Numeric scalar or vector.
|
|
% Defines the upper bound(s) of the data range(s) excluded from
|
|
% background fitting.
|
|
%
|
|
% If energyStart and energyEnd are scalars, the range between
|
|
% energyStart and energyEnd is excluded from the fit (for example,
|
|
% to exclude the main peak).
|
|
%
|
|
% If energyStart is a two-element vector, the range from the first
|
|
% data point to the second value of energyStart is additionally
|
|
% excluded.
|
|
% If energyEnd is a two-element vector, the range from the second
|
|
% value of energyEnd to the end of the dataset is additionally
|
|
% excluded.
|
|
%
|
|
% Specify vector values as comma-separated entries in brackets,
|
|
% for example [395, 407].
|
|
%
|
|
%
|
|
% NAME-VALUE PAIR ARGUMENTS (optional)
|
|
% ------------------------------------
|
|
% 'tableSelection' logical array (default: all rows)
|
|
% Logical index to filter rows in DataTable.
|
|
% Example: tableSelection = DataTable.quality > 0.8
|
|
%
|
|
% 'tableVariable' char vector (default: 'spectrum_countsPerSecond')
|
|
% Spectral variable to process. Must be one of:
|
|
% - 'spectrum_counts'
|
|
% - 'spectrum_countsPerSecond'
|
|
% - 'spectrum_countsNorm'
|
|
%
|
|
% 'mode' char vector (default: 'mean')
|
|
% Background subtraction method:
|
|
% - 'mean': subtract mean intensity in energy range
|
|
% - 'linfit': subtract linear regression (poly1 fit)
|
|
% - 'exponential': subtract exponential fit (exp2)
|
|
% - 'mean, divide': divide by mean (normalization)
|
|
% - 'linfit, divide': divide by linear fit
|
|
% - 'exponential, divide': divide by exponential fit
|
|
%
|
|
% 'energyScale' char vector (default: 'kinetic')
|
|
% Energy axis to use for background calculation:
|
|
% - 'kinetic': uses scientaEkin_eV column
|
|
% - 'binding': uses bindingEnergy_eV column
|
|
% - 'xRay': uses xRayEkin_eV column (for NEXAFS data)
|
|
%
|
|
% 'showFigures' char vector (default: 'on')
|
|
% Display fit visualization:
|
|
% - 'on': show fit plot with data and excluded points
|
|
% - 'off': suppress figure output
|
|
%
|
|
% OUTPUT ARGUMENTS
|
|
% ----------------
|
|
% DataTable table
|
|
% Updated data table with background-corrected spectra. The table
|
|
% variable specified by 'tableVariable' is overwritten with
|
|
% corrected data. If mode contains 'divide', results are written
|
|
% to 'spectrum_countsNorm' instead.
|
|
%
|
|
% Three log columns are added/updated:
|
|
% - logBackgroundMode: fitting method used
|
|
% - logBackgroundRange: [energyStart, energyEnd] for each row
|
|
% - logBackgroundDate: timestamp of processing
|
|
%
|
|
% NOTES
|
|
% -----
|
|
% - Empty or missing spectra in tableVariable will trigger an error.
|
|
%
|
|
% SEE ALSO
|
|
% --------
|
|
% inputParser, fit, fittype, prepareCurveData, excludedata
|
|
|
|
|
|
|
|
% Define conditions and defaults of input
|
|
defaultTableVariable = 'spectrum_countsPerSecond';
|
|
expectedTableVariable = {'spectrum_counts','spectrum_countsPerSecond','spectrum_countsNorm'};
|
|
defaultMode = 'mean';
|
|
expectedMode = {'mean','linfit','exponential','mean, divide','linfit, divide','exponential, divide'};
|
|
defaultTableSelection = true(size(DataTable,1),1);
|
|
defaultShowFigures = 'on';
|
|
expectedShowFigures = {'on','off'};
|
|
defaultEnergyScale = 'kinetic';
|
|
expectedEnergyScale = {'kinetic','binding','xRay'};
|
|
|
|
% Parse input
|
|
p = inputParser;
|
|
addRequired(p,'DataTable',@istable);
|
|
addRequired(p,'energyStart',@isnumeric)
|
|
addRequired(p,'energyEnd',@isnumeric)
|
|
addOptional(p,'tableSelection',defaultTableSelection,@islogical);
|
|
addParameter(p,'tableVariable',defaultTableVariable,...
|
|
@(x) any(validatestring(x,expectedTableVariable)));
|
|
addParameter(p,'mode',defaultMode,...
|
|
@(x) any(validatestring(x,expectedMode)));
|
|
addParameter(p,'showFigures',defaultShowFigures,...
|
|
@(x) any(validatestring(x,expectedShowFigures)));
|
|
addParameter(p,'energyScale',defaultEnergyScale,...
|
|
@(x) any(validatestring(x,expectedEnergyScale)));
|
|
parse(p,DataTable,varargin{:});
|
|
|
|
% check input (all elements of input variable need to have non-zero size)
|
|
if ~all(cell2mat(cellfun(@size,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable),'UniformOutput',false)),'all')
|
|
disp('The table variable appears to be empty. I quit now')
|
|
error('The table variable appears to be empty.')
|
|
end
|
|
|
|
% add empty table variables (if they don't exist)
|
|
existColumn = strcmp('lognappBackgroundRange',DataTable.Properties.VariableNames);
|
|
if ~any(existColumn)
|
|
DataTable.lognappBackgroundRange = NaN(size(DataTable,1),2);
|
|
end
|
|
existColumn = strcmp('lognappBackgroundMode',DataTable.Properties.VariableNames);
|
|
if ~any(existColumn)
|
|
DataTable.lognappBackgroundMode = cellstr(repmat(' ',size(p.Results.DataTable,1),1));
|
|
end
|
|
% add empty variable to table
|
|
existColumn = strcmp('spectrum_countsNorm',DataTable.Properties.VariableNames);
|
|
if ~any(existColumn)
|
|
DataTable.spectrum_countsNorm = cell(size(DataTable,1),1);
|
|
end
|
|
|
|
|
|
% repeat this for all input.
|
|
for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
|
|
|
|
if startsWith(p.Results.mode,'linfit')
|
|
ft = fittype( 'poly1' );
|
|
elseif contains(p.Results.mode,'exponential')
|
|
ft = fittype( 'exp2' );
|
|
end
|
|
|
|
if contains(p.Results.energyScale,'kinetic')
|
|
|
|
if ~startsWith(p.Results.mode,'mean')
|
|
[xData, yData] = prepareCurveData(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii},p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
|
|
excludedPoints = excludedata( xData, yData, 'domain',[p.Results.energyStart(1) p.Results.energyEnd(1)] );
|
|
|
|
if size(p.Results.energyStart,2) ==2 && size(p.Results.energyEnd,2) ==1
|
|
excludedPointsEins = excludedPoints;
|
|
clearvars excludedPoints
|
|
excludedPointsZwei = excludedata( xData, yData, 'domain',[xData(1) p.Results.energyStart(2)] );
|
|
excludedPoints = excludedPointsEins & excludedPointsZwei;
|
|
elseif size(p.Results.energyStart,2) ==2 && size(p.Results.energyEnd,2) ==2
|
|
excludedPointsEins = excludedPoints;
|
|
clearvars excludedPoints
|
|
excludedPointsZwei = excludedata( xData, yData, 'domain',[xData(1) p.Results.energyStart(2)] );
|
|
excludedPointsDrei = excludedata( xData, yData, 'domain',[p.Results.energyEnd(2) xData(end)] );
|
|
excludedPoints = excludedPointsEins & excludedPointsZwei & excludedPointsDrei;
|
|
elseif size(p.Results.energyStart,2) ==1 && size(p.Results.energyEnd,2) ==2
|
|
excludedPointsEins = excludedPoints;
|
|
clearvars excludedPoints
|
|
excludedPointsZwei = excludedata( xData, yData, 'domain',[p.Results.energyEnd(2) xData(end)] );
|
|
excludedPoints = excludedPointsEins & excludedPointsZwei;
|
|
|
|
end
|
|
|
|
opts = fitoptions( 'Method', 'NonlinearLeastSquares','Exclude', ~excludedPoints );
|
|
opts.Display = 'Off';
|
|
|
|
% Fit model to data.
|
|
[fitresult, ~] = fit( xData, yData, ft, opts );
|
|
if contains(p.Results.showFigures,'on')
|
|
% Plot fit with data.
|
|
figure( 'Name', 'untitled fit 1' );
|
|
plot( fitresult, xData, yData, ~excludedPoints );
|
|
% legend( h, 'y vs. x', 'Excluded y vs. x', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
|
|
% Label axes
|
|
xlabel( 'kinetic energy [eV]', 'Interpreter', 'none' );
|
|
ylabel( 'PE intensity', 'Interpreter', 'none' );
|
|
grid on
|
|
title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
|
|
end
|
|
if endsWith(p.Results.mode,'divide')
|
|
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}./fitresult(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii});
|
|
|
|
else
|
|
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii});
|
|
end
|
|
elseif startsWith(p.Results.mode,'mean')
|
|
b_1 = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii} >= p.Results.energyStart,1);
|
|
b_2 = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii} >= p.Results.energyEnd,1);
|
|
|
|
if b_2 == 1
|
|
b_2=10;
|
|
end
|
|
|
|
b_mean = mean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(b_1:b_2),'omitnan');
|
|
tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false); %#ok<*AGROW>
|
|
|
|
|
|
end
|
|
|
|
elseif contains(p.Results.energyScale,'binding')
|
|
if ~startsWith(p.Results.mode,'mean')
|
|
[xData, yData] = prepareCurveData(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii},p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
|
|
|
|
% Set up fittype and options.
|
|
|
|
excludedPoints = ~excludedata( xData, yData, 'domain',[p.Results.energyStart(1) p.Results.energyEnd(1)] );
|
|
if size(p.Results.energyStart,2) ==2 && size(p.Results.energyEnd,2) ==1
|
|
excludedPointsEins = excludedPoints;
|
|
clearvars excludedPoints
|
|
excludedPointsZwei = ~excludedata( xData, yData, 'domain',[xData(end) p.Results.energyStart(2)] );
|
|
excludedPoints = excludedPointsEins | excludedPointsZwei;
|
|
elseif size(p.Results.energyStart,2) ==2 && size(p.Results.energyEnd,2) ==2
|
|
excludedPointsEins = excludedPoints;
|
|
clearvars excludedPoints
|
|
excludedPointsZwei = ~excludedata( xData, yData, 'domain',[xData(end) p.Results.energyStart(2)] );
|
|
excludedPointsDrei = ~excludedata( xData, yData, 'domain',[p.Results.energyEnd(2) xData(1)] );
|
|
excludedPoints = excludedPointsEins | excludedPointsZwei |excludedPointsDrei ;
|
|
elseif size(p.Results.energyStart,2) ==1 && size(p.Results.energyEnd,2) ==2
|
|
excludedPointsEins = excludedPoints;
|
|
clearvars excludedPoints
|
|
excludedPointsZwei = ~excludedata( xData, yData, 'domain',[p.Results.energyEnd(2) xData(1)] );
|
|
excludedPoints = excludedPointsEins | excludedPointsZwei;
|
|
end
|
|
|
|
opts = fitoptions( 'Method', 'NonlinearLeastSquares','Exclude', excludedPoints );
|
|
opts.Display = 'Off';
|
|
|
|
|
|
|
|
% Fit model to data.
|
|
[fitresult, ~] = fit( xData, yData, ft, opts );
|
|
if contains(p.Results.showFigures,'on')
|
|
% Plot fit with data.
|
|
figure( 'Name', 'untitled fit 1' );
|
|
plot( fitresult, xData, yData, excludedPoints );
|
|
% legend( h, 'y vs. x', 'Excluded y vs. x', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
|
|
% Label axes
|
|
xlabel( 'binding energy [eV]', 'Interpreter', 'none' );
|
|
ylabel( 'PE intensity', 'Interpreter', 'none' );
|
|
grid on
|
|
title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
|
|
ax = gca;
|
|
ax.XDir = 'reverse';
|
|
end
|
|
if endsWith(p.Results.mode,'divide')
|
|
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}./fitresult(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii});
|
|
|
|
else
|
|
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii});
|
|
end
|
|
elseif startsWith(p.Results.mode,'mean')
|
|
b_1 = find(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii} <= p.Results.energyStart,1);
|
|
b_2 = find(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii} <= p.Results.energyEnd,1);
|
|
|
|
if b_2 == 1
|
|
b_2=10;
|
|
end
|
|
|
|
b_mean = mean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(b_2:b_1),'omitnan');
|
|
tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false);
|
|
|
|
|
|
end
|
|
elseif contains(p.Results.energyScale,'xRay')
|
|
if contains(p.Results.mode,'linfit') % LINEAR FIT
|
|
|
|
% Fit: 'untitled fit 1'.
|
|
[xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
|
|
p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
|
|
% Set up fittype and options.
|
|
ft = fittype( 'poly1' );
|
|
excludedPoints = excludedata( xData, yData, 'domain',[p.Results.energyStart p.Results.energyEnd] );
|
|
fo = fitoptions('Exclude', excludedPoints);
|
|
% Fit model to data.
|
|
[fitresult, ~] = fit( xData, yData, ft, fo );
|
|
if contains(p.Results.showFigures,'on')
|
|
% Plot fit with data.
|
|
figure( 'Name', 'fit pre-edge' );
|
|
plot( fitresult, xData, yData, excludedPoints );
|
|
xlabel('Photon Energy eV')
|
|
ylabel('partial A-M.E.Y.')
|
|
grid on
|
|
title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
|
|
end
|
|
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii});
|
|
|
|
elseif contains(p.Results.mode,'mean') % mean
|
|
|
|
if p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(1) < p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(end)
|
|
% prepare Data
|
|
|
|
[xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
|
|
p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
|
|
% MEAN
|
|
b_1 = find(xData >= p.Results.energyStart,1);
|
|
b_2 = find(xData >= p.Results.energyEnd,1);
|
|
|
|
b_mean = mean(yData(b_1:b_2));
|
|
elseif p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(1) > p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(end)
|
|
% prepare Data
|
|
[xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
|
|
p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
|
|
% MEAN
|
|
b_1 = find(xData <= p.Results.energyStart,1);
|
|
b_2 = find(xData <= p.Results.energyEnd,1);
|
|
|
|
b_mean = mean(yData(b_2:b_1));
|
|
else
|
|
disp('Warning: Issue with NEXAFS data or energies')
|
|
end
|
|
|
|
|
|
tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false);
|
|
end
|
|
|
|
|
|
end
|
|
end
|
|
|
|
% fill selected rows
|
|
if endsWith(p.Results.mode,'divide')
|
|
DataTable(p.Results.tableSelection,:).spectrum_countsNorm = tmp;
|
|
else
|
|
DataTable(p.Results.tableSelection,:).(p.Results.tableVariable) = tmp;
|
|
end
|
|
% write log to table
|
|
DataTable.lognappBackgroundMode(p.Results.tableSelection,:) = cellstr(repmat(p.Results.mode,size(p.Results.DataTable(p.Results.tableSelection,:),1),1));
|
|
DataTable.lognappBackgroundRange(p.Results.tableSelection,:) = repmat([p.Results.energyStart(1) p.Results.energyEnd(1)],size(p.Results.DataTable(p.Results.tableSelection,:),1),1);
|
|
|
|
|
|
|
|
%--- Define variables for the Universal Logging Block ---
|
|
numRowsTable = height(DataTable);
|
|
processIdx = find(p.Results.tableSelection);
|
|
% -------------------------------------------------------------------
|
|
|
|
% 1. Dynamically get this function's name
|
|
thisFuncName = mfilename;
|
|
colDate = sprintf('log%sDate', thisFuncName);
|
|
colCaller = sprintf('log%sCaller', thisFuncName);
|
|
|
|
|
|
% 2. Probe for column existence, create empty string arrays if missing
|
|
if ~ismember(colDate, DataTable.Properties.VariableNames)
|
|
DataTable.(colDate) = repmat("", numRowsTable, 1);
|
|
end
|
|
if ~ismember(colCaller, DataTable.Properties.VariableNames)
|
|
DataTable.(colCaller) = repmat("", numRowsTable, 1);
|
|
end
|
|
|
|
% 3. Get current time as a string
|
|
currentTimeStr = string(datetime('now', 'Format', 'yyyy-MM-dd HH:mm:ss'));
|
|
|
|
% 4. Determine the calling script using the call stack
|
|
stackInfo = dbstack;
|
|
if length(stackInfo) > 1
|
|
callerScript = string(stackInfo(2).file); % The script that called this function
|
|
else
|
|
callerScript = " "; % Called directly from the command window
|
|
end
|
|
|
|
% 5. Write the log data only to the rows being processed
|
|
DataTable.(colDate)(processIdx) = currentTimeStr;
|
|
DataTable.(colCaller)(processIdx) = callerScript;
|
|
% -------------------------------------------------------------------------
|
|
end |