//BEGIN_HTML <!--
/* -->
<A href="ftp://root.cern.ch/root/Spectrum.doc">Spectrum.doc</A><br>
<A href="ftp://root.cern.ch/root/SpectrumDec.ps.gz">SpectrumDec.ps.gz</A><br>
<A href="ftp://root.cern.ch/root/SpectrumSrc.ps.gz">SpectrumSrc.ps.gz</A><br>
<A href="ftp://root.cern.ch/root/SpectrumBck.ps.gz">SpectrumBck.ps.gz</A><br>
<!--*/
// -->END_HTML
#include "TSpectrum.h"
#include "TPolyMarker.h"
#include "TVirtualPad.h"
#include "TMath.h"
Int_t TSpectrum::fgIterations = 3;
Int_t TSpectrum::fgAverageWindow = 3;
#define PEAK_WINDOW 1024
ClassImp(TSpectrum)
TSpectrum::TSpectrum() :TNamed("Spectrum", "Miroslav Morhac peak finder")
{
Int_t n = 100;
fMaxPeaks = n;
fPosition = new Float_t[n];
fPositionX = new Float_t[n];
fPositionY = new Float_t[n];
fResolution = 1;
fHistogram = 0;
fNPeaks = 0;
}
TSpectrum::TSpectrum(Int_t maxpositions, Float_t resolution) :TNamed("Spectrum", "Miroslav Morhac peak finder")
{
Int_t n = maxpositions;
if (n <= 0) n = 1;
fMaxPeaks = n;
fPosition = new Float_t[n];
fPositionX = new Float_t[n];
fPositionY = new Float_t[n];
fHistogram = 0;
fNPeaks = 0;
SetResolution(resolution);
}
TSpectrum::~TSpectrum()
{
delete [] fPosition;
delete [] fPositionX;
delete [] fPositionY;
delete fHistogram;
}
void TSpectrum::SetAverageWindow(Int_t w)
{
fgAverageWindow = w;
}
void TSpectrum::SetDeconIterations(Int_t n)
{
fgIterations = n;
}
TH1 *TSpectrum::Background(const TH1 * h, int numberIterations, Option_t * option)
{
if (h == 0) return 0;
Int_t dimension = h->GetDimension();
if (dimension > 1) {
Error("Search", "Only implemented for 1-d histograms");
return 0;
}
TString opt = option;
opt.ToLower();
Int_t direction = kBackDecreasingWindow;
if (opt.Contains("backincreasingwindow")) direction = kBackIncreasingWindow;
Int_t filterOrder = kBackOrder2;
if (opt.Contains("backorder4")) filterOrder = kBackOrder4;
if (opt.Contains("backorder6")) filterOrder = kBackOrder6;
if (opt.Contains("backorder8")) filterOrder = kBackOrder8;
Bool_t smoothing = kTRUE;
if (opt.Contains("nosmoothing")) smoothing = kFALSE;
Int_t smoothWindow = kBackSmoothing3;
if (opt.Contains("backsmoothing5")) smoothWindow = kBackSmoothing5;
if (opt.Contains("backsmoothing7")) smoothWindow = kBackSmoothing7;
if (opt.Contains("backsmoothing9")) smoothWindow = kBackSmoothing9;
if (opt.Contains("backsmoothing11")) smoothWindow = kBackSmoothing11;
if (opt.Contains("backsmoothing13")) smoothWindow = kBackSmoothing13;
if (opt.Contains("backsmoothing15")) smoothWindow = kBackSmoothing15;
Bool_t compton = kFALSE;
if (opt.Contains("compton")) compton = kTRUE;
Int_t first = h->GetXaxis()->GetFirst();
Int_t last = h->GetXaxis()->GetLast();
Int_t size = last-first+1;
Int_t i;
Float_t * source = new float[size];
for (i = 0; i < size; i++) source[i] = h->GetBinContent(i + first);
Background(source,size,numberIterations, direction, filterOrder,smoothing,smoothWindow,compton);
Int_t nch = strlen(h->GetName());
char *hbname = new char[nch+20];
sprintf(hbname,"%s_background",h->GetName());
TH1 *hb = (TH1*)h->Clone(hbname);
hb->Reset();
hb->GetListOfFunctions()->Delete();
hb->SetLineColor(2);
for (i=0; i< size; i++) hb->SetBinContent(i+first,source[i]);
hb->SetEntries(size);
if (opt.Contains("same")) {
if (gPad) delete gPad->GetPrimitive(hbname);
hb->Draw("same");
}
delete [] source;
delete [] hbname;
return hb;
}
void TSpectrum::Print(Option_t *) const
{
printf("\nNumber of positions = %d\n",fNPeaks);
for (Int_t i=0;i<fNPeaks;i++) {
printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
}
}
Int_t TSpectrum::Search(const TH1 * hin, Double_t sigma, Option_t * option, Double_t threshold)
{
if (hin == 0) return 0;
Int_t dimension = hin->GetDimension();
if (dimension > 2) {
Error("Search", "Only implemented for 1-d and 2-d histograms");
return 0;
}
if (threshold <=0 || threshold >= 1) {
Warning("Search","threshold must 0<threshold<1, threshol=0.05 assumed");
threshold = 0.05;
}
TString opt = option;
opt.ToLower();
Bool_t background = kTRUE;
if (opt.Contains("nobackground")) {
background = kFALSE;
opt.ReplaceAll("nobackground","");
}
Bool_t markov = kTRUE;
if (opt.Contains("nomarkov")) {
markov = kFALSE;
opt.ReplaceAll("nomarkov","");
}
if (dimension == 1) {
Int_t first = hin->GetXaxis()->GetFirst();
Int_t last = hin->GetXaxis()->GetLast();
Int_t size = last-first+1;
Int_t i, bin, npeaks;
Float_t * source = new float[size];
Float_t * dest = new float[size];
for (i = 0; i < size; i++) source[i] = hin->GetBinContent(i + first);
if (sigma <= 1) {
sigma = size/fMaxPeaks;
if (sigma < 1) sigma = 1;
if (sigma > 8) sigma = 8;
}
npeaks = SearchHighRes(source, dest, size, sigma, 100*threshold, background, fgIterations, markov, fgAverageWindow);
for (i = 0; i < npeaks; i++) {
bin = first + Int_t(fPositionX[i] + 0.5);
fPositionX[i] = hin->GetBinCenter(bin);
fPositionY[i] = hin->GetBinContent(bin);
}
delete [] source;
delete [] dest;
if (opt.Contains("goff"))
return npeaks;
if (!npeaks) return 0;
TPolyMarker * pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
if (pm) {
hin->GetListOfFunctions()->Remove(pm);
delete pm;
}
pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
hin->GetListOfFunctions()->Add(pm);
pm->SetMarkerStyle(23);
pm->SetMarkerColor(kRed);
pm->SetMarkerSize(1.3);
opt.ReplaceAll(" ","");
opt.ReplaceAll(",","");
((TH1*)hin)->Draw(opt.Data());
return npeaks;
}
return 0;
}
void TSpectrum::SetResolution(Float_t resolution)
{
if (resolution > 1)
fResolution = resolution;
else
fResolution = 1;
}
const char *TSpectrum::Background(float *spectrum, int ssize,
int numberIterations,
int direction, int filterOrder,
bool smoothing,int smoothWindow,
bool compton)
{
//Begin_Html <!--
/* -->
<div class=Section1>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:20.0pt'>Background
estimation</span></b></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>Goal:
Separation of useful information (peaks) from useless information (background)</span></i><span
style='font-size:18.0pt'> </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-indent:-18.0pt'><span
style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>method is based on Sensitive
Nonlinear Iterative Peak (SNIP) clipping algorithm</span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-indent:-18.0pt'><span
style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'> new value in the channel i is
calculated</span></p>
<p class=MsoNormal>
<table cellpadding=0 cellspacing=0 align=left>
<tr>
<td width=53 height=23></td>
</tr>
<tr>
<td></td>
<td><img width=486 height=72 src="gif/TSpectrum_Background.gif"></td>
</tr>
</table>
<span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<br clear=ALL>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>where
p = 1, 2,
, numberIterations. In fact it represents second order difference
filter (-1,2,-1).</span></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'>const
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#char" target="_parent">char</a>*
<a
href="http://root.cern.ch/root/html/src/TSpectrum.cxx.html#TSpectrum:Background1General"
target="_parent">Background</a>(<a
href="http://root.cern.ch/root/html/ListOtransTypes.html#float" target="_parent">float</a>
*spectrum, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int"
target="_parent">int</a> ssize, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#int" target="_parent">int</a>
numberIterations, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int"
target="_parent">int</a> direction, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#int" target="_parent">int</a>
filterOrder, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#bool"
target="_parent">bool</a> smoothing, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#int" target="_parent">int</a>
smoothingWindow, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#bool"
target="_parent">bool</a> compton) </span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This
function calculates background spectrum from the source spectrum. The result
is placed in the vector pointed by spectrum pointer. 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.</span><span style='font-size:18.0pt'> </span><span
style='font-size:16.0pt'>On successful completion it returns 0. On error it
returns pointer to the string describing error.</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>spectrum</span></b>-pointer
to the vector of source spectrum </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>ssize</span></b>-length
of the spectrum vector </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>numberIterations</span></b>-maximal
width of clipping window, </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>direction</span></b>-
direction of change of clipping window </p>
<p class=MsoNormal> - possible
values=kBackIncreasingWindow </p>
<p class=MsoNormal> kBackDecreasingWindow
</p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>filterOrder</span></b>-order
of clipping filter, </p>
<p class=MsoNormal> -possible
values=kBackOrder2 </p>
<p class=MsoNormal> kBackOrder4
</p>
<p class=MsoNormal> kBackOrder6
</p>
<p class=MsoNormal> kBackOrder8
</p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>smoothing</span></b>-
logical variable whether the smoothing operation in the estimation of </p>
<p class=MsoNormal> background will be included </p>
<p class=MsoNormal> - possible
values=kFALSE </p>
<p class=MsoNormal> kTRUE
</p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>smoothWindow</span></b>-width
of smoothing window, </p>
<p class=MsoNormal> -possible
values=kBackSmoothing3 </p>
<p class=MsoNormal> kBackSmoothing5
</p>
<p class=MsoNormal> kBackSmoothing7
</p>
<p class=MsoNormal> kBackSmoothing9
</p>
<p class=MsoNormal> kBackSmoothing11
</p>
<p class=MsoNormal> kBackSmoothing13
</p>
<p class=MsoNormal> kBackSmoothing15
</p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>compton</span></b>-
logical variable whether the estimation of Compton edge will be
included </p>
<p class=MsoNormal> - possible
values=kFALSE </p>
<p class=MsoNormal> kTRUE
</p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><i><span style='font-size:18.0pt'>References:</span></i></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[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.</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[2]
</span><span lang=SK style='font-size:16.0pt'> M. Morháč, J. Kliman, V.
Matouek, M. Veselskŭ, I. Turzo</span><span style='font-size:16.0pt'>.:
Background elimination methods for multidimensional gamma-ray spectra. NIM,
A401 (1997) 113-132.</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[3]
D. D. Burgess, R. J. Tervo: Background estimation for gamma-ray spectroscopy.
NIM 214 (1983), 431-434. </span></p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section2>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Example 1 script
Background_incr.c :</span></i></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Background_incr.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
1 Example of the estimation of background for number of iterations=6. Original
spectrum is shown in black color, estimated background in red color.</span></b></p>
<p class=MsoNormal><b><span style='font-size:14.0pt;color:green'> </span></b></p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:green'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate the background estimator (class
TSpectrum).</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Background_incr.C</p>
<p class=MsoNormal>#include <TSpectrum> </p>
<p class=MsoNormal>void Background_incr() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 256;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> TH1F *back = new
TH1F("back","",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d = new
TH1F("d","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> back=(TH1F*) f->Get("back1;1");</p>
<p class=MsoNormal> TCanvas *Background = gROOT->GetListOfCanvases()->FindObject("Background");</p>
<p class=MsoNormal> if (!Background) Background = new
TCanvas("Background","Estimation of background with increasing
window",10,10,1000,700);</p>
<p class=MsoNormal> back->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i
+ 1); </p>
<p class=MsoNormal> s->Background(source,nbins,6,kBackIncreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d->SetLineColor(kRed);</p>
<p class=MsoNormal> d->Draw("SAME L"); </p>
<p class=MsoNormal> }</p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section3>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Example 2 script
Background_decr.c :</span></i></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>In
Figure 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. </span></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Background_decr.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
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.</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate the background estimator (class
TSpectrum).</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Background_decr.C</p>
<p class=MsoNormal>#include <TSpectrum> </p>
<p class=MsoNormal>void Background_decr() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 256;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> TH1F *back = new
TH1F("back","",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d = new
TH1F("d","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> back=(TH1F*) f->Get("back1;1");</p>
<p class=MsoNormal> TCanvas *Background =
gROOT->GetListOfCanvases()->FindObject("Background");</p>
<p class=MsoNormal> if (!Background) Background = new
TCanvas("Background","Estimation of background with decreasing
window",10,10,1000,700);</p>
<p class=MsoNormal> back->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=back->GetBinContent(i + 1); </p>
<p class=MsoNormal> s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d->SetLineColor(kRed);</p>
<p class=MsoNormal> d->Draw("SAME L"); </p>
<p class=MsoNormal> }</p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section4>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Example 3 script
Background_width.c :</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>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
Figure 3.</span></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Background_width.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
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.</span></b></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:18.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><i><span style='font-size:18.0pt'>in general one should set this
parameter so that the value 2*numberIterations+1 was greater than the widths
of preserved objects (peaks).</span></i></p>
<p class=MsoNormal style='text-align:justify'> </p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt;
color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate the influence of the clipping
window width on the estimated background</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Background_width.C</p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal>void Background_width() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 256;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> TH1F *h = new
TH1F("h","",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d1 = new
TH1F("d1","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d2 = new
TH1F("d2","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d3 = new
TH1F("d3","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("back1;1"); </p>
<p class=MsoNormal> TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background");</p>
<p class=MsoNormal> if (!background) background = new
TCanvas("background","Influence of clipping window width on the
estimated background",10,10,1000,700);</p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i
+ 1); </p>
<p class=MsoNormal>
s->Background(source,nbins,4,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d1->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d1->SetLineColor(kRed);</p>
<p class=MsoNormal> d1->Draw("SAME L"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal>
s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d2->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d2->SetLineColor(kBlue);</p>
<p class=MsoNormal> d2->Draw("SAME L"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal>
s->Background(source,nbins,8,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d3->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d3->SetLineColor(kGreen);</p>
<p class=MsoNormal> d3->Draw("SAME L"); </p>
<p class=MsoNormal>}</p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section5>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Example 4 script
Background_width2.c :</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-indent:-18.0pt'><span
style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>another example for very complex
spectrum is given in Figure 4.</span></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Background_width2.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
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.</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate the influence of the clipping
window width on the estimated background</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Background_width2.C</p>
<p class=MsoNormal>#include <TSpectrum> </p>
<p class=MsoNormal>void Background_width2() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 4096;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)4096;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> TH1F *h = new
TH1F("h","",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
</p>
<p class=MsoNormal> TH1F *d2 = new
TH1F("d2","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d3 = new
TH1F("d3","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d4 = new
TH1F("d4","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("back2;1"); </p>
<p class=MsoNormal> TCanvas *background =
gROOT->GetListOfCanvases()->FindObject("background");</p>
<p class=MsoNormal> if (!background) background = new
TCanvas("background","Influence of clipping window width on the
estimated background",10,10,1000,700);</p>
<p class=MsoNormal> h->SetAxisRange(0,1000);</p>
<p class=MsoNormal> h->SetMaximum(20000);</p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1); </p>
<p class=MsoNormal>
s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d1->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d1->SetLineColor(kRed);</p>
<p class=MsoNormal> d1->Draw("SAME L"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal> s->Background(source,nbins,20,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d2->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d2->SetLineColor(kBlue);</p>
<p class=MsoNormal> d2->Draw("SAME L"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal>
s->Background(source,nbins,30,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d3->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d3->SetLineColor(kGreen);</p>
<p class=MsoNormal> d3->Draw("SAME L"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal>
s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d4->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d4->SetLineColor(kMagenta);</p>
<p class=MsoNormal> d4->Draw("SAME L"); </p>
<p class=MsoNormal>}</p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section6>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Example 5 script
Background_order.c :</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>second order difference filter
removes linear (quasi-linear) background and preserves symmetrical peaks.</span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>however if the shape of the
background is more complex one can employ higher-order clipping filters (see
example in Figure 5)</span></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Background_order.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
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.</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate the influence of the clipping
filter difference order on the estimated background</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Background_order.C</p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void Background_order() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 4096;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)4096;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> TH1F *h = new
TH1F("h","",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d1 = new
TH1F("d1","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d2 = new
TH1F("d2","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d3 = new
TH1F("d3","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d4 = new
TH1F("d4","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("back2;1");</p>
<p class=MsoNormal> TCanvas *background =
gROOT->GetListOfCanvases()->FindObject("background");</p>
<p class=MsoNormal> if (!background) background = new
TCanvas("background","Influence of clipping filter difference
order on the estimated background",10,10,1000,700);</p>
<p class=MsoNormal> h->SetAxisRange(1220,1460);</p>
<p class=MsoNormal> h->SetMaximum(11000);</p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1); </p>
<p class=MsoNormal> s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d1->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d1->SetLineColor(kRed);</p>
<p class=MsoNormal> d1->Draw("SAME L"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal> s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder4,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d2->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d2->SetLineColor(kBlue);</p>
<p class=MsoNormal> d2->Draw("SAME L"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal>
s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder6,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d3->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d3->SetLineColor(kGreen);</p>
<p class=MsoNormal> d3->Draw("SAME L"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal> s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder8,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d4->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d4->SetLineColor(kMagenta);</p>
<p class=MsoNormal> d4->Draw("SAME L"); </p>
<p class=MsoNormal>}</p>
</div>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section7>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 6 script
Background_smooth.c :</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>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</span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>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. </span></p>
<p class=MsoNormal style='text-indent:5.7pt'><img width=554 height=104
src="gif/TSpectrum_Background_smooth1.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
7 Principle of background estimation algorithm with simultaneous smoothing</span></b></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Background_smooth2.jpg"></p>
<p class=MsoNormal><b><span style='font-size:16.0pt'>Figure 8 Illustration of
non-smoothing (red line) and smoothing algorithm of background estimation (blue
line).</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate the background estimator (class
TSpectrum) including Compton edges.</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Background_smooth.C</p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void Background_smooth() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 4096;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> TH1F *h = new
TH1F("h","",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d2 = new
TH1F("d2","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("back4;1");</p>
<p class=MsoNormal> TCanvas *background =
gROOT->GetListOfCanvases()->FindObject("background");</p>
<p class=MsoNormal> if (!background) background = new
TCanvas("background","Estimation of background with
noise",10,10,1000,700);</p>
<p class=MsoNormal> h->SetAxisRange(3460,3830); </p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1); </p>
<p class=MsoNormal>
s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d1->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d1->SetLineColor(kRed);</p>
<p class=MsoNormal> d1->Draw("SAME L");</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1); </p>
<p class=MsoNormal>
s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kTRUE,kBackSmoothing3,kFALSE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d2->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d2->SetLineColor(kBlue);</p>
<p class=MsoNormal> d2->Draw("SAME L"); </p>
<p class=MsoNormal>}</p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section8>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 8 script
Background_compton.c :</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>sometimes it is necessary to
include also the Compton edges into the estimate of the background. In Figure 8
we present the example of the synthetic spectrum with Compton edges. </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>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).</span></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Background_compton.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
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).</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate the background estimator (class
TSpectrum) including Compton edges.</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Background_compton.C</p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void Background_compton() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 512;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> TH1F *h = new
TH1F("h","",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d1 = new
TH1F("d1","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("back3;1");</p>
<p class=MsoNormal> TCanvas *background =
gROOT->GetListOfCanvases()->FindObject("background");</p>
<p class=MsoNormal> if (!background) background = new
TCanvas("background","Estimation of background with Compton edges under peaks",10,10,1000,700);</p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1); </p>
<p class=MsoNormal> s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder8,kTRUE,kBackSmoothing5,,kTRUE);
</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d1->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d1->SetLineColor(kRed);</p>
<p class=MsoNormal> d1->Draw("SAME L"); </p>
<p class=MsoNormal>}</p>
</div>
<!-- */
// --> End_Html
int i, j, w, bw, b1, b2, priz;
float 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;
if (ssize <= 0)
return "Wrong Parameters";
if (numberIterations < 1)
return "Width of Clipping Window Must Be Positive";
if (ssize < 2 * numberIterations + 1)
return "Too Large Clipping Window";
if (smoothing == kTRUE && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15)
return "Incorrect width of smoothing window";
float *working_space = new float[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 == kFALSE){
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 == kTRUE){
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 == kFALSE){
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 == kTRUE){
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 == kFALSE){
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 == kTRUE){
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 == kFALSE){
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 == kTRUE){
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 == kTRUE) {
for (i = 0, b2 = 0; i < ssize; i++){
b1 = b2;
a = working_space[i], b = spectrum[i];
j = i;
if (TMath::Abs(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 (TMath::Abs(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;
}
}
}
for (j = 0; j < ssize; j++)
spectrum[j] = working_space[ssize + j];
delete[]working_space;
return 0;
}
const char* TSpectrum::SmoothMarkov(float *source, int ssize, int averWindow)
{
//Begin_Html <!--
/* -->
<div class=Section16>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:20.0pt'>Smoothing</span></b></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Goal: Suppression of
statistical fluctuations</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>the algorithm is based on discrete
Markov chain, which has very simple invariant distribution</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal> <sub><img width=551 height=63
src="gif/TSpectrum_Smoothing1.gif"></sub><span style='font-size:16.0pt;
font-family:Arial'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt;
font-family:Arial'> </span><sub><img width=28 height=36
src="gif/TSpectrum_Smoothing2.gif"></sub><span style='font-size:16.0pt;
font-family:Arial'> being defined from the normalization condition </span><sub><img
width=70 height=52 src="gif/TSpectrum_Smoothing3.gif"></sub></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><span style='font-size:16.0pt;font-family:Arial'> n
is the length of the smoothed spectrum and </span></p>
<p class=MsoNormal>
<table cellpadding=0 cellspacing=0 align=left>
<tr>
<td width=57 height=15></td>
</tr>
<tr>
<td></td>
<td><img width=258 height=60 src="gif/TSpectrum_Smoothing4.gif"></td>
</tr>
</table>
<span style='font-size:16.0pt;font-family:Arial'> </span></p>
<p class=MsoNormal><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<br clear=ALL>
<p class=MsoNormal style='margin-left:34.2pt;text-align:justify'><span
style='font-size:16.0pt;font-family:Arial'>is the probability of the change of
the peak position from channel i to the channel i+1. </span> <sub><img
width=28 height=36 src="gif/TSpectrum_Smoothing5.gif"></sub><span
style='font-size:16.0pt'>is the normalization constant so that </span><sub><img
width=133 height=34 src="gif/TSpectrum_Smoothing6.gif"></sub> <span
style='font-size:16.0pt'>and m is a width of smoothing window. </span></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'>const
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#char" target="_parent">char</a>*
SmoothMarkov(<a href="http://root.cern.ch/root/html/ListOtransTypes.html#float"
target="_parent">float</a> *spectrum, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#int" target="_parent">int</a>
ssize, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int"
target="_parent">int</a> averWindow) </span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This
function calculates smoothed spectrum from the source spectrum based on Markov
chain method. The result is placed in the vector pointed by source pointer. On
successful completion it returns 0. On error it returns pointer to the string
describing error.</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>spectrum</span></b>-pointer
to the vector of source spectrum </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>ssize</span></b>-length
of the spectrum vector </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>averWindow</span></b>-width
of averaging smoothing window </p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><i><span style='font-size:18.0pt'>Reference:</span></i></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[1]
Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
(1996), 45<b>1.</b> </span></p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section17>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 14 script Smoothing.c
:</span></i></p>
<p class=MsoNormal><span style='font-size:16.0pt'> <img width=296 height=182
src="gif/TSpectrum_Smoothing1.jpg"><img width=296 height=182
src="gif/TSpectrum_Smoothing2.jpg"></span></p>
<p class=MsoNormal><b><span style='font-size:16.0pt'>Fig. 23 Original noisy spectrum</span></b><b><span
style='font-size:14.0pt'> </span></b><b><span style='font-size:16.0pt'>Fig.
24 Smoothed spectrum m=3</span></b></p>
<p class=MsoNormal><b><span style='font-size:14.0pt'><img width=299 height=184
src="gif/TSpectrum_Smoothing3.jpg"><img width=299 height=184
src="gif/TSpectrum_Smoothing4.jpg"></span></b></p>
<p class=MsoNormal><b><span style='font-size:16.0pt'>Fig. 25 Smoothed spectrum
m=7 Fig.26 Smoothed spectrum m=10</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate smoothing using Markov algorithm
(class TSpectrum).</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Smoothing.C</p>
<p class=MsoNormal> </p>
<p class=MsoNormal>//#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void Smoothing() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 1024;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> TH1F *h = new TH1F("h","Smoothed spectrum
for m=3",nbins,xmin,xmax);</p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("smooth1;1"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1); </p>
<p class=MsoNormal> TCanvas *Smooth1 =
gROOT->GetListOfCanvases()->FindObject("Smooth1");</p>
<p class=MsoNormal> if (!Smooth1) Smooth1 = new
TCanvas("Smooth1","Smooth1",10,10,1000,700);</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> s->SmoothMarkov(source,1024,3); //3, 7, 10</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) h->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> h->SetAxisRange(330,880);</p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal>}</p>
</div>
<!-- */
// --> End_Html
int xmin, xmax, i, l;
float a, b, maxch;
float nom, nip, nim, sp, sm, area = 0;
if(averWindow <= 0)
return "Averaging Window must be positive";
float *working_space = new float[ssize];
xmin = 0,xmax = ssize - 1;
for(i = 0, maxch = 0; i < ssize; i++){
working_space[i]=0;
if(maxch < source[i])
maxch = source[i];
area += source[i];
}
if(maxch == 0)
return 0 ;
nom = 1;
working_space[xmin] = 1;
for(i = xmin; i < xmax; i++){
nip = source[i] / maxch;
nim = source[i + 1] / maxch;
sp = 0,sm = 0;
for(l = 1; l <= averWindow; l++){
if((i + l) > xmax)
a = source[xmax] / maxch;
else
a = source[i + l] / maxch;
b = a - nip;
if(a + nip <= 0)
a = 1;
else
a = TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
sp = sp + b;
if((i - l + 1) < xmin)
a = source[xmin] / maxch;
else
a = source[i - l + 1] / maxch;
b = a - nim;
if(a + nim <= 0)
a = 1;
else
a = TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
sm = sm + b;
}
a = sp / sm;
a = working_space[i + 1] = working_space[i] * a;
nom = nom + a;
}
for(i = xmin; i <= xmax; i++){
working_space[i] = working_space[i] / nom;
}
for(i = 0; i < ssize; i++)
source[i] = working_space[i] * area;
delete[]working_space;
return 0;
}
const char *TSpectrum::Deconvolution(float *source, const float *response,
int ssize, int numberIterations,
int numberRepetitions, double boost )
{
//Begin_Html <!--
/* -->
<div class=Section9>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:20.0pt'>Deconvolution</span></b></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>Goal:
Improvement of the resolution in spectra, decomposition of multiplets</span></i></p>
<p class=MsoNormal><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:16.0pt'>Mathematical formulation of
the convolution system is</span></p>
<p class=MsoNormal style='margin-left:18.0pt'>
<table cellpadding=0 cellspacing=0 align=left>
<tr>
<td width=0 height=17></td>
</tr>
<tr>
<td></td>
<td><img width=585 height=84 src="gif/TSpectrum_Deconvolution1.gif"></td>
</tr>
</table>
<span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<br clear=ALL>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>where
h(i) is the impulse response function, x, y are input and output vectors, respectively,
N is the length of x and h vectors. In matrix form we have</span></p>
<p class=MsoNormal style='text-align:justify'>
<table cellpadding=0 cellspacing=0 align=left>
<tr>
<td width=4 height=8></td>
</tr>
<tr>
<td></td>
<td><img width=597 height=360 src="gif/TSpectrum_Deconvolution2.gif"></td>
</tr>
</table>
<span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<br clear=ALL>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>let us assume that we know the
response and the output vector (spectrum) of the above given system. </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>the deconvolution represents
solution of the overdetermined system of linear equations, i.e., the
calculation of the vector <b>x.</b></span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>from numerical stability point of
view the operation of deconvolution is extremely critical (ill-posed problem)
as well as time consuming operation. </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>the Gold deconvolution algorithm
proves to work very well, other methods (Fourier, VanCittert etc) oscillate. </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>it is suitable to process positive
definite data (e.g. histograms). </span></p>
<p class=MsoNormal><b><i><span style='font-size:16.0pt'> </span></i></b></p>
<p class=MsoNormal><b><i><span style='font-size:16.0pt'>Gold deconvolution
algorithm</span></i></b></p>
<p class=MsoNormal>
<table cellpadding=0 cellspacing=0 align=left>
<tr>
<td width=46 height=21></td>
</tr>
<tr>
<td></td>
<td><img width=551 height=233 src="gif/TSpectrum_Deconvolution3.gif"></td>
</tr>
</table>
<span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<br clear=ALL>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt;
font-family:Arial'>where L is given number of iterations (numberIterations
parameter).</span></p>
<p class=MsoNormal><b><i><span style='font-size:16.0pt'> </span></i></b></p>
<p class=MsoNormal><span style='position:absolute;z-index:4;margin-left:247px;
margin-top:17px;width:144px;height:36px'><img width=144 height=36
src="gif/TSpectrum_Deconvolution4.gif"></span><b><i><span style='font-size:
16.0pt'>Boosted deconvolution</span></i></b></p>
<p class=MsoNormal style='margin-left:36.0pt;text-indent:-18.0pt'><span
style='font-size:16.0pt'>1.<span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>Set the initial solution </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'>2.<span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>Set required number of repetitions
R and iterations L</span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-indent:-18.0pt'><span
style='font-size:16.0pt'>3.<span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>Set r = 1.</span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-indent:-18.0pt'><span
style='font-size:16.0pt'>4.<span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>Using Gold deconvolution algorithm
for k=1,2,...,L find </span><sub><img width=30 height=24
src="gif/TSpectrum_Deconvolution5.gif"></sub></p>
<p class=MsoNormal style='margin-left:36.0pt;text-indent:-18.0pt'><span
style='font-size:16.0pt'>5.<span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>If r = R stop calculation, else</span></p>
<p class=MsoNormal style='text-indent:36.0pt'><span style='font-size:16.0pt'>a.
apply boosting operation, i.e., set </span><sub><img width=201 height=39
src="gif/TSpectrum_Deconvolution6.gif"></sub></p>
<p class=MsoNormal style='margin-left:36.0pt'><span style='font-size:16.0pt'>
i=0,1,...N-1 and p is boosting coefficient >0.</span></p>
<p class=MsoNormal style='margin-left:36.0pt'><span style='font-size:16.0pt'>b.
r = r + 1</span></p>
<p class=MsoNormal style='margin-left:36.0pt'><span style='font-size:16.0pt'>c.
continue in 4.</span></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'>const
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#char">char</a>* <a
name="TSpectrum:Deconvolution1">Deconvolution</a>(<a
href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a> *source,
const <a href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a>
*respMatrix, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a> ssize,
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a>
numberIterations, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a>
numberRepetitions, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#double">double</a></span></b><b><span
style='font-size:16.0pt'> </span></b><b><span style='font-size:18.0pt'>boost)</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This
function calculates deconvolution from source spectrum according to response
spectrum using Gold deconvolution algorithm. The result is placed in the vector
pointed by source pointer. On successful completion it returns 0. On error it
returns pointer to the string describing error. If desired after every
numberIterations one can apply boosting operation (exponential function with
exponent given by boost coefficient) and repeat it numberRepetitions times.</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>source</span></b>-pointer
to the vector of source spectrum </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>respMatrix</span></b>-pointer
to the vector of response spectrum </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>ssize</span></b>-length
of the spectrum vector </p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>numberIterations</span></b>-number of iterations
(parameter l in the Gold deconvolution </p>
<p class=MsoNormal style='text-align:justify'> algorithm)</p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>numberRepetitions</span></b>-number of repetitions
for boosted deconvolution. It must be </p>
<p class=MsoNormal style='text-align:justify'> greater or equal to one.</p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>boost</span></b>-boosting coefficient, applies only
if numberRepetitions is greater than one. </p>
<p class=MsoNormal style='text-align:justify'> <span style='font-size:
14.0pt;color:fuchsia'>Recommended range <1,2>.</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><i><span style='font-size:18.0pt'>References:</span></i></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[1]
Gold R., ANL-6984, Argonne National Laboratories, Argonne Ill, 1964. </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[2]
Coote G.E., Iterative smoothing and deconvolution of one- and two-dimensional
elemental distribution data, NIM B 130 (1997) 118.</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[3]
</span><span lang=SK style='font-size:16.0pt'>M. Morháč, J. Kliman, V.
Matouek, M. Veselskŭ, I. Turzo</span><span style='font-size:16.0pt'>.:
Efficient one- and two-dimensional Gold deconvolution and its application to
gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[4]
Morháč M., Matouek V., Kliman J., Efficient algorithm of multidimensional
deconvolution and its application to nuclear data processing, Digital Signal
Processing 13 (2003) 144. </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section10>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 8 script
Deconvolution.c :</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>response function (usually peak)
should be shifted left to the first non-zero channel (bin) (see Figure 9)</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><img width=600 height=340
src="gif/TSpectrum_Deconvolution1.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
9 Response spectrum</span></b></p>
<p class=MsoNormal><img width=946 height=407
src="gif/TSpectrum_Deconvolution2.jpg"></p>
<p class=MsoNormal><b><span style='font-size:16.0pt'>Figure 10 Principle how the
response matrix is composed inside of the Deconvolution function</span></b></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Deconvolution3.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
11 Example of Gold deconvolution. The original source spectrum is drawn with
black color, the spectrum after the deconvolution (10000 iterations) with red
color</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate deconvolution function (class
TSpectrum).</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Deconvolution.C</p>
<p class=MsoNormal> </p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void Deconvolution() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 256;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> Float_t * response = new float[nbins]; </p>
<p class=MsoNormal> TH1F *h = new
TH1F("h","Deconvolution",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d = new
TH1F("d","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("decon1;1");</p>
<p class=MsoNormal> TFile *fr = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> d=(TH1F*) fr->Get("decon_response;1"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
response[i]=d->GetBinContent(i + 1); </p>
<p class=MsoNormal> TCanvas *Decon1 =
gROOT->GetListOfCanvases()->FindObject("Decon1");</p>
<p class=MsoNormal> if (!Decon1) Decon1 = new
TCanvas("Decon1","Decon1",10,10,1000,700);</p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> s->Deconvolution(source,response,256,1000,1,1);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d->SetLineColor(kRed);</p>
<p class=MsoNormal> d->Draw("SAME L"); </p>
<p class=MsoNormal>}</p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section11>
<p class=MsoNormal><b><span style='font-size:18.0pt'>Examples of Gold
deconvolution method</span></b></p>
<p class=MsoNormal><b><span style='font-size:18.0pt'> </span></b></p>
<p class=MsoNormal style='margin-left:36.0pt;text-indent:-18.0pt'><span
style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>first let us study the influence
of the number of iterations on the deconvolved spectrum (Figure 12)</span></p>
<p class=MsoNormal><img width=602 height=409
src="gif/TSpectrum_Deconvolution_wide1.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
12 Study of Gold deconvolution algorithm. The original source spectrum is drawn
with black color, spectrum after 100 iterations with red color,</span></b><span
style='font-size:16.0pt'> <b>spectrum after 1000 iterations with blue color,</b>
<b>spectrum after 10000 iterations with green color and </b> <b>spectrum after
100000 iterations with magenta color.</b></span></p>
<p class=MsoNormal style='margin-left:19.95pt;text-align:justify'><span
style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='margin-left:37.05pt;text-align:justify;text-indent:
-17.1pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>for relatively narrow peaks in the
above given example the Gold deconvolution method is able to decompose
overlapping peaks practically to delta - functions.</span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>in the next example we have chosen
a synthetic data (spectrum, 256 channels) consisting of 5 very closely
positioned, relatively wide peaks (sigma =5), with added noise (Figure 13). </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>thin lines represent pure
Gaussians (see Table 1); thick line is a resulting spectrum with additive noise
(10% of the amplitude of small peaks).</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'><img
width=600 height=367 src="gif/TSpectrum_Deconvolution_wide2.jpg"></span></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
13 Testing example of synthetic spectrum composed of 5 Gaussians with added
noise</span></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:18.0pt'> </span></p>
<table class=MsoNormalTable border=0 cellspacing=0 cellpadding=0 width=445
style='width:333.65pt;margin-left:34.95pt'>
<tr style='height:16.6pt'>
<td width=80 valign=top style='width:60.3pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>Peak #</p>
</td>
<td width=132 valign=top style='width:99.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>Position</p>
</td>
<td width=111 valign=top style='width:83.05pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>Height </p>
</td>
<td width=121 valign=top style='width:91.1pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>Area</p>
</td>
</tr>
<tr style='height:16.6pt'>
<td width=80 valign=top style='width:60.3pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>1</p>
</td>
<td width=132 valign=top style='width:99.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>50</p>
</td>
<td width=111 valign=top style='width:83.05pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>500</p>
</td>
<td width=121 valign=top style='width:91.1pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>10159</p>
</td>
</tr>
<tr style='height:16.6pt'>
<td width=80 valign=top style='width:60.3pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>2</p>
</td>
<td width=132 valign=top style='width:99.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>70</p>
</td>
<td width=111 valign=top style='width:83.05pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>3000</p>
</td>
<td width=121 valign=top style='width:91.1pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>60957</p>
</td>
</tr>
<tr style='height:16.6pt'>
<td width=80 valign=top style='width:60.3pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>3</p>
</td>
<td width=132 valign=top style='width:99.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>80</p>
</td>
<td width=111 valign=top style='width:83.05pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>1000</p>
</td>
<td width=121 valign=top style='width:91.1pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>20319</p>
</td>
</tr>
<tr style='height:16.6pt'>
<td width=80 valign=top style='width:60.3pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>4</p>
</td>
<td width=132 valign=top style='width:99.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>100</p>
</td>
<td width=111 valign=top style='width:83.05pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>5000</p>
</td>
<td width=121 valign=top style='width:91.1pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>101596</p>
</td>
</tr>
<tr style='height:16.6pt'>
<td width=80 valign=top style='width:60.3pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>5</p>
</td>
<td width=132 valign=top style='width:99.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>110</p>
</td>
<td width=111 valign=top style='width:83.05pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>500</p>
</td>
<td width=121 valign=top style='width:91.1pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:16.6pt'>
<p class=MsoNormal style='text-align:justify'>10159</p>
</td>
</tr>
</table>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Table
1 Positions, heights and areas of peaks in the spectrum shown in Figure 13</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>in ideal case, we should obtain
the result given in Figure 14. The areas of the Gaussian components of the
spectrum are concentrated completely to delta functions</span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>when solving the overdetermined
system of linear equations with data from Figure 13 in the sense of minimum
least squares criterion without any regularization we obtain the result with
large oscillations (Figure 15). </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>from mathematical point of view,
it is the optimal solution in the unconstrained space of independent variables.
>From physical point of view we are interested only in a meaningful solution. </span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>therefore, we have to employ
regularization techniques (e.g. Gold deconvolution) and/or to confine the space
of allowed solutions to subspace of positive solutions. </span></p>
<p class=MsoNormal style='margin-left:18.0pt;text-align:justify'><span
style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-indent:5.7pt'><span style='font-size:16.0pt'><img
width=589 height=189 src="gif/TSpectrum_Deconvolution_wide3.jpg"></span></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
14 The same spectrum like in Figure 13, outlined bars show the contents of
present components (peaks) </span></b></p>
<p class=MsoNormal style='text-align:justify;text-indent:8.55pt'><span
style='font-size:16.0pt'><img width=585 height=183
src="gif/TSpectrum_Deconvolution_wide4.jpg"></span></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
15 Least squares solution of the system of linear equations without
regularization</span></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 9 script Deconvolution_wide.c
:</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:14.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>when we employ Gold deconvolution
algorithm we obtain the result given in Fig. 16. One can observe that the
resulting spectrum is smooth. On the other hand the method is not able to
decompose completely the peaks in the spectrum.</span></p>
<p class=MsoNormal style='text-align:justify'> <img width=601 height=407
src="gif/TSpectrum_Deconvolution_wide5.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
16 Example of Gold deconvolution for closely positioned wide peaks. The
original source spectrum is drawn with black color, the spectrum after the
deconvolution (10000 iterations) with red color</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate deconvolution function (class
TSpectrum).</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x Deconvolution_wide.C</p>
<p class=MsoNormal> </p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void Deconvolution_wide() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 256;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> Float_t * response = new float[nbins]; </p>
<p class=MsoNormal> TH1F *h = new
TH1F("h","Deconvolution",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d = new
TH1F("d","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("decon3;1");</p>
<p class=MsoNormal> TFile *fr = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> d=(TH1F*)
fr->Get("decon_response_wide;1"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
response[i]=d->GetBinContent(i + 1); </p>
<p class=MsoNormal> TCanvas *Decon1 =
gROOT->GetListOfCanvases()->FindObject("Decon1");</p>
<p class=MsoNormal> if (!Decon1) Decon1 = new
TCanvas("Decon1","Deconvolution of closely positioned
overlapping peaks using Gold deconvolution method",10,10,1000,700);</p>
<p class=MsoNormal> h->SetMaximum(30000);</p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> s->Deconvolution(source,response,256,10000,1,1);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d->SetLineColor(kRed);</p>
<p class=MsoNormal> d->Draw("SAME L"); </p>
<p class=MsoNormal>}</p>
<p class=MsoNormal><i><span style='font-size:16.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 10 script
Deconvolution_wide_boost.c :</span></i></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:14.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>further let us employ boosting
operation into deconvolution (Fig. 17)</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_Deconvolution_wide6.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
17 The original source spectrum is drawn with black color, the spectrum after
the deconvolution with red color. Number of iterations = 200, number of
repetitions = 50 and boosting coefficient = 1.2.</span></b></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'> </span></b></p>
<table class=MsoTableGrid border=1 cellspacing=0 cellpadding=0 width=469
style='width:351.65pt;margin-left:50.55pt;border-collapse:collapse;border:
none'>
<tr style='height:17.15pt'>
<td valign=top style='border:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;
height:17.15pt'>
<p class=MsoNormal style='text-align:justify;text-indent:50.55pt'>Peak #</p>
</td>
<td valign=top style='border:solid windowtext 1.0pt;border-left:none;
padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>Original/Estimated (max)
position</p>
</td>
<td valign=top style='border:solid windowtext 1.0pt;border-left:none;
padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>Original/Estimated area</p>
</td>
</tr>
<tr style='height:17.15pt'>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>1</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>50/49</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>10159/10419</p>
</td>
</tr>
<tr style='height:17.15pt'>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>2</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>70/70</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>60957/58933</p>
</td>
</tr>
<tr style='height:17.15pt'>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>3</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>80/79</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>20319/19935</p>
</td>
</tr>
<tr style='height:17.15pt'>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>4</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>100/100</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:17.15pt'>
<p class=MsoNormal style='text-align:justify'>101596/105413</p>
</td>
</tr>
<tr style='height:18.1pt'>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt;height:18.1pt'>
<p class=MsoNormal style='text-align:justify'>5</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:18.1pt'>
<p class=MsoNormal style='text-align:justify'>110/117</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt;height:18.1pt'>
<p class=MsoNormal style='text-align:justify;page-break-after:avoid'>10159/6676</p>
</td>
</tr>
</table>
<p class=MsoCaption style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoCaption style='text-align:justify'><span style='font-size:16.0pt'>Table
2 Results of the estimation of peaks in spectrum shown in Figure 17</span></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'> </span></b></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-36.0pt'><span style='font-size:14.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>one can observe that peaks are
decomposed practically to delta functions. Number of peaks is correct,
positions of big peaks as well as their areas are relatively well estimated.
However there is a considerable error in the estimation of the position of
small right hand peak.</span></p>
<p class=MsoNormal style='margin-left:18.0pt;text-align:justify'><span
style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal style='text-align:justify'>// Example to illustrate
deconvolution function (class TSpectrum).</p>
<p class=MsoNormal style='text-align:justify'>// To execute this example, do</p>
<p class=MsoNormal style='text-align:justify'>// root > .x
Deconvolution_wide_boost.C</p>
<p class=MsoNormal style='text-align:justify'> </p>
<p class=MsoNormal style='text-align:justify'>#include <TSpectrum></p>
<p class=MsoNormal style='text-align:justify'> </p>
<p class=MsoNormal style='text-align:justify'>void Deconvolution_wide_boost() {</p>
<p class=MsoNormal style='text-align:justify'> Int_t i;</p>
<p class=MsoNormal style='text-align:justify'> Double_t nbins = 256;</p>
<p class=MsoNormal style='text-align:justify'> Double_t xmin = 0;</p>
<p class=MsoNormal style='text-align:justify'> Double_t xmax =
(Double_t)nbins;</p>
<p class=MsoNormal style='text-align:justify'> Float_t * source = new
float[nbins];</p>
<p class=MsoNormal style='text-align:justify'> Float_t * response = new
float[nbins]; </p>
<p class=MsoNormal style='text-align:justify'> TH1F *h = new
TH1F("h","Deconvolution",nbins,xmin,xmax);</p>
<p class=MsoNormal style='text-align:justify'> TH1F *d = new
TH1F("d","",nbins,xmin,xmax); </p>
<p class=MsoNormal style='text-align:justify'> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal style='text-align:justify'> h=(TH1F*)
f->Get("decon3;1");</p>
<p class=MsoNormal style='text-align:justify'> TFile *fr = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal style='text-align:justify'> d=(TH1F*)
fr->Get("decon_response_wide;1"); </p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbins; i++)
response[i]=d->GetBinContent(i + 1); </p>
<p class=MsoNormal style='text-align:justify'> TCanvas *Decon1 =
gROOT->GetListOfCanvases()->FindObject("Decon1");</p>
<p class=MsoNormal style='text-align:justify'> if (!Decon1) Decon1 = new
TCanvas("Decon1","Deconvolution of closely positioned
overlapping peaks using boosted Gold deconvolution
method",10,10,1000,700);</p>
<p class=MsoNormal style='text-align:justify'> h->SetMaximum(110000);</p>
<p class=MsoNormal style='text-align:justify'> h->Draw("L");</p>
<p class=MsoNormal style='text-align:justify'> TSpectrum *s = new
TSpectrum();</p>
<p class=MsoNormal style='text-align:justify'>
s->Deconvolution(source,response,256,200,50,1.2);</p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbins; i++)
d->SetBinContent(i + 1,source[i]); </p>
<p class=MsoNormal style='text-align:justify'> d->SetLineColor(kRed);</p>
<p class=MsoNormal style='text-align:justify'> d->Draw("SAME
L"); </p>
<p class=MsoNormal style='text-align:justify'>}</p>
</div>
<!-- */
// --> End_Html
if (ssize <= 0)
return "Wrong Parameters";
if (numberRepetitions <= 0)
return "Wrong Parameters";
double *working_space = new double[4 * ssize];
int i, j, k, lindex, posit, lh_gold, l, repet;
double lda, ldb, ldc, area, maximum;
area = 0;
lh_gold = -1;
posit = 0;
maximum = 0;
for (i = 0; i < ssize; i++) {
lda = response[i];
if (lda != 0)
lh_gold = i + 1;
working_space[i] = lda;
area += lda;
if (lda > maximum) {
maximum = lda;
posit = i;
}
}
if (lh_gold == -1)
return "ZERO RESPONSE VECTOR";
for (i = 0; i < ssize; i++)
working_space[2 * ssize + i] = source[i];
for (i = 0; i < ssize; i++){
lda = 0;
for (j = 0; j < ssize; j++){
ldb = working_space[j];
k = i + j;
if (k < ssize){
ldc = working_space[k];
lda = lda + ldb * ldc;
}
}
working_space[ssize + i] = lda;
lda = 0;
for (k = 0; k < ssize; k++){
l = k - i;
if (l >= 0){
ldb = working_space[l];
ldc = working_space[2 * ssize + k];
lda = lda + ldb * ldc;
}
}
working_space[3 * ssize + i]=lda;
}
for (i = 0; i < ssize; i++){
working_space[2 * ssize + i] = working_space[3 * ssize + i];
}
for (i = 0; i < ssize; i++)
working_space[i] = 1;
for (repet = 0; repet < numberRepetitions; repet++) {
if (repet != 0) {
for (i = 0; i < ssize; i++)
working_space[i] = TMath::Power(working_space[i], boost);
}
for (lindex = 0; lindex < numberIterations; lindex++) {
for (i = 0; i < ssize; i++) {
if (working_space[2 * ssize + i] > 0.000001
&& working_space[i] > 0.000001) {
lda = 0;
for (j = 0; j < lh_gold; j++) {
ldb = working_space[j + ssize];
if (j != 0){
k = i + j;
ldc = 0;
if (k < ssize)
ldc = working_space[k];
k = i - j;
if (k >= 0)
ldc += working_space[k];
}
else
ldc = working_space[i];
lda = lda + ldb * ldc;
}
ldb = working_space[2 * ssize + i];
if (lda != 0)
lda = ldb / lda;
else
lda = 0;
ldb = working_space[i];
lda = lda * ldb;
working_space[3 * ssize + i] = lda;
}
}
for (i = 0; i < ssize; i++)
working_space[i] = working_space[3 * ssize + i];
}
}
for (i = 0; i < ssize; i++) {
lda = working_space[i];
j = i + posit;
j = j % ssize;
working_space[ssize + j] = lda;
}
for (i = 0; i < ssize; i++)
source[i] = area * working_space[ssize + i];
delete[]working_space;
return 0;
}
const char *TSpectrum::DeconvolutionRL(float *source, const float *response,
int ssize, int numberIterations,
int numberRepetitions, double boost )
{
//Begin_Html <!--
/* -->
<div class=Section12>
<p class=MsoNormal><b><i><span style='font-size:16.0pt'>Richardson-Lucy
deconvolution algorithm</span></i></b></p>
<p class=MsoNormal style='margin-left:32.2pt;text-align:justify;text-indent:
-32.2pt'><span lang=SK style='font-size:14.0pt;font-family:Symbol'>·<span
style='font:7.0pt "Times New Roman"'>
</span></span><span lang=SK style='font-size:16.0pt'>for discrete systems it
has the form</span></p>
<p class=MsoNormal style='margin-left:18.0pt;text-align:justify'><span lang=SK> <sub><img
width=438 height=98 src="gif/TSpectrum_DeconvolutionRL1.gif"></sub>
</span></p>
<p class=MsoNormal style='margin-left:27.0pt;text-align:justify'><span lang=SK> <sub><img
width=124 height=39 src="gif/TSpectrum_DeconvolutionRL2.gif"></sub></span></p>
<p class=MsoNormal style='margin-left:28.5pt;text-align:justify;text-indent:
-28.5pt'><span lang=SK style='font-size:14.0pt;font-family:Symbol'>·<span
style='font:7.0pt "Times New Roman"'>
</span></span><span lang=SK style='font-size:16.0pt'>for positive input data
and response matrix this iterative method forces the deconvoluted spectra to be
non-negative. </span></p>
<p class=MsoNormal style='margin-left:27.0pt;text-align:justify;text-indent:
-27.0pt'><span lang=SK style='font-size:14.0pt;font-family:Symbol'>·<span
style='font:7.0pt "Times New Roman"'>
</span></span><span lang=SK style='font-size:16.0pt'>the Richardson-Lucy
iteration converges to the maximum likelihood solution for Poisson statistics
in the data.</span></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'>const
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#char">char</a>* <a
name="TSpectrum:Deconvolution1">Deconvolution</a>RL(<a
href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a> *source,
const <a href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a>
*respMatrix, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a> ssize,
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a>
numberIterations, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a>
numberRepetitions, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#double">double</a> boost)</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This
function calculates deconvolution from source spectrum according to response
spectrum using Richardson-Lucy deconvolution algorithm. The result is placed in
the vector pointed by source pointer. On successful completion it returns 0. On
error it returns pointer to the string describing error. If desired after every
numberIterations one can apply boosting operation (exponential function with
exponent given by boost coefficient) and repeat it numberRepetitions times
(see Gold deconvolution).</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>source</span></b>-pointer
to the vector of source spectrum </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>respMatrix</span></b>-pointer
to the vector of response spectrum </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>ssize</span></b>-length
of the spectrum vector </p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>numberIterations</span></b>-number of iterations
(parameter l in the Gold deconvolution </p>
<p class=MsoNormal style='text-align:justify'> algorithm)</p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>numberRepetitions</span></b>-number of repetitions
for boosted deconvolution. It must be </p>
<p class=MsoNormal style='text-align:justify'> greater or equal to one.</p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>boost</span></b>-boosting coefficient, applies only
if numberRepetitions is greater than one. </p>
<p class=MsoNormal style='text-align:justify'> <span style='font-size:
14.0pt;color:fuchsia'>Recommended range <1,2>.</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><i><span style='font-size:18.0pt'>References:</span></i></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[1]
Abreu M.C. et al., A four-dimensional deconvolution method to correct NA38
experimental data, NIM A 405 (1998) 139.</span></p>
<p class=MsoNormal style='margin-left:18.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:16.0pt'>[2] Lucy L.B., A.J. 79 (1974) 745.</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[3]
Richardson W.H., J. Opt. Soc. Am. 62 (1972) 55.</span></p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section13>
<p class=MsoNormal><b><span style='font-size:18.0pt'>Examples of Richardson-Lucy
deconvolution method</span></b></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 11 script DeconvolutionRL_wide.c
:</span></i></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:14.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>when we employ Richardson-Lucy
deconvolution algorithm to our data from Fig. 13 we obtain the result given in
Fig. 18. One can observe improvements as compared to the result achieved by
Gold deconvolution.</span></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:14.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>neverthless it is unable to
decompose the multiplet.</span></p>
<p class=MsoNormal style='text-align:justify'> <img width=601 height=407
src="gif/TSpectrum_DeconvolutionRL_wide1.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
18 Example of Richardson-Lucy deconvolution for closely positioned wide peaks.
The original source spectrum is drawn with black color, the spectrum after the
deconvolution (10000 iterations) with red color</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'> </span></b></p>
<p class=MsoNormal>// Example to illustrate deconvolution function (class
TSpectrum).</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x DeconvolutionRL_wide.C</p>
<p class=MsoNormal> </p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void DeconvolutionRL_wide() {</p>
<p class=MsoNormal> Int_t i;</p>
<p class=MsoNormal> Double_t nbins = 256;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> Float_t * response = new float[nbins]; </p>
<p class=MsoNormal> TH1F *h = new
TH1F("h","Deconvolution",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d = new TH1F("d","",nbins,xmin,xmax);
</p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("decon3;1");</p>
<p class=MsoNormal> TFile *fr = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> d=(TH1F*)
fr->Get("decon_response_wide;1"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i
+ 1);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
response[i]=d->GetBinContent(i + 1); </p>
<p class=MsoNormal> TCanvas *Decon1 =
gROOT->GetListOfCanvases()->FindObject("Decon1");</p>
<p class=MsoNormal> if (!Decon1) Decon1 = new
TCanvas("Decon1","Deconvolution of closely positioned
overlapping peaks using Richardson-Lucy deconvolution
method",10,10,1000,700);</p>
<p class=MsoNormal> h->SetMaximum(30000);</p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> s->DeconvolutionRL(source,response,256,10000,1,1);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d->SetBinContent(i +
1,source[i]); </p>
<p class=MsoNormal> d->SetLineColor(kRed);</p>
<p class=MsoNormal> d->Draw("SAME L"); </p>
<p class=MsoNormal>}</p>
<p class=MsoNormal><i><span style='font-size:16.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 12 script
DeconvolutionRL_wide_boost.c :</span></i></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-18.0pt'><span style='font-size:14.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>further let us employ boosting
operation into deconvolution (Fig. 19)</span></p>
<p class=MsoNormal><img width=601 height=407
src="gif/TSpectrum_DeconvolutionRL_wide2.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Figure
19 The original source spectrum is drawn with black color, the spectrum after
the deconvolution with red color. Number of iterations = 200, number of
repetitions = 50 and boosting coefficient = 1.2.</span></b></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'> </span></b></p>
<table class=MsoTableGrid border=1 cellspacing=0 cellpadding=0
style='margin-left:28.5pt;border-collapse:collapse;border:none'>
<tr>
<td valign=top style='border:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>Peak #</p>
</td>
<td valign=top style='border:solid windowtext 1.0pt;border-left:none;
padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>Original/Estimated (max)
position</p>
</td>
<td valign=top style='border:solid windowtext 1.0pt;border-left:none;
padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>Original/Estimated area</p>
</td>
</tr>
<tr>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>1</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>50/51</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>10159/11426</p>
</td>
</tr>
<tr>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>2</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>70/71</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>60957/65003</p>
</td>
</tr>
<tr>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>3</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>80/81</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>20319/12813</p>
</td>
</tr>
<tr>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>4</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>100/100</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>101596/101851</p>
</td>
</tr>
<tr>
<td valign=top style='border:solid windowtext 1.0pt;border-top:none;
padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>5</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify'>110/111</p>
</td>
<td valign=top style='border-top:none;border-left:none;border-bottom:solid windowtext 1.0pt;
border-right:solid windowtext 1.0pt;padding:0mm 5.4pt 0mm 5.4pt'>
<p class=MsoNormal style='text-align:justify;page-break-after:avoid'>10159/8920</p>
</td>
</tr>
</table>
<p class=MsoCaption style='text-align:justify'><span style='font-size:16.0pt'>Table
3 Results of the estimation of peaks in spectrum shown in Figure 19</span></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'> </span></b></p>
<p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
-36.0pt'><span style='font-size:14.0pt'><span style='font:7.0pt "Times New Roman"'>
</span></span><span style='font-size:16.0pt'>one can observe improvements in
the estimation of peak positions as compared to the results achieved by Gold
deconvolution</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal style='text-align:justify'> </p>
<p class=MsoNormal style='text-align:justify'>// Example to illustrate
deconvolution function (class TSpectrum).</p>
<p class=MsoNormal style='text-align:justify'>// To execute this example, do</p>
<p class=MsoNormal style='text-align:justify'>// root > .x
DeconvolutionRL_wide_boost.C</p>
<p class=MsoNormal style='text-align:justify'>#include <TSpectrum></p>
<p class=MsoNormal style='text-align:justify'> </p>
<p class=MsoNormal style='text-align:justify'>void DeconvolutionRL_wide_boost()
{</p>
<p class=MsoNormal style='text-align:justify'> Int_t i;</p>
<p class=MsoNormal style='text-align:justify'> Double_t nbins = 256;</p>
<p class=MsoNormal style='text-align:justify'> Double_t xmin = 0;</p>
<p class=MsoNormal style='text-align:justify'> Double_t xmax =
(Double_t)nbins;</p>
<p class=MsoNormal style='text-align:justify'> Float_t * source = new
float[nbins];</p>
<p class=MsoNormal style='text-align:justify'> Float_t * response = new
float[nbins]; </p>
<p class=MsoNormal style='text-align:justify'> TH1F *h = new
TH1F("h","Deconvolution",nbins,xmin,xmax);</p>
<p class=MsoNormal style='text-align:justify'> TH1F *d = new
TH1F("d","",nbins,xmin,xmax); </p>
<p class=MsoNormal style='text-align:justify'> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal style='text-align:justify'> h=(TH1F*)
f->Get("decon3;1");</p>
<p class=MsoNormal style='text-align:justify'> TFile *fr = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal style='text-align:justify'> d=(TH1F*)
fr->Get("decon_response_wide;1"); </p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbins; i++)
response[i]=d->GetBinContent(i + 1); </p>
<p class=MsoNormal style='text-align:justify'> TCanvas *Decon1 =
gROOT->GetListOfCanvases()->FindObject("Decon1");</p>
<p class=MsoNormal style='text-align:justify'> if (!Decon1) Decon1 = new
TCanvas("Decon1","Deconvolution of closely positioned
overlapping peaks using boosted Richardson-Lucy deconvolution
method",10,10,1000,700);</p>
<p class=MsoNormal style='text-align:justify'> h->SetMaximum(110000);</p>
<p class=MsoNormal style='text-align:justify'> h->Draw("L");</p>
<p class=MsoNormal style='text-align:justify'> TSpectrum *s = new
TSpectrum();</p>
<p class=MsoNormal style='text-align:justify'> s->DeconvolutionRL(source,response,256,200,50,1.2);</p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbins; i++)
d->SetBinContent(i + 1,source[i]); </p>
<p class=MsoNormal style='text-align:justify'> d->SetLineColor(kRed);</p>
<p class=MsoNormal style='text-align:justify'> d->Draw("SAME
L"); </p>
<p class=MsoNormal style='text-align:justify'>}</p>
</div>
<!-- */
// --> End_Html
if (ssize <= 0)
return "Wrong Parameters";
if (numberRepetitions <= 0)
return "Wrong Parameters";
double *working_space = new double[4 * ssize];
int i, j, k, lindex, posit, lh_gold, repet, kmin, kmax;
double lda, ldb, ldc, maximum;
lh_gold = -1;
posit = 0;
maximum = 0;
for (i = 0; i < ssize; i++) {
lda = response[i];
if (lda != 0)
lh_gold = i + 1;
working_space[ssize + i] = lda;
if (lda > maximum) {
maximum = lda;
posit = i;
}
}
if (lh_gold == -1)
return "ZERO RESPONSE VECTOR";
for (i = 0; i < ssize; i++)
working_space[2 * ssize + i] = source[i];
for (i = 0; i < ssize; i++){
if (i <= ssize - lh_gold)
working_space[i] = 1;
else
working_space[i] = 0;
}
for (repet = 0; repet < numberRepetitions; repet++) {
if (repet != 0) {
for (i = 0; i < ssize; i++)
working_space[i] = TMath::Power(working_space[i], boost);
}
for (lindex = 0; lindex < numberIterations; lindex++) {
for (i = 0; i <= ssize - lh_gold; i++){
lda = 0;
if (working_space[i] > 0){
for (j = i; j < i + lh_gold; j++){
ldb = working_space[2 * ssize + j];
if (j < ssize){
if (ldb > 0){
kmax = j;
if (kmax > lh_gold - 1)
kmax = lh_gold - 1;
kmin = j + lh_gold - ssize;
if (kmin < 0)
kmin = 0;
ldc = 0;
for (k = kmax; k >= kmin; k--){
ldc += working_space[ssize + k] * working_space[j - k];
}
if (ldc > 0)
ldb = ldb / ldc;
else
ldb = 0;
}
ldb = ldb * working_space[ssize + j - i];
}
lda += ldb;
}
lda = lda * working_space[i];
}
working_space[3 * ssize + i] = lda;
}
for (i = 0; i < ssize; i++)
working_space[i] = working_space[3 * ssize + i];
}
}
for (i = 0; i < ssize; i++) {
lda = working_space[i];
j = i + posit;
j = j % ssize;
working_space[ssize + j] = lda;
}
for (i = 0; i < ssize; i++)
source[i] = working_space[ssize + i];
delete[]working_space;
return 0;
}
const char *TSpectrum::Unfolding(float *source,
const float **respMatrix,
int ssizex, int ssizey,
int numberIterations,
int numberRepetitions, double boost)
{
//Begin_Html <!--
/* -->
<div class=Section14>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:20.0pt'>Unfolding</span></b></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>Goal:
Decomposition of spectrum to a given set of component spectra</span></i></p>
<p class=MsoNormal><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:16.0pt'>Mathematical formulation of
the discrete linear system is</span></p>
<p class=MsoNormal style='margin-left:18.0pt'>
<table cellpadding=0 cellspacing=0 align=left>
<tr>
<td width=0 height=17></td>
</tr>
<tr>
<td></td>
<td><img width=588 height=89 src="gif/TSpectrum_Unfolding1.gif"></td>
</tr>
</table>
<span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='text-align:justify'>
<table cellpadding=0 cellspacing=0 align=left>
<tr>
<td width=0 height=3></td>
</tr>
<tr>
<td></td>
<td><img width=597 height=228 src="gif/TSpectrum_Unfolding2.gif"></td>
</tr>
</table>
<span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<br clear=ALL>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'>const
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#char">char</a>* Unfolding(<a
href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a> *source,
const <a href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a>
**respMatrix, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a> ssizex,
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a> ssizey, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a>
numberIterations, <a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a>
numberRepetitions, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#double">double</a></span></b><b><span
style='font-size:16.0pt'> </span></b><b><span style='font-size:18.0pt'>boost)</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This
function unfolds source spectrum according to response matrix columns. The
result is placed in the vector pointed by source pointer. The coefficients of
the resulting vector represent contents of the columns (weights) in the input
vector. On successful completion it returns 0. On error it returns pointer to
the string describing error. If desired after every numberIterations one can
apply boosting operation (exponential function with exponent given by boost
coefficient) and repeat it numberRepetitions times. For details we refer to
[1].</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>source</span></b>-pointer
to the vector of source spectrum </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>respMatrix</span></b>-pointer
to the matrix of response spectra </p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>ssizex</span></b>-length
of source spectrum and # of columns of the response matrix</p>
<p class=MsoNormal> <b><span style='font-size:14.0pt'>ssizey</span></b>-length
of destination spectrum and # of rows of the response matrix
</p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>numberIterations</span></b>-number of iterations </p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>numberRepetitions</span></b>-number of repetitions
for boosted deconvolution. It must be </p>
<p class=MsoNormal style='text-align:justify'> greater or equal to one.</p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>boost</span></b>-boosting coefficient, applies only
if numberRepetitions is greater than one. </p>
<p class=MsoNormal style='text-align:justify'> <span style='font-size:
14.0pt;color:fuchsia'> Recommended range <1,2>.</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:14.0pt;
color:fuchsia'>Note!!! sizex must be >= sizey After decomposition the
resulting channels are written back to the first sizey channels of the source
spectrum. </span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><i><span style='font-size:18.0pt'>Reference:</span></i></b></p>
<p class=MsoNormal style='text-align:justify'><span lang=SK style='font-size:
16.0pt'>[1] Jandel M., </span><span style='font-size:16.0pt'>Morháč M., </span><span
lang=SK style='font-size:16.0pt'>Kliman J., Krupa L.,</span><span
style='font-size:16.0pt'> Matou</span><span lang=SK style='font-size:16.0pt'>ek
V., Hamilton J. H., Ramaya A. V.</span><span style='font-size:16.0pt'>:
Decomposition of continuum gamma-ray spectra using synthetized response matrix.
NIM A 516 (2004), 172-183.</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section15>
<p class=MsoNormal><b><span style='font-size:18.0pt'>Example of unfolding</span></b></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'> </span></i></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 13 script Unfolding.c
:</span></i></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'><img width=442 height=648
src="gif/TSpectrum_Unfolding3.gif"></span></i></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt;
font-family:Arial'>Fig. 20 Response matrix composed of neutron spectra of pure
chemical elements</span></b><b><span style='font-size:16.0pt'> </span></b></p>
<p class=MsoNormal><img width=604 height=372
src="gif/TSpectrum_Unfolding2.jpg"></p>
<p class=MsoNormal><b><span style='font-size:16.0pt;font-family:Arial'>Fig. 21
Source neutron spectrum to be decomposed</span></b></p>
<p class=MsoNormal><img width=600 height=360
src="gif/TSpectrum_Unfolding3.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt;
font-family:Arial'>Fig. 22 Spectrum after decomposition, contains 10
coefficients, which correspond to contents of chemical components (dominant
8-th and 10-th components, i.e. O, Si)</span></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal style='text-align:justify'> </p>
<p class=MsoNormal style='text-align:justify'>// Example to illustrate
unfolding function (class TSpectrum).</p>
<p class=MsoNormal style='text-align:justify'>// To execute this example, do</p>
<p class=MsoNormal style='text-align:justify'>// root > .x Unfolding.C</p>
<p class=MsoNormal style='text-align:justify'> </p>
<p class=MsoNormal style='text-align:justify'>#include <TSpectrum></p>
<p class=MsoNormal style='text-align:justify'>void Unfolding() {</p>
<p class=MsoNormal style='text-align:justify'> Int_t i, j;</p>
<p class=MsoNormal style='text-align:justify'> Int_t nbinsx = 2048;</p>
<p class=MsoNormal style='text-align:justify'> Int_t nbinsy = 10; </p>
<p class=MsoNormal style='text-align:justify'> Double_t xmin = 0;</p>
<p class=MsoNormal style='text-align:justify'> Double_t xmax =
(Double_t)nbinsx;</p>
<p class=MsoNormal style='text-align:justify'> Double_t ymin = 0;</p>
<p class=MsoNormal style='text-align:justify'> Double_t ymax =
(Double_t)nbinsy; </p>
<p class=MsoNormal style='text-align:justify'> Float_t * source = new
float[nbinsx];</p>
<p class=MsoNormal style='text-align:justify'> Float_t ** response = new
float *[nbinsy]; </p>
<p class=MsoNormal style='text-align:justify'> for (i=0;i<nbinsy;i++)</p>
<p class=MsoNormal style='text-align:justify'> response[i]=new
float[nbinsx]; </p>
<p class=MsoNormal style='text-align:justify'> TH1F *h = new
TH1F("h","",nbinsx,xmin,xmax);</p>
<p class=MsoNormal style='text-align:justify'> TH1F *d = new
TH1F("d","Decomposition - unfolding",nbinsx,xmin,xmax); </p>
<p class=MsoNormal style='text-align:justify'> TH2F *decon_unf_resp = new
TH2F("decon_unf_resp","Root File",nbinsy,ymin,ymax,nbinsx,xmin,xmax);</p>
<p class=MsoNormal style='text-align:justify'> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal style='text-align:justify'> h=(TH1F*)
f->Get("decon_unf_in;1");</p>
<p class=MsoNormal style='text-align:justify'> TFile *fr = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal style='text-align:justify'> decon_unf_resp = (TH2F*)
fr->Get("decon_unf_resp;1"); </p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbinsx;
i++) source[i] = h->GetBinContent(i + 1);</p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbinsy;
i++){</p>
<p class=MsoNormal style='text-align:justify'> for (j = 0; j< nbinsx;
j++){</p>
<p class=MsoNormal style='text-align:justify'> response[i][j] =
decon_unf_resp->GetBinContent(i + 1, j + 1);</p>
<p class=MsoNormal style='text-align:justify'> }</p>
<p class=MsoNormal style='text-align:justify'> } </p>
<p class=MsoNormal style='text-align:justify'> TCanvas *Decon1 =
gROOT->GetListOfCanvases()->FindObject("Decon1");</p>
<p class=MsoNormal style='text-align:justify'> if (!Decon1) Decon1 = new
TCanvas("Decon1","Decon1",10,10,1000,700);</p>
<p class=MsoNormal style='text-align:justify'> h->Draw("L"); </p>
<p class=MsoNormal style='text-align:justify'> TSpectrum *s = new
TSpectrum();</p>
<p class=MsoNormal style='text-align:justify'> s->Unfolding(source,response,nbinsx,nbinsy,1000,1,1);</p>
<p class=MsoNormal style='text-align:justify'> for (i = 0; i < nbinsy;
i++) d->SetBinContent(i + 1,source[i]); </p>
<p class=MsoNormal style='text-align:justify'> d->SetLineColor(kRed); </p>
<p class=MsoNormal style='text-align:justify'>
d->SetAxisRange(0,nbinsy); </p>
<p class=MsoNormal style='text-align:justify'> d->Draw("");</p>
<p class=MsoNormal style='text-align:justify'>}</p>
</div>
<!-- */
// --> End_Html
int i, j, k, lindex, lhx = 0, repet;
double lda, ldb, ldc, area;
if (ssizex <= 0 || ssizey <= 0)
return "Wrong Parameters";
if (ssizex < ssizey)
return "Sizex must be greater than sizey)";
if (numberIterations <= 0)
return "Number of iterations must be positive";
double *working_space =
new double[ssizex * ssizey + 2 * ssizey * ssizey + 4 * ssizex];
for (j = 0; j < ssizey && lhx != -1; j++) {
area = 0;
lhx = -1;
for (i = 0; i < ssizex; i++) {
lda = respMatrix[j][i];
if (lda != 0) {
lhx = i + 1;
}
working_space[j * ssizex + i] = lda;
area = area + lda;
}
if (lhx != -1) {
for (i = 0; i < ssizex; i++)
working_space[j * ssizex + i] /= area;
}
}
if (lhx == -1)
return ("ZERO COLUMN IN RESPONSE MATRIX");
for (i = 0; i < ssizex; i++)
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
source[i];
for (i = 0; i < ssizey; i++) {
for (j = 0; j < ssizey; j++) {
lda = 0;
for (k = 0; k < ssizex; k++) {
ldb = working_space[ssizex * i + k];
ldc = working_space[ssizex * j + k];
lda = lda + ldb * ldc;
}
working_space[ssizex * ssizey + ssizey * i + j] = lda;
}
lda = 0;
for (k = 0; k < ssizex; k++) {
ldb = working_space[ssizex * i + k];
ldc =
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
k];
lda = lda + ldb * ldc;
}
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
lda;
}
for (i = 0; i < ssizey; i++)
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
for (i = 0; i < ssizey; i++) {
for (j = 0; j < ssizey; j++) {
lda = 0;
for (k = 0; k < ssizey; k++) {
ldb = working_space[ssizex * ssizey + ssizey * i + k];
ldc = working_space[ssizex * ssizey + ssizey * j + k];
lda = lda + ldb * ldc;
}
working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j] =
lda;
}
lda = 0;
for (k = 0; k < ssizey; k++) {
ldb = working_space[ssizex * ssizey + ssizey * i + k];
ldc =
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
k];
lda = lda + ldb * ldc;
}
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
lda;
}
for (i = 0; i < ssizey; i++)
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
for (i = 0; i < ssizey; i++)
working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = 1;
for (repet = 0; repet < numberRepetitions; repet++) {
if (repet != 0) {
for (i = 0; i < ssizey; i++)
working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = TMath::Power(working_space[ssizex * ssizey + 2 * ssizey * ssizey + i], boost);
}
for (lindex = 0; lindex < numberIterations; lindex++) {
for (i = 0; i < ssizey; i++) {
lda = 0;
for (j = 0; j < ssizey; j++) {
ldb =
working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j];
ldc = working_space[ssizex * ssizey + 2 * ssizey * ssizey + j];
lda = lda + ldb * ldc;
}
ldb =
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i];
if (lda != 0) {
lda = ldb / lda;
}
else
lda = 0;
ldb = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
lda = lda * ldb;
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] = lda;
}
for (i = 0; i < ssizey; i++)
working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
}
}
for (i = 0; i < ssizex; i++) {
if (i < ssizey)
source[i] = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
else
source[i] = 0;
}
delete[]working_space;
return 0;
}
Int_t TSpectrum::SearchHighRes(float *source,float *destVector, int ssize,
float sigma, double threshold,
bool backgroundRemove,int deconIterations,
bool markov, int averWindow)
{
//Begin_Html <!--
/* -->
<div class=Section18>
<p class=MsoNormal><b><span style='font-size:20.0pt'>Peaks searching</span></b></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p>
<p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>Goal:
to identify automatically the peaks in spectrum with the presence of the
continuous background and statistical fluctuations - noise.</span></i><span
style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:16.0pt;font-family:Arial'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt;
font-family:Arial'>The common problems connected with correct peak
identification are</span></p>
<ul style='margin-top:0mm' type=disc>
<li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt;
font-family:Arial'>non-sensitivity to noise, i.e., only statistically
relevant peaks should be identified.</span></li>
<li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt;
font-family:Arial'>non-sensitivity of the algorithm to continuous
background</span></li>
<li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt;
font-family:Arial'>ability to identify peaks close to the edges of the
spectrum region. Usually peak finders fail to detect them</span></li>
<li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt;
font-family:Arial'>resolution, decomposition of doublets and multiplets.
The algorithm should be able to recognize close positioned peaks.</span><span
style='font-size:18.0pt'> </span></li>
<li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt;
font-family:Arial'>ability to identify peaks with different sigma</span></li>
</ul>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'><a
href="http://root.cern.ch/root/html/ListOtransTypes.html#Int_t">Int_t</a> <a
name="TSpectrum:SearchHighRes"></a><a
href="http://root.cern.ch/root/html/src/TSpectrum.cxx.html#TSpectrum:SearchHighRes">SearchHighRes</a>(<a
href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a> *source,<a
href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a> *destVector, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a> ssize, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#float">float</a> sigma, <a
href="http://root.cern.ch/root/html/ListOtransTypes.html#double">double</a> threshold,
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#bool">bool</a> backgroundRemove,<a
href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a> deconIterations,
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#bool">bool</a> markov,
<a href="http://root.cern.ch/root/html/ListOtransTypes.html#int">int</a> averWindow)
</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This
function searches for peaks in source spectrum. It is based on deconvolution
method. First the background is removed (if desired), then Markov smoothed spectrum
is calculated (if desired), then the response function is generated according
to given sigma and deconvolution is carried out. The order of peaks is arranged
according to their heights in the spectrum after background elimination. The
highest peak is the first in the list. On success it returns number of found
peaks.</span></p>
<p class=MsoNormal> </p>
<p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>source</span></b>-pointer to the vector of source
spectrum </p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>destVector</span></b>-resulting spectrum after
deconvolution</p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>ssize</span></b>-length of the source and destination
spectra </p>
<p class=MsoNormal style='text-align:justify'> <b><span
style='font-size:14.0pt'>sigma</span></b>-sigma of searched peaks</p>
<p class=MsoNormal style='margin-left:22.8pt;text-align:justify'><b><span
style='font-size:14.0pt'>threshold</span></b>-<span style='font-size:16.0pt'> </span>threshold
value in % for selected peaks, peaks with amplitude less than
threshold*highest_peak/100 are ignored</p>
<p class=MsoNormal style='margin-left:22.8pt;text-align:justify'><b><span
style='font-size:14.0pt'>backgroundRemove</span></b>-<span style='font-size:
16.0pt'> </span>background_remove-logical variable, true if the removal of
background before deconvolution is desired </p>
<p class=MsoNormal style='margin-left:22.8pt;text-align:justify'><b><span
style='font-size:14.0pt'>deconIterations</span></b>-number of iterations in
deconvolution operation</p>
<p class=MsoNormal style='margin-left:22.8pt;text-align:justify'><b><span
style='font-size:14.0pt'>markov</span></b>-logical variable, if it is true,
first the source spectrum is replaced by new spectrum calculated using Markov
chains method </p>
<p class=MsoNormal style='margin-left:19.95pt;text-align:justify;text-indent:
2.85pt'><b><span style='font-size:14.0pt'>averWindow</span></b>-width of
averaging smoothing window </p>
<p class=MsoNormal> </p>
<p class=MsoNormal><img width=600 height=375 src="gif/TSpectrum_Searching1.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig.
27 An example of one-dimensional synthetic spectrum with found peaks denoted by
markers</span></b></p>
<p class=MsoNormal><b><i><span style='font-size:18.0pt'> </span></i></b></p>
<p class=MsoNormal><b><i><span style='font-size:18.0pt'>References:</span></i></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[1]
M.A. Mariscotti: A method for identification of peaks in the presence of
background and its application to spectrum analysis. NIM 50 (1967), 309-320.</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[2]
</span><span lang=SK style='font-size:16.0pt'> M. Morháč, J. Kliman, V.
Matouek, M. Veselskŭ, I. Turzo</span><span style='font-size:16.0pt'>.:Identification
of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
108-125.</span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[3]
Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
(1996), 451.</span></p>
</div>
<!-- */
// --> End_Html
//Begin_Html <!--
/* -->
<div class=Section19>
<p class=MsoNormal><b><span style='font-size:18.0pt'>Examples of peak searching
method</span></b></p>
<p class=MsoNormal><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'><a
href="http://root.cern.ch/root/html/src/TSpectrum.cxx.html#TSpectrum:SearchHighRes"
target="_parent">SearchHighRes</a> function provides users with the possibility
to vary the input parameters and with the access to the output deconvolved data
in the destination spectrum. Based on the output data one can tune the
parameters. </span></p>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 15 script SearchHR1.c:</span></i></p>
<p class=MsoNormal><span style='font-size:16.0pt'><img width=600 height=321
src="gif/TSpectrum_Searching1.jpg"></span></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig.
28 One-dimensional spectrum with found peaks denoted by markers, 3 iterations
steps in the deconvolution</span></b></p>
<p class=MsoNormal><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><span style='font-size:16.0pt'><img width=600 height=323
src="gif/TSpectrum_Searching2.jpg"></span></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig.
29 One-dimensional spectrum with found peaks denoted by markers, 8 iterations
steps in the deconvolution</span></b></p>
<p class=MsoNormal><span style='font-size:16.0pt'> </span></p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal>// Example to illustrate high resolution peak searching
function (class TSpectrum).</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x SearchHR1.C</p>
<p class=MsoNormal> </p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void SearchHR1() {</p>
<p class=MsoNormal> </p>
<p class=MsoNormal> Float_t fPositionX[100];</p>
<p class=MsoNormal> Float_t fPositionY[100]; </p>
<p class=MsoNormal> Int_t fNPeaks = 0; </p>
<p class=MsoNormal> Int_t i,nfound,bin;</p>
<p class=MsoNormal> Double_t nbins = 1024,a;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> Float_t * dest = new float[nbins]; </p>
<p class=MsoNormal> TH1F *h = new TH1F("h","High resolution
peak searching, number of iterations = 3",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d = new
TH1F("d","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("search2;1"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i
+ 1); </p>
<p class=MsoNormal> TCanvas *Search =
gROOT->GetListOfCanvases()->FindObject("Search");</p>
<p class=MsoNormal> if (!Search) Search = new
TCanvas("Search","Search",10,10,1000,700);</p>
<p class=MsoNormal> h->SetMaximum(4000); </p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> nfound = s->SearchHighRes(source, dest, nbins, 8, 2,
kTRUE, 3, kTRUE, 3);</p>
<p class=MsoNormal> Float_t *xpeaks = s->GetPositionX(); </p>
<p class=MsoNormal> for (i = 0; i < nfound; i++) {</p>
<p class=MsoNormal> a=xpeaks[i];</p>
<p class=MsoNormal> bin = 1 + Int_t(a + 0.5);</p>
<p class=MsoNormal> fPositionX[i] = h->GetBinCenter(bin);</p>
<p class=MsoNormal> fPositionY[i] = h->GetBinContent(bin);</p>
<p class=MsoNormal> }</p>
<p class=MsoNormal> TPolyMarker * pm = (TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");</p>
<p class=MsoNormal> if (pm) {</p>
<p class=MsoNormal> h->GetListOfFunctions()->Remove(pm);</p>
<p class=MsoNormal> delete pm;</p>
<p class=MsoNormal> }</p>
<p class=MsoNormal> pm = new TPolyMarker(nfound, fPositionX, fPositionY);</p>
<p class=MsoNormal> h->GetListOfFunctions()->Add(pm);</p>
<p class=MsoNormal> pm->SetMarkerStyle(23);</p>
<p class=MsoNormal> pm->SetMarkerColor(kRed);</p>
<p class=MsoNormal> pm->SetMarkerSize(1.3);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d->SetBinContent(i +
1,dest[i]);</p>
<p class=MsoNormal> d->SetLineColor(kRed); </p>
<p class=MsoNormal> d->Draw("SAME"); </p>
<p class=MsoNormal> printf("Found %d candidate peaks\n",nfound); </p>
<p class=MsoNormal> for(i=0;i<nfound;i++)</p>
<p class=MsoNormal> printf("posx= %d, posy= %d\n",fPositionX[i],
fPositionY[i]); </p>
<p class=MsoNormal> }</p>
<p class=MsoNormal> </p>
<p class=MsoNormal><i><span style='font-size:16.0pt'>Example 16 script SearchHR3.c:</span></i></p>
<table class=MsoNormalTable border=0 cellspacing=0 cellpadding=0 width=131
style='width:97.9pt;margin-left:131.85pt'>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'>Peak #</p>
</td>
</nobr>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>Position</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>Sigma</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>1</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>118</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>26</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>2</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>162</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>41</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>3</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>310</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>4</nobr></p>
</td>
</tr>
<tr style='height:13.5pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:13.5pt'>
<p class=MsoNormal style='text-align:justify'><nobr>4</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:13.5pt'>
<p class=MsoNormal style='text-align:justify'><nobr>330</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:13.5pt'>
<p class=MsoNormal style='text-align:justify'><nobr>8</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>5</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>482</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>22</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>6</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>491</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>26</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>7</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>740</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>21</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>8</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>852</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>15</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>9</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>954</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>12</nobr></p>
</td>
</tr>
<tr style='height:12.75pt'>
<td width=33 valign=top style='width:24.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>10</nobr></p>
</td>
<td width=54 valign=top style='width:40.85pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>989</nobr></p>
</td>
<td width=43 valign=top style='width:32.2pt;border:solid black 1.0pt;
padding:0mm 0mm 0mm 0mm;height:12.75pt'>
<p class=MsoNormal style='text-align:justify'><nobr>13</nobr></p>
</td>
</tr>
</table>
<p class=MsoNormal> </p>
<p class=MsoNormal><b><span style='font-size:14.0pt'>Table 4 Positions and
sigma of peaks in the following examples</span></b></p>
<p class=MsoNormal><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><img width=600 height=328
src="gif/TSpectrum_Searching3.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:14.0pt'>Fig.
30 Influence of number of iterations (3-red, 10-blue, 100- green,
1000-magenta), sigma=8, smoothing width=3</span></b></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:14.0pt'> </span></b></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:14.0pt'><img
width=600 height=321 src="gif/TSpectrum_Searching4.jpg"></span></b></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:14.0pt'>Fig.
31 Influence of sigma (3-red, 8-blue, 20- green, 43-magenta), num. iter.=10,
sm. width=3</span></b></p>
<p class=MsoNormal style='text-align:justify'><span style='font-size:18.0pt'> </span></p>
<p class=MsoNormal><img width=600 height=323
src="gif/TSpectrum_Searching5.jpg"></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:14.0pt'>Fig.
32 Influence smoothing width (0-red, 3-blue, 7- green, 20-magenta), num.
iter.=10, sigma=8</span></b></p>
<p class=MsoNormal style='text-align:justify'><b><span style='font-size:14.0pt'> </span></b></p>
<p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>// Example to illustrate the influence of number of
iterations in deconvolution in high resolution peak searching function (class
TSpectrum).</p>
<p class=MsoNormal>// To execute this example, do</p>
<p class=MsoNormal>// root > .x SearchHR3.C</p>
<p class=MsoNormal> </p>
<p class=MsoNormal>#include <TSpectrum></p>
<p class=MsoNormal> </p>
<p class=MsoNormal>void SearchHR3() {</p>
<p class=MsoNormal> Float_t fPositionX[100];</p>
<p class=MsoNormal> Float_t fPositionY[100]; </p>
<p class=MsoNormal> Int_t fNPeaks = 0; </p>
<p class=MsoNormal> Int_t i,nfound,bin;</p>
<p class=MsoNormal> Double_t nbins = 1024,a;</p>
<p class=MsoNormal> Double_t xmin = 0;</p>
<p class=MsoNormal> Double_t xmax = (Double_t)nbins;</p>
<p class=MsoNormal> Float_t * source = new float[nbins];</p>
<p class=MsoNormal> Float_t * dest = new float[nbins]; </p>
<p class=MsoNormal> TH1F *h = new TH1F("h","Influence of # of
iterations in deconvolution in peak searching",nbins,xmin,xmax);</p>
<p class=MsoNormal> TH1F *d1 = new
TH1F("d1","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d2 = new
TH1F("d2","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d3 = new
TH1F("d3","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TH1F *d4 = new
TH1F("d4","",nbins,xmin,xmax); </p>
<p class=MsoNormal> TFile *f = new
TFile("spectra\\TSpectrum.root");</p>
<p class=MsoNormal> h=(TH1F*) f->Get("search3;1"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1); </p>
<p class=MsoNormal> TCanvas *Search =
gROOT->GetListOfCanvases()->FindObject("Search");</p>
<p class=MsoNormal> if (!Search) Search = new
TCanvas("Search","Search",10,10,1000,700);</p>
<p class=MsoNormal> h->SetMaximum(1300); </p>
<p class=MsoNormal> h->Draw("L");</p>
<p class=MsoNormal> TSpectrum *s = new TSpectrum();</p>
<p class=MsoNormal> nfound = s->SearchHighRes(source, dest, nbins, 8, 2,
kTRUE, 3, kTRUE, 3);</p>
<p class=MsoNormal> Float_t *xpeaks = s->GetPositionX(); </p>
<p class=MsoNormal> for (i = 0; i < nfound; i++) {</p>
<p class=MsoNormal> a=xpeaks[i];</p>
<p class=MsoNormal> bin = 1 + Int_t(a + 0.5);</p>
<p class=MsoNormal> fPositionX[i] = h->GetBinCenter(bin);</p>
<p class=MsoNormal> fPositionY[i] = h->GetBinContent(bin);</p>
<p class=MsoNormal> } </p>
<p class=MsoNormal> TPolyMarker * pm =
(TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");</p>
<p class=MsoNormal> if (pm) {</p>
<p class=MsoNormal> h->GetListOfFunctions()->Remove(pm);</p>
<p class=MsoNormal> delete pm;</p>
<p class=MsoNormal> }</p>
<p class=MsoNormal> pm = new TPolyMarker(nfound, fPositionX, fPositionY);</p>
<p class=MsoNormal> h->GetListOfFunctions()->Add(pm);</p>
<p class=MsoNormal> pm->SetMarkerStyle(23);</p>
<p class=MsoNormal> pm->SetMarkerColor(kRed);</p>
<p class=MsoNormal> pm->SetMarkerSize(1.3);</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d1->SetBinContent(i +
1,dest[i]);</p>
<p class=MsoNormal> h->Draw("");</p>
<p class=MsoNormal> d1->SetLineColor(kRed); </p>
<p class=MsoNormal> d1->Draw("SAME"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal> s->SearchHighRes(source, dest, nbins, 8, 2, kTRUE, 10,
kTRUE, 3); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d2->SetBinContent(i +
1,dest[i]);</p>
<p class=MsoNormal> d2->SetLineColor(kBlue); </p>
<p class=MsoNormal> d2->Draw("SAME");</p>
<p class=MsoNormal> for (i = 0; i < nbins; i++)
source[i]=h->GetBinContent(i + 1);</p>
<p class=MsoNormal> s->SearchHighRes(source, dest, nbins, 8, 2, kTRUE,
100, kTRUE, 3); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d3->SetBinContent(i +
1,dest[i]);</p>
<p class=MsoNormal> d3->SetLineColor(kGreen); </p>
<p class=MsoNormal> d3->Draw("SAME"); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i
+ 1);</p>
<p class=MsoNormal> s->SearchHighRes(source, dest, nbins, 8, 2, kTRUE,
1000, kTRUE, 3); </p>
<p class=MsoNormal> for (i = 0; i < nbins; i++) d4->SetBinContent(i +
1,dest[i]);</p>
<p class=MsoNormal> d4->SetLineColor(kMagenta); </p>
<p class=MsoNormal> d4->Draw("SAME"); </p>
<p class=MsoNormal> printf("Found %d candidate peaks\n",nfound); </p>
<p class=MsoNormal>}</p>
</div>
<!-- */
// --> End_Html
int i, j, numberIterations = (int)(7 * sigma + 0.5);
double a, b, c;
int k, lindex, posit, imin, imax, jmin, jmax, lh_gold, priz;
double lda, ldb, ldc, area, maximum, maximum_decon;
int xmin, xmax, l, peak_index = 0, size_ext = ssize + 2 * numberIterations, shift = numberIterations, bw = 2, w;
double maxch;
double nom, nip, nim, sp, sm, plocha = 0;
double m0low=0,m1low=0,m2low=0,l0low=0,l1low=0,detlow,av,men;
if (sigma < 1) {
Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
return 0;
}
if(threshold<=0 || threshold>=100){
Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
return 0;
}
j = (int) (5.0 * sigma + 0.5);
if (j >= PEAK_WINDOW / 2) {
Error("SearchHighRes", "Too large sigma");
return 0;
}
if (markov == true) {
if (averWindow <= 0) {
Error("SearchHighRes", "Averanging window must be positive");
return 0;
}
}
if(backgroundRemove == true){
if(ssize < 2 * numberIterations + 1){
Error("SearchHighRes", "Too large clipping window");
return 0;
}
}
k = int(2 * sigma+0.5);
if(k >= 2){
for(i = 0;i < k;i++){
a = i,b = source[i];
m0low += 1,m1low += a,m2low += a * a,l0low += b,l1low += a * b;
}
detlow = m0low * m2low - m1low * m1low;
if(detlow != 0)
l1low = (-l0low * m1low + l1low * m0low) / detlow;
else
l1low = 0;
if(l1low > 0)
l1low=0;
}
else{
l1low = 0;
}
i = (int)(7 * sigma + 0.5);
i = 2 * i;
double *working_space = new double [7 * (ssize + i)];
for (j=0;j<7 * (ssize + i);j++) working_space[j] = 0;
for(i = 0; i < size_ext; i++){
if(i < shift){
a = i - shift;
working_space[i + size_ext] = source[0] + l1low * a;
if(working_space[i + size_ext] < 0)
working_space[i + size_ext]=0;
}
else if(i >= ssize + shift){
a = i - (ssize - 1 + shift);
working_space[i + size_ext] = source[ssize - 1];
if(working_space[i + size_ext] < 0)
working_space[i + size_ext]=0;
}
else
working_space[i + size_ext] = source[i - shift];
}
if(backgroundRemove == true){
for(i = 1; i <= numberIterations; i++){
for(j = i; j < size_ext - i; j++){
if(markov == false){
a = working_space[size_ext + j];
b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
if(b < a)
a = b;
working_space[j]=a;
}
else{
a = working_space[size_ext + j];
av = 0;
men = 0;
for (w = j - bw; w <= j + bw; w++){
if ( w >= 0 && w < size_ext){
av += working_space[size_ext + w];
men +=1;
}
}
av = av / men;
b = 0;
men = 0;
for (w = j - i - bw; w <= j - i + bw; w++){
if ( w >= 0 && w < size_ext){
b += working_space[size_ext + w];
men +=1;
}
}
b = b / men;
c = 0;
men = 0;
for (w = j + i - bw; w <= j + i + bw; w++){
if ( w >= 0 && w < size_ext){
c += working_space[size_ext + w];
men +=1;
}
}
c = c / men;
b = (b + c) / 2;
if (b < a)
av = b;
working_space[j]=av;
}
}
for(j = i; j < size_ext - i; j++)
working_space[size_ext + j] = working_space[j];
}
for(j = 0;j < size_ext; j++){
if(j < shift){
a = j - shift;
b = source[0] + l1low * a;
if(b < 0)
b = 0;
working_space[size_ext + j] = b - working_space[size_ext + j];
}
else if(j >= ssize + shift){
a = j - (ssize - 1 + shift);
b = source[ssize - 1];
if(b < 0)
b = 0;
working_space[size_ext + j] = b - working_space[size_ext + j];
}
else{
working_space[size_ext + j] = source[j - shift] - working_space[size_ext + j];
}
}
for(j = 0;j < size_ext; j++){
if(working_space[size_ext + j] < 0)
working_space[size_ext + j] = 0;
}
}
for(i = 0; i < size_ext; i++){
working_space[i + 6*size_ext] = working_space[i + size_ext];
}
if(markov == true){
for(j = 0; j < size_ext; j++)
working_space[2 * size_ext + j] = working_space[size_ext + j];
xmin = 0,xmax = size_ext - 1;
for(i = 0, maxch = 0; i < size_ext; i++){
working_space[i] = 0;
if(maxch < working_space[2 * size_ext + i])
maxch = working_space[2 * size_ext + i];
plocha += working_space[2 * size_ext + i];
}
if(maxch == 0) {
delete [] working_space;
return 0;
}
nom = 1;
working_space[xmin] = 1;
for(i = xmin; i < xmax; i++){
nip = working_space[2 * size_ext + i] / maxch;
nim = working_space[2 * size_ext + i + 1] / maxch;
sp = 0,sm = 0;
for(l = 1; l <= averWindow; l++){
if((i + l) > xmax)
a = working_space[2 * size_ext + xmax] / maxch;
else
a = working_space[2 * size_ext + i + l] / maxch;
b = a - nip;
if(a + nip <= 0)
a=1;
else
a = TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
sp = sp + b;
if((i - l + 1) < xmin)
a = working_space[2 * size_ext + xmin] / maxch;
else
a = working_space[2 * size_ext + i - l + 1] / maxch;
b = a - nim;
if(a + nim <= 0)
a = 1;
else
a = TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
sm = sm + b;
}
a = sp / sm;
a = working_space[i + 1] = working_space[i] * a;
nom = nom + a;
}
for(i = xmin; i <= xmax; i++){
working_space[i] = working_space[i] / nom;
}
for(j = 0; j < size_ext; j++)
working_space[size_ext + j] = working_space[j] * plocha;
for(j = 0; j < size_ext; j++){
working_space[2 * size_ext + j] = working_space[size_ext + j];
}
if(backgroundRemove == true){
for(i = 1; i <= numberIterations; i++){
for(j = i; j < size_ext - i; j++){
a = working_space[size_ext + j];
b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
if(b < a)
a = b;
working_space[j] = a;
}
for(j = i; j < size_ext - i; j++)
working_space[size_ext + j] = working_space[j];
}
for(j = 0; j < size_ext; j++){
working_space[size_ext + j] = working_space[2 * size_ext + j] - working_space[size_ext + j];
}
}
}
area = 0;
lh_gold = -1;
posit = 0;
maximum = 0;
for(i = 0; i < size_ext; i++){
lda = (double)i - 3 * sigma;
lda = lda * lda / (2 * sigma * sigma);
j = (int)(1000 * TMath::Exp(-lda));
lda = j;
if(lda != 0)
lh_gold = i + 1;
working_space[i] = lda;
area = area + lda;
if(lda > maximum){
maximum = lda;
posit = i;
}
}
for(i = 0; i < size_ext; i++)
working_space[2 * size_ext + i] = TMath::Abs(working_space[size_ext + i]);
i = lh_gold - 1;
if(i > size_ext)
i = size_ext;
imin = -i,imax = i;
for(i = imin; i <= imax; i++){
lda = 0;
jmin = 0;
if(i < 0)
jmin = -i;
jmax = lh_gold - 1 - i;
if(jmax > (lh_gold - 1))
jmax = lh_gold - 1;
for(j = jmin;j <= jmax; j++){
ldb = working_space[j];
ldc = working_space[i + j];
lda = lda + ldb * ldc;
}
working_space[size_ext + i - imin] = lda;
}
i = lh_gold - 1;
imin = -i,imax = size_ext + i - 1;
for(i = imin; i <= imax; i++){
lda = 0;
for(j = 0; j <= (lh_gold - 1); j++){
ldb = working_space[j];
k = i + j;
if(k >= 0 && k < size_ext){
ldc = working_space[2 * size_ext + k];
lda = lda + ldb * ldc;
}
}
working_space[4 * size_ext + i - imin] = lda;
}
for(i = imin; i <= imax; i++)
working_space[2 * size_ext + i - imin] = working_space[4 * size_ext + i - imin];
for(i = 0; i < size_ext; i++)
working_space[i] = 1;
for(lindex = 0; lindex < deconIterations; lindex++){
for(i = 0; i < size_ext; i++){
if(TMath::Abs(working_space[2 * size_ext + i]) > 0.00001 && TMath::Abs(working_space[i]) > 0.00001){
lda=0;
jmin = lh_gold - 1;
if(jmin > i)
jmin = i;
jmin = -jmin;
jmax = lh_gold - 1;
if(jmax > (size_ext - 1 - i))
jmax=size_ext-1-i;
for(j = jmin; j <= jmax; j++){
ldb = working_space[j + lh_gold - 1 + size_ext];
ldc = working_space[i + j];
lda = lda + ldb * ldc;
}
ldb = working_space[2 * size_ext + i];
if(lda != 0)
lda = ldb / lda;
else
lda = 0;
ldb = working_space[i];
lda = lda * ldb;
working_space[3 * size_ext + i] = lda;
}
}
for(i = 0; i < size_ext; i++){
working_space[i] = working_space[3 * size_ext + i];
}
}
for(i=0;i<size_ext;i++){
lda = working_space[i];
j = i + posit;
j = j % size_ext;
working_space[size_ext + j] = lda;
}
maximum = 0, maximum_decon = 0;
j = lh_gold - 1;
for(i = 0; i < size_ext - j; i++){
if(i >= shift && i < ssize + shift){
working_space[i] = area * working_space[size_ext + i + j];
if(maximum_decon < working_space[i])
maximum_decon = working_space[i];
if(maximum < working_space[6 * size_ext + i])
maximum = working_space[6 * size_ext + i];
}
else
working_space[i] = 0;
}
lda=1;
if(lda>threshold)
lda=threshold;
lda=lda/100;
for(i = 1; i < size_ext - 1; i++){
if(working_space[i] > working_space[i - 1] && working_space[i] > working_space[i + 1]){
if(i >= shift && i < ssize + shift){
if(working_space[i] > lda*maximum_decon && working_space[6 * size_ext + i] > threshold * maximum / 100.0){
for(j = i - 1, a = 0, b = 0; j <= i + 1; j++){
a += (double)(j - shift) * working_space[j];
b += working_space[j];
}
a = a / b;
if(a < 0)
a = 0;
if(a >= ssize)
a = ssize - 1;
if(peak_index == 0){
fPositionX[0] = a;
peak_index = 1;
}
else{
for(j = 0, priz = 0; j < peak_index && priz == 0; j++){
if(working_space[6 * size_ext + shift + (int)a] > working_space[6 * size_ext + shift + (int)fPositionX[j]])
priz = 1;
}
if(priz == 0){
if(j < fMaxPeaks){
fPositionX[j] = a;
}
}
else{
for(k = peak_index; k >= j; k--){
if(k < fMaxPeaks){
fPositionX[k] = fPositionX[k - 1];
}
}
fPositionX[j - 1] = a;
}
if(peak_index < fMaxPeaks)
peak_index += 1;
}
}
}
}
}
for(i = 0; i < ssize; i++) destVector[i] = working_space[i + shift];
delete [] working_space;
fNPeaks = peak_index;
if(peak_index == fMaxPeaks)
Warning("SearchHighRes", "Peak buffer full");
return fNPeaks;
}
Int_t TSpectrum::Search1HighRes(float *source,float *destVector, int ssize,
float sigma, double threshold,
bool backgroundRemove,int deconIterations,
bool markov, int averWindow)
{
return SearchHighRes(source,destVector,ssize,sigma,threshold,backgroundRemove,
deconIterations,markov,averWindow);
}
Int_t TSpectrum::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
{
TSpectrum s;
return s.Search(hist,sigma,option,threshold);
}
TH1 *TSpectrum::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
{
TSpectrum s;
return s.Background(hist,niter,option);
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.