diff --git a/beamtime/nappBackground.m b/beamtime/nappBackground.m index 11fc1f5..5caeb2e 100644 --- a/beamtime/nappBackground.m +++ b/beamtime/nappBackground.m @@ -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);