Several domains within a spectrum can now be excluded in background fitting

This commit is contained in:
bartels-rausch committed 2026-05-11 11:45:46 +02:00
1 parent 55f972f6ae
commit 2f942c554c
1 file changed
+91 -65
+91 -65
View File
@@ -32,28 +32,42 @@ function DataTable = nappBackground(DataTable,varargin)
% '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.
% Numeric vector of length N (must match length of energyEnd).
% Each element pairs with the corresponding element of energyEnd
% to define an excluded domain:
%
% Position 1 (required):
% Lower bound of the main excluded domain (e.g. the peak).
% Combined with energyEnd(1) as [energyStart(1), energyEnd(1)].
%
% Position 2 (optional, backward compatible):
% Cuts off the near-edge of the spectrum. Excludes from the
% first data point up to energyStart(2) (kinetic), or from
% energyStart(2) to the last data point (binding).
%
% Positions 3+ (optional, new):
% Each pair [energyStart(k), energyEnd(k)] defines an
% additional arbitrary excluded domain.
%
% energyEnd — End energies for excluded ranges
% Numeric scalar or vector.
% Defines the upper bound(s) of the data range(s) excluded from
% background fitting.
% Numeric vector of length N (must match length of energyStart).
%
% If energyStart and energyEnd are scalars, the range between
% energyStart and energyEnd is excluded from the fit (for example,
% to exclude the main peak).
% Position 1 (required):
% Upper bound of the main excluded domain.
%
% 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.
% Position 2 (optional, backward compatible):
% Cuts off the far-edge of the spectrum. Excludes from
% energyEnd(2) to the last data point (kinetic), or from
% the first data point to energyEnd(2) (binding).
%
% Specify vector values as comma-separated entries in brackets,
% for example [395, 407].
% Positions 3+ (optional, new):
% Upper bound of additional excluded domains (pairs with
% energyStart(3), energyStart(4), ...).
%
% Example — main peak + two additional exclusions:
% energyStart = [395, 380, 410]
% energyEnd = [407, NaN, 415]
% (position 2 energyEnd=NaN means no far-edge cutoff)
%
%
% NAME-VALUE PAIR ARGUMENTS (optional)
@@ -138,6 +152,11 @@ addParameter(p,'energyScale',defaultEnergyScale,...
@(x) any(validatestring(x,expectedEnergyScale)));
parse(p,DataTable,varargin{:});
% Pad shorter exclusion vector with NaN so both always have equal length.
% This ensures backward compatibility when only one vector has 2 elements.
energyStart = [p.Results.energyStart, NaN(1, max(0, numel(p.Results.energyEnd) - numel(p.Results.energyStart)))];
energyEnd = [p.Results.energyEnd, NaN(1, max(0, numel(p.Results.energyStart) - numel(p.Results.energyEnd)))];
% 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')
@@ -173,26 +192,29 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
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
% Position 1: main exclusion domain (e.g. the peak) — backward compatible
excludedPoints = excludedata(xData, yData, 'domain', [energyStart(1), energyEnd(1)]);
% Convention: true = point is inside fitting region (kept). Combined with & below.
% Position 2: edge cutoffs — backward compatible
if numel(energyStart) >= 2 && ~isnan(energyStart(2))
excludedPoints = excludedPoints & ...
excludedata(xData, yData, 'domain', [xData(1), energyStart(2)]);
end
if numel(energyEnd) >= 2 && ~isnan(energyEnd(2))
excludedPoints = excludedPoints & ...
excludedata(xData, yData, 'domain', [energyEnd(2), xData(end)]);
end
% Positions 3+: additional arbitrary exclusion domains (new)
nPairs = min(numel(energyStart), numel(energyEnd));
for k = 3:nPairs
if ~isnan(energyStart(k)) && ~isnan(energyEnd(k))
excludedPoints = excludedPoints & ...
excludedata(xData, yData, 'domain', [energyStart(k), energyEnd(k)]);
end
end
opts = fitoptions( 'Method', 'NonlinearLeastSquares','Exclude', ~excludedPoints );
opts.Display = 'Off';
@@ -217,8 +239,8 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
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);
b_1 = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii} >= energyStart,1);
b_2 = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii} >= energyEnd,1);
if b_2 == 1
b_2=10;
@@ -236,24 +258,28 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
% 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
% Position 1: main exclusion domain (e.g. the peak) — backward compatible
excludedPoints = ~excludedata(xData, yData, 'domain', [energyStart(1), energyEnd(1)]);
% Convention: true = point should be excluded from fit. Combined with | below.
% Position 2: edge cutoffs — backward compatible
if numel(energyStart) >= 2 && ~isnan(energyStart(2))
excludedPoints = excludedPoints | ...
~excludedata(xData, yData, 'domain', [xData(end), energyStart(2)]);
end
if numel(energyEnd) >= 2 && ~isnan(energyEnd(2))
excludedPoints = excludedPoints | ...
~excludedata(xData, yData, 'domain', [energyEnd(2), xData(1)]);
end
% Positions 3+: additional arbitrary exclusion domains (new)
nPairs = min(numel(energyStart), numel(energyEnd));
for k = 3:nPairs
if ~isnan(energyStart(k)) && ~isnan(energyEnd(k))
excludedPoints = excludedPoints | ...
~excludedata(xData, yData, 'domain', [energyStart(k), energyEnd(k)]);
end
end
opts = fitoptions( 'Method', 'NonlinearLeastSquares','Exclude', excludedPoints );
opts.Display = 'Off';
@@ -282,8 +308,8 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
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);
b_1 = find(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii} <= energyStart,1);
b_2 = find(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii} <= energyEnd,1);
if b_2 == 1
b_2=10;
@@ -302,7 +328,7 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
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] );
excludedPoints = excludedata( xData, yData, 'domain',[energyStart energyEnd] );
fo = fitoptions('Exclude', excludedPoints);
% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, fo );
@@ -325,8 +351,8 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
[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_1 = find(xData >= energyStart,1);
b_2 = find(xData >= 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)
@@ -334,8 +360,8 @@ for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
[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_1 = find(xData <= energyStart,1);
b_2 = find(xData <= energyEnd,1);
b_mean = mean(yData(b_2:b_1));
else
@@ -358,7 +384,7 @@ else
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);
DataTable.lognappBackgroundRange(p.Results.tableSelection,:) = repmat([energyStart(1) energyEnd(1)],size(p.Results.DataTable(p.Results.tableSelection,:),1),1);