From 4ee8f3022f0076d79a6c20d39c7a809fd1525f35 Mon Sep 17 00:00:00 2001 From: Cristiano Fontana Date: Thu, 24 Nov 2016 17:44:40 +0100 Subject: [PATCH 1/5] Added the one-dimensional background estimator from the ROOT TSpectrum --- fityk/root/background.cpp | 1100 +++++++++++++++++++++++++++++++++++++ fityk/root/background.hpp | 45 ++ 2 files changed, 1145 insertions(+) create mode 100644 fityk/root/background.cpp create mode 100644 fityk/root/background.hpp diff --git a/fityk/root/background.cpp b/fityk/root/background.cpp new file mode 100644 index 00000000..c39459f2 --- /dev/null +++ b/fityk/root/background.cpp @@ -0,0 +1,1100 @@ +// @(#)root/spectrum:$Id$ +// Author: Miroslav Morhac 27/05/99 + +/************************************************************************* + * Copyright (C) 1995-2006, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +// Modified by Cristiano Fontana 17/11/2016 +// Eliminated the dependency on ROOT + +#include +#include + +#include "background.hpp" + +//////////////////////////////////////////////////////////////////////////////// +/// This function calculates background spectrum from source spectrum. +/// The result is placed in the vector pointed by spectrum pointer. +/// The goal is to separate the useful information (peaks) from useless +/// information (background). +/// +/// - method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping +/// algorithm. +/// - new value in the channel "i" is calculated +/// +/// \f[ +/// v_p(i) = min \left\{ v_{p-1}(i)^{\frac{\left[v_{p-1}(i+p)+v_{p-1}(i-p)\right]}{2}} \right\} +/// \f] +/// +/// where p = 1, 2, ..., numberIterations. In fact it represents second order +/// difference filter (-1,2,-1). +/// +/// One can also change the +/// direction of the change of the clipping window, the order of the clipping +/// filter, to include smoothing, to set width of smoothing window and to include +/// the estimation of Compton edges. On successful completion it returns 0. On +/// error it returns pointer to the string describing error. +/// +/// #### Parameters: +/// +/// - spectrum: vector of source spectrum +/// - numberIterations: maximal width of clipping window, +/// - direction: direction of change of clipping window. +/// Possible values: kBackIncreasingWindow, kBackDecreasingWindow +/// - filterOrder: order of clipping filter. +/// Possible values: kBackOrder2, kBackOrder4, kBackOrder6, kBackOrder8 +/// - smoothing: logical variable whether the smoothing operation in the +/// estimation of background will be included. +/// Possible values: kFALSE, kTRUE +/// - smoothWindow: width of smoothing window. +/// Possible values: kBackSmoothing3, kBackSmoothing5, kBackSmoothing7, +/// kBackSmoothing9, kBackSmoothing11, kBackSmoothing13, kBackSmoothing15. +/// - compton: logical variable whether the estimation of Compton edge will be +/// included. Possible values: kFALSE, kTRUE. +/// +/// #### Return: +/// +/// A vector with the estimated background. +/// +/// #### References: +/// +/// 1. C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the +/// quantitative analysis of PIXE spectra in geoscience applications. NIM, B34 +/// (1988), 396-402. +/// +/// 2. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo: +/// Background elimination methods for multidimensional gamma-ray spectra. NIM, +/// A401 (1997) 113-132. +/// +/// 3. D. D. Burgess, R. J. Tervo: Background estimation for gamma-ray +/// spectroscopy. NIM 214 (1983), 431-434. +/// +/// ## Examples from the ROOT documentation +/// +/// ### Example 1 script Background_incr.c: +/// +/// \image html TSpectrum_Background_incr.jpg Fig. 1 Example of the estimation of background for number of iterations=6. Original spectrum is shown in black color, estimated background in red color. +/// +/// #### Script: +/// +/// ~~~ {.cpp} +/// // Example to illustrate the background estimator (class TSpectrum). +/// // To execute this example, do +/// // root > .x Background_incr.C +/// +/// #include +/// +/// void Background_incr() { +/// Int_t i; +/// Double_t nbins = 256; +/// Double_t xmin = 0; +/// Double_t xmax = nbins; +/// Double_t * source = new Double_t[nbins]; +/// TH1F *back = new TH1F("back","",nbins,xmin,xmax); +/// TH1F *d = new TH1F("d","",nbins,xmin,xmax); +/// TFile *f = new TFile("spectra/TSpectrum.root"); +/// back=(TH1F*) f->Get("back1;1"); +/// TCanvas *Background = gROOT->GetListOfCanvases()->FindObject("Background"); +/// if (!Background) Background = +/// new TCanvas("Background", +/// "Estimation of background with increasing window", +/// 10,10,1000,700); +/// back->Draw("L"); +/// TSpectrum *s = new TSpectrum(); +/// for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1); +/// s->Background(source,nbins,6,kBackIncreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]); +/// d->SetLineColor(kRed); +/// d->Draw("SAME L"); +/// } +/// ~~~ +/// +/// ### Example 2 script Background_decr.c: +/// +/// In Fig. 1. one can notice that at the edges of the peaks the estimated +/// background goes under the peaks. An alternative approach is to decrease the +/// clipping window from a given value numberIterations to the value of one, which +/// is presented in this example. +/// +/// \image html TSpectrum_Background_decr.jpg Fig. 2 Example of the estimation of background for numberIterations=6 using decreasing clipping window algorithm. Original spectrum is shown in black color, estimated background in red color. +/// +/// #### Script: +/// +/// ~~~ {.cpp} +/// // Example to illustrate the background estimator (class TSpectrum). +/// // To execute this example, do +/// // root > .x Background_decr.C +/// +/// #include +/// +/// void Background_decr() { +/// Int_t i; +/// Double_t nbins = 256; +/// Double_t xmin = 0; +/// Double_t xmax = nbins; +/// Double_t * source = new Double_t[nbins]; +/// TH1F *back = new TH1F("back","",nbins,xmin,xmax); +/// TH1F *d = new TH1F("d","",nbins,xmin,xmax); +/// TFile *f = new TFile("spectra/TSpectrum.root"); +/// back=(TH1F*) f->Get("back1;1"); +/// TCanvas *Background = gROOT->GetListOfCanvases()->FindObject("Background"); +/// if (!Background) Background = +/// new TCanvas("Background","Estimation of background with decreasing window", +/// 10,10,1000,700); +/// back->Draw("L"); +/// TSpectrum *s = new TSpectrum(); +/// for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1); +/// s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]); +/// d->SetLineColor(kRed); +/// d->Draw("SAME L"); +/// } +/// ~~~ +/// +/// ### Example 3 script Background_width.c: +/// +/// The question is how to choose the width of the clipping window, i.e., +/// numberIterations parameter. The influence of this parameter on the estimated +/// background is illustrated in Fig. 3. +/// +/// \image html TSpectrum_Background_width.jpg Fig. 3 Example of the influence of clipping window width on the estimated background for numberIterations=4 (red line), 6 (blue line) 8 (green line) using decreasing clipping window algorithm. +/// +/// in general one should set this parameter so that the value +/// 2*numberIterations+1 was greater than the widths of preserved objects (peaks). +/// +/// #### Script: +/// +/// ~~~ {.cpp} +/// // Example to illustrate the influence of the clipping window width on the +/// // estimated background. To execute this example, do: +/// // root > .x Background_width.C +/// +/// #include +/// +/// void Background_width() { +/// Int_t i; +/// Double_t nbins = 256; +/// Double_t xmin = 0; +/// Double_t xmax = nbins; +/// Double_t * source = new Double_t[nbins]; +/// TH1F *h = new TH1F("h","",nbins,xmin,xmax); +/// TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax); +/// TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax); +/// TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax); +/// TFile *f = new TFile("spectra/TSpectrum.root"); +/// h=(TH1F*) f->Get("back1;1"); +/// TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background"); +/// if (!background) background = new TCanvas("background", +/// "Influence of clipping window width on the estimated background", +/// 10,10,1000,700); +/// h->Draw("L"); +/// TSpectrum *s = new TSpectrum(); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,4,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]); +/// d1->SetLineColor(kRed); +/// d1->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]); +/// d2->SetLineColor(kBlue); +/// d2->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,8,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]); +/// d3->SetLineColor(kGreen); +/// d3->Draw("SAME L"); +/// } +/// ~~~ +/// +/// ### Example 4 script Background_width2.c: +/// +/// another example for very complex spectrum is given in Fig. 4. +/// +/// \image html TSpectrum_Background_width2.jpg Fig. 4 Example of the influence of clipping window width on the estimated background for numberIterations=10 (red line), 20 (blue line), 30 (green line) and 40 (magenta line) using decreasing clipping window algorithm. +/// +/// #### Script: +/// +/// ~~~ {.cpp} +/// // Example to illustrate the influence of the clipping window width on the +/// // estimated background. To execute this example, do: +/// // root > .x Background_width2.C +/// +/// #include +/// +/// void Background_width2() { +/// Int_t i; +/// Double_t nbins = 4096; +/// Double_t xmin = 0; +/// Double_t xmax = 4096; +/// Double_t * source = new Double_t[nbins]; +/// TH1F *h = new TH1F("h","",nbins,xmin,xmax); +/// TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax); +/// TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax); +/// TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax); +/// TH1F *d4 = new TH1F("d4","",nbins,xmin,xmax); +/// TFile *f = new TFile("spectra/TSpectrum.root"); +/// h=(TH1F*) f->Get("back2;1"); +/// TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background"); +/// if (!background) background = new TCanvas("background", +/// "Influence of clipping window width on the estimated background", +/// 10,10,1000,700); +/// h->SetAxisRange(0,1000); +/// h->SetMaximum(20000); +/// h->Draw("L"); +/// TSpectrum *s = new TSpectrum(); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]); +/// d1->SetLineColor(kRed); +/// d1->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,20,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]); +/// d2->SetLineColor(kBlue); +/// d2->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,30,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]); +/// d3->SetLineColor(kGreen); +/// d3->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d4->SetBinContent(i + 1,source[i]); +/// d4->SetLineColor(kMagenta); +/// d4->Draw("SAME L"); +/// } +/// ~~~ +/// +/// ### Example 5 script Background_order.c: +/// +/// Second order difference filter removes linear (quasi-linear) background and +/// preserves symmetrical peaks. However if the shape of the background is more +/// complex one can employ higher-order clipping filters (see example in Fig. 5) +/// +/// \image html TSpectrum_Background_order.jpg Fig. 5 Example of the influence of clipping filter difference order on the estimated background for fNnumberIterations=40, 2-nd order red line, 4-th order blue line, 6-th order green line and 8-th order magenta line, and using decreasing clipping window algorithm. +/// +/// #### Script: +/// +/// ~~~ {.cpp} +/// // Example to illustrate the influence of the clipping filter difference order +/// // on the estimated background. To execute this example, do +/// // root > .x Background_order.C +/// +/// #include +/// +/// void Background_order() { +/// Int_t i; +/// Double_t nbins = 4096; +/// Double_t xmin = 0; +/// Double_t xmax = 4096; +/// Double_t * source = new Double_t[nbins]; +/// TH1F *h = new TH1F("h","",nbins,xmin,xmax); +/// TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax); +/// TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax); +/// TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax); +/// TH1F *d4 = new TH1F("d4","",nbins,xmin,xmax); +/// TFile *f = new TFile("spectra/TSpectrum.root"); +/// h=(TH1F*) f->Get("back2;1"); +/// TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background"); +/// if (!background) background = new TCanvas("background", +/// "Influence of clipping filter difference order on the estimated background", +/// 10,10,1000,700); +/// h->SetAxisRange(1220,1460); +/// h->SetMaximum(11000); +/// h->Draw("L"); +/// TSpectrum *s = new TSpectrum(); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]); +/// d1->SetLineColor(kRed); +/// d1->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder4,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]); +/// d2->SetLineColor(kBlue); +/// d2->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder6,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]); +/// d3->SetLineColor(kGreen); +/// d3->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder8,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d4->SetBinContent(i + 1,source[i]); +/// d4->SetLineColor(kMagenta); +/// d4->Draw("SAME L"); +/// } +/// ~~~ +/// +/// ### Example 6 script Background_smooth.c: +/// +/// The estimate of the background can be influenced by noise present in the +/// spectrum. We proposed the algorithm of the background estimate with +/// simultaneous smoothing. In the original algorithm without smoothing, the +/// estimated background snatches the lower spikes in the noise. Consequently, +/// the areas of peaks are biased by this error. +/// +/// \image html TSpectrum_Background_smooth1.jpg Fig. 7 Principle of background estimation algorithm with simultaneous smoothing. +/// \image html TSpectrum_Background_smooth2.jpg Fig. 8 Illustration of non-smoothing (red line) and smoothing algorithm of background estimation (blue line). +/// +/// #### Script: +/// +/// ~~~ {.cpp} +/// // Example to illustrate the background estimator (class TSpectrum) including +/// // Compton edges. To execute this example, do: +/// // root > .x Background_smooth.C +/// +/// #include +/// +/// void Background_smooth() { +/// Int_t i; +/// Double_t nbins = 4096; +/// Double_t xmin = 0; +/// Double_t xmax = nbins; +/// Double_t * source = new Double_t[nbins]; +/// TH1F *h = new TH1F("h","",nbins,xmin,xmax); +/// TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax); +/// TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax); +/// TFile *f = new TFile("spectra/TSpectrum.root"); +/// h=(TH1F*) f->Get("back4;1"); +/// TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background"); +/// if (!background) background = new TCanvas("background", +/// "Estimation of background with noise",10,10,1000,700); +/// h->SetAxisRange(3460,3830); +/// h->Draw("L"); +/// TSpectrum *s = new TSpectrum(); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]); +/// d1->SetLineColor(kRed); +/// d1->Draw("SAME L"); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kTRUE, +/// kBackSmoothing3,kFALSE); +/// for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]); +/// d2->SetLineColor(kBlue); +/// d2->Draw("SAME L"); +/// } +/// ~~~ +/// +/// ### Example 8 script Background_compton.c: +/// +/// Sometimes it is necessary to include also the Compton edges into the estimate of +/// the background. In Fig. 8 we present the example of the synthetic spectrum +/// with Compton edges. The background was estimated using the 8-th order filter +/// with the estimation of the Compton edges using decreasing +/// clipping window algorithm (numberIterations=10) with smoothing +/// (smoothingWindow=5). +/// +/// \image html TSpectrum_Background_compton.jpg Fig. 8 Example of the estimate of the background with Compton edges (red line) for numberIterations=10, 8-th order difference filter, using decreasing clipping window algorithm and smoothing (smoothingWindow=5). +/// +/// #### Script: +/// +/// ~~~ {.cpp} +/// // Example to illustrate the background estimator (class TSpectrum) including +/// // Compton edges. To execute this example, do: +/// // root > .x Background_compton.C +/// +/// #include +/// +/// void Background_compton() { +/// Int_t i; +/// Double_t nbins = 512; +/// Double_t xmin = 0; +/// Double_t xmax = nbins; +/// Double_t * source = new Double_t[nbins]; +/// TH1F *h = new TH1F("h","",nbins,xmin,xmax); +/// TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax); +/// TFile *f = new TFile("spectra/TSpectrum.root"); +/// h=(TH1F*) f->Get("back3;1"); +/// TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background"); +/// if (!background) background = new TCanvas("background", +/// "Estimation of background with Compton edges under peaks",10,10,1000,700); +/// h->Draw("L"); +/// TSpectrum *s = new TSpectrum(); +/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1); +/// s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder8,kTRUE, +/// kBackSmoothing5,,kTRUE); +/// for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]); +/// d1->SetLineColor(kRed); +/// d1->Draw("SAME L"); +/// } +/// ~~~ + +std::vector background::background(std::vector spectrum, + int numberIterations, + int direction, + int filterOrder, + bool smoothing, + int smoothWindow, + bool compton) +{ + int i, j, w, bw, b1, b2, priz; + double a, b, c, d, e, yb1, yb2, ai, av, men, b4, c4, d4, e4, b6, c6, d6, e6, f6, g6, b8, c8, d8, e8, f8, g8, h8, i8; + + const unsigned int ssize = spectrum.size(); + + if (ssize <= 0) + throw "Wrong Parameters"; + + if (numberIterations < 1) + throw "Width of Clipping Window Must Be Positive"; + + if (ssize < 2 * numberIterations + 1) + throw "Too Large Clipping Window"; + + if (smoothing == true && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15) + throw "Incorrect width of smoothing window"; + + std::vector working_space(2 * ssize); + + for (i = 0; i < ssize; i++) + { + working_space[i] = spectrum[i]; + working_space[i + ssize] = spectrum[i]; + } + + bw=(smoothWindow-1)/2; + + if (direction == kBackIncreasingWindow) + i = 1; + else if(direction == kBackDecreasingWindow) + i = numberIterations; + + if (filterOrder == kBackOrder2) { + do{ + for (j = i; j < ssize - i; j++) { + if (smoothing == false){ + a = working_space[ssize + j]; + b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0; + if (b < a) + a = b; + working_space[j] = a; + } + + else if (smoothing == true){ + a = working_space[ssize + j]; + av = 0; + men = 0; + for (w = j - bw; w <= j + bw; w++){ + if ( w >= 0 && w < ssize){ + av += working_space[ssize + w]; + men +=1; + } + } + av = av / men; + b = 0; + men = 0; + for (w = j - i - bw; w <= j - i + bw; w++){ + if ( w >= 0 && w < ssize){ + b += working_space[ssize + w]; + men +=1; + } + } + b = b / men; + c = 0; + men = 0; + for (w = j + i - bw; w <= j + i + bw; w++){ + if ( w >= 0 && w < ssize){ + c += working_space[ssize + w]; + men +=1; + } + } + c = c / men; + b = (b + c) / 2; + if (b < a) + av = b; + working_space[j]=av; + } + } + for (j = i; j < ssize - i; j++) + working_space[ssize + j] = working_space[j]; + if (direction == kBackIncreasingWindow) + i+=1; + else if(direction == kBackDecreasingWindow) + i-=1; + }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1)); + } + + else if (filterOrder == kBackOrder4) { + do{ + for (j = i; j < ssize - i; j++) { + if (smoothing == false){ + a = working_space[ssize + j]; + b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0; + c = 0; + ai = i / 2; + c -= working_space[ssize + j - (int) (2 * ai)] / 6; + c += 4 * working_space[ssize + j - (int) ai] / 6; + c += 4 * working_space[ssize + j + (int) ai] / 6; + c -= working_space[ssize + j + (int) (2 * ai)] / 6; + if (b < c) + b = c; + if (b < a) + a = b; + working_space[j] = a; + } + + else if (smoothing == true){ + a = working_space[ssize + j]; + av = 0; + men = 0; + for (w = j - bw; w <= j + bw; w++){ + if ( w >= 0 && w < ssize){ + av += working_space[ssize + w]; + men +=1; + } + } + av = av / men; + b = 0; + men = 0; + for (w = j - i - bw; w <= j - i + bw; w++){ + if ( w >= 0 && w < ssize){ + b += working_space[ssize + w]; + men +=1; + } + } + b = b / men; + c = 0; + men = 0; + for (w = j + i - bw; w <= j + i + bw; w++){ + if ( w >= 0 && w < ssize){ + c += working_space[ssize + w]; + men +=1; + } + } + c = c / men; + b = (b + c) / 2; + ai = i / 2; + b4 = 0, men = 0; + for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + b4 += working_space[ssize + w]; + men +=1; + } + } + b4 = b4 / men; + c4 = 0, men = 0; + for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + c4 += working_space[ssize + w]; + men +=1; + } + } + c4 = c4 / men; + d4 = 0, men = 0; + for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + d4 += working_space[ssize + w]; + men +=1; + } + } + d4 = d4 / men; + e4 = 0, men = 0; + for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + e4 += working_space[ssize + w]; + men +=1; + } + } + e4 = e4 / men; + b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6; + if (b < b4) + b = b4; + if (b < a) + av = b; + working_space[j]=av; + } + } + for (j = i; j < ssize - i; j++) + working_space[ssize + j] = working_space[j]; + if (direction == kBackIncreasingWindow) + i+=1; + else if(direction == kBackDecreasingWindow) + i-=1; + }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1)); + } + + else if (filterOrder == kBackOrder6) { + do{ + for (j = i; j < ssize - i; j++) { + if (smoothing == false){ + a = working_space[ssize + j]; + b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0; + c = 0; + ai = i / 2; + c -= working_space[ssize + j - (int) (2 * ai)] / 6; + c += 4 * working_space[ssize + j - (int) ai] / 6; + c += 4 * working_space[ssize + j + (int) ai] / 6; + c -= working_space[ssize + j + (int) (2 * ai)] / 6; + d = 0; + ai = i / 3; + d += working_space[ssize + j - (int) (3 * ai)] / 20; + d -= 6 * working_space[ssize + j - (int) (2 * ai)] / 20; + d += 15 * working_space[ssize + j - (int) ai] / 20; + d += 15 * working_space[ssize + j + (int) ai] / 20; + d -= 6 * working_space[ssize + j + (int) (2 * ai)] / 20; + d += working_space[ssize + j + (int) (3 * ai)] / 20; + if (b < d) + b = d; + if (b < c) + b = c; + if (b < a) + a = b; + working_space[j] = a; + } + + else if (smoothing == true){ + a = working_space[ssize + j]; + av = 0; + men = 0; + for (w = j - bw; w <= j + bw; w++){ + if ( w >= 0 && w < ssize){ + av += working_space[ssize + w]; + men +=1; + } + } + av = av / men; + b = 0; + men = 0; + for (w = j - i - bw; w <= j - i + bw; w++){ + if ( w >= 0 && w < ssize){ + b += working_space[ssize + w]; + men +=1; + } + } + b = b / men; + c = 0; + men = 0; + for (w = j + i - bw; w <= j + i + bw; w++){ + if ( w >= 0 && w < ssize){ + c += working_space[ssize + w]; + men +=1; + } + } + c = c / men; + b = (b + c) / 2; + ai = i / 2; + b4 = 0, men = 0; + for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + b4 += working_space[ssize + w]; + men +=1; + } + } + b4 = b4 / men; + c4 = 0, men = 0; + for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + c4 += working_space[ssize + w]; + men +=1; + } + } + c4 = c4 / men; + d4 = 0, men = 0; + for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + d4 += working_space[ssize + w]; + men +=1; + } + } + d4 = d4 / men; + e4 = 0, men = 0; + for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + e4 += working_space[ssize + w]; + men +=1; + } + } + e4 = e4 / men; + b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6; + ai = i / 3; + b6 = 0, men = 0; + for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + b6 += working_space[ssize + w]; + men +=1; + } + } + b6 = b6 / men; + c6 = 0, men = 0; + for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + c6 += working_space[ssize + w]; + men +=1; + } + } + c6 = c6 / men; + d6 = 0, men = 0; + for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + d6 += working_space[ssize + w]; + men +=1; + } + } + d6 = d6 / men; + e6 = 0, men = 0; + for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + e6 += working_space[ssize + w]; + men +=1; + } + } + e6 = e6 / men; + f6 = 0, men = 0; + for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + f6 += working_space[ssize + w]; + men +=1; + } + } + f6 = f6 / men; + g6 = 0, men = 0; + for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + g6 += working_space[ssize + w]; + men +=1; + } + } + g6 = g6 / men; + b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20; + if (b < b6) + b = b6; + if (b < b4) + b = b4; + if (b < a) + av = b; + working_space[j]=av; + } + } + for (j = i; j < ssize - i; j++) + working_space[ssize + j] = working_space[j]; + if (direction == kBackIncreasingWindow) + i+=1; + else if(direction == kBackDecreasingWindow) + i-=1; + }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1)); + } + + else if (filterOrder == kBackOrder8) { + do{ + for (j = i; j < ssize - i; j++) { + if (smoothing == false){ + a = working_space[ssize + j]; + b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0; + c = 0; + ai = i / 2; + c -= working_space[ssize + j - (int) (2 * ai)] / 6; + c += 4 * working_space[ssize + j - (int) ai] / 6; + c += 4 * working_space[ssize + j + (int) ai] / 6; + c -= working_space[ssize + j + (int) (2 * ai)] / 6; + d = 0; + ai = i / 3; + d += working_space[ssize + j - (int) (3 * ai)] / 20; + d -= 6 * working_space[ssize + j - (int) (2 * ai)] / 20; + d += 15 * working_space[ssize + j - (int) ai] / 20; + d += 15 * working_space[ssize + j + (int) ai] / 20; + d -= 6 * working_space[ssize + j + (int) (2 * ai)] / 20; + d += working_space[ssize + j + (int) (3 * ai)] / 20; + e = 0; + ai = i / 4; + e -= working_space[ssize + j - (int) (4 * ai)] / 70; + e += 8 * working_space[ssize + j - (int) (3 * ai)] / 70; + e -= 28 * working_space[ssize + j - (int) (2 * ai)] / 70; + e += 56 * working_space[ssize + j - (int) ai] / 70; + e += 56 * working_space[ssize + j + (int) ai] / 70; + e -= 28 * working_space[ssize + j + (int) (2 * ai)] / 70; + e += 8 * working_space[ssize + j + (int) (3 * ai)] / 70; + e -= working_space[ssize + j + (int) (4 * ai)] / 70; + if (b < e) + b = e; + if (b < d) + b = d; + if (b < c) + b = c; + if (b < a) + a = b; + working_space[j] = a; + } + + else if (smoothing == true){ + a = working_space[ssize + j]; + av = 0; + men = 0; + for (w = j - bw; w <= j + bw; w++){ + if ( w >= 0 && w < ssize){ + av += working_space[ssize + w]; + men +=1; + } + } + av = av / men; + b = 0; + men = 0; + for (w = j - i - bw; w <= j - i + bw; w++){ + if ( w >= 0 && w < ssize){ + b += working_space[ssize + w]; + men +=1; + } + } + b = b / men; + c = 0; + men = 0; + for (w = j + i - bw; w <= j + i + bw; w++){ + if ( w >= 0 && w < ssize){ + c += working_space[ssize + w]; + men +=1; + } + } + c = c / men; + b = (b + c) / 2; + ai = i / 2; + b4 = 0, men = 0; + for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + b4 += working_space[ssize + w]; + men +=1; + } + } + b4 = b4 / men; + c4 = 0, men = 0; + for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + c4 += working_space[ssize + w]; + men +=1; + } + } + c4 = c4 / men; + d4 = 0, men = 0; + for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + d4 += working_space[ssize + w]; + men +=1; + } + } + d4 = d4 / men; + e4 = 0, men = 0; + for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + e4 += working_space[ssize + w]; + men +=1; + } + } + e4 = e4 / men; + b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6; + ai = i / 3; + b6 = 0, men = 0; + for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + b6 += working_space[ssize + w]; + men +=1; + } + } + b6 = b6 / men; + c6 = 0, men = 0; + for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + c6 += working_space[ssize + w]; + men +=1; + } + } + c6 = c6 / men; + d6 = 0, men = 0; + for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + d6 += working_space[ssize + w]; + men +=1; + } + } + d6 = d6 / men; + e6 = 0, men = 0; + for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + e6 += working_space[ssize + w]; + men +=1; + } + } + e6 = e6 / men; + f6 = 0, men = 0; + for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + f6 += working_space[ssize + w]; + men +=1; + } + } + f6 = f6 / men; + g6 = 0, men = 0; + for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + g6 += working_space[ssize + w]; + men +=1; + } + } + g6 = g6 / men; + b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20; + ai = i / 4; + b8 = 0, men = 0; + for (w = j - (int)(4 * ai) - bw; w <= j - (int)(4 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + b8 += working_space[ssize + w]; + men +=1; + } + } + b8 = b8 / men; + c8 = 0, men = 0; + for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + c8 += working_space[ssize + w]; + men +=1; + } + } + c8 = c8 / men; + d8 = 0, men = 0; + for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + d8 += working_space[ssize + w]; + men +=1; + } + } + d8 = d8 / men; + e8 = 0, men = 0; + for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + e8 += working_space[ssize + w]; + men +=1; + } + } + e8 = e8 / men; + f8 = 0, men = 0; + for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ + if (w >= 0 && w < ssize){ + f8 += working_space[ssize + w]; + men +=1; + } + } + f8 = f8 / men; + g8 = 0, men = 0; + for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + g8 += working_space[ssize + w]; + men +=1; + } + } + g8 = g8 / men; + h8 = 0, men = 0; + for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + h8 += working_space[ssize + w]; + men +=1; + } + } + h8 = h8 / men; + i8 = 0, men = 0; + for (w = j + (int)(4 * ai) - bw; w <= j + (int)(4 * ai) + bw; w++){ + if (w >= 0 && w < ssize){ + i8 += working_space[ssize + w]; + men +=1; + } + } + i8 = i8 / men; + b8 = ( -b8 + 8 * c8 - 28 * d8 + 56 * e8 - 56 * f8 - 28 * g8 + 8 * h8 - i8)/70; + if (b < b8) + b = b8; + if (b < b6) + b = b6; + if (b < b4) + b = b4; + if (b < a) + av = b; + working_space[j]=av; + } + } + for (j = i; j < ssize - i; j++) + working_space[ssize + j] = working_space[j]; + if (direction == kBackIncreasingWindow) + i += 1; + else if(direction == kBackDecreasingWindow) + i -= 1; + }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1)); + } + + if (compton == true) { + for (i = 0, b2 = 0; i < ssize; i++){ + b1 = b2; + a = working_space[i], b = spectrum[i]; + j = i; + if (fabs(a - b) >= 1) { + b1 = i - 1; + if (b1 < 0) + b1 = 0; + yb1 = working_space[b1]; + for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){ + a = working_space[b2], b = spectrum[b2]; + c = c + b - yb1; + if (fabs(a - b) < 1) { + priz = 1; + yb2 = b; + } + } + if (b2 == ssize) + b2 -= 1; + yb2 = working_space[b2]; + if (yb1 <= yb2){ + for (j = b1, c = 0; j <= b2; j++){ + b = spectrum[j]; + c = c + b - yb1; + } + if (c > 1){ + c = (yb2 - yb1) / c; + for (j = b1, d = 0; j <= b2 && j < ssize; j++){ + b = spectrum[j]; + d = d + b - yb1; + a = c * d + yb1; + working_space[ssize + j] = a; + } + } + } + + else{ + for (j = b2, c = 0; j >= b1; j--){ + b = spectrum[j]; + c = c + b - yb2; + } + if (c > 1){ + c = (yb1 - yb2) / c; + for (j = b2, d = 0;j >= b1 && j >= 0; j--){ + b = spectrum[j]; + d = d + b - yb2; + a = c * d + yb2; + working_space[ssize + j] = a; + } + } + } + i=b2; + } + } + } + + working_space.resize(ssize); + + return working_space; +} diff --git a/fityk/root/background.hpp b/fityk/root/background.hpp new file mode 100644 index 00000000..efe016b8 --- /dev/null +++ b/fityk/root/background.hpp @@ -0,0 +1,45 @@ +// @(#)root/spectrum:$Id$ +// Author: Miroslav Morhac 27/05/99 + +/************************************************************************* + * Copyright (C) 1995-2006, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +// Modified by Cristiano Fontana 17/11/2016 +// Eliminated the dependency on ROOT + +#ifndef __BACKGROUND_HPP__ +#define __BACKGROUND_HPP__ + +namespace background { + + enum { + kBackIncreasingWindow =0, + kBackDecreasingWindow =1, + kBackOrder2 =2, + kBackOrder4 =4, + kBackOrder6 =6, + kBackOrder8 =8, + kBackSmoothing3 =3, + kBackSmoothing5 =5, + kBackSmoothing7 =7, + kBackSmoothing9 =9, + kBackSmoothing11 =11, + kBackSmoothing13 =13, + kBackSmoothing15 =15 + }; + + std::vector background(std::vector spectrum, + int numberIterations, + int direction, + int filterOrder, + bool smoothing, + int smoothWindow, + bool compton); +} + +#endif From ce63a5588bd448e495c502541fe6a49369f44105 Mon Sep 17 00:00:00 2001 From: Cristiano Fontana Date: Fri, 28 Dec 2018 08:52:10 +0100 Subject: [PATCH 2/5] Modified to play nicely with fityk datatypes --- fityk/root/background.cpp | 236 +++++++++++++++++++------------------- fityk/root/background.hpp | 46 ++++---- 2 files changed, 145 insertions(+), 137 deletions(-) diff --git a/fityk/root/background.cpp b/fityk/root/background.cpp index c39459f2..f18bb323 100644 --- a/fityk/root/background.cpp +++ b/fityk/root/background.cpp @@ -15,6 +15,7 @@ #include #include +#include "fityk.h" #include "background.hpp" //////////////////////////////////////////////////////////////////////////////// @@ -441,7 +442,7 @@ /// } /// ~~~ -std::vector background::background(std::vector spectrum, +std::vector ROOT::background(const std::vector spectrum, int numberIterations, int direction, int filterOrder, @@ -455,18 +456,21 @@ std::vector background::background(std::vector spectrum, const unsigned int ssize = spectrum.size(); if (ssize <= 0) - throw "Wrong Parameters"; + throw fityk::ExecuteError("Wrong vector size"); if (numberIterations < 1) - throw "Width of Clipping Window Must Be Positive"; + throw fityk::ExecuteError("Width of Clipping Window Must Be Positive"); if (ssize < 2 * numberIterations + 1) - throw "Too Large Clipping Window"; + throw fityk::ExecuteError("Too Large Clipping Window"); + + if (filterOrder != kBackOrder2 && filterOrder != kBackOrder4 && filterOrder != kBackOrder6 && filterOrder != kBackOrder8) + throw fityk::ExecuteError("Incorrect order of clipping filter, possible values are: 2, 4, 6 or 8"); if (smoothing == true && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15) - throw "Incorrect width of smoothing window"; + throw fityk::ExecuteError("Incorrect width of smoothing window"); - std::vector working_space(2 * ssize); + std::vector working_space(2 * ssize); for (i = 0; i < ssize; i++) { @@ -485,20 +489,20 @@ std::vector background::background(std::vector spectrum, do{ for (j = i; j < ssize - i; j++) { if (smoothing == false){ - a = working_space[ssize + j]; - b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0; + a = working_space[ssize + j].y; + b = (working_space[ssize + j - i].y + working_space[ssize + j + i].y) / 2.0; if (b < a) a = b; - working_space[j] = a; + working_space[j].y = a; } else if (smoothing == true){ - a = working_space[ssize + j]; + a = working_space[ssize + j].y; av = 0; men = 0; for (w = j - bw; w <= j + bw; w++){ if ( w >= 0 && w < ssize){ - av += working_space[ssize + w]; + av += working_space[ssize + w].y; men +=1; } } @@ -507,7 +511,7 @@ std::vector background::background(std::vector spectrum, men = 0; for (w = j - i - bw; w <= j - i + bw; w++){ if ( w >= 0 && w < ssize){ - b += working_space[ssize + w]; + b += working_space[ssize + w].y; men +=1; } } @@ -516,7 +520,7 @@ std::vector background::background(std::vector spectrum, men = 0; for (w = j + i - bw; w <= j + i + bw; w++){ if ( w >= 0 && w < ssize){ - c += working_space[ssize + w]; + c += working_space[ssize + w].y; men +=1; } } @@ -524,11 +528,11 @@ std::vector background::background(std::vector spectrum, b = (b + c) / 2; if (b < a) av = b; - working_space[j]=av; + working_space[j].y=av; } } for (j = i; j < ssize - i; j++) - working_space[ssize + j] = working_space[j]; + working_space[ssize + j].y = working_space[j].y; if (direction == kBackIncreasingWindow) i+=1; else if(direction == kBackDecreasingWindow) @@ -540,28 +544,28 @@ std::vector background::background(std::vector spectrum, do{ for (j = i; j < ssize - i; j++) { if (smoothing == false){ - a = working_space[ssize + j]; - b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0; + a = working_space[ssize + j].y; + b = (working_space[ssize + j - i].y + working_space[ssize + j + i].y) / 2.0; c = 0; ai = i / 2; - c -= working_space[ssize + j - (int) (2 * ai)] / 6; - c += 4 * working_space[ssize + j - (int) ai] / 6; - c += 4 * working_space[ssize + j + (int) ai] / 6; - c -= working_space[ssize + j + (int) (2 * ai)] / 6; + c -= working_space[ssize + j - (int) (2 * ai)].y / 6; + c += 4 * working_space[ssize + j - (int) ai].y / 6; + c += 4 * working_space[ssize + j + (int) ai].y / 6; + c -= working_space[ssize + j + (int) (2 * ai)].y / 6; if (b < c) b = c; if (b < a) a = b; - working_space[j] = a; + working_space[j].y = a; } else if (smoothing == true){ - a = working_space[ssize + j]; + a = working_space[ssize + j].y; av = 0; men = 0; for (w = j - bw; w <= j + bw; w++){ if ( w >= 0 && w < ssize){ - av += working_space[ssize + w]; + av += working_space[ssize + w].y; men +=1; } } @@ -570,7 +574,7 @@ std::vector background::background(std::vector spectrum, men = 0; for (w = j - i - bw; w <= j - i + bw; w++){ if ( w >= 0 && w < ssize){ - b += working_space[ssize + w]; + b += working_space[ssize + w].y; men +=1; } } @@ -579,7 +583,7 @@ std::vector background::background(std::vector spectrum, men = 0; for (w = j + i - bw; w <= j + i + bw; w++){ if ( w >= 0 && w < ssize){ - c += working_space[ssize + w]; + c += working_space[ssize + w].y; men +=1; } } @@ -589,7 +593,7 @@ std::vector background::background(std::vector spectrum, b4 = 0, men = 0; for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - b4 += working_space[ssize + w]; + b4 += working_space[ssize + w].y; men +=1; } } @@ -597,7 +601,7 @@ std::vector background::background(std::vector spectrum, c4 = 0, men = 0; for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - c4 += working_space[ssize + w]; + c4 += working_space[ssize + w].y; men +=1; } } @@ -605,7 +609,7 @@ std::vector background::background(std::vector spectrum, d4 = 0, men = 0; for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - d4 += working_space[ssize + w]; + d4 += working_space[ssize + w].y; men +=1; } } @@ -613,7 +617,7 @@ std::vector background::background(std::vector spectrum, e4 = 0, men = 0; for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - e4 += working_space[ssize + w]; + e4 += working_space[ssize + w].y; men +=1; } } @@ -623,11 +627,11 @@ std::vector background::background(std::vector spectrum, b = b4; if (b < a) av = b; - working_space[j]=av; + working_space[j].y=av; } } for (j = i; j < ssize - i; j++) - working_space[ssize + j] = working_space[j]; + working_space[ssize + j].y = working_space[j].y; if (direction == kBackIncreasingWindow) i+=1; else if(direction == kBackDecreasingWindow) @@ -639,38 +643,38 @@ std::vector background::background(std::vector spectrum, do{ for (j = i; j < ssize - i; j++) { if (smoothing == false){ - a = working_space[ssize + j]; - b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0; + a = working_space[ssize + j].y; + b = (working_space[ssize + j - i].y + working_space[ssize + j + i].y) / 2.0; c = 0; ai = i / 2; - c -= working_space[ssize + j - (int) (2 * ai)] / 6; - c += 4 * working_space[ssize + j - (int) ai] / 6; - c += 4 * working_space[ssize + j + (int) ai] / 6; - c -= working_space[ssize + j + (int) (2 * ai)] / 6; + c -= working_space[ssize + j - (int) (2 * ai)].y / 6; + c += 4 * working_space[ssize + j - (int) ai].y / 6; + c += 4 * working_space[ssize + j + (int) ai].y / 6; + c -= working_space[ssize + j + (int) (2 * ai)].y / 6; d = 0; ai = i / 3; - d += working_space[ssize + j - (int) (3 * ai)] / 20; - d -= 6 * working_space[ssize + j - (int) (2 * ai)] / 20; - d += 15 * working_space[ssize + j - (int) ai] / 20; - d += 15 * working_space[ssize + j + (int) ai] / 20; - d -= 6 * working_space[ssize + j + (int) (2 * ai)] / 20; - d += working_space[ssize + j + (int) (3 * ai)] / 20; + d += working_space[ssize + j - (int) (3 * ai)].y / 20; + d -= 6 * working_space[ssize + j - (int) (2 * ai)].y / 20; + d += 15 * working_space[ssize + j - (int) ai].y / 20; + d += 15 * working_space[ssize + j + (int) ai].y / 20; + d -= 6 * working_space[ssize + j + (int) (2 * ai)].y / 20; + d += working_space[ssize + j + (int) (3 * ai)].y / 20; if (b < d) b = d; if (b < c) b = c; if (b < a) a = b; - working_space[j] = a; + working_space[j].y = a; } else if (smoothing == true){ - a = working_space[ssize + j]; + a = working_space[ssize + j].y; av = 0; men = 0; for (w = j - bw; w <= j + bw; w++){ if ( w >= 0 && w < ssize){ - av += working_space[ssize + w]; + av += working_space[ssize + w].y; men +=1; } } @@ -679,7 +683,7 @@ std::vector background::background(std::vector spectrum, men = 0; for (w = j - i - bw; w <= j - i + bw; w++){ if ( w >= 0 && w < ssize){ - b += working_space[ssize + w]; + b += working_space[ssize + w].y; men +=1; } } @@ -688,7 +692,7 @@ std::vector background::background(std::vector spectrum, men = 0; for (w = j + i - bw; w <= j + i + bw; w++){ if ( w >= 0 && w < ssize){ - c += working_space[ssize + w]; + c += working_space[ssize + w].y; men +=1; } } @@ -698,7 +702,7 @@ std::vector background::background(std::vector spectrum, b4 = 0, men = 0; for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - b4 += working_space[ssize + w]; + b4 += working_space[ssize + w].y; men +=1; } } @@ -706,7 +710,7 @@ std::vector background::background(std::vector spectrum, c4 = 0, men = 0; for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - c4 += working_space[ssize + w]; + c4 += working_space[ssize + w].y; men +=1; } } @@ -714,7 +718,7 @@ std::vector background::background(std::vector spectrum, d4 = 0, men = 0; for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - d4 += working_space[ssize + w]; + d4 += working_space[ssize + w].y; men +=1; } } @@ -722,7 +726,7 @@ std::vector background::background(std::vector spectrum, e4 = 0, men = 0; for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - e4 += working_space[ssize + w]; + e4 += working_space[ssize + w].y; men +=1; } } @@ -732,7 +736,7 @@ std::vector background::background(std::vector spectrum, b6 = 0, men = 0; for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - b6 += working_space[ssize + w]; + b6 += working_space[ssize + w].y; men +=1; } } @@ -740,7 +744,7 @@ std::vector background::background(std::vector spectrum, c6 = 0, men = 0; for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - c6 += working_space[ssize + w]; + c6 += working_space[ssize + w].y; men +=1; } } @@ -748,7 +752,7 @@ std::vector background::background(std::vector spectrum, d6 = 0, men = 0; for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - d6 += working_space[ssize + w]; + d6 += working_space[ssize + w].y; men +=1; } } @@ -756,7 +760,7 @@ std::vector background::background(std::vector spectrum, e6 = 0, men = 0; for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - e6 += working_space[ssize + w]; + e6 += working_space[ssize + w].y; men +=1; } } @@ -764,7 +768,7 @@ std::vector background::background(std::vector spectrum, f6 = 0, men = 0; for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - f6 += working_space[ssize + w]; + f6 += working_space[ssize + w].y; men +=1; } } @@ -772,7 +776,7 @@ std::vector background::background(std::vector spectrum, g6 = 0, men = 0; for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - g6 += working_space[ssize + w]; + g6 += working_space[ssize + w].y; men +=1; } } @@ -784,11 +788,11 @@ std::vector background::background(std::vector spectrum, b = b4; if (b < a) av = b; - working_space[j]=av; + working_space[j].y=av; } } for (j = i; j < ssize - i; j++) - working_space[ssize + j] = working_space[j]; + working_space[ssize + j].y = working_space[j].y; if (direction == kBackIncreasingWindow) i+=1; else if(direction == kBackDecreasingWindow) @@ -800,32 +804,32 @@ std::vector background::background(std::vector spectrum, do{ for (j = i; j < ssize - i; j++) { if (smoothing == false){ - a = working_space[ssize + j]; - b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0; + a = working_space[ssize + j].y; + b = (working_space[ssize + j - i].y + working_space[ssize + j + i].y) / 2.0; c = 0; ai = i / 2; - c -= working_space[ssize + j - (int) (2 * ai)] / 6; - c += 4 * working_space[ssize + j - (int) ai] / 6; - c += 4 * working_space[ssize + j + (int) ai] / 6; - c -= working_space[ssize + j + (int) (2 * ai)] / 6; + c -= working_space[ssize + j - (int) (2 * ai)].y / 6; + c += 4 * working_space[ssize + j - (int) ai].y / 6; + c += 4 * working_space[ssize + j + (int) ai].y / 6; + c -= working_space[ssize + j + (int) (2 * ai)].y / 6; d = 0; ai = i / 3; - d += working_space[ssize + j - (int) (3 * ai)] / 20; - d -= 6 * working_space[ssize + j - (int) (2 * ai)] / 20; - d += 15 * working_space[ssize + j - (int) ai] / 20; - d += 15 * working_space[ssize + j + (int) ai] / 20; - d -= 6 * working_space[ssize + j + (int) (2 * ai)] / 20; - d += working_space[ssize + j + (int) (3 * ai)] / 20; + d += working_space[ssize + j - (int) (3 * ai)].y / 20; + d -= 6 * working_space[ssize + j - (int) (2 * ai)].y / 20; + d += 15 * working_space[ssize + j - (int) ai].y / 20; + d += 15 * working_space[ssize + j + (int) ai].y / 20; + d -= 6 * working_space[ssize + j + (int) (2 * ai)].y / 20; + d += working_space[ssize + j + (int) (3 * ai)].y / 20; e = 0; ai = i / 4; - e -= working_space[ssize + j - (int) (4 * ai)] / 70; - e += 8 * working_space[ssize + j - (int) (3 * ai)] / 70; - e -= 28 * working_space[ssize + j - (int) (2 * ai)] / 70; - e += 56 * working_space[ssize + j - (int) ai] / 70; - e += 56 * working_space[ssize + j + (int) ai] / 70; - e -= 28 * working_space[ssize + j + (int) (2 * ai)] / 70; - e += 8 * working_space[ssize + j + (int) (3 * ai)] / 70; - e -= working_space[ssize + j + (int) (4 * ai)] / 70; + e -= working_space[ssize + j - (int) (4 * ai)].y / 70; + e += 8 * working_space[ssize + j - (int) (3 * ai)].y / 70; + e -= 28 * working_space[ssize + j - (int) (2 * ai)].y / 70; + e += 56 * working_space[ssize + j - (int) ai].y / 70; + e += 56 * working_space[ssize + j + (int) ai].y / 70; + e -= 28 * working_space[ssize + j + (int) (2 * ai)].y / 70; + e += 8 * working_space[ssize + j + (int) (3 * ai)].y / 70; + e -= working_space[ssize + j + (int) (4 * ai)].y / 70; if (b < e) b = e; if (b < d) @@ -834,16 +838,16 @@ std::vector background::background(std::vector spectrum, b = c; if (b < a) a = b; - working_space[j] = a; + working_space[j].y = a; } else if (smoothing == true){ - a = working_space[ssize + j]; + a = working_space[ssize + j].y; av = 0; men = 0; for (w = j - bw; w <= j + bw; w++){ if ( w >= 0 && w < ssize){ - av += working_space[ssize + w]; + av += working_space[ssize + w].y; men +=1; } } @@ -852,7 +856,7 @@ std::vector background::background(std::vector spectrum, men = 0; for (w = j - i - bw; w <= j - i + bw; w++){ if ( w >= 0 && w < ssize){ - b += working_space[ssize + w]; + b += working_space[ssize + w].y; men +=1; } } @@ -861,7 +865,7 @@ std::vector background::background(std::vector spectrum, men = 0; for (w = j + i - bw; w <= j + i + bw; w++){ if ( w >= 0 && w < ssize){ - c += working_space[ssize + w]; + c += working_space[ssize + w].y; men +=1; } } @@ -871,7 +875,7 @@ std::vector background::background(std::vector spectrum, b4 = 0, men = 0; for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - b4 += working_space[ssize + w]; + b4 += working_space[ssize + w].y; men +=1; } } @@ -879,7 +883,7 @@ std::vector background::background(std::vector spectrum, c4 = 0, men = 0; for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - c4 += working_space[ssize + w]; + c4 += working_space[ssize + w].y; men +=1; } } @@ -887,7 +891,7 @@ std::vector background::background(std::vector spectrum, d4 = 0, men = 0; for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - d4 += working_space[ssize + w]; + d4 += working_space[ssize + w].y; men +=1; } } @@ -895,7 +899,7 @@ std::vector background::background(std::vector spectrum, e4 = 0, men = 0; for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - e4 += working_space[ssize + w]; + e4 += working_space[ssize + w].y; men +=1; } } @@ -905,7 +909,7 @@ std::vector background::background(std::vector spectrum, b6 = 0, men = 0; for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - b6 += working_space[ssize + w]; + b6 += working_space[ssize + w].y; men +=1; } } @@ -913,7 +917,7 @@ std::vector background::background(std::vector spectrum, c6 = 0, men = 0; for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - c6 += working_space[ssize + w]; + c6 += working_space[ssize + w].y; men +=1; } } @@ -921,7 +925,7 @@ std::vector background::background(std::vector spectrum, d6 = 0, men = 0; for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - d6 += working_space[ssize + w]; + d6 += working_space[ssize + w].y; men +=1; } } @@ -929,7 +933,7 @@ std::vector background::background(std::vector spectrum, e6 = 0, men = 0; for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - e6 += working_space[ssize + w]; + e6 += working_space[ssize + w].y; men +=1; } } @@ -937,7 +941,7 @@ std::vector background::background(std::vector spectrum, f6 = 0, men = 0; for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - f6 += working_space[ssize + w]; + f6 += working_space[ssize + w].y; men +=1; } } @@ -945,7 +949,7 @@ std::vector background::background(std::vector spectrum, g6 = 0, men = 0; for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - g6 += working_space[ssize + w]; + g6 += working_space[ssize + w].y; men +=1; } } @@ -955,7 +959,7 @@ std::vector background::background(std::vector spectrum, b8 = 0, men = 0; for (w = j - (int)(4 * ai) - bw; w <= j - (int)(4 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - b8 += working_space[ssize + w]; + b8 += working_space[ssize + w].y; men +=1; } } @@ -963,7 +967,7 @@ std::vector background::background(std::vector spectrum, c8 = 0, men = 0; for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - c8 += working_space[ssize + w]; + c8 += working_space[ssize + w].y; men +=1; } } @@ -971,7 +975,7 @@ std::vector background::background(std::vector spectrum, d8 = 0, men = 0; for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - d8 += working_space[ssize + w]; + d8 += working_space[ssize + w].y; men +=1; } } @@ -979,7 +983,7 @@ std::vector background::background(std::vector spectrum, e8 = 0, men = 0; for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - e8 += working_space[ssize + w]; + e8 += working_space[ssize + w].y; men +=1; } } @@ -987,7 +991,7 @@ std::vector background::background(std::vector spectrum, f8 = 0, men = 0; for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){ if (w >= 0 && w < ssize){ - f8 += working_space[ssize + w]; + f8 += working_space[ssize + w].y; men +=1; } } @@ -995,7 +999,7 @@ std::vector background::background(std::vector spectrum, g8 = 0, men = 0; for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - g8 += working_space[ssize + w]; + g8 += working_space[ssize + w].y; men +=1; } } @@ -1003,7 +1007,7 @@ std::vector background::background(std::vector spectrum, h8 = 0, men = 0; for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - h8 += working_space[ssize + w]; + h8 += working_space[ssize + w].y; men +=1; } } @@ -1011,7 +1015,7 @@ std::vector background::background(std::vector spectrum, i8 = 0, men = 0; for (w = j + (int)(4 * ai) - bw; w <= j + (int)(4 * ai) + bw; w++){ if (w >= 0 && w < ssize){ - i8 += working_space[ssize + w]; + i8 += working_space[ssize + w].y; men +=1; } } @@ -1025,11 +1029,11 @@ std::vector background::background(std::vector spectrum, b = b4; if (b < a) av = b; - working_space[j]=av; + working_space[j].y = av; } } for (j = i; j < ssize - i; j++) - working_space[ssize + j] = working_space[j]; + working_space[ssize + j].y = working_space[j].y; if (direction == kBackIncreasingWindow) i += 1; else if(direction == kBackDecreasingWindow) @@ -1040,15 +1044,15 @@ std::vector background::background(std::vector spectrum, if (compton == true) { for (i = 0, b2 = 0; i < ssize; i++){ b1 = b2; - a = working_space[i], b = spectrum[i]; + a = working_space[i].y, b = spectrum[i].y; j = i; if (fabs(a - b) >= 1) { b1 = i - 1; if (b1 < 0) b1 = 0; - yb1 = working_space[b1]; + yb1 = working_space[b1].y; for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){ - a = working_space[b2], b = spectrum[b2]; + a = working_space[b2].y, b = spectrum[b2].y; c = c + b - yb1; if (fabs(a - b) < 1) { priz = 1; @@ -1057,35 +1061,35 @@ std::vector background::background(std::vector spectrum, } if (b2 == ssize) b2 -= 1; - yb2 = working_space[b2]; + yb2 = working_space[b2].y; if (yb1 <= yb2){ for (j = b1, c = 0; j <= b2; j++){ - b = spectrum[j]; + b = spectrum[j].y; c = c + b - yb1; } if (c > 1){ c = (yb2 - yb1) / c; for (j = b1, d = 0; j <= b2 && j < ssize; j++){ - b = spectrum[j]; + b = spectrum[j].y; d = d + b - yb1; a = c * d + yb1; - working_space[ssize + j] = a; + working_space[ssize + j].y = a; } } } else{ for (j = b2, c = 0; j >= b1; j--){ - b = spectrum[j]; + b = spectrum[j].y; c = c + b - yb2; } if (c > 1){ c = (yb1 - yb2) / c; for (j = b2, d = 0;j >= b1 && j >= 0; j--){ - b = spectrum[j]; + b = spectrum[j].y; d = d + b - yb2; a = c * d + yb2; - working_space[ssize + j] = a; + working_space[ssize + j].y = a; } } } diff --git a/fityk/root/background.hpp b/fityk/root/background.hpp index efe016b8..44862f1a 100644 --- a/fityk/root/background.hpp +++ b/fityk/root/background.hpp @@ -15,31 +15,35 @@ #ifndef __BACKGROUND_HPP__ #define __BACKGROUND_HPP__ -namespace background { +#include + +#include "fityk.h" + +namespace ROOT { enum { - kBackIncreasingWindow =0, - kBackDecreasingWindow =1, - kBackOrder2 =2, - kBackOrder4 =4, - kBackOrder6 =6, - kBackOrder8 =8, - kBackSmoothing3 =3, - kBackSmoothing5 =5, - kBackSmoothing7 =7, - kBackSmoothing9 =9, - kBackSmoothing11 =11, - kBackSmoothing13 =13, - kBackSmoothing15 =15 + kBackIncreasingWindow = 1, + kBackDecreasingWindow = -1, + kBackOrder2 = 2, + kBackOrder4 = 4, + kBackOrder6 = 6, + kBackOrder8 = 8, + kBackSmoothing3 = 3, + kBackSmoothing5 = 5, + kBackSmoothing7 = 7, + kBackSmoothing9 = 9, + kBackSmoothing11 = 11, + kBackSmoothing13 = 13, + kBackSmoothing15 = 15 }; - std::vector background(std::vector spectrum, - int numberIterations, - int direction, - int filterOrder, - bool smoothing, - int smoothWindow, - bool compton); + std::vector background(const std::vector spectrum, + int numberIterations, + int direction, + int filterOrder, + bool smoothing, + int smoothWindow, + bool compton); } #endif From 8db857da5a38bca2e642416bb0c23d8b1955ec0b Mon Sep 17 00:00:00 2001 From: Cristiano Fontana Date: Fri, 28 Dec 2018 08:53:17 +0100 Subject: [PATCH 3/5] Added the ROOT source codes to the default compilation --- fityk/Makefile.am | 2 ++ 1 file changed, 2 insertions(+) diff --git a/fityk/Makefile.am b/fityk/Makefile.am index 1eb1ab43..b9ba04bf 100644 --- a/fityk/Makefile.am +++ b/fityk/Makefile.am @@ -13,6 +13,7 @@ libfityk_la_SOURCES = logic.cpp view.cpp lexer.cpp eparser.cpp cparser.cpp \ runner.cpp info.cpp common.cpp data.cpp var.cpp mgr.cpp \ tplate.cpp func.cpp udf.cpp bfunc.cpp f_fcjasym.cpp ast.cpp \ vm.cpp transform.cpp settings.cpp ui.cpp ui_api.cpp \ + root/background.cpp \ luabridge.cpp GAfit.cpp LMfit.cpp guess.cpp NMfit.cpp \ model.cpp fit.cpp voigt.cpp numfuncs.cpp fityk.cpp \ \ @@ -20,6 +21,7 @@ libfityk_la_SOURCES = logic.cpp view.cpp lexer.cpp eparser.cpp cparser.cpp \ runner.h info.h common.h data.h var.h mgr.h \ tplate.h func.h udf.h bfunc.h f_fcjasym.h ast.h \ vm.h transform.h settings.h ui.h luabridge.h \ + root/background.hpp \ GAfit.h LMfit.h guess.h NMfit.h \ model.h fit.h voigt.h numfuncs.h \ swig/fityk_lua.cpp swig/luarun.h \ From a54a6f9f224e43d83defffd35fefa31a8f403342 Mon Sep 17 00:00:00 2001 From: Cristiano Fontana Date: Fri, 28 Dec 2018 09:24:14 +0100 Subject: [PATCH 4/5] Added the SNIP background calculation as a function on datasets --- fityk/eparser.cpp | 5 ++ fityk/transform.cpp | 115 ++++++++++++++++++++++++++++++++++++++++++++ fityk/vm.cpp | 2 +- fityk/vm.h | 1 + 4 files changed, 122 insertions(+), 1 deletion(-) diff --git a/fityk/eparser.cpp b/fityk/eparser.cpp index b21a7db2..99d9786e 100644 --- a/fityk/eparser.cpp +++ b/fityk/eparser.cpp @@ -79,6 +79,7 @@ const char* function_name(int op) case OP_DT_SUM_SAME_X: return "sum_same_x"; case OP_DT_AVG_SAME_X: return "avg_same_x"; case OP_DT_SHIRLEY_BG: return "shirley_bg"; + case OP_DT_SNIP_BG: return "snip_bg"; // 2-args functions case OP_MOD: return "mod"; case OP_MIN2: return "min2"; @@ -138,6 +139,8 @@ int get_function_narg(int op) case OP_RANDNORM: case OP_RANDU: return 2; + case OP_DT_SNIP_BG: + return 5; // Fityk functions case OP_FUNC: case OP_SUM_F: @@ -819,6 +822,8 @@ void ExpressionParser::parse_expr(Lexer& lex, int default_ds, put_function(OP_DT_AVG_SAME_X); else if (mode == kDatasetTrMode && word == "shirley_bg") put_function(OP_DT_SHIRLEY_BG); + else if (mode == kDatasetTrMode && word == "snip_bg") + put_function(OP_DT_SNIP_BG); else lex.throw_syntax_error("unknown function: " + word); diff --git a/fityk/transform.cpp b/fityk/transform.cpp index 0e675713..ce4f0ce7 100644 --- a/fityk/transform.cpp +++ b/fityk/transform.cpp @@ -6,6 +6,8 @@ #include "transform.h" #include "logic.h" #include "data.h" +#include "root/background.hpp" +#include using namespace std; @@ -95,6 +97,94 @@ void shirley_bg(vector &pp) pp[i].y = B[i]; } +template +struct vector_slice { + typename vector::iterator begin; + typename vector::iterator end; +}; + +// Determines the first active slice of the data points. +// It reads from the vector slice determined by the given +// iterators. +// It returns the iterators to the begin of the active +// slice (first active datapoint) and to the end of the +// slice (datapoint *after* the last active one). +struct vector_slice first_active_slice(const vector &pp, + vector::iterator begin, + vector::iterator end) +{ + vector::iterator active_start = begin; + + for (vector::iterator it = begin; it != end; it++) { + if (it->is_active) { + active_start = it; + break; + } else { + active_start += 1; + } + } + + vector::iterator active_end = active_start; + + for (vector::iterator it = active_start; it != end; it++) { + if (it->is_active) { + active_end = it + 1; + } else { + break; + } + } + + return {active_start, active_end}; +} + +// 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. +void snip_bg_slice(vector &pp, + vector::iterator begin, + vector::iterator end, + int window_width, + int direction, + int filter_order, + bool smoothing, + bool smooth_window, + bool estimate_compton) +{ + vector_slice slice = first_active_slice(pp, begin, end); + + if (slice.begin == slice.end) { + return; + } else { + vector bg_input(slice.begin, slice.end); + + 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 = slice.begin; it != slice.end; it++) { + const unsigned int index = it - slice.begin; + + it->y = bg_output[index].y; + } + } + + snip_bg_slice(pp, slice.end, pp.end(), + window_width, + direction, + filter_order, + smoothing, + smooth_window, + estimate_compton); + } +} + } // anonymous namespace namespace fityk { @@ -202,6 +292,31 @@ void run_data_transform(const DataKeeper& dk, const VMData& vm, Data* data_out) shirley_bg(stackPtr->points); break; + case OP_DT_SNIP_BG: + stackPtr -= 4; + + if ((stackPtr)->is_num) { + throw ExecuteError(op2str(*i) + " is defined only for @n"); + } else { + const int window_width = (stackPtr+1)->num; + const int direction = (stackPtr+2)->num >= 0 ? ROOT::kBackIncreasingWindow : ROOT::kBackDecreasingWindow; + const int filter_order = (stackPtr+3)->num; + const bool smoothing = false; + const int smooth_window = ROOT::kBackSmoothing3; + const bool estimate_compton = (stackPtr+4)->num > 0 ? true : false; + + snip_bg_slice(stackPtr->points, + stackPtr->points.begin(), + stackPtr->points.end(), + window_width, + direction, + filter_order, + smoothing, + smooth_window, + estimate_compton); + } + break; + case OP_AND: // do nothing break; diff --git a/fityk/vm.cpp b/fityk/vm.cpp index 6e19b734..a6e467da 100644 --- a/fityk/vm.cpp +++ b/fityk/vm.cpp @@ -109,7 +109,7 @@ string op2str(int op) OP_(FUNC) OP_(SUM_F) OP_(SUM_Z) OP_(NUMAREA) OP_(FINDX) OP_(FIND_EXTR) OP_(TILDE) - OP_(DATASET) OP_(DT_SUM_SAME_X) OP_(DT_AVG_SAME_X) OP_(DT_SHIRLEY_BG) + OP_(DATASET) OP_(DT_SUM_SAME_X) OP_(DT_AVG_SAME_X) OP_(DT_SHIRLEY_BG) OP_(DT_SNIP_BG) OP_(OPEN_ROUND) OP_(OPEN_SQUARE) } return S(op); // unreachable (if all OPs are listed above) diff --git a/fityk/vm.h b/fityk/vm.h index 0aa3437c..7532b011 100644 --- a/fityk/vm.h +++ b/fityk/vm.h @@ -101,6 +101,7 @@ enum Op OP_DT_SUM_SAME_X, OP_DT_AVG_SAME_X, OP_DT_SHIRLEY_BG, + OP_DT_SNIP_BG, // these two are not VM operators, but are handy to have here, // they and are used in implementation of shunting yard algorithm From 3a75228cca674f83669fa97667969afd4d3dbf1a Mon Sep 17 00:00:00 2001 From: Cristiano Fontana Date: Fri, 28 Dec 2018 09:24:34 +0100 Subject: [PATCH 5/5] Added some documentation for the SNIP background --- doc/data.rst | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/doc/data.rst b/doc/data.rst index 053ed1e4..35ebd947 100644 --- a/doc/data.rst +++ b/doc/data.rst @@ -637,11 +637,24 @@ and functions: Calculates Shirley background (useful in X-ray photoelectron spectroscopy). +``snip_bg(@n, window_width, direction, filter_order, estimate_compton)`` + Calculates the spectrum background applying a Sensitive Nonlinear Iterative Peak clipping algorithm (SNIP) (useful in gamma-ray spectroscopy). + This SNIP algorithm is taken from the `TSpectrum `_ class of `ROOT `_. + It requires the parameters: + 1. ``window_width``: maximal width of clipping window; + 2. ``direction``: direction of change of clipping window, if it is a positive number the window will be increasing; + 3. ``filter_order``: order of clipping filter, possible values are 2, 4, 6 or 8; + 4. ``estimate_compton``: if a positive number the algorithm will try to estimate the Compton edge. + Active points are modified according to the SNIP algorithm. + Inactive points are left as they are, assuming that they are already considered background. + Subtracting the calculated background from the initial spectrum gives zero in all inactive points. + Examples:: @+ = @0 # duplicate the dataset @+ = @0 and @1 # create a new dataset from @0 and @1 @0 = @0 - shirley_bg(@0) # remove Shirley background + @0 = @0 - snip_bg(@0, 30, 1, 2, 0) # remove gamma background @0 = @0 - @1 # subtract @1 from @0 @0 = @0 - 0.28*@1 # subtract scaled dataset @1 from @0