ROOT  6.06/09
Reference Guide
TSpectrum.cxx
Go to the documentation of this file.
1 // @(#)root/spectrum:$Id$
2 // Author: Miroslav Morhac 27/05/99
3 
4 #include "TSpectrum.h"
5 #include "TPolyMarker.h"
6 #include "TVirtualPad.h"
7 #include "TList.h"
8 #include "TH1.h"
9 #include "TMath.h"
10 
11 
12 ////////////////////////////////////////////////////////////////////////////////
13 
14 ////////////////////////////////////////////////////////////////////////////////
15 
16 /** \class TSpectrum
17  \ingroup Hist
18  \brief Advanced Spectra Processing
19 
20  ## Advanced spectra processing
21 
22  This class contains advanced spectra processing functions for:
23 
24  - One-dimensional background estimation
25  - One-dimensional smoothing
26  - One-dimensional deconvolution
27  - One-dimensional peak search
28 
29  Author:
30 
31  Miroslav Morhac
32  Institute of Physics
33  Slovak Academy of Sciences
34  Dubravska cesta 9, 842 28 BRATISLAVA
35  SLOVAKIA
36  email:fyzimiro@savba.sk, fax:+421 7 54772479
37 
38  The original code in C has been repackaged as a C++ class by R.Brun.
39 
40  The algorithms in this class have been published in the following references:
41 
42  1. M.Morhac et al.: Background elimination methods for multidimensional coincidence gamma-ray spectra. Nuclear Instruments and Methods in Physics Research A 401 (1997) 113-132.
43  2. M.Morhac et al.: Efficient one- and two-dimensional Gold deconvolution and its application to gamma-ray spectra decomposition. Nuclear Instruments and Methods in Physics Research A 401 (1997) 385-408.
44  3. M.Morhac et al.: Identification of peaks in multidimensional coincidence gamma-ray spectra. Nuclear Instruments and Methods in Research Physics A 443(2000), 108-125.
45 
46  These NIM papers are also available as doc or ps files from:
47 
48  - [Spectrum.doc](ftp://root.cern.ch/root/Spectrum.doc)
49 
50  - [SpectrumDec.ps.gz](ftp://root.cern.ch/root/SpectrumDec.ps.gz)
51 
52  - [SpectrumSrc.ps.gz](ftp://root.cern.ch/root/SpectrumSrc.ps.gz)
53 
54  - [SpectrumBck.ps.gz](ftp://root.cern.ch/root/SpectrumBck.ps.gz)
55 */
56 
57 Int_t TSpectrum::fgIterations = 3;
58 Int_t TSpectrum::fgAverageWindow = 3;
59 
60 #define PEAK_WINDOW 1024
61 ClassImp(TSpectrum)
62 
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 
66 TSpectrum::TSpectrum() :TNamed("Spectrum", "Miroslav Morhac peak finder")
67 {
68  /* Begin_Html
69  Constructor.
70  End_Html */
71 
72  Int_t n = 100;
73  fMaxPeaks = n;
74  fPosition = new Double_t[n];
75  fPositionX = new Double_t[n];
76  fPositionY = new Double_t[n];
77  fResolution = 1;
78  fHistogram = 0;
79  fNPeaks = 0;
80 }
81 
82 
83 ////////////////////////////////////////////////////////////////////////////////
84 
85 TSpectrum::TSpectrum(Int_t maxpositions, Double_t resolution)
86  :TNamed("Spectrum", "Miroslav Morhac peak finder")
87 {
88  /* Begin_Html
89  <ul>
90  <li> maxpositions: maximum number of peaks
91  <li> resolution: determines resolution of the neighboring peaks
92  default value is 1 correspond to 3 sigma distance
93  between peaks. Higher values allow higher resolution
94  (smaller distance between peaks.
95  May be set later through SetResolution.
96  </ul>
97  End_Html */
98 
99  Int_t n = maxpositions;
100  if (n <= 0) n = 1;
101  fMaxPeaks = n;
102  fPosition = new Double_t[n];
103  fPositionX = new Double_t[n];
104  fPositionY = new Double_t[n];
105  fHistogram = 0;
106  fNPeaks = 0;
107  SetResolution(resolution);
108 }
109 
110 
111 ////////////////////////////////////////////////////////////////////////////////
112 
113 TSpectrum::~TSpectrum()
114 {
115  /* Begin_Html
116  Destructor.
117  End_Html */
118 
119  delete [] fPosition;
120  delete [] fPositionX;
121  delete [] fPositionY;
122  delete fHistogram;
123 }
124 
125 
126 ////////////////////////////////////////////////////////////////////////////////
127 
128 void TSpectrum::SetAverageWindow(Int_t w)
129 {
130  /* Begin_Html
131  Static function: Set average window of searched peaks
132  (see TSpectrum::SearchHighRes).
133  End_Html */
134 
135  fgAverageWindow = w;
136 }
137 
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 
141 void TSpectrum::SetDeconIterations(Int_t n)
142 {
143  /* Begin_Html
144  Static function: Set max number of decon iterations in deconvolution
145  operation (see TSpectrum::SearchHighRes).
146  End_Html */
147 
148  fgIterations = n;
149 }
150 
151 
152 ////////////////////////////////////////////////////////////////////////////////
153 
154 TH1 *TSpectrum::Background(const TH1 * h, int numberIterations,
155  Option_t * option)
156 {
157  /* Begin_Html
158  <b>One-dimensional background estimation function.</b>
159  <p>
160  This function calculates the background spectrum in the input histogram h.
161  The background is returned as a histogram.
162  <p>
163  Function parameters:
164  <ul>
165  <li> h: input 1-d histogram
166  <li> numberIterations, (default value = 20)
167  Increasing numberIterations make the result smoother and lower.
168  <li> option: may contain one of the following options:
169  <ul>
170  <li> to set the direction parameter
171  "BackIncreasingWindow". By default the direction is BackDecreasingWindow
172  <li> filterOrder-order of clipping filter, (default "BackOrder2")
173  -possible values= "BackOrder4"
174  "BackOrder6"
175  "BackOrder8"
176  <li> "nosmoothing"- if selected, the background is not smoothed
177  By default the background is smoothed.
178  <li> smoothWindow-width of smoothing window, (default is "BackSmoothing3")
179  -possible values= "BackSmoothing5"
180  "BackSmoothing7"
181  "BackSmoothing9"
182  "BackSmoothing11"
183  "BackSmoothing13"
184  "BackSmoothing15"
185  <li> "Compton" if selected the estimation of Compton edge
186  will be included.
187  <li> "same" : if this option is specified, the resulting background
188  histogram is superimposed on the picture in the current pad.
189  </ul>
190  </ul>
191  NOTE that the background is only evaluated in the current range of h.
192  ie, if h has a bin range (set via h->GetXaxis()->SetRange(binmin,binmax),
193  the returned histogram will be created with the same number of bins
194  as the input histogram h, but only bins from binmin to binmax will be filled
195  with the estimated background.
196  End_Html */
197 
198  if (h == 0) return 0;
199  Int_t dimension = h->GetDimension();
200  if (dimension > 1) {
201  Error("Search", "Only implemented for 1-d histograms");
202  return 0;
203  }
204  TString opt = option;
205  opt.ToLower();
206 
207  //set options
208  Int_t direction = kBackDecreasingWindow;
209  if (opt.Contains("backincreasingwindow")) direction = kBackIncreasingWindow;
210  Int_t filterOrder = kBackOrder2;
211  if (opt.Contains("backorder4")) filterOrder = kBackOrder4;
212  if (opt.Contains("backorder6")) filterOrder = kBackOrder6;
213  if (opt.Contains("backorder8")) filterOrder = kBackOrder8;
214  Bool_t smoothing = kTRUE;
215  if (opt.Contains("nosmoothing")) smoothing = kFALSE;
216  Int_t smoothWindow = kBackSmoothing3;
217  if (opt.Contains("backsmoothing5")) smoothWindow = kBackSmoothing5;
218  if (opt.Contains("backsmoothing7")) smoothWindow = kBackSmoothing7;
219  if (opt.Contains("backsmoothing9")) smoothWindow = kBackSmoothing9;
220  if (opt.Contains("backsmoothing11")) smoothWindow = kBackSmoothing11;
221  if (opt.Contains("backsmoothing13")) smoothWindow = kBackSmoothing13;
222  if (opt.Contains("backsmoothing15")) smoothWindow = kBackSmoothing15;
223  Bool_t compton = kFALSE;
224  if (opt.Contains("compton")) compton = kTRUE;
225 
226  Int_t first = h->GetXaxis()->GetFirst();
227  Int_t last = h->GetXaxis()->GetLast();
228  Int_t size = last-first+1;
229  Int_t i;
230  Double_t * source = new Double_t[size];
231  for (i = 0; i < size; i++) source[i] = h->GetBinContent(i + first);
232 
233  //find background (source is input and in output contains the background
234  Background(source,size,numberIterations, direction, filterOrder,smoothing,
235  smoothWindow,compton);
236 
237  //create output histogram containing backgound
238  //only bins in the range of the input histogram are filled
239  Int_t nch = strlen(h->GetName());
240  char *hbname = new char[nch+20];
241  snprintf(hbname,nch+20,"%s_background",h->GetName());
242  TH1 *hb = (TH1*)h->Clone(hbname);
243  hb->Reset();
244  hb->GetListOfFunctions()->Delete();
245  hb->SetLineColor(2);
246  for (i=0; i< size; i++) hb->SetBinContent(i+first,source[i]);
247  hb->SetEntries(size);
248 
249  //if option "same is specified, draw the result in the pad
250  if (opt.Contains("same")) {
251  if (gPad) delete gPad->GetPrimitive(hbname);
252  hb->Draw("same");
253  }
254  delete [] source;
255  delete [] hbname;
256  return hb;
257 }
258 
259 
260 ////////////////////////////////////////////////////////////////////////////////
261 
262 void TSpectrum::Print(Option_t *) const
263 {
264  /* Begin_Html
265  Print the array of positions.
266  End_Html */
267 
268  printf("\nNumber of positions = %d\n",fNPeaks);
269  for (Int_t i=0;i<fNPeaks;i++) {
270  printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
271  }
272 }
273 
274 
275 ////////////////////////////////////////////////////////////////////////////////
276 
277 Int_t TSpectrum::Search(const TH1 * hin, Double_t sigma, Option_t * option,
278  Double_t threshold)
279 {
280  /* Begin_Html
281  <b>One-dimensional peak search function</b>
282  <p>
283  This function searches for peaks in source spectrum in hin
284  The number of found peaks and their positions are written into
285  the members fNpeaks and fPositionX.
286  The search is performed in the current histogram range.
287  <p>
288  Function parameters:
289  <ul>
290  <li> hin: pointer to the histogram of source spectrum
291  <li> sigma: sigma of searched peaks, for details we refer to manual
292  <li> threshold: (default=0.05) peaks with amplitude less than
293  threshold*highest_peak are discarded. 0<threshold<1
294  </ul>
295  By default, the background is removed before deconvolution.
296  Specify the option "nobackground" to not remove the background.
297  <p>
298  By default the "Markov" chain algorithm is used.
299  Specify the option "noMarkov" to disable this algorithm
300  Note that by default the source spectrum is replaced by a new spectrum
301  <p>
302  By default a polymarker object is created and added to the list of
303  functions of the histogram. The histogram is drawn with the specified
304  option and the polymarker object drawn on top of the histogram.
305  The polymarker coordinates correspond to the npeaks peaks found in
306  the histogram.
307  <p>
308  A pointer to the polymarker object can be retrieved later via:
309  <pre>
310  TList *functions = hin->GetListOfFunctions();
311  TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
312  </pre>
313  Specify the option "goff" to disable the storage and drawing of the
314  polymarker.
315  <p>
316  To disable the final drawing of the histogram with the search results (in case
317  you want to draw it yourself) specify "nodraw" in the options parameter.
318  End_Html */
319 
320  if (hin == 0) return 0;
321  Int_t dimension = hin->GetDimension();
322  if (dimension > 2) {
323  Error("Search", "Only implemented for 1-d and 2-d histograms");
324  return 0;
325  }
326  if (threshold <=0 || threshold >= 1) {
327  Warning("Search","threshold must 0<threshold<1, threshol=0.05 assumed");
328  threshold = 0.05;
329  }
330  TString opt = option;
331  opt.ToLower();
332  Bool_t background = kTRUE;
333  if (opt.Contains("nobackground")) {
334  background = kFALSE;
335  opt.ReplaceAll("nobackground","");
336  }
337  Bool_t markov = kTRUE;
338  if (opt.Contains("nomarkov")) {
339  markov = kFALSE;
340  opt.ReplaceAll("nomarkov","");
341  }
342  Bool_t draw = kTRUE;
343  if (opt.Contains("nodraw")) {
344  draw = kFALSE;
345  opt.ReplaceAll("nodraw","");
346  }
347  if (dimension == 1) {
348  Int_t first = hin->GetXaxis()->GetFirst();
349  Int_t last = hin->GetXaxis()->GetLast();
350  Int_t size = last-first+1;
351  Int_t i, bin, npeaks;
352  Double_t * source = new Double_t[size];
353  Double_t * dest = new Double_t[size];
354  for (i = 0; i < size; i++) source[i] = hin->GetBinContent(i + first);
355  if (sigma < 1) {
356  sigma = size/fMaxPeaks;
357  if (sigma < 1) sigma = 1;
358  if (sigma > 8) sigma = 8;
359  }
360  npeaks = SearchHighRes(source, dest, size, sigma, 100*threshold,
361  background, fgIterations, markov, fgAverageWindow);
362 
363  for (i = 0; i < npeaks; i++) {
364  bin = first + Int_t(fPositionX[i] + 0.5);
365  fPositionX[i] = hin->GetBinCenter(bin);
366  fPositionY[i] = hin->GetBinContent(bin);
367  }
368  delete [] source;
369  delete [] dest;
370 
371  if (opt.Contains("goff"))
372  return npeaks;
373  if (!npeaks) return 0;
374  TPolyMarker * pm =
375  (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
376  if (pm) {
377  hin->GetListOfFunctions()->Remove(pm);
378  delete pm;
379  }
380  pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
381  hin->GetListOfFunctions()->Add(pm);
382  pm->SetMarkerStyle(23);
383  pm->SetMarkerColor(kRed);
384  pm->SetMarkerSize(1.3);
385  opt.ReplaceAll(" ","");
386  opt.ReplaceAll(",","");
387  if (draw)
388  ((TH1*)hin)->Draw(opt.Data());
389  return npeaks;
390  }
391  return 0;
392 }
393 
394 
395 ////////////////////////////////////////////////////////////////////////////////
396 
397 void TSpectrum::SetResolution(Double_t resolution)
398 {
399  /* Begin_Html
400  resolution: determines resolution of the neighboring peaks
401  default value is 1 correspond to 3 sigma distance
402  between peaks. Higher values allow higher resolution
403  (smaller distance between peaks.
404  May be set later through SetResolution.
405  End_Html */
406 
407  if (resolution > 1)
408  fResolution = resolution;
409  else
410  fResolution = 1;
411 }
412 
413 
414 ////////////////////////////////////////////////////////////////////////////////
415 
416 const char *TSpectrum::Background(Double_t *spectrum, int ssize,
417  int numberIterations,
418  int direction, int filterOrder,
419  bool smoothing,int smoothWindow,
420  bool compton)
421 {
422 /* Begin_Html
423 This function calculates background spectrum from source spectrum.
424 The result is placed in the vector pointed by spe1945ctrum pointer.
425 The goal is to separate the useful information (peaks) from useless
426 information (background).
427 
428 <ul>
429 <li> method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping
430  algorithm.
431 <li> new value in the channel "i" is calculated
432 </ul>
433 
434 <img width=486 height=72 src="gif/TSpectrum_Background.gif">
435 
436 where p = 1, 2, ..., numberIterations. In fact it represents second order
437 difference filter (-1,2,-1).
438 
439 One can also change the
440 direction of the change of the clipping window, the order of the clipping
441 filter, to include smoothing, to set width of smoothing window and to include
442 the estimation of Compton edges. On successful completion it returns 0. On
443 error it returns pointer to the string describing error.
444 
445 <h4>Parameters:</h4>
446 <ul>
447 <li> spectrum: pointer to the vector of source spectrum
448 <li> ssize: length of the spectrum vector
449 <li> numberIterations: maximal width of clipping window,
450 <li> direction: direction of change of clipping window.
451  Possible values: kBackIncreasingWindow, kBackDecreasingWindow
452 <li> filterOrder: order of clipping filter.
453  Possible values: kBackOrder2, kBackOrder4, kBackOrder6, kBackOrder8
454 <li> smoothing: logical variable whether the smoothing operation in the
455  estimation of background will be included.
456  Possible values: kFALSE, kTRUE
457 <li> smoothWindow: width of smoothing window.
458  Possible values: kBackSmoothing3, kBackSmoothing5, kBackSmoothing7,
459  kBackSmoothing9, kBackSmoothing11, kBackSmoothing13, kBackSmoothing15.
460 <li> compton: logical variable whether the estimation of Compton edge will be
461  included. Possible values: kFALSE, kTRUE.
462 </ul>
463 
464 
465 <h4>References:</h4>
466 <ol>
467 <li> C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the
468 quantitative analysis of PIXE spectra in geoscience applications. NIM, B34
469 (1988), 396-402.
470 
471 <li> M. Morháč, J. Kliman, V. Matoušek, M. Veselský, I. Turzo:
472 Background elimination methods for multidimensional gamma-ray spectra. NIM,
473 A401 (1997) 113-132.
474 
475 <li> D. D. Burgess, R. J. Tervo: Background estimation for gamma-ray
476 spectroscopy. NIM 214 (1983), 431-434.
477 </ol>
478 
479 Example 1 script Background_incr.c:
480 <p>
481 <img width=601 height=407 src="gif/TSpectrum_Background_incr.jpg">
482 <p>
483 Figure 1 Example of the estimation of background for number of iterations=6.
484 Original spectrum is shown in black color, estimated background in red color.
485 <p>
486 Script:
487 <pre>
488 // Example to illustrate the background estimator (class TSpectrum).
489 // To execute this example, do
490 // root > .x Background_incr.C
491 
492 #include <TSpectrum>
493 
494 void Background_incr() {
495  Int_t i;
496  Double_t nbins = 256;
497  Double_t xmin = 0;
498  Double_t xmax = nbins;
499  Double_t * source = new Double_t[nbins];
500  TH1F *back = new TH1F("back","",nbins,xmin,xmax);
501  TH1F *d = new TH1F("d","",nbins,xmin,xmax);
502  TFile *f = new TFile("spectra\\TSpectrum.root");
503  back=(TH1F*) f->Get("back1;1");
504  TCanvas *Background = gROOT->GetListOfCanvases()->FindObject("Background");
505  if (!Background) Background =
506  new TCanvas("Background",
507  "Estimation of background with increasing window",
508  10,10,1000,700);
509  back->Draw("L");
510  TSpectrum *s = new TSpectrum();
511  for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
512  s->Background(source,nbins,6,kBackIncreasingWindow,kBackOrder2,kFALSE,
513  kBackSmoothing3,kFALSE);
514  for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
515  d->SetLineColor(kRed);
516  d->Draw("SAME L");
517 }
518 </pre>
519 
520 Example 2 script Background_decr.c:
521 <p>
522 In Figure 1. one can notice that at the edges of the peaks the estimated
523 background goes under the peaks. An alternative approach is to decrease the
524 clipping window from a given value numberIterations to the value of one, which
525 is presented in this example.
526 <p>
527 <img width=601 height=407 src="gif/TSpectrum_Background_decr.jpg">
528 <p>
529 Figure 2 Example of the estimation of background for numberIterations=6 using
530 decreasing clipping window algorithm. Original spectrum is shown in black
531 color, estimated background in red color.
532 <p>
533 Script:
534 
535 <pre>
536 // Example to illustrate the background estimator (class TSpectrum).
537 // To execute this example, do
538 // root > .x Background_decr.C
539 
540 #include <TSpectrum>
541 
542 void Background_decr() {
543  Int_t i;
544  Double_t nbins = 256;
545  Double_t xmin = 0;
546  Double_t xmax = nbins;
547  Double_t * source = new Double_t[nbins];
548  TH1F *back = new TH1F("back","",nbins,xmin,xmax);
549  TH1F *d = new TH1F("d","",nbins,xmin,xmax);
550  TFile *f = new TFile("spectra\\TSpectrum.root");
551  back=(TH1F*) f->Get("back1;1");
552  TCanvas *Background = gROOT->GetListOfCanvases()->FindObject("Background");
553  if (!Background) Background =
554  new TCanvas("Background","Estimation of background with decreasing window",
555  10,10,1000,700);
556  back->Draw("L");
557  TSpectrum *s = new TSpectrum();
558  for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
559  s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE,
560  kBackSmoothing3,kFALSE);
561  for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
562  d->SetLineColor(kRed);
563  d->Draw("SAME L");
564 }
565 </pre>
566 
567 Example 3 script Background_width.c:
568 <p>
569 The question is how to choose the width of the clipping window, i.e.,
570 numberIterations parameter. The influence of this parameter on the estimated
571 background is illustrated in Figure 3.
572 <p>
573 <img width=601 height=407 src="gif/TSpectrum_Background_width.jpg">
574 <p>
575 Figure 3 Example of the influence of clipping window width on the estimated
576 background for numberIterations=4 (red line), 6 (blue line) 8 (green line) using
577 decreasing clipping window algorithm.
578 
579 <p>
580 in general one should set this parameter so that the value
581 2*numberIterations+1 was greater than the widths of preserved objects (peaks).
582 
583 <p>
584 Script:
585 
586 <pre>
587 // Example to illustrate the influence of the clipping window width on the
588 // estimated background. To execute this example, do:
589 // root > .x Background_width.C
590 
591 #include <TSpectrum>
592 
593 void Background_width() {
594  Int_t i;
595  Double_t nbins = 256;
596  Double_t xmin = 0;
597  Double_t xmax = nbins;
598  Double_t * source = new Double_t[nbins];
599  TH1F *h = new TH1F("h","",nbins,xmin,xmax);
600  TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
601  TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax);
602  TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax);
603  TFile *f = new TFile("spectra\\TSpectrum.root");
604  h=(TH1F*) f->Get("back1;1");
605  TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background");
606  if (!background) background = new TCanvas("background",
607  "Influence of clipping window width on the estimated background",
608  10,10,1000,700);
609  h->Draw("L");
610  TSpectrum *s = new TSpectrum();
611  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
612  s->Background(source,nbins,4,kBackDecreasingWindow,kBackOrder2,kFALSE,
613  kBackSmoothing3,kFALSE);
614  for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]);
615  d1->SetLineColor(kRed);
616  d1->Draw("SAME L");
617  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
618  s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE,
619  kBackSmoothing3,kFALSE);
620  for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]);
621  d2->SetLineColor(kBlue);
622  d2->Draw("SAME L");
623  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
624  s->Background(source,nbins,8,kBackDecreasingWindow,kBackOrder2,kFALSE,
625  kBackSmoothing3,kFALSE);
626  for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]);
627  d3->SetLineColor(kGreen);
628  d3->Draw("SAME L");
629 }
630 </pre>
631 
632 Example 4 script Background_width2.c:
633 <p>
634 another example for very complex spectrum is given in Figure 4.
635 <p>
636 <img width=601 height=407 src="gif/TSpectrum_Background_width2.jpg">
637 <p>
638 Figure 4 Example of the influence of clipping window width on the estimated
639 background for numberIterations=10 (red line), 20 (blue line), 30 (green line)
640 and 40 (magenta line) using decreasing clipping window algorithm.
641 
642 <p>
643 Script:
644 
645 <pre>
646 // Example to illustrate the influence of the clipping window width on the
647 // estimated background. To execute this example, do:
648 // root > .x Background_width2.C
649 
650 #include <TSpectrum>
651 
652 void Background_width2() {
653  Int_t i;
654  Double_t nbins = 4096;
655  Double_t xmin = 0;
656  Double_t xmax = 4096;
657  Double_t * source = new Double_t[nbins];
658  TH1F *h = new TH1F("h","",nbins,xmin,xmax);
659  TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
660  TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax);
661  TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax);
662  TH1F *d4 = new TH1F("d4","",nbins,xmin,xmax);
663  TFile *f = new TFile("spectra\\TSpectrum.root");
664  h=(TH1F*) f->Get("back2;1");
665  TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background");
666  if (!background) background = new TCanvas("background",
667  "Influence of clipping window width on the estimated background",
668  10,10,1000,700);
669  h->SetAxisRange(0,1000);
670  h->SetMaximum(20000);
671  h->Draw("L");
672  TSpectrum *s = new TSpectrum();
673  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
674  s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder2,kFALSE,
675  kBackSmoothing3,kFALSE);
676  for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]);
677  d1->SetLineColor(kRed);
678  d1->Draw("SAME L");
679  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
680  s->Background(source,nbins,20,kBackDecreasingWindow,kBackOrder2,kFALSE,
681  kBackSmoothing3,kFALSE);
682  for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]);
683  d2->SetLineColor(kBlue);
684  d2->Draw("SAME L");
685  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
686  s->Background(source,nbins,30,kBackDecreasingWindow,kBackOrder2,kFALSE,
687  kBackSmoothing3,kFALSE);
688  for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]);
689  d3->SetLineColor(kGreen);
690  d3->Draw("SAME L");
691  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
692  s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder2,kFALSE,
693  kBackSmoothing3,kFALSE);
694  for (i = 0; i < nbins; i++) d4->SetBinContent(i + 1,source[i]);
695  d4->SetLineColor(kMagenta);
696  d4->Draw("SAME L");
697 }
698 </pre>
699 
700 Example 5 script Background_order.c:
701 <p>
702 Second order difference filter removes linear (quasi-linear) background and
703 preserves symmetrical peaks. However if the shape of the background is more
704 complex one can employ higher-order clipping filters (see example in Figure 5)
705 <p>
706 <img width=601 height=407 src="gif/TSpectrum_Background_order.jpg">
707 <p>
708 Figure 5 Example of the influence of clipping filter difference order on the
709 estimated background for fNnumberIterations=40, 2-nd order red line, 4-th order
710 blue line, 6-th order green line and 8-th order magenta line, and using
711 decreasing clipping window algorithm.
712 <p>
713 Script:
714 <pre>
715 // Example to illustrate the influence of the clipping filter difference order
716 // on the estimated background. To execute this example, do
717 // root > .x Background_order.C
718 
719 #include <TSpectrum>
720 
721 void Background_order() {
722  Int_t i;
723  Double_t nbins = 4096;
724  Double_t xmin = 0;
725  Double_t xmax = 4096;
726  Double_t * source = new Double_t[nbins];
727  TH1F *h = new TH1F("h","",nbins,xmin,xmax);
728  TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
729  TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax);
730  TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax);
731  TH1F *d4 = new TH1F("d4","",nbins,xmin,xmax);
732  TFile *f = new TFile("spectra\\TSpectrum.root");
733  h=(TH1F*) f->Get("back2;1");
734  TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background");
735  if (!background) background = new TCanvas("background",
736  "Influence of clipping filter difference order on the estimated background",
737  10,10,1000,700);
738  h->SetAxisRange(1220,1460);
739  h->SetMaximum(11000);
740  h->Draw("L");
741  TSpectrum *s = new TSpectrum();
742  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
743  s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder2,kFALSE,
744  kBackSmoothing3,kFALSE);
745  for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]);
746  d1->SetLineColor(kRed);
747  d1->Draw("SAME L");
748  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
749  s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder4,kFALSE,
750  kBackSmoothing3,kFALSE);
751  for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]);
752  d2->SetLineColor(kBlue);
753  d2->Draw("SAME L");
754  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
755  s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder6,kFALSE,
756  kBackSmoothing3,kFALSE);
757  for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]);
758  d3->SetLineColor(kGreen);
759  d3->Draw("SAME L");
760  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
761  s->Background(source,nbins,40,kBackDecreasingWindow,kBackOrder8,kFALSE,
762  kBackSmoothing3,kFALSE);
763  for (i = 0; i < nbins; i++) d4->SetBinContent(i + 1,source[i]);
764  d4->SetLineColor(kMagenta);
765  d4->Draw("SAME L");
766 }
767 </pre>
768 
769 Example 6 script Background_smooth.c:
770 <p>
771 The estimate of the background can be influenced by noise present in the
772 spectrum. We proposed the algorithm of the background estimate with
773 simultaneous smoothing. In the original algorithm without smoothing, the
774 estimated background snatches the lower spikes in the noise. Consequently,
775 the areas of peaks are biased by this error.
776 <p>
777 <img width=554 height=104 src="gif/TSpectrum_Background_smooth1.jpg">
778 <p>
779 Figure 7 Principle of background estimation algorithm with simultaneous
780 smoothing.
781 <p>
782 <img width=601 height=407 src="gif/TSpectrum_Background_smooth2.jpg">
783 <p>
784 Figure 8 Illustration of non-smoothing (red line) and smoothing algorithm of
785 background estimation (blue line).
786 
787 <p>
788 
789 Script:
790 
791 <pre>
792 // Example to illustrate the background estimator (class TSpectrum) including
793 // Compton edges. To execute this example, do:
794 // root > .x Background_smooth.C
795 
796 #include <TSpectrum>
797 
798 void Background_smooth() {
799  Int_t i;
800  Double_t nbins = 4096;
801  Double_t xmin = 0;
802  Double_t xmax = nbins;
803  Double_t * source = new Double_t[nbins];
804  TH1F *h = new TH1F("h","",nbins,xmin,xmax);
805  TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
806  TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax);
807  TFile *f = new TFile("spectra\\TSpectrum.root");
808  h=(TH1F*) f->Get("back4;1");
809  TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background");
810  if (!background) background = new TCanvas("background",
811  "Estimation of background with noise",10,10,1000,700);
812  h->SetAxisRange(3460,3830);
813  h->Draw("L");
814  TSpectrum *s = new TSpectrum();
815  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
816  s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kFALSE,
817  kBackSmoothing3,kFALSE);
818  for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]);
819  d1->SetLineColor(kRed);
820  d1->Draw("SAME L");
821  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
822  s->Background(source,nbins,6,kBackDecreasingWindow,kBackOrder2,kTRUE,
823  kBackSmoothing3,kFALSE);
824  for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]);
825  d2->SetLineColor(kBlue);
826  d2->Draw("SAME L");
827 }
828 </pre>
829 
830 Example 8 script Background_compton.c:
831 <p>
832 Sometimes it is necessary to include also the Compton edges into the estimate of
833 the background. In Figure 8 we present the example of the synthetic spectrum
834 with Compton edges. The background was estimated using the 8-th order filter
835 with the estimation of the Compton edges using decreasing
836 clipping window algorithm (numberIterations=10) with smoothing
837 (smoothingWindow=5).
838 <p>
839 <img width=601 height=407 src="gif/TSpectrum_Background_compton.jpg">
840 <p>
841 Figure 8 Example of the estimate of the background with Compton edges (red
842 line) for numberIterations=10, 8-th order difference filter, using decreasing
843 clipping window algorithm and smoothing (smoothingWindow=5).
844 <p>
845 Script:
846 
847 <pre>
848 // Example to illustrate the background estimator (class TSpectrum) including
849 // Compton edges. To execute this example, do:
850 // root > .x Background_compton.C
851 
852 #include <TSpectrum>
853 
854 void Background_compton() {
855  Int_t i;
856  Double_t nbins = 512;
857  Double_t xmin = 0;
858  Double_t xmax = nbins;
859  Double_t * source = new Double_t[nbins];
860  TH1F *h = new TH1F("h","",nbins,xmin,xmax);
861  TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
862  TFile *f = new TFile("spectra\\TSpectrum.root");
863  h=(TH1F*) f->Get("back3;1");
864  TCanvas *background = gROOT->GetListOfCanvases()->FindObject("background");
865  if (!background) background = new TCanvas("background",
866  "Estimation of background with Compton edges under peaks",10,10,1000,700);
867  h->Draw("L");
868  TSpectrum *s = new TSpectrum();
869  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
870  s->Background(source,nbins,10,kBackDecreasingWindow,kBackOrder8,kTRUE,
871  kBackSmoothing5,,kTRUE);
872  for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]);
873  d1->SetLineColor(kRed);
874  d1->Draw("SAME L");
875 }
876 </pre>
877 
878 End_Html */
879 
880  int i, j, w, bw, b1, b2, priz;
881  Double_t 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;
882  if (ssize <= 0)
883  return "Wrong Parameters";
884  if (numberIterations < 1)
885  return "Width of Clipping Window Must Be Positive";
886  if (ssize < 2 * numberIterations + 1)
887  return "Too Large Clipping Window";
888  if (smoothing == kTRUE && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15)
889  return "Incorrect width of smoothing window";
890  Double_t *working_space = new Double_t[2 * ssize];
891  for (i = 0; i < ssize; i++){
892  working_space[i] = spectrum[i];
893  working_space[i + ssize] = spectrum[i];
894  }
895  bw=(smoothWindow-1)/2;
896  if (direction == kBackIncreasingWindow)
897  i = 1;
898  else if(direction == kBackDecreasingWindow)
899  i = numberIterations;
900  if (filterOrder == kBackOrder2) {
901  do{
902  for (j = i; j < ssize - i; j++) {
903  if (smoothing == kFALSE){
904  a = working_space[ssize + j];
905  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
906  if (b < a)
907  a = b;
908  working_space[j] = a;
909  }
910 
911  else if (smoothing == kTRUE){
912  a = working_space[ssize + j];
913  av = 0;
914  men = 0;
915  for (w = j - bw; w <= j + bw; w++){
916  if ( w >= 0 && w < ssize){
917  av += working_space[ssize + w];
918  men +=1;
919  }
920  }
921  av = av / men;
922  b = 0;
923  men = 0;
924  for (w = j - i - bw; w <= j - i + bw; w++){
925  if ( w >= 0 && w < ssize){
926  b += working_space[ssize + w];
927  men +=1;
928  }
929  }
930  b = b / men;
931  c = 0;
932  men = 0;
933  for (w = j + i - bw; w <= j + i + bw; w++){
934  if ( w >= 0 && w < ssize){
935  c += working_space[ssize + w];
936  men +=1;
937  }
938  }
939  c = c / men;
940  b = (b + c) / 2;
941  if (b < a)
942  av = b;
943  working_space[j]=av;
944  }
945  }
946  for (j = i; j < ssize - i; j++)
947  working_space[ssize + j] = working_space[j];
948  if (direction == kBackIncreasingWindow)
949  i+=1;
950  else if(direction == kBackDecreasingWindow)
951  i-=1;
952  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
953  }
954 
955  else if (filterOrder == kBackOrder4) {
956  do{
957  for (j = i; j < ssize - i; j++) {
958  if (smoothing == kFALSE){
959  a = working_space[ssize + j];
960  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
961  c = 0;
962  ai = i / 2;
963  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
964  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
965  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
966  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
967  if (b < c)
968  b = c;
969  if (b < a)
970  a = b;
971  working_space[j] = a;
972  }
973 
974  else if (smoothing == kTRUE){
975  a = working_space[ssize + j];
976  av = 0;
977  men = 0;
978  for (w = j - bw; w <= j + bw; w++){
979  if ( w >= 0 && w < ssize){
980  av += working_space[ssize + w];
981  men +=1;
982  }
983  }
984  av = av / men;
985  b = 0;
986  men = 0;
987  for (w = j - i - bw; w <= j - i + bw; w++){
988  if ( w >= 0 && w < ssize){
989  b += working_space[ssize + w];
990  men +=1;
991  }
992  }
993  b = b / men;
994  c = 0;
995  men = 0;
996  for (w = j + i - bw; w <= j + i + bw; w++){
997  if ( w >= 0 && w < ssize){
998  c += working_space[ssize + w];
999  men +=1;
1000  }
1001  }
1002  c = c / men;
1003  b = (b + c) / 2;
1004  ai = i / 2;
1005  b4 = 0, men = 0;
1006  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1007  if (w >= 0 && w < ssize){
1008  b4 += working_space[ssize + w];
1009  men +=1;
1010  }
1011  }
1012  b4 = b4 / men;
1013  c4 = 0, men = 0;
1014  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1015  if (w >= 0 && w < ssize){
1016  c4 += working_space[ssize + w];
1017  men +=1;
1018  }
1019  }
1020  c4 = c4 / men;
1021  d4 = 0, men = 0;
1022  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1023  if (w >= 0 && w < ssize){
1024  d4 += working_space[ssize + w];
1025  men +=1;
1026  }
1027  }
1028  d4 = d4 / men;
1029  e4 = 0, men = 0;
1030  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1031  if (w >= 0 && w < ssize){
1032  e4 += working_space[ssize + w];
1033  men +=1;
1034  }
1035  }
1036  e4 = e4 / men;
1037  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1038  if (b < b4)
1039  b = b4;
1040  if (b < a)
1041  av = b;
1042  working_space[j]=av;
1043  }
1044  }
1045  for (j = i; j < ssize - i; j++)
1046  working_space[ssize + j] = working_space[j];
1047  if (direction == kBackIncreasingWindow)
1048  i+=1;
1049  else if(direction == kBackDecreasingWindow)
1050  i-=1;
1051  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1052  }
1053 
1054  else if (filterOrder == kBackOrder6) {
1055  do{
1056  for (j = i; j < ssize - i; j++) {
1057  if (smoothing == kFALSE){
1058  a = working_space[ssize + j];
1059  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
1060  c = 0;
1061  ai = i / 2;
1062  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
1063  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
1064  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
1065  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
1066  d = 0;
1067  ai = i / 3;
1068  d += working_space[ssize + j - (Int_t) (3 * ai)] / 20;
1069  d -= 6 * working_space[ssize + j - (Int_t) (2 * ai)] / 20;
1070  d += 15 * working_space[ssize + j - (Int_t) ai] / 20;
1071  d += 15 * working_space[ssize + j + (Int_t) ai] / 20;
1072  d -= 6 * working_space[ssize + j + (Int_t) (2 * ai)] / 20;
1073  d += working_space[ssize + j + (Int_t) (3 * ai)] / 20;
1074  if (b < d)
1075  b = d;
1076  if (b < c)
1077  b = c;
1078  if (b < a)
1079  a = b;
1080  working_space[j] = a;
1081  }
1082 
1083  else if (smoothing == kTRUE){
1084  a = working_space[ssize + j];
1085  av = 0;
1086  men = 0;
1087  for (w = j - bw; w <= j + bw; w++){
1088  if ( w >= 0 && w < ssize){
1089  av += working_space[ssize + w];
1090  men +=1;
1091  }
1092  }
1093  av = av / men;
1094  b = 0;
1095  men = 0;
1096  for (w = j - i - bw; w <= j - i + bw; w++){
1097  if ( w >= 0 && w < ssize){
1098  b += working_space[ssize + w];
1099  men +=1;
1100  }
1101  }
1102  b = b / men;
1103  c = 0;
1104  men = 0;
1105  for (w = j + i - bw; w <= j + i + bw; w++){
1106  if ( w >= 0 && w < ssize){
1107  c += working_space[ssize + w];
1108  men +=1;
1109  }
1110  }
1111  c = c / men;
1112  b = (b + c) / 2;
1113  ai = i / 2;
1114  b4 = 0, men = 0;
1115  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1116  if (w >= 0 && w < ssize){
1117  b4 += working_space[ssize + w];
1118  men +=1;
1119  }
1120  }
1121  b4 = b4 / men;
1122  c4 = 0, men = 0;
1123  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1124  if (w >= 0 && w < ssize){
1125  c4 += working_space[ssize + w];
1126  men +=1;
1127  }
1128  }
1129  c4 = c4 / men;
1130  d4 = 0, men = 0;
1131  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1132  if (w >= 0 && w < ssize){
1133  d4 += working_space[ssize + w];
1134  men +=1;
1135  }
1136  }
1137  d4 = d4 / men;
1138  e4 = 0, men = 0;
1139  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1140  if (w >= 0 && w < ssize){
1141  e4 += working_space[ssize + w];
1142  men +=1;
1143  }
1144  }
1145  e4 = e4 / men;
1146  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1147  ai = i / 3;
1148  b6 = 0, men = 0;
1149  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
1150  if (w >= 0 && w < ssize){
1151  b6 += working_space[ssize + w];
1152  men +=1;
1153  }
1154  }
1155  b6 = b6 / men;
1156  c6 = 0, men = 0;
1157  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1158  if (w >= 0 && w < ssize){
1159  c6 += working_space[ssize + w];
1160  men +=1;
1161  }
1162  }
1163  c6 = c6 / men;
1164  d6 = 0, men = 0;
1165  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1166  if (w >= 0 && w < ssize){
1167  d6 += working_space[ssize + w];
1168  men +=1;
1169  }
1170  }
1171  d6 = d6 / men;
1172  e6 = 0, men = 0;
1173  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1174  if (w >= 0 && w < ssize){
1175  e6 += working_space[ssize + w];
1176  men +=1;
1177  }
1178  }
1179  e6 = e6 / men;
1180  f6 = 0, men = 0;
1181  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1182  if (w >= 0 && w < ssize){
1183  f6 += working_space[ssize + w];
1184  men +=1;
1185  }
1186  }
1187  f6 = f6 / men;
1188  g6 = 0, men = 0;
1189  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1190  if (w >= 0 && w < ssize){
1191  g6 += working_space[ssize + w];
1192  men +=1;
1193  }
1194  }
1195  g6 = g6 / men;
1196  b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1197  if (b < b6)
1198  b = b6;
1199  if (b < b4)
1200  b = b4;
1201  if (b < a)
1202  av = b;
1203  working_space[j]=av;
1204  }
1205  }
1206  for (j = i; j < ssize - i; j++)
1207  working_space[ssize + j] = working_space[j];
1208  if (direction == kBackIncreasingWindow)
1209  i+=1;
1210  else if(direction == kBackDecreasingWindow)
1211  i-=1;
1212  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1213  }
1214 
1215  else if (filterOrder == kBackOrder8) {
1216  do{
1217  for (j = i; j < ssize - i; j++) {
1218  if (smoothing == kFALSE){
1219  a = working_space[ssize + j];
1220  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
1221  c = 0;
1222  ai = i / 2;
1223  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
1224  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
1225  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
1226  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
1227  d = 0;
1228  ai = i / 3;
1229  d += working_space[ssize + j - (Int_t) (3 * ai)] / 20;
1230  d -= 6 * working_space[ssize + j - (Int_t) (2 * ai)] / 20;
1231  d += 15 * working_space[ssize + j - (Int_t) ai] / 20;
1232  d += 15 * working_space[ssize + j + (Int_t) ai] / 20;
1233  d -= 6 * working_space[ssize + j + (Int_t) (2 * ai)] / 20;
1234  d += working_space[ssize + j + (Int_t) (3 * ai)] / 20;
1235  e = 0;
1236  ai = i / 4;
1237  e -= working_space[ssize + j - (Int_t) (4 * ai)] / 70;
1238  e += 8 * working_space[ssize + j - (Int_t) (3 * ai)] / 70;
1239  e -= 28 * working_space[ssize + j - (Int_t) (2 * ai)] / 70;
1240  e += 56 * working_space[ssize + j - (Int_t) ai] / 70;
1241  e += 56 * working_space[ssize + j + (Int_t) ai] / 70;
1242  e -= 28 * working_space[ssize + j + (Int_t) (2 * ai)] / 70;
1243  e += 8 * working_space[ssize + j + (Int_t) (3 * ai)] / 70;
1244  e -= working_space[ssize + j + (Int_t) (4 * ai)] / 70;
1245  if (b < e)
1246  b = e;
1247  if (b < d)
1248  b = d;
1249  if (b < c)
1250  b = c;
1251  if (b < a)
1252  a = b;
1253  working_space[j] = a;
1254  }
1255 
1256  else if (smoothing == kTRUE){
1257  a = working_space[ssize + j];
1258  av = 0;
1259  men = 0;
1260  for (w = j - bw; w <= j + bw; w++){
1261  if ( w >= 0 && w < ssize){
1262  av += working_space[ssize + w];
1263  men +=1;
1264  }
1265  }
1266  av = av / men;
1267  b = 0;
1268  men = 0;
1269  for (w = j - i - bw; w <= j - i + bw; w++){
1270  if ( w >= 0 && w < ssize){
1271  b += working_space[ssize + w];
1272  men +=1;
1273  }
1274  }
1275  b = b / men;
1276  c = 0;
1277  men = 0;
1278  for (w = j + i - bw; w <= j + i + bw; w++){
1279  if ( w >= 0 && w < ssize){
1280  c += working_space[ssize + w];
1281  men +=1;
1282  }
1283  }
1284  c = c / men;
1285  b = (b + c) / 2;
1286  ai = i / 2;
1287  b4 = 0, men = 0;
1288  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1289  if (w >= 0 && w < ssize){
1290  b4 += working_space[ssize + w];
1291  men +=1;
1292  }
1293  }
1294  b4 = b4 / men;
1295  c4 = 0, men = 0;
1296  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1297  if (w >= 0 && w < ssize){
1298  c4 += working_space[ssize + w];
1299  men +=1;
1300  }
1301  }
1302  c4 = c4 / men;
1303  d4 = 0, men = 0;
1304  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1305  if (w >= 0 && w < ssize){
1306  d4 += working_space[ssize + w];
1307  men +=1;
1308  }
1309  }
1310  d4 = d4 / men;
1311  e4 = 0, men = 0;
1312  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1313  if (w >= 0 && w < ssize){
1314  e4 += working_space[ssize + w];
1315  men +=1;
1316  }
1317  }
1318  e4 = e4 / men;
1319  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
1320  ai = i / 3;
1321  b6 = 0, men = 0;
1322  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
1323  if (w >= 0 && w < ssize){
1324  b6 += working_space[ssize + w];
1325  men +=1;
1326  }
1327  }
1328  b6 = b6 / men;
1329  c6 = 0, men = 0;
1330  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1331  if (w >= 0 && w < ssize){
1332  c6 += working_space[ssize + w];
1333  men +=1;
1334  }
1335  }
1336  c6 = c6 / men;
1337  d6 = 0, men = 0;
1338  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1339  if (w >= 0 && w < ssize){
1340  d6 += working_space[ssize + w];
1341  men +=1;
1342  }
1343  }
1344  d6 = d6 / men;
1345  e6 = 0, men = 0;
1346  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1347  if (w >= 0 && w < ssize){
1348  e6 += working_space[ssize + w];
1349  men +=1;
1350  }
1351  }
1352  e6 = e6 / men;
1353  f6 = 0, men = 0;
1354  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1355  if (w >= 0 && w < ssize){
1356  f6 += working_space[ssize + w];
1357  men +=1;
1358  }
1359  }
1360  f6 = f6 / men;
1361  g6 = 0, men = 0;
1362  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1363  if (w >= 0 && w < ssize){
1364  g6 += working_space[ssize + w];
1365  men +=1;
1366  }
1367  }
1368  g6 = g6 / men;
1369  b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1370  ai = i / 4;
1371  b8 = 0, men = 0;
1372  for (w = j - (Int_t)(4 * ai) - bw; w <= j - (Int_t)(4 * ai) + bw; w++){
1373  if (w >= 0 && w < ssize){
1374  b8 += working_space[ssize + w];
1375  men +=1;
1376  }
1377  }
1378  b8 = b8 / men;
1379  c8 = 0, men = 0;
1380  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
1381  if (w >= 0 && w < ssize){
1382  c8 += working_space[ssize + w];
1383  men +=1;
1384  }
1385  }
1386  c8 = c8 / men;
1387  d8 = 0, men = 0;
1388  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1389  if (w >= 0 && w < ssize){
1390  d8 += working_space[ssize + w];
1391  men +=1;
1392  }
1393  }
1394  d8 = d8 / men;
1395  e8 = 0, men = 0;
1396  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1397  if (w >= 0 && w < ssize){
1398  e8 += working_space[ssize + w];
1399  men +=1;
1400  }
1401  }
1402  e8 = e8 / men;
1403  f8 = 0, men = 0;
1404  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1405  if (w >= 0 && w < ssize){
1406  f8 += working_space[ssize + w];
1407  men +=1;
1408  }
1409  }
1410  f8 = f8 / men;
1411  g8 = 0, men = 0;
1412  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1413  if (w >= 0 && w < ssize){
1414  g8 += working_space[ssize + w];
1415  men +=1;
1416  }
1417  }
1418  g8 = g8 / men;
1419  h8 = 0, men = 0;
1420  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1421  if (w >= 0 && w < ssize){
1422  h8 += working_space[ssize + w];
1423  men +=1;
1424  }
1425  }
1426  h8 = h8 / men;
1427  i8 = 0, men = 0;
1428  for (w = j + (Int_t)(4 * ai) - bw; w <= j + (Int_t)(4 * ai) + bw; w++){
1429  if (w >= 0 && w < ssize){
1430  i8 += working_space[ssize + w];
1431  men +=1;
1432  }
1433  }
1434  i8 = i8 / men;
1435  b8 = ( -b8 + 8 * c8 - 28 * d8 + 56 * e8 - 56 * f8 - 28 * g8 + 8 * h8 - i8)/70;
1436  if (b < b8)
1437  b = b8;
1438  if (b < b6)
1439  b = b6;
1440  if (b < b4)
1441  b = b4;
1442  if (b < a)
1443  av = b;
1444  working_space[j]=av;
1445  }
1446  }
1447  for (j = i; j < ssize - i; j++)
1448  working_space[ssize + j] = working_space[j];
1449  if (direction == kBackIncreasingWindow)
1450  i += 1;
1451  else if(direction == kBackDecreasingWindow)
1452  i -= 1;
1453  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1454  }
1455 
1456  if (compton == kTRUE) {
1457  for (i = 0, b2 = 0; i < ssize; i++){
1458  b1 = b2;
1459  a = working_space[i], b = spectrum[i];
1460  j = i;
1461  if (TMath::Abs(a - b) >= 1) {
1462  b1 = i - 1;
1463  if (b1 < 0)
1464  b1 = 0;
1465  yb1 = working_space[b1];
1466  for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){
1467  a = working_space[b2], b = spectrum[b2];
1468  c = c + b - yb1;
1469  if (TMath::Abs(a - b) < 1) {
1470  priz = 1;
1471  yb2 = b;
1472  }
1473  }
1474  if (b2 == ssize)
1475  b2 -= 1;
1476  yb2 = working_space[b2];
1477  if (yb1 <= yb2){
1478  for (j = b1, c = 0; j <= b2; j++){
1479  b = spectrum[j];
1480  c = c + b - yb1;
1481  }
1482  if (c > 1){
1483  c = (yb2 - yb1) / c;
1484  for (j = b1, d = 0; j <= b2 && j < ssize; j++){
1485  b = spectrum[j];
1486  d = d + b - yb1;
1487  a = c * d + yb1;
1488  working_space[ssize + j] = a;
1489  }
1490  }
1491  }
1492 
1493  else{
1494  for (j = b2, c = 0; j >= b1; j--){
1495  b = spectrum[j];
1496  c = c + b - yb2;
1497  }
1498  if (c > 1){
1499  c = (yb1 - yb2) / c;
1500  for (j = b2, d = 0;j >= b1 && j >= 0; j--){
1501  b = spectrum[j];
1502  d = d + b - yb2;
1503  a = c * d + yb2;
1504  working_space[ssize + j] = a;
1505  }
1506  }
1507  }
1508  i=b2;
1509  }
1510  }
1511  }
1512 
1513  for (j = 0; j < ssize; j++)
1514  spectrum[j] = working_space[ssize + j];
1515  delete[]working_space;
1516  return 0;
1517 }
1518 
1519 
1520 ////////////////////////////////////////////////////////////////////////////////
1521 
1522 const char* TSpectrum::SmoothMarkov(Double_t *source, int ssize, int averWindow)
1523 {
1524  /* Begin_Html
1525  <b>One-dimensional markov spectrum smoothing function</b>
1526  <p>
1527  This function calculates smoothed spectrum from source spectrum based on
1528  Markov chain method. The result is placed in the array pointed by source
1529  pointer. On successful completion it returns 0. On error it returns pointer
1530  to the string describing error.
1531  <p>
1532  Function parameters:
1533  <ul>
1534  <li> source: pointer to the array of source spectrum
1535  <li> ssize: length of source array
1536  <li> averWindow: width of averaging smoothing window
1537  </ul>
1538  The goal of this function is the suppression of the statistical fluctuations.
1539  The algorithm is based on discrete Markov chain, which has very simple
1540  invariant distribution:
1541  <img width=551 height=63 src="gif/TSpectrum_Smoothing1.gif">
1542  <p>
1543  <img width=28 height=36 src="gif/TSpectrum_Smoothing2.gif"> being defined
1544  from the normalization condition
1545  <img width=70 height=52 src="gif/TSpectrum_Smoothing3.gif">.
1546  n is the length of the smoothed spectrum and
1547  <img width=258 height=60 src="gif/TSpectrum_Smoothing4.gif">
1548  <p>
1549  Reference:
1550  <ol>
1551  <li> Z.K. Silagadze, A new algorithm for automatic photopeak searches.
1552  NIM A 376 (1996), 451.
1553  </ol>
1554  <p>
1555  Example 14 - script Smoothing.c
1556  <p>
1557  <img width=296 height=182 src="gif/TSpectrum_Smoothing1.jpg">
1558  Fig. 23 Original noisy spectrum
1559  <p>
1560  <img width=296 height=182 src="gif/TSpectrum_Smoothing2.jpg">
1561  Fig. 24 Smoothed spectrum m=3
1562  <p>
1563  <img width=299 height=184 src="gif/TSpectrum_Smoothing3.jpg">
1564  Fig. 25 Smoothed spectrum
1565  <p>
1566  <img width=299 height=184 src="gif/TSpectrum_Smoothing4.jpg">
1567  Fig.26 Smoothed spectrum m=10
1568  <p>
1569  Script:
1570  <pre>
1571  // Example to illustrate smoothing using Markov algorithm (class TSpectrum).
1572  // To execute this example, do
1573  // root > .x Smoothing.C
1574 
1575  void Smoothing() {
1576  Int_t i;
1577  Double_t nbins = 1024;
1578  Double_t xmin = 0;
1579  Double_t xmax = nbins;
1580  Double_t * source = new Double_t[nbins];
1581  TH1F *h = new TH1F("h","Smoothed spectrum for m=3",nbins,xmin,xmax);
1582  TFile *f = new TFile("spectra\\TSpectrum.root");
1583  h=(TH1F*) f->Get("smooth1;1");
1584  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
1585  TCanvas *Smooth1 = gROOT->GetListOfCanvases()->FindObject("Smooth1");
1586  if (!Smooth1) Smooth1 = new TCanvas("Smooth1","Smooth1",10,10,1000,700);
1587  TSpectrum *s = new TSpectrum();
1588  s->SmoothMarkov(source,1024,3); //3, 7, 10
1589  for (i = 0; i < nbins; i++) h->SetBinContent(i + 1,source[i]);
1590  h->SetAxisRange(330,880);
1591  h->Draw("L");
1592  }
1593  </pre>
1594  End_Html */
1595 
1596  int xmin, xmax, i, l;
1597  Double_t a, b, maxch;
1598  Double_t nom, nip, nim, sp, sm, area = 0;
1599  if(averWindow <= 0)
1600  return "Averaging Window must be positive";
1601  Double_t *working_space = new Double_t[ssize];
1602  xmin = 0,xmax = ssize - 1;
1603  for(i = 0, maxch = 0; i < ssize; i++){
1604  working_space[i]=0;
1605  if(maxch < source[i])
1606  maxch = source[i];
1607 
1608  area += source[i];
1609  }
1610  if(maxch == 0) {
1611  delete [] working_space;
1612  return 0 ;
1613  }
1614 
1615  nom = 1;
1616  working_space[xmin] = 1;
1617  for(i = xmin; i < xmax; i++){
1618  nip = source[i] / maxch;
1619  nim = source[i + 1] / maxch;
1620  sp = 0,sm = 0;
1621  for(l = 1; l <= averWindow; l++){
1622  if((i + l) > xmax)
1623  a = source[xmax] / maxch;
1624 
1625  else
1626  a = source[i + l] / maxch;
1627  b = a - nip;
1628  if(a + nip <= 0)
1629  a = 1;
1630 
1631  else
1632  a = TMath::Sqrt(a + nip);
1633  b = b / a;
1634  b = TMath::Exp(b);
1635  sp = sp + b;
1636  if((i - l + 1) < xmin)
1637  a = source[xmin] / maxch;
1638 
1639  else
1640  a = source[i - l + 1] / maxch;
1641  b = a - nim;
1642  if(a + nim <= 0)
1643  a = 1;
1644  else
1645  a = TMath::Sqrt(a + nim);
1646  b = b / a;
1647  b = TMath::Exp(b);
1648  sm = sm + b;
1649  }
1650  a = sp / sm;
1651  a = working_space[i + 1] = working_space[i] * a;
1652  nom = nom + a;
1653  }
1654  for(i = xmin; i <= xmax; i++){
1655  working_space[i] = working_space[i] / nom;
1656  }
1657  for(i = 0; i < ssize; i++)
1658  source[i] = working_space[i] * area;
1659  delete[]working_space;
1660  return 0;
1661 }
1662 
1663 
1664 ////////////////////////////////////////////////////////////////////////////////
1665 
1666 const char *TSpectrum::Deconvolution(Double_t *source, const Double_t *response,
1667  int ssize, int numberIterations,
1668  int numberRepetitions, Double_t boost )
1669 {
1670  /* Begin_Html
1671  <b>One-dimensional deconvolution function</b>
1672  <p>
1673  This function calculates deconvolution from source spectrum according to
1674  response spectrum using Gold deconvolution algorithm. The result is placed
1675  in the vector pointed by source pointer. On successful completion it
1676  returns 0. On error it returns pointer to the string describing error. If
1677  desired after every numberIterations one can apply boosting operation
1678  (exponential function with exponent given by boost coefficient) and repeat
1679  it numberRepetitions times.
1680  <p>
1681  Function parameters:
1682  <ul>
1683  <li>source: pointer to the vector of source spectrum
1684  <li>response: pointer to the vector of response spectrum
1685  <li>ssize: length of source and response spectra
1686  numberIterations, for details we refer to the reference given below
1687  numberRepetitions, for repeated boosted deconvolution
1688  boost, boosting coefficient
1689  </ul>
1690  The goal of this function is the improvement of the resolution in spectra,
1691  decomposition of multiplets. The mathematical formulation of
1692  the convolution system is:
1693  <p>
1694  <img width=585 height=84 src="gif/TSpectrum_Deconvolution1.gif">
1695  <p>
1696  where h(i) is the impulse response function, x, y are input and output
1697  vectors, respectively, N is the length of x and h vectors. In matrix form
1698  we have:
1699  <p>
1700  <img width=597 height=360 src="gif/TSpectrum_Deconvolution2.gif">
1701  <p>
1702  Let us assume that we know the response and the output vector (spectrum) of
1703  the above given system. The deconvolution represents solution of the
1704  overdetermined system of linear equations, i.e., the calculation of the
1705  vector <b>x</b>. From numerical stability point of view the operation of
1706  deconvolution is extremely critical (ill-posed problem) as well as time
1707  consuming operation. The Gold deconvolution algorithm proves to work very
1708  well, other methods (Fourier, VanCittert etc) oscillate. It is suitable to
1709  process positive definite data (e.g. histograms).
1710  <p>
1711  <b>Gold deconvolution algorithm:</b>
1712  <p>
1713  <img width=551 height=233 src="gif/TSpectrum_Deconvolution3.gif">
1714  <p>
1715  Where L is given number of iterations (numberIterations parameter).
1716  <p>
1717  <b>Boosted deconvolution:</b>
1718  <ol>
1719  <li> Set the initial solution:
1720  End_Html Begin_Latex x^{(0)} = [1,1,...,1]^{T} End_Latex Begin_Html
1721  <li> Set required number of repetitions R and iterations L.
1722  <li> Set r = 1.
1723  <li>Using Gold deconvolution algorithm for k=1,2,...,L find
1724  End_Html Begin_Latex x^{(L)} End_Latex Begin_Html
1725  <li> If r = R stop calculation, else
1726  <ol>
1727  <li> Apply boosting operation, i.e., set
1728  End_Html Begin_Latex x^{(0)}(i) = [x^{(L)}(i)]^{p} End_Latex Begin_Html
1729  i=0,1,...N-1 and p is boosting coefficient &gt;0.
1730  <li> r = r + 1
1731  <li> continue in 4.
1732  </ol>
1733  </ol>
1734  <p>
1735  <b>References:</b>
1736  <ol>
1737  <li> Gold R., ANL-6984, Argonne National Laboratories, Argonne Ill, 1964.
1738  <li> Coote G.E., Iterative smoothing and deconvolution of one- and two-dimensional
1739  elemental distribution data, NIM B 130 (1997) 118.
1740  <li> M. Morháč, J. Kliman, V. Matoušek, M. Veselský,
1741  I. Turzo: Efficient one- and two-dimensional Gold deconvolution and
1742  its application to gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
1743  <li> Morháč M., Matoušek V., Kliman J., Efficient algorithm of multidimensional
1744  deconvolution and its application to nuclear data processing, Digital Signal
1745  Processing 13 (2003) 144.
1746  </ol>
1747  <p>
1748  <i>Example 8 - script Deconvolution.c :</i>
1749  <p>
1750  response function (usually peak) should be shifted left to the first
1751  non-zero channel (bin) (see Figure 9)
1752  <p>
1753  <img width=600 height=340 src="gif/TSpectrum_Deconvolution1.jpg">
1754  <p>
1755  Figure 9 Response spectrum.
1756  <p>
1757  <img width=946 height=407 src="gif/TSpectrum_Deconvolution2.jpg">
1758  <p>
1759  Figure 10 Principle how the response matrix is composed inside of the
1760  Deconvolution function.
1761  <img width=601 height=407 src="gif/TSpectrum_Deconvolution3.jpg">
1762  <p>
1763  Figure 11 Example of Gold deconvolution. The original source spectrum is
1764  drawn with black color, the spectrum after the deconvolution (10000
1765  iterations) with red color.
1766  <p>
1767  Script:
1768  <p>
1769  <pre>
1770  // Example to illustrate deconvolution function (class TSpectrum).
1771  // To execute this example, do
1772  // root > .x Deconvolution.C
1773 
1774  #include <TSpectrum>
1775 
1776  void Deconvolution() {
1777  Int_t i;
1778  Double_t nbins = 256;
1779  Double_t xmin = 0;
1780  Double_t xmax = nbins;
1781  Double_t * source = new Double_t[nbins];
1782  Double_t * response = new Double_t[nbins];
1783  TH1F *h = new TH1F("h","Deconvolution",nbins,xmin,xmax);
1784  TH1F *d = new TH1F("d","",nbins,xmin,xmax);
1785  TFile *f = new TFile("spectra\\TSpectrum.root");
1786  h=(TH1F*) f->Get("decon1;1");
1787  TFile *fr = new TFile("spectra\\TSpectrum.root");
1788  d=(TH1F*) fr->Get("decon_response;1");
1789  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
1790  for (i = 0; i < nbins; i++) response[i]=d->GetBinContent(i + 1);
1791  TCanvas *Decon1 = gROOT->GetListOfCanvases()->FindObject("Decon1");
1792  if (!Decon1) Decon1 = new TCanvas("Decon1","Decon1",10,10,1000,700);
1793  h->Draw("L");
1794  TSpectrum *s = new TSpectrum();
1795  s->Deconvolution(source,response,256,1000,1,1);
1796  for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
1797  d->SetLineColor(kRed);
1798  d->Draw("SAME L");
1799  }
1800  </pre>
1801  <p>
1802  <b>Examples of Gold deconvolution method:</b>
1803  <p>
1804  First let us study the influence of the number of iterations on the
1805  deconvolved spectrum (Figure 12).
1806  <p>
1807  <img width=602 height=409 src="gif/TSpectrum_Deconvolution_wide1.jpg">
1808  <p>
1809  Figure 12 Study of Gold deconvolution algorithm.The original source spectrum
1810  is drawn with black color, spectrum after 100 iterations with red color,
1811  spectrum after 1000 iterations with blue color, spectrum after 10000
1812  iterations with green color and spectrum after 100000 iterations with
1813  magenta color.
1814  <p>
1815  For relatively narrow peaks in the above given example the Gold
1816  deconvolution method is able to decompose overlapping peaks practically to
1817  delta - functions. In the next example we have chosen a synthetic data
1818  (spectrum, 256 channels) consisting of 5 very closely positioned, relatively
1819  wide peaks (sigma =5), with added noise (Figure 13). Thin lines represent
1820  pure Gaussians (see Table 1); thick line is a resulting spectrum with
1821  additive noise (10% of the amplitude of small peaks).
1822  <p>
1823  <img width=600 height=367 src="gif/TSpectrum_Deconvolution_wide2.jpg">
1824  <p>
1825  Figure 13 Testing example of synthetic spectrum composed of 5 Gaussians with
1826  added noise.
1827  <p>
1828  <table border=solid><tr>
1829  <td> Peak # </td><td> Position </td><td> Height </td><td> Area </td>
1830  </tr><tr>
1831  <td> 1 </td><td> 50 </td><td> 500 </td><td> 10159 </td>
1832  </tr><tr>
1833  <td> 2 </td><td> 70 </td><td> 3000 </td><td> 60957 </td>
1834  </tr><tr>
1835  <td> 3 </td><td> 80 </td><td> 1000 </td><td> 20319 </td>
1836  </tr><tr>
1837  <td> 4 </td><td> 100 </td><td> 5000 </td><td> 101596 </td>
1838  </tr><tr>
1839  <td> 5 </td><td> 110 </td><td> 500 </td><td> 10159 </td>
1840  </tr></table>
1841  <p>
1842  Table 1 Positions, heights and areas of peaks in the spectrum shown in
1843  Figure 13.
1844  <p>
1845  In ideal case, we should obtain the result given in Figure 14. The areas of
1846  the Gaussian components of the spectrum are concentrated completely to
1847  delta-functions. When solving the overdetermined system of linear equations
1848  with data from Figure 13 in the sense of minimum least squares criterion
1849  without any regularization we obtain the result with large oscillations
1850  (Figure 15). From mathematical point of view, it is the optimal solution in
1851  the unconstrained space of independent variables. From physical point of
1852  view we are interested only in a meaningful solution. Therefore, we have to
1853  employ regularization techniques (e.g. Gold deconvolution) and/or to
1854  confine the space of allowed solutions to subspace of positive solutions.
1855  <p>
1856  <img width=589 height=189 src="gif/TSpectrum_Deconvolution_wide3.jpg">
1857  <p>
1858  Figure 14 The same spectrum like in Figure 13, outlined bars show the
1859  contents of present components (peaks).
1860  <img width=585 height=183 src="gif/TSpectrum_Deconvolution_wide4.jpg">
1861  <p>
1862  Figure 15 Least squares solution of the system of linear equations without
1863  regularization.
1864  <p>
1865  <i>Example 9 - script Deconvolution_wide.c</i>
1866  <p>
1867  When we employ Gold deconvolution algorithm we obtain the result given in
1868  Fig. 16. One can observe that the resulting spectrum is smooth. On the
1869  other hand the method is not able to decompose completely the peaks in the
1870  spectrum.
1871  <p>
1872  <img width=601 height=407 src="gif/TSpectrum_Deconvolution_wide5.jpg">
1873  Figure 16 Example of Gold deconvolution for closely positioned wide peaks.
1874  The original source spectrum is drawn with black color, the spectrum after
1875  the deconvolution (10000 iterations) with red color.
1876  <p>
1877  Script:
1878  <p>
1879  <pre>
1880  // Example to illustrate deconvolution function (class TSpectrum).
1881  // To execute this example, do
1882  // root > .x Deconvolution_wide.C
1883 
1884  #include <TSpectrum>
1885 
1886  void Deconvolution_wide() {
1887  Int_t i;
1888  Double_t nbins = 256;
1889  Double_t xmin = 0;
1890  Double_t xmax = nbins;
1891  Double_t * source = new Double_t[nbins];
1892  Double_t * response = new Double_t[nbins];
1893  TH1F *h = new TH1F("h","Deconvolution",nbins,xmin,xmax);
1894  TH1F *d = new TH1F("d","",nbins,xmin,xmax);
1895  TFile *f = new TFile("spectra\\TSpectrum.root");
1896  h=(TH1F*) f->Get("decon3;1");
1897  TFile *fr = new TFile("spectra\\TSpectrum.root");
1898  d=(TH1F*) fr->Get("decon_response_wide;1");
1899  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
1900  for (i = 0; i < nbins; i++) response[i]=d->GetBinContent(i + 1);
1901  TCanvas *Decon1 = gROOT->GetListOfCanvases()->FindObject("Decon1");
1902  if (!Decon1) Decon1 = new TCanvas("Decon1",
1903  "Deconvolution of closely positioned overlapping peaks using Gold deconvolution method",10,10,1000,700);
1904  h->SetMaximum(30000);
1905  h->Draw("L");
1906  TSpectrum *s = new TSpectrum();
1907  s->Deconvolution(source,response,256,10000,1,1);
1908  for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
1909  d->SetLineColor(kRed);
1910  d->Draw("SAME L");
1911  }
1912  </pre>
1913  <p>
1914  <i>Example 10 - script Deconvolution_wide_boost.c :</i>
1915  <p>
1916  Further let us employ boosting operation into deconvolution (Fig. 17).
1917  <p>
1918  <img width=601 height=407 src="gif/TSpectrum_Deconvolution_wide6.jpg">
1919  <p>
1920  Figure 17 The original source spectrum is drawn with black color, the
1921  spectrum after the deconvolution with red color. Number of iterations = 200,
1922  number of repetitions = 50 and boosting coefficient = 1.2.
1923  <p>
1924  <table border=solid><tr>
1925  <td> Peak # </td> <td> Original/Estimated (max) position </td> <td> Original/Estimated area </td>
1926  </tr> <tr>
1927  <td> 1 </td> <td> 50/49 </td> <td> 10159/10419 </td>
1928  </tr> <tr>
1929  <td> 2 </td> <td> 70/70 </td> <td> 60957/58933 </td>
1930  </tr> <tr>
1931  <td> 3 </td> <td> 80/79 </td> <td> 20319/19935 </td>
1932  </tr> <tr>
1933  <td> 4 </td> <td> 100/100 </td> <td> 101596/105413 </td>
1934  </tr> <tr>
1935  <td> 5 </td> <td> 110/117 </td> <td> 10159/6676 </td>
1936  </tr> </table>
1937  <p>
1938  Table 2 Results of the estimation of peaks in spectrum shown in Figure 17.
1939  <p>
1940  One can observe that peaks are decomposed practically to delta functions.
1941  Number of peaks is correct, positions of big peaks as well as their areas
1942  are relatively well estimated. However there is a considerable error in
1943  the estimation of the position of small right hand peak.
1944  <p>
1945  Script:
1946  <p>
1947  <pre>
1948  // Example to illustrate deconvolution function (class TSpectrum).
1949  // To execute this example, do
1950  // root > .x Deconvolution_wide_boost.C
1951 
1952  #include <TSpectrum>
1953 
1954  void Deconvolution_wide_boost() {
1955  Int_t i;
1956  Double_t nbins = 256;
1957  Double_t xmin = 0;
1958  Double_t xmax = nbins;
1959  Double_t * source = new Double_t[nbins];
1960  Double_t * response = new Double_t[nbins];
1961  TH1F *h = new TH1F("h","Deconvolution",nbins,xmin,xmax);
1962  TH1F *d = new TH1F("d","",nbins,xmin,xmax);
1963  TFile *f = new TFile("spectra\\TSpectrum.root");
1964  h=(TH1F*) f->Get("decon3;1");
1965  TFile *fr = new TFile("spectra\\TSpectrum.root");
1966  d=(TH1F*) fr->Get("decon_response_wide;1");
1967  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
1968  for (i = 0; i < nbins; i++) response[i]=d->GetBinContent(i + 1);
1969  TCanvas *Decon1 = gROOT->GetListOfCanvases()->FindObject("Decon1");
1970  if (!Decon1) Decon1 = new TCanvas("Decon1",
1971  "Deconvolution of closely positioned overlapping peaks using boosted Gold deconvolution method",10,10,1000,700);
1972  h->SetMaximum(110000);
1973  h->Draw("L");
1974  TSpectrum *s = new TSpectrum();
1975  s->Deconvolution(source,response,256,200,50,1.2);
1976  for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
1977  d->SetLineColor(kRed);
1978  d->Draw("SAME L");
1979  }
1980  </pre>
1981  End_Html */
1982 
1983  if (ssize <= 0)
1984  return "Wrong Parameters";
1985 
1986  if (numberRepetitions <= 0)
1987  return "Wrong Parameters";
1988 
1989  // working_space-pointer to the working vector
1990  // (its size must be 4*ssize of source spectrum)
1991  Double_t *working_space = new Double_t[4 * ssize];
1992  int i, j, k, lindex, posit, lh_gold, l, repet;
1993  Double_t lda, ldb, ldc, area, maximum;
1994  area = 0;
1995  lh_gold = -1;
1996  posit = 0;
1997  maximum = 0;
1998 
1999 //read response vector
2000  for (i = 0; i < ssize; i++) {
2001  lda = response[i];
2002  if (lda != 0)
2003  lh_gold = i + 1;
2004  working_space[i] = lda;
2005  area += lda;
2006  if (lda > maximum) {
2007  maximum = lda;
2008  posit = i;
2009  }
2010  }
2011  if (lh_gold == -1) {
2012  delete [] working_space;
2013  return "ZERO RESPONSE VECTOR";
2014  }
2015 
2016 //read source vector
2017  for (i = 0; i < ssize; i++)
2018  working_space[2 * ssize + i] = source[i];
2019 
2020 // create matrix at*a and vector at*y
2021  for (i = 0; i < ssize; i++){
2022  lda = 0;
2023  for (j = 0; j < ssize; j++){
2024  ldb = working_space[j];
2025  k = i + j;
2026  if (k < ssize){
2027  ldc = working_space[k];
2028  lda = lda + ldb * ldc;
2029  }
2030  }
2031  working_space[ssize + i] = lda;
2032  lda = 0;
2033  for (k = 0; k < ssize; k++){
2034  l = k - i;
2035  if (l >= 0){
2036  ldb = working_space[l];
2037  ldc = working_space[2 * ssize + k];
2038  lda = lda + ldb * ldc;
2039  }
2040  }
2041  working_space[3 * ssize + i]=lda;
2042  }
2043 
2044 // move vector at*y
2045  for (i = 0; i < ssize; i++){
2046  working_space[2 * ssize + i] = working_space[3 * ssize + i];
2047  }
2048 
2049 //initialization of resulting vector
2050  for (i = 0; i < ssize; i++)
2051  working_space[i] = 1;
2052 
2053  //**START OF ITERATIONS**
2054  for (repet = 0; repet < numberRepetitions; repet++) {
2055  if (repet != 0) {
2056  for (i = 0; i < ssize; i++)
2057  working_space[i] = TMath::Power(working_space[i], boost);
2058  }
2059  for (lindex = 0; lindex < numberIterations; lindex++) {
2060  for (i = 0; i < ssize; i++) {
2061  if (working_space[2 * ssize + i] > 0.000001
2062  && working_space[i] > 0.000001) {
2063  lda = 0;
2064  for (j = 0; j < lh_gold; j++) {
2065  ldb = working_space[j + ssize];
2066  if (j != 0){
2067  k = i + j;
2068  ldc = 0;
2069  if (k < ssize)
2070  ldc = working_space[k];
2071  k = i - j;
2072  if (k >= 0)
2073  ldc += working_space[k];
2074  }
2075 
2076  else
2077  ldc = working_space[i];
2078  lda = lda + ldb * ldc;
2079  }
2080  ldb = working_space[2 * ssize + i];
2081  if (lda != 0)
2082  lda = ldb / lda;
2083 
2084  else
2085  lda = 0;
2086  ldb = working_space[i];
2087  lda = lda * ldb;
2088  working_space[3 * ssize + i] = lda;
2089  }
2090  }
2091  for (i = 0; i < ssize; i++)
2092  working_space[i] = working_space[3 * ssize + i];
2093  }
2094  }
2095 
2096 //shift resulting spectrum
2097  for (i = 0; i < ssize; i++) {
2098  lda = working_space[i];
2099  j = i + posit;
2100  j = j % ssize;
2101  working_space[ssize + j] = lda;
2102  }
2103 
2104 //write back resulting spectrum
2105  for (i = 0; i < ssize; i++)
2106  source[i] = area * working_space[ssize + i];
2107  delete[]working_space;
2108  return 0;
2109 }
2110 
2111 
2112 ////////////////////////////////////////////////////////////////////////////////
2113 
2114 const char *TSpectrum::DeconvolutionRL(Double_t *source, const Double_t *response,
2115  int ssize, int numberIterations,
2116  int numberRepetitions, Double_t boost )
2117 {
2118  /* Begin_Html
2119  <b>One-dimensional deconvolution function.</b>
2120  <p>
2121  This function calculates deconvolution from source spectrum according to
2122  response spectrum using Richardson-Lucy deconvolution algorithm. The result
2123  is placed in the vector pointed by source pointer. On successful completion
2124  it returns 0. On error it returns pointer to the string describing error.
2125  If desired after every numberIterations one can apply boosting operation
2126  (exponential function with exponent given by boost coefficient) and repeat
2127  it numberRepetitions times (see Gold deconvolution).
2128  <p>
2129  Function parameters:
2130  <ul>
2131  <li> source: pointer to the vector of source spectrum
2132  <li> response: pointer to the vector of response spectrum
2133  <li> ssize: length of source and response spectra
2134  numberIterations, for details we refer to the reference given above
2135  numberRepetitions, for repeated boosted deconvolution
2136  boost, boosting coefficient
2137  </ul>
2138  <p>
2139  <b>Richardson-Lucy deconvolution algorithm:</b>
2140  <p>
2141  For discrete systems it has the form:
2142  <p>
2143  <img width=438 height=98 src="gif/TSpectrum_DeconvolutionRL1.gif">
2144  <p>
2145  <img width=124 height=39 src="gif/TSpectrum_DeconvolutionRL2.gif">
2146  <p>
2147  for positive input data and response matrix this iterative method forces
2148  the deconvoluted spectra to be non-negative. The Richardson-Lucy
2149  iteration converges to the maximum likelihood solution for Poisson statistics
2150  in the data.
2151  <p>
2152  <b>References:</b>
2153  <ol>
2154  <li> Abreu M.C. et al., A four-dimensional deconvolution method to correct NA38
2155  experimental data, NIM A 405 (1998) 139.
2156  <li> Lucy L.B., A.J. 79 (1974) 745.
2157  <li> Richardson W.H., J. Opt. Soc. Am. 62 (1972) 55.
2158  </ol>
2159  <p>
2160  <b>Examples of Richardson-Lucy deconvolution method:</b>
2161  <p>
2162  <i>Example 11 - script DeconvolutionRL_wide.c :</i>
2163  <p>
2164  When we employ Richardson-Lucy deconvolution algorithm to our data from
2165  Fig. 13 we obtain the result given in Fig. 18. One can observe improvements
2166  as compared to the result achieved by Gold deconvolution. Neverthless it is
2167  unable to decompose the multiplet.
2168  <p>
2169  <img width=601 height=407 src="gif/TSpectrum_DeconvolutionRL_wide1.jpg">
2170  Figure 18 Example of Richardson-Lucy deconvolution for closely positioned
2171  wide peaks. The original source spectrum is drawn with black color, the
2172  spectrum after the deconvolution (10000 iterations) with red color.
2173  <p>
2174  Script:
2175  <p>
2176  <pre>
2177  // Example to illustrate deconvolution function (class TSpectrum).
2178  // To execute this example, do
2179  // root > .x DeconvolutionRL_wide.C
2180 
2181  #include <TSpectrum>
2182 
2183  void DeconvolutionRL_wide() {
2184  Int_t i;
2185  Double_t nbins = 256;
2186  Double_t xmin = 0;
2187  Double_t xmax = nbins;
2188  Double_t * source = new Double_t[nbins];
2189  Double_t * response = new Double_t[nbins];
2190  TH1F *h = new TH1F("h","Deconvolution",nbins,xmin,xmax);
2191  TH1F *d = new TH1F("d","",nbins,xmin,xmax);
2192  TFile *f = new TFile("spectra\\TSpectrum.root");
2193  h=(TH1F*) f->Get("decon3;1");
2194  TFile *fr = new TFile("spectra\\TSpectrum.root");
2195  d=(TH1F*) fr->Get("decon_response_wide;1");
2196  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
2197  for (i = 0; i < nbins; i++) response[i]=d->GetBinContent(i + 1);
2198  TCanvas *Decon1 = gROOT->GetListOfCanvases()->FindObject("Decon1");
2199  if (!Decon1) Decon1 = new TCanvas("Decon1",
2200  "Deconvolution of closely positioned overlapping peaks using Richardson-Lucy deconvolution method",
2201  10,10,1000,700);
2202  h->SetMaximum(30000);
2203  h->Draw("L");
2204  TSpectrum *s = new TSpectrum();
2205  s->DeconvolutionRL(source,response,256,10000,1,1);
2206  for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
2207  d->SetLineColor(kRed);
2208  d->Draw("SAME L");
2209  }
2210  </pre>
2211  <p>
2212  <i>Example 12 - script DeconvolutionRL_wide_boost.c :</i>
2213  <p>
2214  Further let us employ boosting operation into deconvolution (Fig. 19).
2215  <img width=601 height=407 src="gif/TSpectrum_DeconvolutionRL_wide2.jpg">
2216  <p>
2217  Figure 19 The original source spectrum is drawn with black color, the
2218  spectrum after the deconvolution with red color. Number of iterations = 200,
2219  number of repetitions = 50 and boosting coefficient = 1.2.
2220  <p>
2221  <table border=solid>
2222  <tr><td> Peak # </td><td> Original/Estimated (max) position </td><td> Original/Estimated area </td></tr>
2223  <tr><td> 1 </td><td> 50/51 </td><td> 10159/11426 </td></tr>
2224  <tr><td> 2 </td><td> 70/71 </td><td> 60957/65003 </td></tr>
2225  <tr><td> 3 </td><td> 80/81 </td><td> 20319/12813 </td></tr>
2226  <tr><td> 4 </td><td> 100/100 </td><td> 101596/101851 </td></tr>
2227  <tr><td> 5 </td><td> 110/111 </td><td> 10159/8920 </td></tr>
2228  </table>
2229  <p>
2230  Table 3 Results of the estimation of peaks in spectrum shown in Figure 19.
2231  <p>
2232  One can observe improvements in the estimation of peak positions as compared
2233  to the results achieved by Gold deconvolution.
2234  <p>
2235  Script:
2236  <pre>
2237  // Example to illustrate deconvolution function (class TSpectrum).
2238  // To execute this example, do
2239  // root > .x DeconvolutionRL_wide_boost.C
2240 
2241  #include <TSpectrum>
2242 
2243  void DeconvolutionRL_wide_boost() {
2244  Int_t i;
2245  Double_t nbins = 256;
2246  Double_t xmin = 0;
2247  Double_t xmax = nbins;
2248  Double_t * source = new Double_t[nbins];
2249  Double_t * response = new Double_t[nbins];
2250  TH1F *h = new TH1F("h","Deconvolution",nbins,xmin,xmax);
2251  TH1F *d = new TH1F("d","",nbins,xmin,xmax);
2252  TFile *f = new TFile("spectra\\TSpectrum.root");
2253  h=(TH1F*) f->Get("decon3;1");
2254  TFile *fr = new TFile("spectra\\TSpectrum.root");
2255  d=(TH1F*) fr->Get("decon_response_wide;1");
2256  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
2257  for (i = 0; i < nbins; i++) response[i]=d->GetBinContent(i + 1);
2258  TCanvas *Decon1 = gROOT->GetListOfCanvases()->FindObject("Decon1");
2259  if (!Decon1) Decon1 = new TCanvas("Decon1",
2260  "Deconvolution of closely positioned overlapping peaks using boosted Richardson-Lucy deconvolution method",
2261  10,10,1000,700);
2262  h->SetMaximum(110000);
2263  h->Draw("L");
2264  TSpectrum *s = new TSpectrum();
2265  s->DeconvolutionRL(source,response,256,200,50,1.2);
2266  for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
2267  d->SetLineColor(kRed);
2268  d->Draw("SAME L");
2269  }
2270  </pre>
2271  End_Html */
2272 
2273  if (ssize <= 0)
2274  return "Wrong Parameters";
2275 
2276  if (numberRepetitions <= 0)
2277  return "Wrong Parameters";
2278 
2279  // working_space-pointer to the working vector
2280  // (its size must be 4*ssize of source spectrum)
2281  Double_t *working_space = new Double_t[4 * ssize];
2282  int i, j, k, lindex, posit, lh_gold, repet, kmin, kmax;
2283  Double_t lda, ldb, ldc, maximum;
2284  lh_gold = -1;
2285  posit = 0;
2286  maximum = 0;
2287 
2288 //read response vector
2289  for (i = 0; i < ssize; i++) {
2290  lda = response[i];
2291  if (lda != 0)
2292  lh_gold = i + 1;
2293  working_space[ssize + i] = lda;
2294  if (lda > maximum) {
2295  maximum = lda;
2296  posit = i;
2297  }
2298  }
2299  if (lh_gold == -1) {
2300  delete [] working_space;
2301  return "ZERO RESPONSE VECTOR";
2302  }
2303 
2304 //read source vector
2305  for (i = 0; i < ssize; i++)
2306  working_space[2 * ssize + i] = source[i];
2307 
2308 //initialization of resulting vector
2309  for (i = 0; i < ssize; i++){
2310  if (i <= ssize - lh_gold)
2311  working_space[i] = 1;
2312 
2313  else
2314  working_space[i] = 0;
2315 
2316  }
2317  //**START OF ITERATIONS**
2318  for (repet = 0; repet < numberRepetitions; repet++) {
2319  if (repet != 0) {
2320  for (i = 0; i < ssize; i++)
2321  working_space[i] = TMath::Power(working_space[i], boost);
2322  }
2323  for (lindex = 0; lindex < numberIterations; lindex++) {
2324  for (i = 0; i <= ssize - lh_gold; i++){
2325  lda = 0;
2326  if (working_space[i] > 0){//x[i]
2327  for (j = i; j < i + lh_gold; j++){
2328  ldb = working_space[2 * ssize + j];//y[j]
2329  if (j < ssize){
2330  if (ldb > 0){//y[j]
2331  kmax = j;
2332  if (kmax > lh_gold - 1)
2333  kmax = lh_gold - 1;
2334  kmin = j + lh_gold - ssize;
2335  if (kmin < 0)
2336  kmin = 0;
2337  ldc = 0;
2338  for (k = kmax; k >= kmin; k--){
2339  ldc += working_space[ssize + k] * working_space[j - k];//h[k]*x[j-k]
2340  }
2341  if (ldc > 0)
2342  ldb = ldb / ldc;
2343 
2344  else
2345  ldb = 0;
2346  }
2347  ldb = ldb * working_space[ssize + j - i];//y[j]*h[j-i]/suma(h[j][k]x[k])
2348  }
2349  lda += ldb;
2350  }
2351  lda = lda * working_space[i];
2352  }
2353  working_space[3 * ssize + i] = lda;
2354  }
2355  for (i = 0; i < ssize; i++)
2356  working_space[i] = working_space[3 * ssize + i];
2357  }
2358  }
2359 
2360 //shift resulting spectrum
2361  for (i = 0; i < ssize; i++) {
2362  lda = working_space[i];
2363  j = i + posit;
2364  j = j % ssize;
2365  working_space[ssize + j] = lda;
2366  }
2367 
2368 //write back resulting spectrum
2369  for (i = 0; i < ssize; i++)
2370  source[i] = working_space[ssize + i];
2371  delete[]working_space;
2372  return 0;
2373 }
2374 
2375 
2376 ////////////////////////////////////////////////////////////////////////////////
2377 
2378 const char *TSpectrum::Unfolding(Double_t *source,
2379  const Double_t **respMatrix,
2380  int ssizex, int ssizey,
2381  int numberIterations,
2382  int numberRepetitions, Double_t boost)
2383 {
2384  /* Begin_Html
2385  <b>One-dimensional unfolding function</b>
2386  <p>
2387  This function unfolds source spectrum according to response matrix columns.
2388  The result is placed in the vector pointed by source pointer.
2389  The coefficients of the resulting vector represent contents of the columns
2390  (weights) in the input vector. On successful completion it returns 0. On
2391  error it returns pointer to the string describing error. If desired after
2392  every numberIterations one can apply boosting operation (exponential
2393  function with exponent given by boost coefficient) and repeat it
2394  numberRepetitions times. For details we refer to [1].
2395  <p>
2396  Function parameters:
2397  <ul>
2398  <li> source: pointer to the vector of source spectrum
2399  <li> respMatrix: pointer to the matrix of response spectra
2400  <li> ssizex: length of source spectrum and # of columns of the response
2401  matrix. ssizex must be >= ssizey.
2402  <li> ssizey: length of destination spectrum and # of rows of the response
2403  matrix.
2404  <li> numberIterations: number of iterations
2405  <li> numberRepetitions: number of repetitions for boosted deconvolution.
2406  It must be greater or equal to one.
2407  <li> boost: boosting coefficient, applies only if numberRepetitions is
2408  greater than one.
2409  </ul>
2410  <p>
2411  <b>Unfolding:</b>
2412  <p>
2413  The goal is the decomposition of spectrum to a given set of component
2414  spectra.
2415  <p>
2416  The mathematical formulation of the discrete linear system is:
2417  <p>
2418  <img width=588 height=89 src="gif/TSpectrum_Unfolding1.gif">
2419  <p>
2420  <img width=597 height=228 src="gif/TSpectrum_Unfolding2.gif">
2421  <p>
2422  <b>References:</b>
2423  <ol>
2424  <li> Jandel M., Morháč M., Kliman J., Krupa L., Matoušek
2425  V., Hamilton J. H., Ramaya A. V.:
2426  Decomposition of continuum gamma-ray spectra using synthetized response matrix.
2427  NIM A 516 (2004), 172-183.
2428  </ol>
2429  <p>
2430  <b>Example of unfolding:</b>
2431  <p>
2432  <i>Example 13 - script Unfolding.c:</i>
2433  <p>
2434  <img width=442 height=648 src="gif/TSpectrum_Unfolding3.gif">
2435  <p>
2436  Fig. 20 Response matrix composed of neutron spectra of pure
2437  chemical elements.
2438  <img width=604 height=372 src="gif/TSpectrum_Unfolding2.jpg">
2439  <p>
2440  Fig. 21 Source neutron spectrum to be decomposed
2441  <P>
2442  <img width=600 height=360 src="gif/TSpectrum_Unfolding3.jpg">
2443  <p>
2444  Fig. 22 Spectrum after decomposition, contains 10 coefficients, which
2445  correspond to contents of chemical components (dominant 8-th and 10-th
2446  components, i.e. O, Si)
2447  <p>
2448  Script:
2449  <pre>
2450  // Example to illustrate unfolding function (class TSpectrum).
2451  // To execute this example, do
2452  // root > .x Unfolding.C
2453 
2454  #include <TSpectrum>
2455 
2456  void Unfolding() {
2457  Int_t i, j;
2458  Int_t nbinsx = 2048;
2459  Int_t nbinsy = 10;
2460  Double_t xmin = 0;
2461  Double_t xmax = nbinsx;
2462  Double_t ymin = 0;
2463  Double_t ymax = nbinsy;
2464  Double_t * source = new Double_t[nbinsx];
2465  Double_t ** response = new Double_t *[nbinsy];
2466  for (i=0;i<nbinsy;i++) response[i]=new Double_t[nbinsx];
2467  TH1F *h = new TH1F("h","",nbinsx,xmin,xmax);
2468  TH1F *d = new TH1F("d","Decomposition - unfolding",nbinsx,xmin,xmax);
2469  TH2F *decon_unf_resp = new TH2F("decon_unf_resp","Root File",nbinsy,ymin,ymax,nbinsx,xmin,xmax);
2470  TFile *f = new TFile("spectra\\TSpectrum.root");
2471  h=(TH1F*) f->Get("decon_unf_in;1");
2472  TFile *fr = new TFile("spectra\\TSpectrum.root");
2473  decon_unf_resp = (TH2F*) fr->Get("decon_unf_resp;1");
2474  for (i = 0; i < nbinsx; i++) source[i] = h->GetBinContent(i + 1);
2475  for (i = 0; i < nbinsy; i++){
2476  for (j = 0; j< nbinsx; j++){
2477  response[i][j] = decon_unf_resp->GetBinContent(i + 1, j + 1);
2478  }
2479  }
2480  TCanvas *Decon1 = gROOT->GetListOfCanvases()->FindObject("Decon1");
2481  if (!Decon1) Decon1 = new TCanvas("Decon1","Decon1",10,10,1000,700);
2482  h->Draw("L");
2483  TSpectrum *s = new TSpectrum();
2484  s->Unfolding(source,response,nbinsx,nbinsy,1000,1,1);
2485  for (i = 0; i < nbinsy; i++) d->SetBinContent(i + 1,source[i]);
2486  d->SetLineColor(kRed);
2487  d->SetAxisRange(0,nbinsy);
2488  d->Draw("");
2489  }
2490  </pre>
2491  End_Html */
2492 
2493  int i, j, k, lindex, lhx = 0, repet;
2494  Double_t lda, ldb, ldc, area;
2495  if (ssizex <= 0 || ssizey <= 0)
2496  return "Wrong Parameters";
2497  if (ssizex < ssizey)
2498  return "Sizex must be greater than sizey)";
2499  if (numberIterations <= 0)
2500  return "Number of iterations must be positive";
2501  Double_t *working_space =
2502  new Double_t[ssizex * ssizey + 2 * ssizey * ssizey + 4 * ssizex];
2503 
2504 /*read response matrix*/
2505  for (j = 0; j < ssizey && lhx != -1; j++) {
2506  area = 0;
2507  lhx = -1;
2508  for (i = 0; i < ssizex; i++) {
2509  lda = respMatrix[j][i];
2510  if (lda != 0) {
2511  lhx = i + 1;
2512  }
2513  working_space[j * ssizex + i] = lda;
2514  area = area + lda;
2515  }
2516  if (lhx != -1) {
2517  for (i = 0; i < ssizex; i++)
2518  working_space[j * ssizex + i] /= area;
2519  }
2520  }
2521  if (lhx == -1) {
2522  delete [] working_space;
2523  return ("ZERO COLUMN IN RESPONSE MATRIX");
2524  }
2525 
2526 /*read source vector*/
2527  for (i = 0; i < ssizex; i++)
2528  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2529  source[i];
2530 
2531 /*create matrix at*a + at*y */
2532  for (i = 0; i < ssizey; i++) {
2533  for (j = 0; j < ssizey; j++) {
2534  lda = 0;
2535  for (k = 0; k < ssizex; k++) {
2536  ldb = working_space[ssizex * i + k];
2537  ldc = working_space[ssizex * j + k];
2538  lda = lda + ldb * ldc;
2539  }
2540  working_space[ssizex * ssizey + ssizey * i + j] = lda;
2541  }
2542  lda = 0;
2543  for (k = 0; k < ssizex; k++) {
2544  ldb = working_space[ssizex * i + k];
2545  ldc =
2546  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
2547  k];
2548  lda = lda + ldb * ldc;
2549  }
2550  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
2551  lda;
2552  }
2553 
2554 /*move vector at*y*/
2555  for (i = 0; i < ssizey; i++)
2556  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2557  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2558 
2559 /*create matrix at*a*at*a + vector at*a*at*y */
2560  for (i = 0; i < ssizey; i++) {
2561  for (j = 0; j < ssizey; j++) {
2562  lda = 0;
2563  for (k = 0; k < ssizey; k++) {
2564  ldb = working_space[ssizex * ssizey + ssizey * i + k];
2565  ldc = working_space[ssizex * ssizey + ssizey * j + k];
2566  lda = lda + ldb * ldc;
2567  }
2568  working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j] =
2569  lda;
2570  }
2571  lda = 0;
2572  for (k = 0; k < ssizey; k++) {
2573  ldb = working_space[ssizex * ssizey + ssizey * i + k];
2574  ldc =
2575  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
2576  k];
2577  lda = lda + ldb * ldc;
2578  }
2579  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
2580  lda;
2581  }
2582 
2583 /*move at*a*at*y*/
2584  for (i = 0; i < ssizey; i++)
2585  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
2586  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2587 
2588 /*initialization in resulting vector */
2589  for (i = 0; i < ssizey; i++)
2590  working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = 1;
2591 
2592  /***START OF ITERATIONS***/
2593  for (repet = 0; repet < numberRepetitions; repet++) {
2594  if (repet != 0) {
2595  for (i = 0; i < ssizey; i++)
2596  working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = TMath::Power(working_space[ssizex * ssizey + 2 * ssizey * ssizey + i], boost);
2597  }
2598  for (lindex = 0; lindex < numberIterations; lindex++) {
2599  for (i = 0; i < ssizey; i++) {
2600  lda = 0;
2601  for (j = 0; j < ssizey; j++) {
2602  ldb =
2603  working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j];
2604  ldc = working_space[ssizex * ssizey + 2 * ssizey * ssizey + j];
2605  lda = lda + ldb * ldc;
2606  }
2607  ldb =
2608  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i];
2609  if (lda != 0) {
2610  lda = ldb / lda;
2611  }
2612 
2613  else
2614  lda = 0;
2615  ldb = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2616  lda = lda * ldb;
2617  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] = lda;
2618  }
2619  for (i = 0; i < ssizey; i++)
2620  working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
2621  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2622  }
2623  }
2624 
2625 /*write back resulting spectrum*/
2626  for (i = 0; i < ssizex; i++) {
2627  if (i < ssizey)
2628  source[i] = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2629 
2630  else
2631  source[i] = 0;
2632  }
2633  delete[]working_space;
2634  return 0;
2635 }
2636 
2637 
2638 ////////////////////////////////////////////////////////////////////////////////
2639 
2640 Int_t TSpectrum::SearchHighRes(Double_t *source,Double_t *destVector, int ssize,
2641  Double_t sigma, Double_t threshold,
2642  bool backgroundRemove,int deconIterations,
2643  bool markov, int averWindow)
2644 {
2645  /* Begin_Html
2646  <b>One-dimensional high-resolution peak search function</b>
2647  <p>
2648  This function searches for peaks in source spectrum. It is based on
2649  deconvolution method. First the background is removed (if desired), then
2650  Markov smoothed spectrum is calculated (if desired), then the response
2651  function is generated according to given sigma and deconvolution is
2652  carried out. The order of peaks is arranged according to their heights in
2653  the spectrum after background elimination. The highest peak is the first in
2654  the list. On success it returns number of found peaks.
2655  <p>
2656  <b>Function parameters:</b>
2657  <ul>
2658  <li> source: pointer to the vector of source spectrum.
2659  <li> destVector: pointer to the vector of resulting deconvolved spectrum.
2660  <li> ssize: length of source spectrum.
2661  <li> sigma: sigma of searched peaks, for details we refer to manual.
2662  <li> threshold: threshold value in % for selected peaks, peaks with
2663  amplitude less than threshold*highest_peak/100
2664  are ignored, see manual.
2665  <li> backgroundRemove: logical variable, set if the removal of
2666  background before deconvolution is desired.
2667  <li> deconIterations-number of iterations in deconvolution operation.
2668  <li> markov: logical variable, if it is true, first the source spectrum
2669  is replaced by new spectrum calculated using Markov
2670  chains method.
2671  <li> averWindow: averanging window of searched peaks, for details
2672  we refer to manual (applies only for Markov method).
2673  </ul>
2674  <p>
2675  <b>Peaks searching:</b>
2676  <p>
2677  The goal of this function is to identify automatically the peaks in spectrum
2678  with the presence of the continuous background and statistical
2679  fluctuations - noise.
2680  <p>
2681  The common problems connected with correct peak identification are:
2682  <ul>
2683  <li> non-sensitivity to noise, i.e., only statistically
2684  relevant peaks should be identified.
2685  <li> non-sensitivity of the algorithm to continuous
2686  background.
2687  <li> ability to identify peaks close to the edges of the
2688  spectrum region. Usually peak finders fail to detect them.
2689  <li> resolution, decomposition of Double_tts and multiplets.
2690  The algorithm should be able to recognize close positioned peaks.
2691  <li> ability to identify peaks with different sigma.
2692  </ul>
2693  <img width=600 height=375 src="gif/TSpectrum_Searching1.jpg">
2694  <p>
2695  Fig. 27 An example of one-dimensional synthetic spectrum with found peaks
2696  denoted by markers.
2697  <p>
2698  <b>References:</b>
2699  <ol>
2700  <li> M.A. Mariscotti: A method for identification of peaks in the presence of
2701  background and its application to spectrum analysis. NIM 50 (1967),
2702  309-320.
2703  <li> M. Morháč, J. Kliman, V. Matoušek, M. Veselský,
2704  I. Turzo.:Identification of peaks in
2705  multidimensional coincidence gamma-ray spectra. NIM, A443 (2000) 108-125.
2706  <li> Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM
2707  A 376 (1996), 451.
2708  </ol>
2709  <p>
2710  <b>Examples of peak searching method:</b>
2711  <p>
2712  The SearchHighRes function provides users with the possibility to vary the
2713  input parameters and with the access to the output deconvolved data in the
2714  destination spectrum. Based on the output data one can tune the parameters.
2715  <p>
2716  Example 15 - script SearchHR1.c:
2717  <img width=600 height=321 src="gif/TSpectrum_Searching1.jpg">
2718  <p>
2719  Fig. 28 One-dimensional spectrum with found peaks denoted by markers, 3
2720  iterations steps in the deconvolution.
2721  <p>
2722  <img width=600 height=323 src="gif/TSpectrum_Searching2.jpg">
2723  Fig. 29 One-dimensional spectrum with found peaks denoted by markers, 8
2724  iterations steps in the deconvolution.
2725  <p>
2726  Script:
2727  <pre>
2728  // Example to illustrate high resolution peak searching function (class TSpectrum).
2729  // To execute this example, do
2730  // root > .x SearchHR1.C
2731 
2732  #include <TSpectrum>
2733 
2734  void SearchHR1() {
2735  Double_t fPositionX[100];
2736  Double_t fPositionY[100];
2737  Int_t fNPeaks = 0;
2738  Int_t i,nfound,bin;
2739  Double_t nbins = 1024,a;
2740  Double_t xmin = 0;
2741  Double_t xmax = nbins;
2742  Double_t * source = new Double_t[nbins];
2743  Double_t * dest = new Double_t[nbins];
2744  TH1F *h = new TH1F("h","High resolution peak searching, number of iterations = 3",nbins,xmin,xmax);
2745  TH1F *d = new TH1F("d","",nbins,xmin,xmax);
2746  TFile *f = new TFile("spectra\\TSpectrum.root");
2747  h=(TH1F*) f->Get("search2;1");
2748  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
2749  TCanvas *Search = gROOT->GetListOfCanvases()->FindObject("Search");
2750  if (!Search) Search = new TCanvas("Search","Search",10,10,1000,700);
2751  h->SetMaximum(4000);
2752  h->Draw("L");
2753  TSpectrum *s = new TSpectrum();
2754  nfound = s->SearchHighRes(source, dest, nbins, 8, 2, kTRUE, 3, kTRUE, 3);
2755  Double_t *xpeaks = s->GetPositionX();
2756  for (i = 0; i < nfound; i++) {
2757  a=xpeaks[i];
2758  bin = 1 + Int_t(a + 0.5);
2759  fPositionX[i] = h->GetBinCenter(bin);
2760  fPositionY[i] = h->GetBinContent(bin);
2761  }
2762  TPolyMarker * pm = (TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");
2763  if (pm) {
2764  h->GetListOfFunctions()->Remove(pm);
2765  delete pm;
2766  }
2767  pm = new TPolyMarker(nfound, fPositionX, fPositionY);
2768  h->GetListOfFunctions()->Add(pm);
2769  pm->SetMarkerStyle(23);
2770  pm->SetMarkerColor(kRed);
2771  pm->SetMarkerSize(1.3);
2772  for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,dest[i]);
2773  d->SetLineColor(kRed);
2774  d->Draw("SAME");
2775  printf("Found %d candidate peaks\n",nfound);
2776  for(i=0;i<nfound;i++)
2777  printf("posx= %d, posy= %d\n",fPositionX[i], fPositionY[i]);
2778  }
2779  </pre>
2780  <p>
2781  Example 16 - script SearchHR3.c:
2782  <p>
2783  <table border=solid>
2784  <tr><td> Peak # </td><td> Position </td><td> Sigma </td></tr>
2785  <tr><td> 1 </td><td> 118 </td><td> 26 </td></tr>
2786  <tr><td> 2 </td><td> 162 </td><td> 41 </td></tr>
2787  <tr><td> 3 </td><td> 310 </td><td> 4 </td></tr>
2788  <tr><td> 4 </td><td> 330 </td><td> 8 </td></tr>
2789  <tr><td> 5 </td><td> 482 </td><td> 22 </td></tr>
2790  <tr><td> 6 </td><td> 491 </td><td> 26 </td></tr>
2791  <tr><td> 7 </td><td> 740 </td><td> 21 </td></tr>
2792  <tr><td> 8 </td><td> 852 </td><td> 15 </td></tr>
2793  <tr><td> 9 </td><td> 954 </td><td> 12 </td></tr>
2794  <tr><td> 10 </td><td> 989 </td><td> 13 </td></tr>
2795  </table>
2796  <p>
2797  Table 4 Positions and sigma of peaks in the following examples.
2798  <p>
2799  <img width=600 height=328 src="gif/TSpectrum_Searching3.jpg">
2800  <p>
2801  Fig. 30 Influence of number of iterations (3-red, 10-blue, 100- green,
2802  1000-magenta), sigma=8, smoothing width=3.
2803  <p>
2804  <img width=600 height=321 src="gif/TSpectrum_Searching4.jpg">
2805  <p>
2806  Fig. 31 Influence of sigma (3-red, 8-blue, 20- green, 43-magenta),
2807  num. iter.=10, sm. width=3.
2808  <p>
2809  <img width=600 height=323 src="gif/TSpectrum_Searching5.jpg"></p>
2810  <p>
2811  Fig. 32 Influence smoothing width (0-red, 3-blue, 7- green, 20-magenta), num.
2812  iter.=10, sigma=8.
2813  <p>
2814  Script:
2815  <pre>
2816  // Example to illustrate the influence of number of iterations in deconvolution in high resolution peak searching function (class TSpectrum).
2817  // To execute this example, do
2818  // root > .x SearchHR3.C
2819 
2820  #include <TSpectrum>
2821 
2822  void SearchHR3() {
2823  Double_t fPositionX[100];
2824  Double_t fPositionY[100];
2825  Int_t fNPeaks = 0;
2826  Int_t i,nfound,bin;
2827  Double_t nbins = 1024,a;
2828  Double_t xmin = 0;
2829  Double_t xmax = nbins;
2830  Double_t * source = new Double_t[nbins];
2831  Double_t * dest = new Double_t[nbins];
2832  TH1F *h = new TH1F("h","Influence of # of iterations in deconvolution in peak searching",nbins,xmin,xmax);
2833  TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
2834  TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax);
2835  TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax);
2836  TH1F *d4 = new TH1F("d4","",nbins,xmin,xmax);
2837  TFile *f = new TFile("spectra\\TSpectrum.root");
2838  h=(TH1F*) f->Get("search3;1");
2839  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
2840  TCanvas *Search = gROOT->GetListOfCanvases()->FindObject("Search");
2841  if (!Search) Search = new TCanvas("Search","Search",10,10,1000,700);
2842  h->SetMaximum(1300);
2843  h->Draw("L");
2844  TSpectrum *s = new TSpectrum();
2845  nfound = s->SearchHighRes(source, dest, nbins, 8, 2, kTRUE, 3, kTRUE, 3);
2846  Double_t *xpeaks = s->GetPositionX();
2847  for (i = 0; i < nfound; i++) {
2848  a=xpeaks[i];
2849  bin = 1 + Int_t(a + 0.5);
2850  fPositionX[i] = h->GetBinCenter(bin);
2851  fPositionY[i] = h->GetBinContent(bin);
2852  }
2853  TPolyMarker * pm = (TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");
2854  if (pm) {
2855  h->GetListOfFunctions()->Remove(pm);
2856  delete pm;
2857  }
2858  pm = new TPolyMarker(nfound, fPositionX, fPositionY);
2859  h->GetListOfFunctions()->Add(pm);
2860  pm->SetMarkerStyle(23);
2861  pm->SetMarkerColor(kRed);
2862  pm->SetMarkerSize(1.3);
2863  for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,dest[i]);
2864  h->Draw("");
2865  d1->SetLineColor(kRed);
2866  d1->Draw("SAME");
2867  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
2868  s->SearchHighRes(source, dest, nbins, 8, 2, kTRUE, 10, kTRUE, 3);
2869  for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,dest[i]);
2870  d2->SetLineColor(kBlue);
2871  d2->Draw("SAME");
2872  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
2873  s->SearchHighRes(source, dest, nbins, 8, 2, kTRUE, 100, kTRUE, 3);
2874  for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,dest[i]);
2875  d3->SetLineColor(kGreen);
2876  d3->Draw("SAME");
2877  for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
2878  s->SearchHighRes(source, dest, nbins, 8, 2, kTRUE, 1000, kTRUE, 3);
2879  for (i = 0; i < nbins; i++) d4->SetBinContent(i + 1,dest[i]);
2880  d4->SetLineColor(kMagenta);
2881  d4->Draw("SAME");
2882  printf("Found %d candidate peaks\n",nfound);
2883  }
2884  </pre>
2885  End_Html */
2886 
2887  int i, j, numberIterations = (Int_t)(7 * sigma + 0.5);
2888  Double_t a, b, c;
2889  int k, lindex, posit, imin, imax, jmin, jmax, lh_gold, priz;
2890  Double_t lda, ldb, ldc, area, maximum, maximum_decon;
2891  int xmin, xmax, l, peak_index = 0, size_ext = ssize + 2 * numberIterations, shift = numberIterations, bw = 2, w;
2892  Double_t maxch;
2893  Double_t nom, nip, nim, sp, sm, plocha = 0;
2894  Double_t m0low=0,m1low=0,m2low=0,l0low=0,l1low=0,detlow,av,men;
2895  if (sigma < 1) {
2896  Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
2897  return 0;
2898  }
2899 
2900  if(threshold<=0 || threshold>=100){
2901  Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
2902  return 0;
2903  }
2904 
2905  j = (Int_t) (5.0 * sigma + 0.5);
2906  if (j >= PEAK_WINDOW / 2) {
2907  Error("SearchHighRes", "Too large sigma");
2908  return 0;
2909  }
2910 
2911  if (markov == true) {
2912  if (averWindow <= 0) {
2913  Error("SearchHighRes", "Averanging window must be positive");
2914  return 0;
2915  }
2916  }
2917 
2918  if(backgroundRemove == true){
2919  if(ssize < 2 * numberIterations + 1){
2920  Error("SearchHighRes", "Too large clipping window");
2921  return 0;
2922  }
2923  }
2924 
2925  k = int(2 * sigma+0.5);
2926  if(k >= 2){
2927  for(i = 0;i < k;i++){
2928  a = i,b = source[i];
2929  m0low += 1,m1low += a,m2low += a * a,l0low += b,l1low += a * b;
2930  }
2931  detlow = m0low * m2low - m1low * m1low;
2932  if(detlow != 0)
2933  l1low = (-l0low * m1low + l1low * m0low) / detlow;
2934 
2935  else
2936  l1low = 0;
2937  if(l1low > 0)
2938  l1low=0;
2939  }
2940 
2941  else{
2942  l1low = 0;
2943  }
2944 
2945  i = (Int_t)(7 * sigma + 0.5);
2946  i = 2 * i;
2947  Double_t *working_space = new Double_t [7 * (ssize + i)];
2948  for (j=0;j<7 * (ssize + i);j++) working_space[j] = 0;
2949  for(i = 0; i < size_ext; i++){
2950  if(i < shift){
2951  a = i - shift;
2952  working_space[i + size_ext] = source[0] + l1low * a;
2953  if(working_space[i + size_ext] < 0)
2954  working_space[i + size_ext]=0;
2955  }
2956 
2957  else if(i >= ssize + shift){
2958  a = i - (ssize - 1 + shift);
2959  working_space[i + size_ext] = source[ssize - 1];
2960  if(working_space[i + size_ext] < 0)
2961  working_space[i + size_ext]=0;
2962  }
2963 
2964  else
2965  working_space[i + size_ext] = source[i - shift];
2966  }
2967 
2968  if(backgroundRemove == true){
2969  for(i = 1; i <= numberIterations; i++){
2970  for(j = i; j < size_ext - i; j++){
2971  if(markov == false){
2972  a = working_space[size_ext + j];
2973  b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2974  if(b < a)
2975  a = b;
2976 
2977  working_space[j]=a;
2978  }
2979 
2980  else{
2981  a = working_space[size_ext + j];
2982  av = 0;
2983  men = 0;
2984  for (w = j - bw; w <= j + bw; w++){
2985  if ( w >= 0 && w < size_ext){
2986  av += working_space[size_ext + w];
2987  men +=1;
2988  }
2989  }
2990  av = av / men;
2991  b = 0;
2992  men = 0;
2993  for (w = j - i - bw; w <= j - i + bw; w++){
2994  if ( w >= 0 && w < size_ext){
2995  b += working_space[size_ext + w];
2996  men +=1;
2997  }
2998  }
2999  b = b / men;
3000  c = 0;
3001  men = 0;
3002  for (w = j + i - bw; w <= j + i + bw; w++){
3003  if ( w >= 0 && w < size_ext){
3004  c += working_space[size_ext + w];
3005  men +=1;
3006  }
3007  }
3008  c = c / men;
3009  b = (b + c) / 2;
3010  if (b < a)
3011  av = b;
3012  working_space[j]=av;
3013  }
3014  }
3015  for(j = i; j < size_ext - i; j++)
3016  working_space[size_ext + j] = working_space[j];
3017  }
3018  for(j = 0;j < size_ext; j++){
3019  if(j < shift){
3020  a = j - shift;
3021  b = source[0] + l1low * a;
3022  if (b < 0) b = 0;
3023  working_space[size_ext + j] = b - working_space[size_ext + j];
3024  }
3025 
3026  else if(j >= ssize + shift){
3027  a = j - (ssize - 1 + shift);
3028  b = source[ssize - 1];
3029  if (b < 0) b = 0;
3030  working_space[size_ext + j] = b - working_space[size_ext + j];
3031  }
3032 
3033  else{
3034  working_space[size_ext + j] = source[j - shift] - working_space[size_ext + j];
3035  }
3036  }
3037  for(j = 0;j < size_ext; j++){
3038  if(working_space[size_ext + j] < 0) working_space[size_ext + j] = 0;
3039  }
3040  }
3041 
3042  for(i = 0; i < size_ext; i++){
3043  working_space[i + 6*size_ext] = working_space[i + size_ext];
3044  }
3045 
3046  if(markov == true){
3047  for(j = 0; j < size_ext; j++)
3048  working_space[2 * size_ext + j] = working_space[size_ext + j];
3049  xmin = 0,xmax = size_ext - 1;
3050  for(i = 0, maxch = 0; i < size_ext; i++){
3051  working_space[i] = 0;
3052  if(maxch < working_space[2 * size_ext + i])
3053  maxch = working_space[2 * size_ext + i];
3054  plocha += working_space[2 * size_ext + i];
3055  }
3056  if(maxch == 0) {
3057  delete [] working_space;
3058  return 0;
3059  }
3060 
3061  nom = 1;
3062  working_space[xmin] = 1;
3063  for(i = xmin; i < xmax; i++){
3064  nip = working_space[2 * size_ext + i] / maxch;
3065  nim = working_space[2 * size_ext + i + 1] / maxch;
3066  sp = 0,sm = 0;
3067  for(l = 1; l <= averWindow; l++){
3068  if((i + l) > xmax)
3069  a = working_space[2 * size_ext + xmax] / maxch;
3070 
3071  else
3072  a = working_space[2 * size_ext + i + l] / maxch;
3073 
3074  b = a - nip;
3075  if(a + nip <= 0)
3076  a=1;
3077 
3078  else
3079  a = TMath::Sqrt(a + nip);
3080 
3081  b = b / a;
3082  b = TMath::Exp(b);
3083  sp = sp + b;
3084  if((i - l + 1) < xmin)
3085  a = working_space[2 * size_ext + xmin] / maxch;
3086 
3087  else
3088  a = working_space[2 * size_ext + i - l + 1] / maxch;
3089 
3090  b = a - nim;
3091  if(a + nim <= 0)
3092  a = 1;
3093 
3094  else
3095  a = TMath::Sqrt(a + nim);
3096 
3097  b = b / a;
3098  b = TMath::Exp(b);
3099  sm = sm + b;
3100  }
3101  a = sp / sm;
3102  a = working_space[i + 1] = working_space[i] * a;
3103  nom = nom + a;
3104  }
3105  for(i = xmin; i <= xmax; i++){
3106  working_space[i] = working_space[i] / nom;
3107  }
3108  for(j = 0; j < size_ext; j++)
3109  working_space[size_ext + j] = working_space[j] * plocha;
3110  for(j = 0; j < size_ext; j++){
3111  working_space[2 * size_ext + j] = working_space[size_ext + j];
3112  }
3113  if(backgroundRemove == true){
3114  for(i = 1; i <= numberIterations; i++){
3115  for(j = i; j < size_ext - i; j++){
3116  a = working_space[size_ext + j];
3117  b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
3118  if(b < a)
3119  a = b;
3120  working_space[j] = a;
3121  }
3122  for(j = i; j < size_ext - i; j++)
3123  working_space[size_ext + j] = working_space[j];
3124  }
3125  for(j = 0; j < size_ext; j++){
3126  working_space[size_ext + j] = working_space[2 * size_ext + j] - working_space[size_ext + j];
3127  }
3128  }
3129  }
3130 //deconvolution starts
3131  area = 0;
3132  lh_gold = -1;
3133  posit = 0;
3134  maximum = 0;
3135 //generate response vector
3136  for(i = 0; i < size_ext; i++){
3137  lda = (Double_t)i - 3 * sigma;
3138  lda = lda * lda / (2 * sigma * sigma);
3139  j = (Int_t)(1000 * TMath::Exp(-lda));
3140  lda = j;
3141  if(lda != 0)
3142  lh_gold = i + 1;
3143 
3144  working_space[i] = lda;
3145  area = area + lda;
3146  if(lda > maximum){
3147  maximum = lda;
3148  posit = i;
3149  }
3150  }
3151 //read source vector
3152  for(i = 0; i < size_ext; i++)
3153  working_space[2 * size_ext + i] = TMath::Abs(working_space[size_ext + i]);
3154 //create matrix at*a(vector b)
3155  i = lh_gold - 1;
3156  if(i > size_ext)
3157  i = size_ext;
3158 
3159  imin = -i,imax = i;
3160  for(i = imin; i <= imax; i++){
3161  lda = 0;
3162  jmin = 0;
3163  if(i < 0)
3164  jmin = -i;
3165  jmax = lh_gold - 1 - i;
3166  if(jmax > (lh_gold - 1))
3167  jmax = lh_gold - 1;
3168 
3169  for(j = jmin;j <= jmax; j++){
3170  ldb = working_space[j];
3171  ldc = working_space[i + j];
3172  lda = lda + ldb * ldc;
3173  }
3174  working_space[size_ext + i - imin] = lda;
3175  }
3176 //create vector p
3177  i = lh_gold - 1;
3178  imin = -i,imax = size_ext + i - 1;
3179  for(i = imin; i <= imax; i++){
3180  lda = 0;
3181  for(j = 0; j <= (lh_gold - 1); j++){
3182  ldb = working_space[j];
3183  k = i + j;
3184  if(k >= 0 && k < size_ext){
3185  ldc = working_space[2 * size_ext + k];
3186  lda = lda + ldb * ldc;
3187  }
3188 
3189  }
3190  working_space[4 * size_ext + i - imin] = lda;
3191  }
3192 //move vector p
3193  for(i = imin; i <= imax; i++)
3194  working_space[2 * size_ext + i - imin] = working_space[4 * size_ext + i - imin];
3195 //initialization of resulting vector
3196  for(i = 0; i < size_ext; i++)
3197  working_space[i] = 1;
3198 //START OF ITERATIONS
3199  for(lindex = 0; lindex < deconIterations; lindex++){
3200  for(i = 0; i < size_ext; i++){
3201  if(TMath::Abs(working_space[2 * size_ext + i]) > 0.00001 && TMath::Abs(working_space[i]) > 0.00001){
3202  lda=0;
3203  jmin = lh_gold - 1;
3204  if(jmin > i)
3205  jmin = i;
3206 
3207  jmin = -jmin;
3208  jmax = lh_gold - 1;
3209  if(jmax > (size_ext - 1 - i))
3210  jmax=size_ext-1-i;
3211 
3212  for(j = jmin; j <= jmax; j++){
3213  ldb = working_space[j + lh_gold - 1 + size_ext];
3214  ldc = working_space[i + j];
3215  lda = lda + ldb * ldc;
3216  }
3217  ldb = working_space[2 * size_ext + i];
3218  if(lda != 0)
3219  lda = ldb / lda;
3220 
3221  else
3222  lda = 0;
3223 
3224  ldb = working_space[i];
3225  lda = lda * ldb;
3226  working_space[3 * size_ext + i] = lda;
3227  }
3228  }
3229  for(i = 0; i < size_ext; i++){
3230  working_space[i] = working_space[3 * size_ext + i];
3231  }
3232  }
3233 //shift resulting spectrum
3234  for(i=0;i<size_ext;i++){
3235  lda = working_space[i];
3236  j = i + posit;
3237  j = j % size_ext;
3238  working_space[size_ext + j] = lda;
3239  }
3240 //write back resulting spectrum
3241  maximum = 0, maximum_decon = 0;
3242  j = lh_gold - 1;
3243  for(i = 0; i < size_ext - j; i++){
3244  if(i >= shift && i < ssize + shift){
3245  working_space[i] = area * working_space[size_ext + i + j];
3246  if(maximum_decon < working_space[i])
3247  maximum_decon = working_space[i];
3248  if(maximum < working_space[6 * size_ext + i])
3249  maximum = working_space[6 * size_ext + i];
3250  }
3251 
3252  else
3253  working_space[i] = 0;
3254  }
3255  lda=1;
3256  if(lda>threshold)
3257  lda=threshold;
3258  lda=lda/100;
3259 
3260 //searching for peaks in deconvolved spectrum
3261  for(i = 1; i < size_ext - 1; i++){
3262  if(working_space[i] > working_space[i - 1] && working_space[i] > working_space[i + 1]){
3263  if(i >= shift && i < ssize + shift){
3264  if(working_space[i] > lda*maximum_decon && working_space[6 * size_ext + i] > threshold * maximum / 100.0){
3265  for(j = i - 1, a = 0, b = 0; j <= i + 1; j++){
3266  a += (Double_t)(j - shift) * working_space[j];
3267  b += working_space[j];
3268  }
3269  a = a / b;
3270  if(a < 0)
3271  a = 0;
3272 
3273  if(a >= ssize)
3274  a = ssize - 1;
3275  if(peak_index == 0){
3276  fPositionX[0] = a;
3277  peak_index = 1;
3278  }
3279 
3280  else{
3281  for(j = 0, priz = 0; j < peak_index && priz == 0; j++){
3282  if(working_space[6 * size_ext + shift + (Int_t)a] > working_space[6 * size_ext + shift + (Int_t)fPositionX[j]])
3283  priz = 1;
3284  }
3285  if(priz == 0){
3286  if(j < fMaxPeaks){
3287  fPositionX[j] = a;
3288  }
3289  }
3290 
3291  else{
3292  for(k = peak_index; k >= j; k--){
3293  if(k < fMaxPeaks){
3294  fPositionX[k] = fPositionX[k - 1];
3295  }
3296  }
3297  fPositionX[j - 1] = a;
3298  }
3299  if(peak_index < fMaxPeaks)
3300  peak_index += 1;
3301  }
3302  }
3303  }
3304  }
3305  }
3306 
3307  for(i = 0; i < ssize; i++) destVector[i] = working_space[i + shift];
3308  delete [] working_space;
3309  fNPeaks = peak_index;
3310  if(peak_index == fMaxPeaks)
3311  Warning("SearchHighRes", "Peak buffer full");
3312  return fNPeaks;
3313 }
3314 
3315 
3316 ////////////////////////////////////////////////////////////////////////////////
3317 
3318 Int_t TSpectrum::Search1HighRes(Double_t *source,Double_t *destVector, int ssize,
3319  Double_t sigma, Double_t threshold,
3320  bool backgroundRemove,int deconIterations,
3321  bool markov, int averWindow)
3322 {
3323  /* Begin_Html
3324  Old name of SearcHighRes introduced for back compatibility.
3325  This function will be removed after the June 2006 release
3326  End_Html */
3327 
3328  return SearchHighRes(source,destVector,ssize,sigma,threshold,backgroundRemove,
3329  deconIterations,markov,averWindow);
3330 }
3331 
3332 
3333 ////////////////////////////////////////////////////////////////////////////////
3334 
3335 Int_t TSpectrum::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
3336 {
3337  /* Begin_Html
3338  Static function, interface to TSpectrum::Search.
3339  End_Html */
3340 
3341  TSpectrum s;
3342  return s.Search(hist,sigma,option,threshold);
3343 }
3344 
3345 
3346 ////////////////////////////////////////////////////////////////////////////////
3347 
3348 TH1 *TSpectrum::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
3349 {
3350  /* Begin_Html
3351  Static function, interface to TSpectrum::Background.
3352  End_Html */
3353 
3354  TSpectrum s;
3355  return s.Background(hist,niter,option);
3356 }
TAxis::GetFirst
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:429
TH1::GetListOfFunctions
TList * GetListOfFunctions() const
Definition: TH1.h:244
TSpectrum::Search1HighRes
Int_t Search1HighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
Definition: TSpectrum.cxx:3318
xmin
float xmin
Definition: THbookFile.cxx:93
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4629
Option_t
const char Option_t
Definition: RtypesCore.h:62
kRed
Definition: Rtypes.h:61
TString::ReplaceAll
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:635
TH1::GetDimension
virtual Int_t GetDimension() const
Definition: TH1.h:283
h
TH1 * h
Definition: legend2.C:5
TSpectrum::StaticBackground
static TH1 * StaticBackground(const TH1 *hist, Int_t niter=20, Option_t *option="")
Definition: TSpectrum.cxx:3348
TSpectrum::DeconvolutionRL
const char * DeconvolutionRL(Double_t *source, const Double_t *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
Definition: TSpectrum.cxx:2114
background
Double_t background(Double_t *x, Double_t *par)
Definition: GAMinTutorial.cxx:23
TString
Basic string class.
Definition: TString.h:137
ClassImp
ClassImp(TSpectrum) TSpectrum
Definition: TSpectrum.cxx:61
TString::ToLower
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1088
Int_t
int Int_t
Definition: RtypesCore.h:41
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
a
TArc * a
Definition: textangle.C:12
kFALSE
const Bool_t kFALSE
Definition: Rtypes.h:92
TList::FindObject
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:496
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
TH1::Clone
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2565
TString::Data
const char * Data() const
Definition: TString.h:349
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:33
TSpectrum::fPositionX
Double_t * fPositionX
Definition: TSpectrum.h:30
TSpectrum::SearchHighRes
Int_t SearchHighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
Definition: TSpectrum.cxx:2640
TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
Definition: TAttMarker.h:51
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
TObject::Delete
virtual void Delete(Option_t *option="")
Delete this object.
Definition: TObject.cxx:228
TH1::GetBinCenter
virtual Double_t GetBinCenter(Int_t bin) const
return bin center for 1D historam Better to use h1.GetXaxis().GetBinCenter(bin)
Definition: TH1.cxx:8470
TSpectrum::Search
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05)
Definition: TSpectrum.cxx:277
TSpectrum::~TSpectrum
virtual ~TSpectrum()
Definition: TSpectrum.cxx:113
TList::Remove
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:674
ROOT::Math::VectorUtil::boost
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only re...
Definition: VectorUtil.h:329
TSpectrum::fHistogram
TH1 * fHistogram
Definition: TSpectrum.h:33
TSpectrum::fPositionY
Double_t * fPositionY
Definition: TSpectrum.h:31
TSpectrum::Deconvolution
const char * Deconvolution(Double_t *source, const Double_t *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
Definition: TSpectrum.cxx:1666
l
TLine * l
Definition: textangle.C:4
TSpectrum::SetAverageWindow
static void SetAverageWindow(Int_t w=3)
Definition: TSpectrum.cxx:128
TSpectrum::fgAverageWindow
static Int_t fgAverageWindow
Definition: TSpectrum.h:34
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
Definition: TAttMarker.h:53
xmax
float xmax
Definition: THbookFile.cxx:93
PEAK_WINDOW
#define PEAK_WINDOW
Definition: TSpectrum.cxx:60
TSpectrum::SetDeconIterations
static void SetDeconIterations(Int_t n=3)
Definition: TSpectrum.cxx:141
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
Definition: TAttMarker.h:54
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:495
TPolyMarker
A PolyMarker is defined by an array on N points in a 2-D space.
Definition: TPolyMarker.h:37
Double_t
double Double_t
Definition: RtypesCore.h:55
TSpectrum::Background
virtual TH1 * Background(const TH1 *hist, Int_t niter=20, Option_t *option="")
Definition: TSpectrum.cxx:154
TSpectrum::fgIterations
static Int_t fgIterations
Definition: TSpectrum.h:35
printf
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
TSpectrum::fMaxPeaks
Int_t fMaxPeaks
Definition: TSpectrum.h:27
TH1
The TH1 histogram class.
Definition: TH1.h:80
TSpectrum::SetResolution
void SetResolution(Double_t resolution=1)
Definition: TSpectrum.cxx:397
TSpectrum
Advanced Spectra Processing.
Definition: TSpectrum.h:20
TAxis::GetLast
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:440
hbname
#define hbname
Definition: THbookFile.cxx:123
TSpectrum::fPosition
Double_t * fPosition
Definition: TSpectrum.h:29
sample_config.c
tuple c
Definition: sample_config.py:3
TSpectrum::Print
virtual void Print(Option_t *option="") const
Print TNamed name and title.
Definition: TSpectrum.cxx:262
TList::Add
virtual void Add(TObject *obj)
Definition: TList.h:81
dest
#define dest(otri, vertexptr)
Definition: triangle.c:1040
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
TSpectrum::StaticSearch
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
Definition: TSpectrum.cxx:3335
gPad
#define gPad
Definition: TVirtualPad.h:288
TSpectrum::fNPeaks
Int_t fNPeaks
Definition: TSpectrum.h:28
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
kTRUE
const Bool_t kTRUE
Definition: Rtypes.h:91
TSpectrum::SmoothMarkov
const char * SmoothMarkov(Double_t *source, Int_t ssize, Int_t averWindow)
Definition: TSpectrum.cxx:1522
n
const Int_t n
Definition: legend1.C:16
TSpectrum::Unfolding
const char * Unfolding(Double_t *source, const Double_t **respMatrix, Int_t ssizex, Int_t ssizey, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
Definition: TSpectrum.cxx:2378
TSpectrum::fResolution
Double_t fResolution
Definition: TSpectrum.h:32
TH1::GetXaxis
TAxis * GetXaxis()
Definition: TH1.h:319
TObject::Warning
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:904