From f6c9581fb040786a82d6e4133efae4fcc61d573c Mon Sep 17 00:00:00 2001 From: Marc Hahn Date: Wed, 20 May 2020 19:37:58 +0200 Subject: [PATCH 1/8] Add variable region for Shirley background substraction for XPS analysis. --- fityk/transform.cpp | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/fityk/transform.cpp b/fityk/transform.cpp index 0e675713..3c68764e 100644 --- a/fityk/transform.cpp +++ b/fityk/transform.cpp @@ -72,8 +72,30 @@ void shirley_bg(vector &pp) const int max_iter = 50; const double max_rdiff = 1e-6; const int n = pp.size(); - double ya = pp[0].y; // lowest bg - double yb = pp[n-1].y; // highest bg + + /* Determine the range for the shirley background from the active range in the dataset + * the first/last point is determined by the first/last active datapoint */ + unsigned long startXindex = 0; + unsigned long endXindex =pp.size()-1; + + for (unsigned long i = 0; i < pp.size(); i++) { + if( pp[i].is_active==true){ + startXindex = i; + i = pp.size(); + } + } + for (unsigned long i = pp.size(); i == 0; --i) { + if( pp[i].is_active==true){ + endXindex = i; + i = 0; + } + } + + //double ya = pp[0].y; // lowest bg + //double yb = pp[n-1].y; // highest bg + double ya = pp[startXindex].y; // lowest bg + double yb = pp[endXindex].y; // highest bg + double dy = yb - ya; vector B(n, ya); vector PA(n, 0.); @@ -95,6 +117,8 @@ void shirley_bg(vector &pp) pp[i].y = B[i]; } + + } // anonymous namespace namespace fityk { From 9e45ca9370f85085c2abab8b864911605ab49947 Mon Sep 17 00:00:00 2001 From: Marc Benjamin Hahn Date: Fri, 22 May 2020 21:14:32 +0200 Subject: [PATCH 2/8] Reimplement Shirley background substraction and extend it to enable selection of an active region. --- fityk/transform.cpp | 96 ++++++++++++++++++++++++++++----------------- 1 file changed, 60 insertions(+), 36 deletions(-) diff --git a/fityk/transform.cpp b/fityk/transform.cpp index 3c68764e..a9178bf1 100644 --- a/fityk/transform.cpp +++ b/fityk/transform.cpp @@ -67,54 +67,78 @@ void merge_same_x(vector &pp, bool avg) } } + void shirley_bg(vector &pp) { - const int max_iter = 50; - const double max_rdiff = 1e-6; - const int n = pp.size(); - - /* Determine the range for the shirley background from the active range in the dataset - * the first/last point is determined by the first/last active datapoint */ - unsigned long startXindex = 0; - unsigned long endXindex =pp.size()-1; - - for (unsigned long i = 0; i < pp.size(); i++) { + /* Set criteria to stop itertions */ + const int max_iter = 50; + const double max_rdiff = 1e-6; // impement: check for differences + const unsigned long n = pp.size(); + + /* The user has to set an appropriate region around the peak since the outcome is sensitive to that */ + /* Determine the range for the shirley background calculation from the active range in the dataset + * the shirley BG will be calculated between the first and last active datapoint */ + unsigned long startXindex = 0; + unsigned long endXindex = n-1; + + for (unsigned long i = 0; i < n; i++) { if( pp[i].is_active==true){ startXindex = i; - i = pp.size(); + break; } } - for (unsigned long i = pp.size(); i == 0; --i) { + + for (unsigned long i = endXindex; i > 0; i--) { if( pp[i].is_active==true){ endXindex = i; - i = 0; + break; } } - //double ya = pp[0].y; // lowest bg - //double yb = pp[n-1].y; // highest bg - double ya = pp[startXindex].y; // lowest bg - double yb = pp[endXindex].y; // highest bg - - double dy = yb - ya; - vector B(n, ya); - vector PA(n, 0.); - double old_A = 0; - for (int iter = 0; iter < max_iter; ++iter) { - vector Y(n); - for (int i = 0; i < n; ++i) - Y[i] = pp[i].y - B[i]; - for (int i = 1; i < n; ++i) - PA[i] = PA[i-1] + (Y[i] + Y[i-1]) / 2 * (pp[i].x - pp[i-1].x); - double rel_diff = old_A != 0. ? fabs(PA[n-1] - old_A) / old_A : 1.; - if (rel_diff < max_rdiff) - break; - old_A = PA[n-1]; - for (int i = 0; i < n; ++i) - B[i] = ya + dy / PA[n-1] * PA[i]; + /* Here we assume that we are in BE representation and Alow belongs to smaller BE and Ahigh to higher BE */ + double ymin = pp[startXindex].y; + double ymax = pp[endXindex].y; + + /* Determine minimum element from active data. + * In theory this should be ymin. + * but this might not be the case for a bad SNR */ + double yminElement = ymin; + for (unsigned long i = 0; i < n; i++) { + if( pp[i].y < yminElement){ + yminElement = pp[i].y; + } + } + /* create starting background */ + std::vector BG(n, yminElement); + + /* recursive determination of the shirley background */ + for (int iter = 0; iter < max_iter; iter++) { + + /* Determine all background BG(Energy) values in dependence of the Energy */ + for (unsigned long EnergyIdx = 0; EnergyIdx < n ; EnergyIdx++){ + + double Alow = 0.0; + double Ahigh = 0.0; + + /* calculate Areas A1 and A2 for each energy value */ + for (unsigned long x = startXindex+1; x < EnergyIdx ; x++){ + Alow = Alow + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]; + } + + for (unsigned long x = EnergyIdx; x < endXindex ; x++){ + Ahigh = Ahigh + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]; + } + + /* The formula is implemented as given in the CasaXPS manual (2006) "Peak fitting in XPS" chapter */ + BG[EnergyIdx] = ymin + (ymax-ymin) * (Alow / (Ahigh+Alow)); + + } + } + + /* copy the resulting background */ + for (unsigned long i = 0; i < n; i++){ + pp[i].y = BG[i]; } - for (int i = 0; i < n; ++i) - pp[i].y = B[i]; } From cd331d41086bdf34574db52c0997897cd1ef36b0 Mon Sep 17 00:00:00 2001 From: MarcBHahn <60821053+MarcBHahn@users.noreply.github.com> Date: Tue, 26 May 2020 15:36:56 +0200 Subject: [PATCH 3/8] Improved ShirleyBG substraction Improves Shirley BG substraction by calculation of ymin and ymax parameter from the average of up to five surrounding datapoints. This improves the Shirley BG estimation for XPS data with bad signal to noise ratio. --- fityk/transform.cpp | 37 +++++++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 6 deletions(-) diff --git a/fityk/transform.cpp b/fityk/transform.cpp index a9178bf1..b50e7cc4 100644 --- a/fityk/transform.cpp +++ b/fityk/transform.cpp @@ -95,13 +95,38 @@ void shirley_bg(vector &pp) } } - /* Here we assume that we are in BE representation and Alow belongs to smaller BE and Ahigh to higher BE */ - double ymin = pp[startXindex].y; - double ymax = pp[endXindex].y; + /* To set ymin and ymax, we assume that we are in BE representation and Alow belongs to smaller BE and Ahigh to higher BE */ - /* Determine minimum element from active data. - * In theory this should be ymin. - * but this might not be the case for a bad SNR */ + /* Straight forward determination of ymin and ymax, but this can result to poor results when SNR is high. */ + //double ymin = pp[startXindex].y; + //double ymax = pp[endXindex].y; + + /* To improve behaviour for bad SNR we average up to 5 y values around the x-index*/ + double ymin = 0.00001; + double ymax = 0.00001; + + int count_start = 0; + int count_end = 0; + + for(unsigned long i=0; i<5; i++){ + /* check if argument is within the range of the array indices */ + if( (int(startXindex+i)-2 >= 0) && ( startXindex+i-2 < n) ){ + ymin = ymin + pp[startXindex+i-2].y; + count_start++; + } + + /* check if argument is within the range of the array indices */ + if( (int(endXindex+i)-2 >= 0) && ( endXindex+i-2 < n) ){ + ymax = ymax + pp[endXindex+i-2].y; + count_end++; + } + } + /* divide by number of summations */ + ymin=ymin/double(count_start); + ymax=ymax/double(count_end); + + + /* Determine minimum element from data to obtain a reasonable starting value for the fit */ double yminElement = ymin; for (unsigned long i = 0; i < n; i++) { if( pp[i].y < yminElement){ From 2f14c76f7243abfa82eada45f9ba6ba16508ae92 Mon Sep 17 00:00:00 2001 From: MarcBHahn <60821053+MarcBHahn@users.noreply.github.com> Date: Mon, 18 Oct 2021 19:54:44 +0200 Subject: [PATCH 4/8] Update transform.cpp Correct comments. --- fityk/transform.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/fityk/transform.cpp b/fityk/transform.cpp index b50e7cc4..0e5b72ee 100644 --- a/fityk/transform.cpp +++ b/fityk/transform.cpp @@ -72,7 +72,7 @@ void shirley_bg(vector &pp) { /* Set criteria to stop itertions */ const int max_iter = 50; - const double max_rdiff = 1e-6; // impement: check for differences + const double max_rdiff = 1e-6; // implement: check for differences const unsigned long n = pp.size(); /* The user has to set an appropriate region around the peak since the outcome is sensitive to that */ @@ -97,11 +97,11 @@ void shirley_bg(vector &pp) /* To set ymin and ymax, we assume that we are in BE representation and Alow belongs to smaller BE and Ahigh to higher BE */ - /* Straight forward determination of ymin and ymax, but this can result to poor results when SNR is high. */ + /* Straight forward determination of ymin and ymax, but this can result in poor results when the SNR is high. */ //double ymin = pp[startXindex].y; //double ymax = pp[endXindex].y; - /* To improve behaviour for bad SNR we average up to 5 y values around the x-index*/ + /* To improve the behaviour for bad SNR we average up to 5 y values around the x-index*/ double ymin = 0.00001; double ymax = 0.00001; @@ -154,7 +154,7 @@ void shirley_bg(vector &pp) Ahigh = Ahigh + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]; } - /* The formula is implemented as given in the CasaXPS manual (2006) "Peak fitting in XPS" chapter */ + /* The formula is implemented as given in the CasaXPS manual (2006) "Peak fitting in XPS" chapter */ BG[EnergyIdx] = ymin + (ymax-ymin) * (Alow / (Ahigh+Alow)); } From c67e7953f6cda690486c959a7e13538393dfd622 Mon Sep 17 00:00:00 2001 From: MarcBHahn <60821053+MarcBHahn@users.noreply.github.com> Date: Mon, 18 Oct 2021 20:28:56 +0200 Subject: [PATCH 5/8] Update transform.cpp Adjust indices of Shirley background removal procedure for the analysis of XPS spectra. --- fityk/transform.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fityk/transform.cpp b/fityk/transform.cpp index 0e5b72ee..799a8384 100644 --- a/fityk/transform.cpp +++ b/fityk/transform.cpp @@ -146,11 +146,11 @@ void shirley_bg(vector &pp) double Ahigh = 0.0; /* calculate Areas A1 and A2 for each energy value */ - for (unsigned long x = startXindex+1; x < EnergyIdx ; x++){ + for (unsigned long x = startXindex+1; x <= EnergyIdx ; x++){ Alow = Alow + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]; } - for (unsigned long x = EnergyIdx; x < endXindex ; x++){ + for (unsigned long x = EnergyIdx+1; x < endXindex ; x++){ Ahigh = Ahigh + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]; } From e5e75ab5fdc6d6e136a17329bacc10605a82bcc9 Mon Sep 17 00:00:00 2001 From: MarcBHahn <60821053+MarcBHahn@users.noreply.github.com> Date: Mon, 18 Oct 2021 20:57:15 +0200 Subject: [PATCH 6/8] Prepare for merge Add the changes from main branch to prepare a new version fpr the pull request --- fityk/transform.cpp | 46 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/fityk/transform.cpp b/fityk/transform.cpp index 799a8384..6f9aa6e9 100644 --- a/fityk/transform.cpp +++ b/fityk/transform.cpp @@ -166,6 +166,52 @@ void shirley_bg(vector &pp) } } +// Calculates the SNIP background iteratively in the +// given slice of the vector of points. +// The active points are modified in-place, inactive +// points are left as they are, assuming that they +// are already considered background. +vector::iterator snip_bg_slice(vector::iterator begin, + vector::iterator end, + int window_width, + int direction, + int filter_order, + bool smoothing, + bool smooth_window, + bool estimate_compton) +{ + // Determine the first active slice of the data points. + while (begin != end && !begin->is_active) + ++begin; + if (begin == end) + return end; + for (vector::iterator it = begin + 1; it != end; ++it) + if (!it->is_active) { + end = it; + break; + } + if (end <= begin) // just in case + return end; + + vector bg_input(end - begin); + for (vector::iterator it = begin; it != end; it++) + bg_input[it - begin] = it->y; + + vector bg_output = ROOT::background(bg_input, + window_width, + direction, + filter_order, + smoothing, + smooth_window, + estimate_compton); + + if (bg_output.size() == bg_input.size()) { + for (vector::iterator it = begin; it != end; it++) + it->y = bg_output[it - begin]; + } + + return end; +} } // anonymous namespace From 1dfb964b24855e615a280165b0b2d4e53c35ed46 Mon Sep 17 00:00:00 2001 From: MarcBHahn <60821053+MarcBHahn@users.noreply.github.com> Date: Wed, 15 Jun 2022 14:41:04 +0200 Subject: [PATCH 7/8] Improve initial guess for shirleyBG determination Improve initial guess for shirleyBG determination to improve the starting values for the iterative procedure. --- fityk/transform.cpp | 49 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 39 insertions(+), 10 deletions(-) diff --git a/fityk/transform.cpp b/fityk/transform.cpp index b6277e51..7309973e 100644 --- a/fityk/transform.cpp +++ b/fityk/transform.cpp @@ -72,7 +72,7 @@ void merge_same_x(vector &pp, bool avg) void shirley_bg(vector &pp) { /* Set criteria to stop itertions */ - const int max_iter = 50; + const int max_iter = 100; const double max_rdiff = 1e-6; // implement: check for differences const unsigned long n = pp.size(); @@ -123,9 +123,10 @@ void shirley_bg(vector &pp) } } /* divide by number of summations */ - ymin=ymin/double(count_start); - ymax=ymax/double(count_end); - + if(count_start>0 && count_end>0){ + ymin=ymin/double(count_start); + ymax=ymax/double(count_end); + } /* Determine minimum element from data to obtain a reasonable starting value for the fit */ double yminElement = ymin; @@ -137,26 +138,54 @@ void shirley_bg(vector &pp) /* create starting background */ std::vector BG(n, yminElement); + /* set initial BG values by assuming a linear background of the form y = a* x + y0*/ + double divisor = 1.0; + if( (endXindex - startXindex) > 0){ + divisor = (double)(endXindex - startXindex); + } + + double initial_slope = (pp[endXindex].y - pp[startXindex].y) / divisor; + + for (unsigned long EnergyIdx = 0; EnergyIdx < n ; EnergyIdx++){ + + if (EnergyIdx < startXindex){ + BG[EnergyIdx] = pp[startXindex].y; + } + + if (EnergyIdx >= startXindex && EnergyIdx <= endXindex){ + BG[EnergyIdx] = initial_slope * (EnergyIdx-startXindex) + pp[startXindex].y; + } + + if (EnergyIdx > endXindex){ + BG[EnergyIdx] = pp[endXindex].y; + } + + } + /* recursive determination of the shirley background */ for (int iter = 0; iter < max_iter; iter++) { - /* Determine all background BG(Energy) values in dependence of the Energy */ + /* Determine all background BG(Energy) values in dependence of the binding energy */ for (unsigned long EnergyIdx = 0; EnergyIdx < n ; EnergyIdx++){ double Alow = 0.0; double Ahigh = 0.0; /* calculate Areas A1 and A2 for each energy value */ - for (unsigned long x = startXindex+1; x <= EnergyIdx ; x++){ - Alow = Alow + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]; + for (unsigned long x = startXindex+1; x <= EnergyIdx; x++){ + + Alow = Alow + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]*(pp[x].x-pp[x-1].x); + } - for (unsigned long x = EnergyIdx+1; x < endXindex ; x++){ - Ahigh = Ahigh + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]; + for (unsigned long x = EnergyIdx+1; x < endXindex; x++){ + + Ahigh = Ahigh + ((pp[x].x-pp[x-1].x)*(pp[x].y+pp[x-1].y)/2.0) - BG[x]*(pp[x].x-pp[x-1].x); + } /* The formula is implemented as given in the CasaXPS manual (2006) "Peak fitting in XPS" chapter */ - BG[EnergyIdx] = ymin + (ymax-ymin) * (Alow / (Ahigh+Alow)); + BG[EnergyIdx] = ymin + (ymax-ymin) * (Alow/(Ahigh+Alow)); } } From df4554b3f30bd1e5e63e7137d5344a0fda006fce Mon Sep 17 00:00:00 2001 From: MarcBHahn <60821053+MarcBHahn@users.noreply.github.com> Date: Wed, 15 Jun 2022 14:42:42 +0200 Subject: [PATCH 8/8] Update INSTALL Add installation information about autogen.sh --- INSTALL | 1 + 1 file changed, 1 insertion(+) diff --git a/INSTALL b/INSTALL index ff78a40a..de12d69a 100644 --- a/INSTALL +++ b/INSTALL @@ -14,6 +14,7 @@ Prerequisites Installation ============ +$ ./autogen.sh $ ./configure (optionally followed by options, see: Configuration Options) $ make # make install