Logo ROOT   6.10/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 /** \class TSpectrum
12  \ingroup Spectrum
13  \brief Advanced Spectra Processing
14  \author Miroslav Morhac
15 
16  This class contains advanced spectra processing functions for:
17 
18  - One-dimensional background estimation
19  - One-dimensional smoothing
20  - One-dimensional deconvolution
21  - One-dimensional peak search
22 
23  The algorithms in this class have been published in the following references:
24 
25  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.
26  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.
27  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.
28 
29  These NIM papers are also available as doc or ps files from:
30 
31  - [Spectrum.doc](https://root.cern.ch/download/Spectrum.doc)
32  - [SpectrumDec.ps.gz](https://root.cern.ch/download/SpectrumDec.ps.gz)
33  - [SpectrumSrc.ps.gz](https://root.cern.ch/download/SpectrumSrc.ps.gz)
34  - [SpectrumBck.ps.gz](https://root.cern.ch/download/SpectrumBck.ps.gz)
35 
36  See also the
37  [online documentation](https://root.cern.ch/guides/tspectrum-manual) and
38  [tutorials](https://root.cern.ch/doc/master/group__tutorial__spectrum.html).
39 */
40 
43 
44 #define PEAK_WINDOW 1024
46 
47 ////////////////////////////////////////////////////////////////////////////////
48 /// Constructor.
49 
50 TSpectrum::TSpectrum() :TNamed("Spectrum", "Miroslav Morhac peak finder")
51 {
52  Int_t n = 100;
53  fMaxPeaks = n;
54  fPosition = new Double_t[n];
55  fPositionX = new Double_t[n];
56  fPositionY = new Double_t[n];
57  fResolution = 1;
58  fHistogram = 0;
59  fNPeaks = 0;
60 }
61 
62 ////////////////////////////////////////////////////////////////////////////////
63 /// - maxpositions: maximum number of peaks
64 /// - resolution: *NOT USED* determines resolution of the neighbouring peaks
65 /// default value is 1 correspond to 3 sigma distance
66 /// between peaks. Higher values allow higher resolution
67 /// (smaller distance between peaks.
68 /// May be set later through SetResolution.
69 
70 TSpectrum::TSpectrum(Int_t maxpositions, Double_t resolution)
71  :TNamed("Spectrum", "Miroslav Morhac peak finder")
72 {
73  Int_t n = maxpositions;
74  if (n <= 0) n = 1;
75  fMaxPeaks = n;
76  fPosition = new Double_t[n];
77  fPositionX = new Double_t[n];
78  fPositionY = new Double_t[n];
79  fHistogram = 0;
80  fNPeaks = 0;
81  SetResolution(resolution);
82 }
83 
84 ////////////////////////////////////////////////////////////////////////////////
85 /// Destructor.
86 
88 {
89  delete [] fPosition;
90  delete [] fPositionX;
91  delete [] fPositionY;
92  delete fHistogram;
93 }
94 
95 ////////////////////////////////////////////////////////////////////////////////
96 /// Static function: Set average window of searched peaks
97 /// (see TSpectrum::SearchHighRes).
98 
100 {
101  fgAverageWindow = w;
102 }
103 
104 ////////////////////////////////////////////////////////////////////////////////
105 /// Static function: Set max number of decon iterations in deconvolution
106 /// operation (see TSpectrum::SearchHighRes).
107 
109 {
110  fgIterations = n;
111 }
112 
113 ////////////////////////////////////////////////////////////////////////////////
114 /// One-dimensional background estimation function.
115 ///
116 /// This function calculates the background spectrum in the input histogram h.
117 /// The background is returned as a histogram.
118 ///
119 /// #### Parameters:
120 ///
121 /// - h: input 1-d histogram
122 /// - numberIterations, (default value = 20)
123 /// Increasing numberIterations make the result smoother and lower.
124 /// - option: may contain one of the following options:
125 ///
126 /// - to set the direction parameter
127 /// "BackIncreasingWindow". By default the direction is BackDecreasingWindow
128 /// - filterOrder-order of clipping filter, (default "BackOrder2")
129 /// -possible values= "BackOrder4"
130 /// "BackOrder6"
131 /// "BackOrder8"
132 /// - "nosmoothing"- if selected, the background is not smoothed
133 /// By default the background is smoothed.
134 /// - smoothWindow-width of smoothing window, (default is "BackSmoothing3")
135 /// -possible values= "BackSmoothing5"
136 /// "BackSmoothing7"
137 /// "BackSmoothing9"
138 /// "BackSmoothing11"
139 /// "BackSmoothing13"
140 /// "BackSmoothing15"
141 /// - "Compton" if selected the estimation of Compton edge
142 /// will be included.
143 /// - "same" : if this option is specified, the resulting background
144 /// histogram is superimposed on the picture in the current pad.
145 ///
146 /// NOTE that the background is only evaluated in the current range of h.
147 /// ie, if h has a bin range (set via `h->GetXaxis()->SetRange(binmin,binmax)`,
148 /// the returned histogram will be created with the same number of bins
149 /// as the input histogram h, but only bins from `binmin` to `binmax` will be filled
150 /// with the estimated background.
151 
152 TH1 *TSpectrum::Background(const TH1 * h, int numberIterations,
153  Option_t * option)
154 {
155  if (h == 0) return 0;
156  Int_t dimension = h->GetDimension();
157  if (dimension > 1) {
158  Error("Search", "Only implemented for 1-d histograms");
159  return 0;
160  }
161  TString opt = option;
162  opt.ToLower();
163 
164  //set options
165  Int_t direction = kBackDecreasingWindow;
166  if (opt.Contains("backincreasingwindow")) direction = kBackIncreasingWindow;
167  Int_t filterOrder = kBackOrder2;
168  if (opt.Contains("backorder4")) filterOrder = kBackOrder4;
169  if (opt.Contains("backorder6")) filterOrder = kBackOrder6;
170  if (opt.Contains("backorder8")) filterOrder = kBackOrder8;
171  Bool_t smoothing = kTRUE;
172  if (opt.Contains("nosmoothing")) smoothing = kFALSE;
173  Int_t smoothWindow = kBackSmoothing3;
174  if (opt.Contains("backsmoothing5")) smoothWindow = kBackSmoothing5;
175  if (opt.Contains("backsmoothing7")) smoothWindow = kBackSmoothing7;
176  if (opt.Contains("backsmoothing9")) smoothWindow = kBackSmoothing9;
177  if (opt.Contains("backsmoothing11")) smoothWindow = kBackSmoothing11;
178  if (opt.Contains("backsmoothing13")) smoothWindow = kBackSmoothing13;
179  if (opt.Contains("backsmoothing15")) smoothWindow = kBackSmoothing15;
180  Bool_t compton = kFALSE;
181  if (opt.Contains("compton")) compton = kTRUE;
182 
183  Int_t first = h->GetXaxis()->GetFirst();
184  Int_t last = h->GetXaxis()->GetLast();
185  Int_t size = last-first+1;
186  Int_t i;
187  Double_t * source = new Double_t[size];
188  for (i = 0; i < size; i++) source[i] = h->GetBinContent(i + first);
189 
190  //find background (source is input and in output contains the background
191  Background(source,size,numberIterations, direction, filterOrder,smoothing,
192  smoothWindow,compton);
193 
194  //create output histogram containing background
195  //only bins in the range of the input histogram are filled
196  Int_t nch = strlen(h->GetName());
197  char *hbname = new char[nch+20];
198  snprintf(hbname,nch+20,"%s_background",h->GetName());
199  TH1 *hb = (TH1*)h->Clone(hbname);
200  hb->Reset();
201  hb->GetListOfFunctions()->Delete();
202  hb->SetLineColor(2);
203  for (i=0; i< size; i++) hb->SetBinContent(i+first,source[i]);
204  hb->SetEntries(size);
205 
206  //if option "same is specified, draw the result in the pad
207  if (opt.Contains("same")) {
208  if (gPad) delete gPad->GetPrimitive(hbname);
209  hb->Draw("same");
210  }
211  delete [] source;
212  delete [] hbname;
213  return hb;
214 }
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 /// Print the array of positions.
218 
220 {
221  printf("\nNumber of positions = %d\n",fNPeaks);
222  for (Int_t i=0;i<fNPeaks;i++) {
223  printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
224  }
225 }
226 
227 ////////////////////////////////////////////////////////////////////////////////
228 /// One-dimensional peak search function
229 ///
230 /// This function searches for peaks in source spectrum in hin
231 /// The number of found peaks and their positions are written into
232 /// the members fNpeaks and fPositionX.
233 /// The search is performed in the current histogram range.
234 ///
235 /// #### Parameters:
236 ///
237 /// - hin: pointer to the histogram of source spectrum
238 /// - sigma: sigma of searched peaks, for details we refer to manual
239 /// - threshold: (default=0.05) peaks with amplitude less than
240 /// threshold*highest_peak are discarded. 0<threshold<1
241 ///
242 /// By default, the background is removed before deconvolution.
243 /// Specify the option "nobackground" to not remove the background.
244 ///
245 /// By default the "Markov" chain algorithm is used.
246 /// Specify the option "noMarkov" to disable this algorithm
247 /// Note that by default the source spectrum is replaced by a new spectrum
248 ///
249 /// By default a polymarker object is created and added to the list of
250 /// functions of the histogram. The histogram is drawn with the specified
251 /// option and the polymarker object drawn on top of the histogram.
252 /// The polymarker coordinates correspond to the npeaks peaks found in
253 /// the histogram.
254 ///
255 /// A pointer to the polymarker object can be retrieved later via:
256 /// ~~~ {.cpp}
257 /// TList *functions = hin->GetListOfFunctions();
258 /// TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
259 /// ~~~
260 /// Specify the option "goff" to disable the storage and drawing of the
261 /// polymarker.
262 ///
263 /// To disable the final drawing of the histogram with the search results (in case
264 /// you want to draw it yourself) specify "nodraw" in the options parameter.
265 
267  Double_t threshold)
268 {
269  if (hin == 0) return 0;
270  Int_t dimension = hin->GetDimension();
271  if (dimension > 2) {
272  Error("Search", "Only implemented for 1-d and 2-d histograms");
273  return 0;
274  }
275  if (threshold <=0 || threshold >= 1) {
276  Warning("Search","threshold must 0<threshold<1, threshold=0.05 assumed");
277  threshold = 0.05;
278  }
279  TString opt = option;
280  opt.ToLower();
282  if (opt.Contains("nobackground")) {
283  background = kFALSE;
284  opt.ReplaceAll("nobackground","");
285  }
286  Bool_t markov = kTRUE;
287  if (opt.Contains("nomarkov")) {
288  markov = kFALSE;
289  opt.ReplaceAll("nomarkov","");
290  }
291  Bool_t draw = kTRUE;
292  if (opt.Contains("nodraw")) {
293  draw = kFALSE;
294  opt.ReplaceAll("nodraw","");
295  }
296  if (dimension == 1) {
297  Int_t first = hin->GetXaxis()->GetFirst();
298  Int_t last = hin->GetXaxis()->GetLast();
299  Int_t size = last-first+1;
300  Int_t i, bin, npeaks;
301  Double_t * source = new Double_t[size];
302  Double_t * dest = new Double_t[size];
303  for (i = 0; i < size; i++) source[i] = hin->GetBinContent(i + first);
304  if (sigma < 1) {
305  sigma = size/fMaxPeaks;
306  if (sigma < 1) sigma = 1;
307  if (sigma > 8) sigma = 8;
308  }
309  npeaks = SearchHighRes(source, dest, size, sigma, 100*threshold,
310  background, fgIterations, markov, fgAverageWindow);
311 
312  for (i = 0; i < npeaks; i++) {
313  bin = first + Int_t(fPositionX[i] + 0.5);
314  fPositionX[i] = hin->GetBinCenter(bin);
315  fPositionY[i] = hin->GetBinContent(bin);
316  }
317  delete [] source;
318  delete [] dest;
319 
320  if (opt.Contains("goff"))
321  return npeaks;
322  if (!npeaks) return 0;
323  TPolyMarker * pm =
324  (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
325  if (pm) {
326  hin->GetListOfFunctions()->Remove(pm);
327  delete pm;
328  }
329  pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
330  hin->GetListOfFunctions()->Add(pm);
331  pm->SetMarkerStyle(23);
332  pm->SetMarkerColor(kRed);
333  pm->SetMarkerSize(1.3);
334  opt.ReplaceAll(" ","");
335  opt.ReplaceAll(",","");
336  if (draw)
337  ((TH1*)hin)->Draw(opt.Data());
338  return npeaks;
339  }
340  return 0;
341 }
342 
343 ////////////////////////////////////////////////////////////////////////////////
344 /// *NOT USED*
345 /// resolution: determines resolution of the neighbouring peaks
346 /// default value is 1 correspond to 3 sigma distance
347 /// between peaks. Higher values allow higher resolution
348 /// (smaller distance between peaks.
349 /// May be set later through SetResolution.
350 
352 {
353  if (resolution > 1)
354  fResolution = resolution;
355  else
356  fResolution = 1;
357 }
358 
359 ////////////////////////////////////////////////////////////////////////////////
360 /// This function calculates background spectrum from source spectrum.
361 /// The result is placed in the vector pointed by spe1945ctrum pointer.
362 /// The goal is to separate the useful information (peaks) from useless
363 /// information (background).
364 ///
365 /// - method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping
366 /// algorithm.
367 /// - new value in the channel "i" is calculated
368 ///
369 /// \f[
370 /// v_p(i) = min \left\{ v_{p-1}(i)^{\frac{\left[v_{p-1}(i+p)+v_{p-1}(i-p)\right]}{2}} \right\}
371 /// \f]
372 ///
373 /// where p = 1, 2, ..., numberIterations. In fact it represents second order
374 /// difference filter (-1,2,-1).
375 ///
376 /// One can also change the
377 /// direction of the change of the clipping window, the order of the clipping
378 /// filter, to include smoothing, to set width of smoothing window and to include
379 /// the estimation of Compton edges. On successful completion it returns 0. On
380 /// error it returns pointer to the string describing error.
381 ///
382 /// #### Parameters:
383 ///
384 /// - spectrum: pointer to the vector of source spectrum
385 /// - ssize: length of the spectrum vector
386 /// - numberIterations: maximal width of clipping window,
387 /// - direction: direction of change of clipping window.
388 /// Possible values: kBackIncreasingWindow, kBackDecreasingWindow
389 /// - filterOrder: order of clipping filter.
390 /// Possible values: kBackOrder2, kBackOrder4, kBackOrder6, kBackOrder8
391 /// - smoothing: logical variable whether the smoothing operation in the
392 /// estimation of background will be included.
393 /// Possible values: kFALSE, kTRUE
394 /// - smoothWindow: width of smoothing window.
395 /// Possible values: kBackSmoothing3, kBackSmoothing5, kBackSmoothing7,
396 /// kBackSmoothing9, kBackSmoothing11, kBackSmoothing13, kBackSmoothing15.
397 /// - compton: logical variable whether the estimation of Compton edge will be
398 /// included. Possible values: kFALSE, kTRUE.
399 ///
400 /// #### References:
401 ///
402 /// 1. C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the
403 /// quantitative analysis of PIXE spectra in geoscience applications. NIM, B34
404 /// (1988), 396-402.
405 ///
406 /// 2. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo:
407 /// Background elimination methods for multidimensional gamma-ray spectra. NIM,
408 /// A401 (1997) 113-132.
409 ///
410 /// 3. D. D. Burgess, R. J. Tervo: Background estimation for gamma-ray
411 /// spectroscopy. NIM 214 (1983), 431-434.
412 ///
413 /// ### Example 1 script Background_incr.C:
414 ///
415 /// Example of the estimation of background for number of iterations=6.
416 /// Original spectrum is shown in black color, estimated background in red color.
417 ///
418 /// Begin_Macro(source)
419 /// ../../../tutorials/spectrum/Background_incr.C
420 /// End_Macro
421 ///
422 /// ### Example 2 script Background_decr.C:
423 ///
424 /// In Example 1. one can notice that at the edges of the peaks the estimated
425 /// background goes under the peaks. An alternative approach is to decrease the
426 /// clipping window from a given value numberIterations to the value of one, which
427 /// is presented in this example.
428 ///
429 /// Example of the estimation of background for numberIterations=6 using
430 /// decreasing clipping window algorithm. Original spectrum is shown in black
431 /// color, estimated background in red color.
432 ///
433 /// Begin_Macro(source)
434 /// ../../../tutorials/spectrum/Background_decr.C
435 /// End_Macro
436 ///
437 /// ### Example 3 script Background_width.C:
438 ///
439 /// The question is how to choose the width of the clipping window, i.e.,
440 /// numberIterations parameter. The influence of this parameter on the estimated
441 /// background is illustrated in Example 3.
442 ///
443 /// Example of the influence of clipping window width on the estimated background
444 /// for numberIterations=4 (red line), 6 (orange line) 8 (green line) using decreasing
445 /// clipping window algorithm.
446 ///
447 /// in general one should set this parameter so that the value
448 /// 2*numberIterations+1 was greater than the widths of preserved objects (peaks).
449 ///
450 /// Begin_Macro(source)
451 /// ../../../tutorials/spectrum/Background_width.C
452 /// End_Macro
453 ///
454 /// ### Example 4 script Background_width2.C:
455 ///
456 /// another example for very complex spectrum is given here.
457 ///
458 /// Example of the influence of clipping window width on the estimated background
459 /// for numberIterations=10 (red line), 20 (blue line), 30 (green line) and
460 /// 40 (magenta line) using decreasing clipping window algorithm.
461 ///
462 /// Begin_Macro(source)
463 /// ../../../tutorials/spectrum/Background_width2.C
464 /// End_Macro
465 ///
466 /// ### Example 5 script Background_order.C:
467 ///
468 /// Second order difference filter removes linear (quasi-linear) background and
469 /// preserves symmetrical peaks. However if the shape of the background is more
470 /// complex one can employ higher-order clipping filters.
471 ///
472 /// Example of the influence of clipping filter difference order on the estimated
473 /// background for fNnumberIterations=40, 2-nd order red line, 4-th order blue line,
474 /// 6-th order green line and 8-th order magenta line, and using decreasing
475 /// clipping window algorithm.
476 ///
477 /// Begin_Macro(source)
478 /// ../../../tutorials/spectrum/Background_order.C
479 /// End_Macro
480 ///
481 /// ### Example 6 script Background_smooth.C:
482 ///
483 /// The estimate of the background can be influenced by noise present in the
484 /// spectrum. We proposed the algorithm of the background estimate with
485 /// simultaneous smoothing. In the original algorithm without smoothing, the
486 /// estimated background snatches the lower spikes in the noise. Consequently,
487 /// the areas of peaks are biased by this error.
488 ///
489 /// \image html TSpectrum_Background_smooth1.jpg Principle of background estimation algorithm with simultaneous smoothing.
490 ///
491 /// Begin_Macro(source)
492 /// ../../../tutorials/spectrum/Background_smooth.C
493 /// End_Macro
494 ///
495 /// ### Example 8 script Background_compton.C:
496 ///
497 /// Sometimes it is necessary to include also the Compton edges into the estimate of
498 /// the background. This example presents the synthetic spectrum
499 /// with Compton edges. The background was estimated using the 8-th order filter
500 /// with the estimation of the Compton edges using decreasing
501 /// clipping window algorithm (numberIterations=10) with smoothing
502 /// (smoothingWindow=5).
503 ///
504 /// Example of the estimate of the background with Compton edges (red line) for
505 /// numberIterations=10, 8-th order difference filter, using decreasing clipping
506 /// window algorithm and smoothing (smoothingWindow=5).
507 ///
508 /// Begin_Macro(source)
509 /// ../../../tutorials/spectrum/Background_compton.C
510 /// End_Macro
511 
512 const char *TSpectrum::Background(Double_t *spectrum, int ssize,
513  int numberIterations,
514  int direction, int filterOrder,
515  bool smoothing,int smoothWindow,
516  bool compton)
517 {
518  int i, j, w, bw, b1, b2, priz;
519  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;
520  if (ssize <= 0)
521  return "Wrong Parameters";
522  if (numberIterations < 1)
523  return "Width of Clipping Window Must Be Positive";
524  if (ssize < 2 * numberIterations + 1)
525  return "Too Large Clipping Window";
526  if (smoothing == kTRUE && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15)
527  return "Incorrect width of smoothing window";
528  Double_t *working_space = new Double_t[2 * ssize];
529  for (i = 0; i < ssize; i++){
530  working_space[i] = spectrum[i];
531  working_space[i + ssize] = spectrum[i];
532  }
533  bw=(smoothWindow-1)/2;
534  if (direction == kBackIncreasingWindow)
535  i = 1;
536  else if(direction == kBackDecreasingWindow)
537  i = numberIterations;
538  if (filterOrder == kBackOrder2) {
539  do{
540  for (j = i; j < ssize - i; j++) {
541  if (smoothing == kFALSE){
542  a = working_space[ssize + j];
543  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
544  if (b < a)
545  a = b;
546  working_space[j] = a;
547  }
548 
549  else if (smoothing == kTRUE){
550  a = working_space[ssize + j];
551  av = 0;
552  men = 0;
553  for (w = j - bw; w <= j + bw; w++){
554  if ( w >= 0 && w < ssize){
555  av += working_space[ssize + w];
556  men +=1;
557  }
558  }
559  av = av / men;
560  b = 0;
561  men = 0;
562  for (w = j - i - bw; w <= j - i + bw; w++){
563  if ( w >= 0 && w < ssize){
564  b += working_space[ssize + w];
565  men +=1;
566  }
567  }
568  b = b / men;
569  c = 0;
570  men = 0;
571  for (w = j + i - bw; w <= j + i + bw; w++){
572  if ( w >= 0 && w < ssize){
573  c += working_space[ssize + w];
574  men +=1;
575  }
576  }
577  c = c / men;
578  b = (b + c) / 2;
579  if (b < a)
580  av = b;
581  working_space[j]=av;
582  }
583  }
584  for (j = i; j < ssize - i; j++)
585  working_space[ssize + j] = working_space[j];
586  if (direction == kBackIncreasingWindow)
587  i+=1;
588  else if(direction == kBackDecreasingWindow)
589  i-=1;
590  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
591  }
592 
593  else if (filterOrder == kBackOrder4) {
594  do{
595  for (j = i; j < ssize - i; j++) {
596  if (smoothing == kFALSE){
597  a = working_space[ssize + j];
598  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
599  c = 0;
600  ai = i / 2;
601  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
602  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
603  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
604  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
605  if (b < c)
606  b = c;
607  if (b < a)
608  a = b;
609  working_space[j] = a;
610  }
611 
612  else if (smoothing == kTRUE){
613  a = working_space[ssize + j];
614  av = 0;
615  men = 0;
616  for (w = j - bw; w <= j + bw; w++){
617  if ( w >= 0 && w < ssize){
618  av += working_space[ssize + w];
619  men +=1;
620  }
621  }
622  av = av / men;
623  b = 0;
624  men = 0;
625  for (w = j - i - bw; w <= j - i + bw; w++){
626  if ( w >= 0 && w < ssize){
627  b += working_space[ssize + w];
628  men +=1;
629  }
630  }
631  b = b / men;
632  c = 0;
633  men = 0;
634  for (w = j + i - bw; w <= j + i + bw; w++){
635  if ( w >= 0 && w < ssize){
636  c += working_space[ssize + w];
637  men +=1;
638  }
639  }
640  c = c / men;
641  b = (b + c) / 2;
642  ai = i / 2;
643  b4 = 0, men = 0;
644  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
645  if (w >= 0 && w < ssize){
646  b4 += working_space[ssize + w];
647  men +=1;
648  }
649  }
650  b4 = b4 / men;
651  c4 = 0, men = 0;
652  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
653  if (w >= 0 && w < ssize){
654  c4 += working_space[ssize + w];
655  men +=1;
656  }
657  }
658  c4 = c4 / men;
659  d4 = 0, men = 0;
660  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
661  if (w >= 0 && w < ssize){
662  d4 += working_space[ssize + w];
663  men +=1;
664  }
665  }
666  d4 = d4 / men;
667  e4 = 0, men = 0;
668  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
669  if (w >= 0 && w < ssize){
670  e4 += working_space[ssize + w];
671  men +=1;
672  }
673  }
674  e4 = e4 / men;
675  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
676  if (b < b4)
677  b = b4;
678  if (b < a)
679  av = b;
680  working_space[j]=av;
681  }
682  }
683  for (j = i; j < ssize - i; j++)
684  working_space[ssize + j] = working_space[j];
685  if (direction == kBackIncreasingWindow)
686  i+=1;
687  else if(direction == kBackDecreasingWindow)
688  i-=1;
689  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
690  }
691 
692  else if (filterOrder == kBackOrder6) {
693  do{
694  for (j = i; j < ssize - i; j++) {
695  if (smoothing == kFALSE){
696  a = working_space[ssize + j];
697  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
698  c = 0;
699  ai = i / 2;
700  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
701  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
702  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
703  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
704  d = 0;
705  ai = i / 3;
706  d += working_space[ssize + j - (Int_t) (3 * ai)] / 20;
707  d -= 6 * working_space[ssize + j - (Int_t) (2 * ai)] / 20;
708  d += 15 * working_space[ssize + j - (Int_t) ai] / 20;
709  d += 15 * working_space[ssize + j + (Int_t) ai] / 20;
710  d -= 6 * working_space[ssize + j + (Int_t) (2 * ai)] / 20;
711  d += working_space[ssize + j + (Int_t) (3 * ai)] / 20;
712  if (b < d)
713  b = d;
714  if (b < c)
715  b = c;
716  if (b < a)
717  a = b;
718  working_space[j] = a;
719  }
720 
721  else if (smoothing == kTRUE){
722  a = working_space[ssize + j];
723  av = 0;
724  men = 0;
725  for (w = j - bw; w <= j + bw; w++){
726  if ( w >= 0 && w < ssize){
727  av += working_space[ssize + w];
728  men +=1;
729  }
730  }
731  av = av / men;
732  b = 0;
733  men = 0;
734  for (w = j - i - bw; w <= j - i + bw; w++){
735  if ( w >= 0 && w < ssize){
736  b += working_space[ssize + w];
737  men +=1;
738  }
739  }
740  b = b / men;
741  c = 0;
742  men = 0;
743  for (w = j + i - bw; w <= j + i + bw; w++){
744  if ( w >= 0 && w < ssize){
745  c += working_space[ssize + w];
746  men +=1;
747  }
748  }
749  c = c / men;
750  b = (b + c) / 2;
751  ai = i / 2;
752  b4 = 0, men = 0;
753  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
754  if (w >= 0 && w < ssize){
755  b4 += working_space[ssize + w];
756  men +=1;
757  }
758  }
759  b4 = b4 / men;
760  c4 = 0, men = 0;
761  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
762  if (w >= 0 && w < ssize){
763  c4 += working_space[ssize + w];
764  men +=1;
765  }
766  }
767  c4 = c4 / men;
768  d4 = 0, men = 0;
769  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
770  if (w >= 0 && w < ssize){
771  d4 += working_space[ssize + w];
772  men +=1;
773  }
774  }
775  d4 = d4 / men;
776  e4 = 0, men = 0;
777  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
778  if (w >= 0 && w < ssize){
779  e4 += working_space[ssize + w];
780  men +=1;
781  }
782  }
783  e4 = e4 / men;
784  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
785  ai = i / 3;
786  b6 = 0, men = 0;
787  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
788  if (w >= 0 && w < ssize){
789  b6 += working_space[ssize + w];
790  men +=1;
791  }
792  }
793  b6 = b6 / men;
794  c6 = 0, men = 0;
795  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
796  if (w >= 0 && w < ssize){
797  c6 += working_space[ssize + w];
798  men +=1;
799  }
800  }
801  c6 = c6 / men;
802  d6 = 0, men = 0;
803  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
804  if (w >= 0 && w < ssize){
805  d6 += working_space[ssize + w];
806  men +=1;
807  }
808  }
809  d6 = d6 / men;
810  e6 = 0, men = 0;
811  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
812  if (w >= 0 && w < ssize){
813  e6 += working_space[ssize + w];
814  men +=1;
815  }
816  }
817  e6 = e6 / men;
818  f6 = 0, men = 0;
819  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
820  if (w >= 0 && w < ssize){
821  f6 += working_space[ssize + w];
822  men +=1;
823  }
824  }
825  f6 = f6 / men;
826  g6 = 0, men = 0;
827  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
828  if (w >= 0 && w < ssize){
829  g6 += working_space[ssize + w];
830  men +=1;
831  }
832  }
833  g6 = g6 / men;
834  b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
835  if (b < b6)
836  b = b6;
837  if (b < b4)
838  b = b4;
839  if (b < a)
840  av = b;
841  working_space[j]=av;
842  }
843  }
844  for (j = i; j < ssize - i; j++)
845  working_space[ssize + j] = working_space[j];
846  if (direction == kBackIncreasingWindow)
847  i+=1;
848  else if(direction == kBackDecreasingWindow)
849  i-=1;
850  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
851  }
852 
853  else if (filterOrder == kBackOrder8) {
854  do{
855  for (j = i; j < ssize - i; j++) {
856  if (smoothing == kFALSE){
857  a = working_space[ssize + j];
858  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
859  c = 0;
860  ai = i / 2;
861  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
862  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
863  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
864  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
865  d = 0;
866  ai = i / 3;
867  d += working_space[ssize + j - (Int_t) (3 * ai)] / 20;
868  d -= 6 * working_space[ssize + j - (Int_t) (2 * ai)] / 20;
869  d += 15 * working_space[ssize + j - (Int_t) ai] / 20;
870  d += 15 * working_space[ssize + j + (Int_t) ai] / 20;
871  d -= 6 * working_space[ssize + j + (Int_t) (2 * ai)] / 20;
872  d += working_space[ssize + j + (Int_t) (3 * ai)] / 20;
873  e = 0;
874  ai = i / 4;
875  e -= working_space[ssize + j - (Int_t) (4 * ai)] / 70;
876  e += 8 * working_space[ssize + j - (Int_t) (3 * ai)] / 70;
877  e -= 28 * working_space[ssize + j - (Int_t) (2 * ai)] / 70;
878  e += 56 * working_space[ssize + j - (Int_t) ai] / 70;
879  e += 56 * working_space[ssize + j + (Int_t) ai] / 70;
880  e -= 28 * working_space[ssize + j + (Int_t) (2 * ai)] / 70;
881  e += 8 * working_space[ssize + j + (Int_t) (3 * ai)] / 70;
882  e -= working_space[ssize + j + (Int_t) (4 * ai)] / 70;
883  if (b < e)
884  b = e;
885  if (b < d)
886  b = d;
887  if (b < c)
888  b = c;
889  if (b < a)
890  a = b;
891  working_space[j] = a;
892  }
893 
894  else if (smoothing == kTRUE){
895  a = working_space[ssize + j];
896  av = 0;
897  men = 0;
898  for (w = j - bw; w <= j + bw; w++){
899  if ( w >= 0 && w < ssize){
900  av += working_space[ssize + w];
901  men +=1;
902  }
903  }
904  av = av / men;
905  b = 0;
906  men = 0;
907  for (w = j - i - bw; w <= j - i + bw; w++){
908  if ( w >= 0 && w < ssize){
909  b += working_space[ssize + w];
910  men +=1;
911  }
912  }
913  b = b / men;
914  c = 0;
915  men = 0;
916  for (w = j + i - bw; w <= j + i + bw; w++){
917  if ( w >= 0 && w < ssize){
918  c += working_space[ssize + w];
919  men +=1;
920  }
921  }
922  c = c / men;
923  b = (b + c) / 2;
924  ai = i / 2;
925  b4 = 0, men = 0;
926  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
927  if (w >= 0 && w < ssize){
928  b4 += working_space[ssize + w];
929  men +=1;
930  }
931  }
932  b4 = b4 / men;
933  c4 = 0, men = 0;
934  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
935  if (w >= 0 && w < ssize){
936  c4 += working_space[ssize + w];
937  men +=1;
938  }
939  }
940  c4 = c4 / men;
941  d4 = 0, men = 0;
942  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
943  if (w >= 0 && w < ssize){
944  d4 += working_space[ssize + w];
945  men +=1;
946  }
947  }
948  d4 = d4 / men;
949  e4 = 0, men = 0;
950  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
951  if (w >= 0 && w < ssize){
952  e4 += working_space[ssize + w];
953  men +=1;
954  }
955  }
956  e4 = e4 / men;
957  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
958  ai = i / 3;
959  b6 = 0, men = 0;
960  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
961  if (w >= 0 && w < ssize){
962  b6 += working_space[ssize + w];
963  men +=1;
964  }
965  }
966  b6 = b6 / men;
967  c6 = 0, men = 0;
968  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
969  if (w >= 0 && w < ssize){
970  c6 += working_space[ssize + w];
971  men +=1;
972  }
973  }
974  c6 = c6 / men;
975  d6 = 0, men = 0;
976  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
977  if (w >= 0 && w < ssize){
978  d6 += working_space[ssize + w];
979  men +=1;
980  }
981  }
982  d6 = d6 / men;
983  e6 = 0, men = 0;
984  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
985  if (w >= 0 && w < ssize){
986  e6 += working_space[ssize + w];
987  men +=1;
988  }
989  }
990  e6 = e6 / men;
991  f6 = 0, men = 0;
992  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
993  if (w >= 0 && w < ssize){
994  f6 += working_space[ssize + w];
995  men +=1;
996  }
997  }
998  f6 = f6 / men;
999  g6 = 0, men = 0;
1000  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1001  if (w >= 0 && w < ssize){
1002  g6 += working_space[ssize + w];
1003  men +=1;
1004  }
1005  }
1006  g6 = g6 / men;
1007  b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1008  ai = i / 4;
1009  b8 = 0, men = 0;
1010  for (w = j - (Int_t)(4 * ai) - bw; w <= j - (Int_t)(4 * ai) + bw; w++){
1011  if (w >= 0 && w < ssize){
1012  b8 += working_space[ssize + w];
1013  men +=1;
1014  }
1015  }
1016  b8 = b8 / men;
1017  c8 = 0, men = 0;
1018  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
1019  if (w >= 0 && w < ssize){
1020  c8 += working_space[ssize + w];
1021  men +=1;
1022  }
1023  }
1024  c8 = c8 / men;
1025  d8 = 0, men = 0;
1026  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1027  if (w >= 0 && w < ssize){
1028  d8 += working_space[ssize + w];
1029  men +=1;
1030  }
1031  }
1032  d8 = d8 / men;
1033  e8 = 0, men = 0;
1034  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1035  if (w >= 0 && w < ssize){
1036  e8 += working_space[ssize + w];
1037  men +=1;
1038  }
1039  }
1040  e8 = e8 / men;
1041  f8 = 0, men = 0;
1042  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1043  if (w >= 0 && w < ssize){
1044  f8 += working_space[ssize + w];
1045  men +=1;
1046  }
1047  }
1048  f8 = f8 / men;
1049  g8 = 0, men = 0;
1050  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1051  if (w >= 0 && w < ssize){
1052  g8 += working_space[ssize + w];
1053  men +=1;
1054  }
1055  }
1056  g8 = g8 / men;
1057  h8 = 0, men = 0;
1058  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1059  if (w >= 0 && w < ssize){
1060  h8 += working_space[ssize + w];
1061  men +=1;
1062  }
1063  }
1064  h8 = h8 / men;
1065  i8 = 0, men = 0;
1066  for (w = j + (Int_t)(4 * ai) - bw; w <= j + (Int_t)(4 * ai) + bw; w++){
1067  if (w >= 0 && w < ssize){
1068  i8 += working_space[ssize + w];
1069  men +=1;
1070  }
1071  }
1072  i8 = i8 / men;
1073  b8 = ( -b8 + 8 * c8 - 28 * d8 + 56 * e8 - 56 * f8 - 28 * g8 + 8 * h8 - i8)/70;
1074  if (b < b8)
1075  b = b8;
1076  if (b < b6)
1077  b = b6;
1078  if (b < b4)
1079  b = b4;
1080  if (b < a)
1081  av = b;
1082  working_space[j]=av;
1083  }
1084  }
1085  for (j = i; j < ssize - i; j++)
1086  working_space[ssize + j] = working_space[j];
1087  if (direction == kBackIncreasingWindow)
1088  i += 1;
1089  else if(direction == kBackDecreasingWindow)
1090  i -= 1;
1091  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1092  }
1093 
1094  if (compton == kTRUE) {
1095  for (i = 0, b2 = 0; i < ssize; i++){
1096  b1 = b2;
1097  a = working_space[i], b = spectrum[i];
1098  j = i;
1099  if (TMath::Abs(a - b) >= 1) {
1100  b1 = i - 1;
1101  if (b1 < 0)
1102  b1 = 0;
1103  yb1 = working_space[b1];
1104  for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){
1105  a = working_space[b2], b = spectrum[b2];
1106  c = c + b - yb1;
1107  if (TMath::Abs(a - b) < 1) {
1108  priz = 1;
1109  yb2 = b;
1110  }
1111  }
1112  if (b2 == ssize)
1113  b2 -= 1;
1114  yb2 = working_space[b2];
1115  if (yb1 <= yb2){
1116  for (j = b1, c = 0; j <= b2; j++){
1117  b = spectrum[j];
1118  c = c + b - yb1;
1119  }
1120  if (c > 1){
1121  c = (yb2 - yb1) / c;
1122  for (j = b1, d = 0; j <= b2 && j < ssize; j++){
1123  b = spectrum[j];
1124  d = d + b - yb1;
1125  a = c * d + yb1;
1126  working_space[ssize + j] = a;
1127  }
1128  }
1129  }
1130 
1131  else{
1132  for (j = b2, c = 0; j >= b1; j--){
1133  b = spectrum[j];
1134  c = c + b - yb2;
1135  }
1136  if (c > 1){
1137  c = (yb1 - yb2) / c;
1138  for (j = b2, d = 0;j >= b1 && j >= 0; j--){
1139  b = spectrum[j];
1140  d = d + b - yb2;
1141  a = c * d + yb2;
1142  working_space[ssize + j] = a;
1143  }
1144  }
1145  }
1146  i=b2;
1147  }
1148  }
1149  }
1150 
1151  for (j = 0; j < ssize; j++)
1152  spectrum[j] = working_space[ssize + j];
1153  delete[]working_space;
1154  return 0;
1155 }
1156 
1157 ////////////////////////////////////////////////////////////////////////////////
1158 /// One-dimensional markov spectrum smoothing function
1159 ///
1160 /// This function calculates smoothed spectrum from source spectrum based on
1161 /// Markov chain method. The result is placed in the array pointed by source
1162 /// pointer. On successful completion it returns 0. On error it returns pointer
1163 /// to the string describing error.
1164 ///
1165 /// #### Parameters:
1166 ///
1167 /// - source: pointer to the array of source spectrum
1168 /// - ssize: length of source array
1169 /// - averWindow: width of averaging smoothing window
1170 ///
1171 /// The goal of this function is the suppression of the statistical fluctuations.
1172 /// The algorithm is based on discrete Markov chain, which has very simple
1173 /// invariant distribution:
1174 ///
1175 /// \f[
1176 /// U_2 = \frac{p_{1,2}}{p_{2,1}}U_1, U_3 = \frac{p_{2,3}}{p_{3,2}}U_2U_1, ... , U_n = \frac{p_{n-1,n}}{p_{n,n-1}}U_{n-1}...U_2U_1
1177 /// \f]
1178 /// \f$ U_1\f$ being defined from the normalization condition
1179 /// \f$ \sum_{i=1}^{n} U_i=1\f$. \f$ n \f$ is the length of the smoothed spectrum and
1180 /// \f[
1181 /// p_{i,i\pm 1} = A_i\sum_{k=1}^{m} exp\left[ \frac{y(i\pm k)-y(i)}{y(i\pm k)+y(i)}\right]
1182 /// \f]
1183 ///
1184 /// #### Reference:
1185 ///
1186 /// 1. Z.K. Silagadze, A new algorithm for automatic photopeak searches.
1187 /// NIM A 376 (1996), 451.
1188 ///
1189 /// ### Example 14 - script Smoothing.C
1190 ///
1191 /// Begin_Macro(source)
1192 /// ../../../tutorials/spectrum/Smoothing.C
1193 /// End_Macro
1194 
1195 const char* TSpectrum::SmoothMarkov(Double_t *source, int ssize, int averWindow)
1196 {
1197  int xmin, xmax, i, l;
1198  Double_t a, b, maxch;
1199  Double_t nom, nip, nim, sp, sm, area = 0;
1200  if(averWindow <= 0)
1201  return "Averaging Window must be positive";
1202  Double_t *working_space = new Double_t[ssize];
1203  xmin = 0,xmax = ssize - 1;
1204  for(i = 0, maxch = 0; i < ssize; i++){
1205  working_space[i]=0;
1206  if(maxch < source[i])
1207  maxch = source[i];
1208 
1209  area += source[i];
1210  }
1211  if(maxch == 0) {
1212  delete [] working_space;
1213  return 0 ;
1214  }
1215 
1216  nom = 1;
1217  working_space[xmin] = 1;
1218  for(i = xmin; i < xmax; i++){
1219  nip = source[i] / maxch;
1220  nim = source[i + 1] / maxch;
1221  sp = 0,sm = 0;
1222  for(l = 1; l <= averWindow; l++){
1223  if((i + l) > xmax)
1224  a = source[xmax] / maxch;
1225 
1226  else
1227  a = source[i + l] / maxch;
1228  b = a - nip;
1229  if(a + nip <= 0)
1230  a = 1;
1231 
1232  else
1233  a = TMath::Sqrt(a + nip);
1234  b = b / a;
1235  b = TMath::Exp(b);
1236  sp = sp + b;
1237  if((i - l + 1) < xmin)
1238  a = source[xmin] / maxch;
1239 
1240  else
1241  a = source[i - l + 1] / maxch;
1242  b = a - nim;
1243  if(a + nim <= 0)
1244  a = 1;
1245  else
1246  a = TMath::Sqrt(a + nim);
1247  b = b / a;
1248  b = TMath::Exp(b);
1249  sm = sm + b;
1250  }
1251  a = sp / sm;
1252  a = working_space[i + 1] = working_space[i] * a;
1253  nom = nom + a;
1254  }
1255  for(i = xmin; i <= xmax; i++){
1256  working_space[i] = working_space[i] / nom;
1257  }
1258  for(i = 0; i < ssize; i++)
1259  source[i] = working_space[i] * area;
1260  delete[]working_space;
1261  return 0;
1262 }
1263 
1264 ////////////////////////////////////////////////////////////////////////////////
1265 /// One-dimensional deconvolution function
1266 ///
1267 /// This function calculates deconvolution from source spectrum according to
1268 /// response spectrum using Gold deconvolution algorithm. The result is placed
1269 /// in the vector pointed by source pointer. On successful completion it
1270 /// returns 0. On error it returns pointer to the string describing error. If
1271 /// desired after every numberIterations one can apply boosting operation
1272 /// (exponential function with exponent given by boost coefficient) and repeat
1273 /// it numberRepetitions times.
1274 ///
1275 /// #### Parameters:
1276 ///
1277 /// - source: pointer to the vector of source spectrum
1278 /// - response: pointer to the vector of response spectrum
1279 /// - ssize: length of source and response spectra
1280 /// - numberIterations, for details we refer to the reference given below
1281 /// - numberRepetitions, for repeated boosted deconvolution
1282 /// - boost, boosting coefficient
1283 ///
1284 /// The goal of this function is the improvement of the resolution in spectra,
1285 /// decomposition of multiplets. The mathematical formulation of
1286 /// the convolution system is:
1287 ///
1288 /// \f[
1289 /// y(i) = \sum_{k=0}^{N-1} h(i-k)x(k), i=0,1,2,...,N-1
1290 /// \f]
1291 ///
1292 /// where h(i) is the impulse response function, x, y are input and output
1293 /// vectors, respectively, N is the length of x and h vectors. In matrix form
1294 /// we have:
1295 /**
1296  \f[
1297  \begin{bmatrix}
1298  y(0) \\
1299  y(1) \\
1300  \dots \\
1301  y(2N-2) \\
1302  y(2N-1)
1303  \end{bmatrix}
1304  =
1305  \begin{bmatrix}
1306  h(0) & 0 & 0 & \dots & 0 \\
1307  h(1) & h(0) & 0 & \dots & \dots \\
1308  \dots & h(1) & h(0) & \dots & \dots \\
1309  \dots & \dots & h(1) & \dots & \dots \\
1310  \dots & \dots & \dots & \dots & \dots \\
1311  h(N-1) & \dots & \dots &\dots & 0 \\
1312  0 & h(N-1) & \dots & \dots & h(0) \\
1313  0 & 0 & h(N-1) & \dots & h(1) \\
1314  \dots & \dots & \dots & \dots & \dots \\
1315  0 & 0 & 0 & \dots & h(N-1)
1316  \end{bmatrix}
1317  \begin{bmatrix}
1318  x(0) \\
1319  x(1) \\
1320  \dots \\
1321  x(N-1)
1322  \end{bmatrix}
1323  \f]
1324 */
1325 /// Let us assume that we know the response and the output vector (spectrum) of
1326 /// the above given system. The deconvolution represents solution of the
1327 /// overdetermined system of linear equations, i.e., the calculation of the
1328 /// vector x. From numerical stability point of view the operation of
1329 /// deconvolution is extremely critical (ill-posed problem) as well as time
1330 /// consuming operation. The Gold deconvolution algorithm proves to work very
1331 /// well, other methods (Fourier, VanCittert etc) oscillate. It is suitable to
1332 /// process positive definite data (e.g. histograms).
1333 ///
1334 /// #### Gold deconvolution algorithm:
1335 /**
1336  \f[
1337  y = Hx \\
1338  H^{T}y = H^{T}Hx \\
1339  y^{'} = H^{'}x \\
1340  x_{i}^{(k+1)} = \frac{y_{i}^{'}}{\sum_{m=0}^{N-1}H_{im}^{'}x_{m}{(k)}}x_{i}{(k)}, i=0,1,2,...,N-1 \\
1341  k = 1,2,3,...,L
1342  x^{0} = [1,1, ..., 1]^T
1343  \f]
1344 */
1345 /// Where L is given number of iterations (numberIterations parameter).
1346 ///
1347 /// #### Boosted deconvolution:
1348 ///
1349 /// 1. Set the initial solution:
1350 /// \f$ x^{(0)} = [1,1,...,1]^{T} \f$
1351 /// 2. Set required number of repetitions R and iterations L.
1352 /// 3. Set r = 1.
1353 /// 4. Using Gold deconvolution algorithm for k=1,2,...,L find
1354 /// \f$ x^{(L)} \f$
1355 /// 5. If r = R stop calculation, else
1356 ///
1357 /// 1. Apply boosting operation, i.e., set
1358 /// \f$ x^{(0)}(i) = [x^{(L)}(i)]^{p} \f$
1359 /// i=0,1,...N-1 and p is boosting coefficient >0.
1360 /// 2. r = r + 1
1361 /// 3. continue in 4.
1362 ///
1363 /// #### References:
1364 ///
1365 /// 1. Gold R., ANL-6984, Argonne National Laboratories, Argonne Ill, 1964.
1366 /// 2. Coote G.E., Iterative smoothing and deconvolution of one- and two-dimensional
1367 /// elemental distribution data, NIM B 130 (1997) 118.
1368 /// 3. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky,
1369 /// I. Turzo: Efficient one- and two-dimensional Gold deconvolution and
1370 /// its application to gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
1371 /// 4. Morhac; M., Matouoek V., Kliman J., Efficient algorithm of multidimensional
1372 /// deconvolution and its application to nuclear data processing, Digital Signal
1373 /// Processing 13 (2003) 144.
1374 ///
1375 /// ### Example 8 - script Deconvolution.C :
1376 ///
1377 /// response function (usually peak) should be shifted left to the first
1378 /// non-zero channel (bin).
1379 ///
1380 /// \image html TSpectrum_Deconvolution2.jpg Principle how the response matrix is composed inside of the Deconvolution function.
1381 ///
1382 /// Begin_Macro(source)
1383 /// ../../../tutorials/spectrum/Deconvolution.C
1384 /// End_Macro
1385 ///
1386 /// ### Examples of Gold deconvolution method:
1387 ///
1388 /// First let us study the influence of the number of iterations on the
1389 /// deconvolved spectrum (Fig. 12).
1390 ///
1391 /// \image html TSpectrum_Deconvolution_wide1.jpg Fig. 12 Study of Gold deconvolution algorithm.The original source spectrum is drawn with black color, spectrum after 100 iterations with red color, spectrum after 1000 iterations with blue color, spectrum after 10000 iterations with green color and spectrum after 100000 iterations with magenta color.
1392 ///
1393 /// For relatively narrow peaks in the above given example the Gold
1394 /// deconvolution method is able to decompose overlapping peaks practically to
1395 /// delta - functions. In the next example we have chosen a synthetic data
1396 /// (spectrum, 256 channels) consisting of 5 very closely positioned, relatively
1397 /// wide peaks (sigma =5), with added noise (Fig. 13). Thin lines represent
1398 /// pure Gaussians (see Table 1); thick line is a resulting spectrum with
1399 /// additive noise (10% of the amplitude of small peaks).
1400 ///
1401 /// \image html TSpectrum_Deconvolution_wide2.jpg Fig. 13 Testing example of synthetic spectrum composed of 5 Gaussians with added noise.
1402 ///
1403 /// | Peak # | Position | Height | Area |
1404 /// |----------|----------|--------|--------|
1405 /// | 1 | 50 | 500 | 10159 |
1406 /// | 2 | 70 | 3000 | 60957 |
1407 /// | 3 | 80 | 1000 | 20319 |
1408 /// | 4 | 100 | 5000 | 101596 |
1409 /// | 5 | 110 | 500 | 10159 |
1410 ///
1411 /// Table 1 Positions, heights and areas of peaks in the spectrum shown in Fig. 13.
1412 ///
1413 /// In ideal case, we should obtain the result given in Fig. 14. The areas of
1414 /// the Gaussian components of the spectrum are concentrated completely to
1415 /// delta-functions. When solving the overdetermined system of linear equations
1416 /// with data from Fig. 13 in the sense of minimum least squares criterion
1417 /// without any regularisation we obtain the result with large oscillations
1418 /// (Fig. 15). From mathematical point of view, it is the optimal solution in
1419 /// the unconstrained space of independent variables. From physical point of
1420 /// view we are interested only in a meaningful solution. Therefore, we have to
1421 /// employ regularisation techniques (e.g. Gold deconvolution) and/or to
1422 /// confine the space of allowed solutions to subspace of positive solutions.
1423 ///
1424 /// \image html TSpectrum_Deconvolution_wide3.jpg Fig. 14 The same spectrum like in Fig. 13, outlined bars show the contents of present components (peaks).
1425 /// \image html TSpectrum_Deconvolution_wide4.jpg Fig. 15 Least squares solution of the system of linear equations without regularisation.
1426 ///
1427 /// ### Example 9 - script Deconvolution_wide.C
1428 ///
1429 /// When we employ Gold deconvolution algorithm we obtain the result given in
1430 /// Fig. 16. One can observe that the resulting spectrum is smooth. On the
1431 /// other hand the method is not able to decompose completely the peaks in the
1432 /// spectrum.
1433 ///
1434 /// Example of Gold deconvolution for closely positioned wide peaks. The original
1435 /// source spectrum is drawn with black color, the spectrum after the deconvolution
1436 /// (10000 iterations) with red color.
1437 ///
1438 /// Begin_Macro(source)
1439 /// ../../../tutorials/spectrum/Deconvolution_wide.C
1440 /// End_Macro
1441 ///
1442 /// ### Example 10 - script Deconvolution_wide_boost.C :
1443 ///
1444 /// Further let us employ boosting operation into deconvolution (Fig. 17).
1445 ///
1446 /// The original source spectrum is drawn with black color, the spectrum after
1447 /// the deconvolution with red color. Number of iterations = 200, number of
1448 /// repetitions = 50 and boosting coefficient = 1.2.
1449 ///
1450 /// One can observe that peaks are decomposed practically to delta functions.
1451 /// Number of peaks is correct, positions of big peaks as well as their areas
1452 /// are relatively well estimated. However there is a considerable error in
1453 /// the estimation of the position of small right hand peak.
1454 ///
1455 /// Begin_Macro(source)
1456 /// ../../../tutorials/spectrum/Deconvolution_wide_boost.C
1457 /// End_Macro
1458 
1459 const char *TSpectrum::Deconvolution(Double_t *source, const Double_t *response,
1460  int ssize, int numberIterations,
1461  int numberRepetitions, Double_t boost )
1462 {
1463  if (ssize <= 0)
1464  return "Wrong Parameters";
1465 
1466  if (numberRepetitions <= 0)
1467  return "Wrong Parameters";
1468 
1469  // working_space-pointer to the working vector
1470  // (its size must be 4*ssize of source spectrum)
1471  Double_t *working_space = new Double_t[4 * ssize];
1472  int i, j, k, lindex, posit, lh_gold, l, repet;
1473  Double_t lda, ldb, ldc, area, maximum;
1474  area = 0;
1475  lh_gold = -1;
1476  posit = 0;
1477  maximum = 0;
1478 
1479 //read response vector
1480  for (i = 0; i < ssize; i++) {
1481  lda = response[i];
1482  if (lda != 0)
1483  lh_gold = i + 1;
1484  working_space[i] = lda;
1485  area += lda;
1486  if (lda > maximum) {
1487  maximum = lda;
1488  posit = i;
1489  }
1490  }
1491  if (lh_gold == -1) {
1492  delete [] working_space;
1493  return "ZERO RESPONSE VECTOR";
1494  }
1495 
1496 //read source vector
1497  for (i = 0; i < ssize; i++)
1498  working_space[2 * ssize + i] = source[i];
1499 
1500 // create matrix at*a and vector at*y
1501  for (i = 0; i < ssize; i++){
1502  lda = 0;
1503  for (j = 0; j < ssize; j++){
1504  ldb = working_space[j];
1505  k = i + j;
1506  if (k < ssize){
1507  ldc = working_space[k];
1508  lda = lda + ldb * ldc;
1509  }
1510  }
1511  working_space[ssize + i] = lda;
1512  lda = 0;
1513  for (k = 0; k < ssize; k++){
1514  l = k - i;
1515  if (l >= 0){
1516  ldb = working_space[l];
1517  ldc = working_space[2 * ssize + k];
1518  lda = lda + ldb * ldc;
1519  }
1520  }
1521  working_space[3 * ssize + i]=lda;
1522  }
1523 
1524 // move vector at*y
1525  for (i = 0; i < ssize; i++){
1526  working_space[2 * ssize + i] = working_space[3 * ssize + i];
1527  }
1528 
1529 //initialization of resulting vector
1530  for (i = 0; i < ssize; i++)
1531  working_space[i] = 1;
1532 
1533  //**START OF ITERATIONS**
1534  for (repet = 0; repet < numberRepetitions; repet++) {
1535  if (repet != 0) {
1536  for (i = 0; i < ssize; i++)
1537  working_space[i] = TMath::Power(working_space[i], boost);
1538  }
1539  for (lindex = 0; lindex < numberIterations; lindex++) {
1540  for (i = 0; i < ssize; i++) {
1541  if (working_space[2 * ssize + i] > 0.000001
1542  && working_space[i] > 0.000001) {
1543  lda = 0;
1544  for (j = 0; j < lh_gold; j++) {
1545  ldb = working_space[j + ssize];
1546  if (j != 0){
1547  k = i + j;
1548  ldc = 0;
1549  if (k < ssize)
1550  ldc = working_space[k];
1551  k = i - j;
1552  if (k >= 0)
1553  ldc += working_space[k];
1554  }
1555 
1556  else
1557  ldc = working_space[i];
1558  lda = lda + ldb * ldc;
1559  }
1560  ldb = working_space[2 * ssize + i];
1561  if (lda != 0)
1562  lda = ldb / lda;
1563 
1564  else
1565  lda = 0;
1566  ldb = working_space[i];
1567  lda = lda * ldb;
1568  working_space[3 * ssize + i] = lda;
1569  }
1570  }
1571  for (i = 0; i < ssize; i++)
1572  working_space[i] = working_space[3 * ssize + i];
1573  }
1574  }
1575 
1576 //shift resulting spectrum
1577  for (i = 0; i < ssize; i++) {
1578  lda = working_space[i];
1579  j = i + posit;
1580  j = j % ssize;
1581  working_space[ssize + j] = lda;
1582  }
1583 
1584 //write back resulting spectrum
1585  for (i = 0; i < ssize; i++)
1586  source[i] = area * working_space[ssize + i];
1587  delete[]working_space;
1588  return 0;
1589 }
1590 
1591 
1592 ////////////////////////////////////////////////////////////////////////////////
1593 /// One-dimensional deconvolution function.
1594 ///
1595 /// This function calculates deconvolution from source spectrum according to
1596 /// response spectrum using Richardson-Lucy deconvolution algorithm. The result
1597 /// is placed in the vector pointed by source pointer. On successful completion
1598 /// it returns 0. On error it returns pointer to the string describing error.
1599 /// If desired after every numberIterations one can apply boosting operation
1600 /// (exponential function with exponent given by boost coefficient) and repeat
1601 /// it numberRepetitions times (see Gold deconvolution).
1602 ///
1603 /// #### Parameters:
1604 ///
1605 /// - source: pointer to the vector of source spectrum
1606 /// - response: pointer to the vector of response spectrum
1607 /// - ssize: length of source and response spectra
1608 /// - numberIterations, for details we refer to the reference given above
1609 /// - numberRepetitions, for repeated boosted deconvolution
1610 /// - boost, boosting coefficient
1611 ///
1612 /// ### Richardson-Lucy deconvolution algorithm:
1613 ///
1614 /// For discrete systems it has the form:
1615 /**
1616  \f[
1617  x^{(n)}(i) = x^{(n-1)}(i) \sum_{j=0}^{N-1}h(i,j)\frac{y(j)}{\sum_{k=0}^{M-1}h(j,k)x^{(n-1)}(k)} \\
1618  i \in \left<0,M-1\right>
1619  \f]
1620 */
1621 /// for positive input data and response matrix this iterative method forces
1622 /// the deconvoluted spectra to be non-negative. The Richardson-Lucy
1623 /// iteration converges to the maximum likelihood solution for Poisson statistics
1624 /// in the data.
1625 ///
1626 /// #### References:
1627 ///
1628 /// 1. Abreu M.C. et al., A four-dimensional deconvolution method to correct NA38
1629 /// experimental data, NIM A 405 (1998) 139.
1630 /// 2. Lucy L.B., A.J. 79 (1974) 745.
1631 /// 3. Richardson W.H., J. Opt. Soc. Am. 62 (1972) 55.
1632 ///
1633 /// ### Examples of Richardson-Lucy deconvolution method:
1634 ///
1635 /// ### Example 11 - script DeconvolutionRL_wide.C :
1636 ///
1637 /// When we employ Richardson-Lucy deconvolution algorithm to our data from
1638 /// Fig. 13 we obtain the following result. One can observe improvements
1639 /// as compared to the result achieved by Gold deconvolution. Nevertheless it is
1640 /// unable to decompose the multiplet.
1641 ///
1642 /// Example of Richardson-Lucy deconvolution for closely positioned wide peaks.
1643 /// The original source spectrum is drawn with black color, the spectrum after
1644 /// the deconvolution (10000 iterations) with red color.
1645 ///
1646 /// Begin_Macro(source)
1647 /// ../../../tutorials/spectrum/DeconvolutionRL_wide.C
1648 /// End_Macro
1649 ///
1650 /// ### Example 12 - script DeconvolutionRL_wide_boost.C :
1651 ///
1652 /// Further let us employ boosting operation into deconvolution.
1653 ///
1654 /// The original source spectrum is drawn with black color, the spectrum after
1655 /// the deconvolution with red color. Number of iterations = 200, number of
1656 /// repetitions = 50 and boosting coefficient = 1.2.
1657 ///
1658 /// One can observe improvements in the estimation of peak positions as compared
1659 /// to the results achieved by Gold deconvolution.
1660 ///
1661 /// Begin_Macro(source)
1662 /// ../../../tutorials/spectrum/DeconvolutionRL_wide_boost.C
1663 /// End_Macro
1664 
1665 const char *TSpectrum::DeconvolutionRL(Double_t *source, const Double_t *response,
1666  int ssize, int numberIterations,
1667  int numberRepetitions, Double_t boost )
1668 {
1669  if (ssize <= 0)
1670  return "Wrong Parameters";
1671 
1672  if (numberRepetitions <= 0)
1673  return "Wrong Parameters";
1674 
1675  // working_space-pointer to the working vector
1676  // (its size must be 4*ssize of source spectrum)
1677  Double_t *working_space = new Double_t[4 * ssize];
1678  int i, j, k, lindex, posit, lh_gold, repet, kmin, kmax;
1679  Double_t lda, ldb, ldc, maximum;
1680  lh_gold = -1;
1681  posit = 0;
1682  maximum = 0;
1683 
1684 //read response vector
1685  for (i = 0; i < ssize; i++) {
1686  lda = response[i];
1687  if (lda != 0)
1688  lh_gold = i + 1;
1689  working_space[ssize + i] = lda;
1690  if (lda > maximum) {
1691  maximum = lda;
1692  posit = i;
1693  }
1694  }
1695  if (lh_gold == -1) {
1696  delete [] working_space;
1697  return "ZERO RESPONSE VECTOR";
1698  }
1699 
1700 //read source vector
1701  for (i = 0; i < ssize; i++)
1702  working_space[2 * ssize + i] = source[i];
1703 
1704 //initialization of resulting vector
1705  for (i = 0; i < ssize; i++){
1706  if (i <= ssize - lh_gold)
1707  working_space[i] = 1;
1708 
1709  else
1710  working_space[i] = 0;
1711 
1712  }
1713  //**START OF ITERATIONS**
1714  for (repet = 0; repet < numberRepetitions; repet++) {
1715  if (repet != 0) {
1716  for (i = 0; i < ssize; i++)
1717  working_space[i] = TMath::Power(working_space[i], boost);
1718  }
1719  for (lindex = 0; lindex < numberIterations; lindex++) {
1720  for (i = 0; i <= ssize - lh_gold; i++){
1721  lda = 0;
1722  if (working_space[i] > 0){//x[i]
1723  for (j = i; j < i + lh_gold; j++){
1724  ldb = working_space[2 * ssize + j];//y[j]
1725  if (j < ssize){
1726  if (ldb > 0){//y[j]
1727  kmax = j;
1728  if (kmax > lh_gold - 1)
1729  kmax = lh_gold - 1;
1730  kmin = j + lh_gold - ssize;
1731  if (kmin < 0)
1732  kmin = 0;
1733  ldc = 0;
1734  for (k = kmax; k >= kmin; k--){
1735  ldc += working_space[ssize + k] * working_space[j - k];//h[k]*x[j-k]
1736  }
1737  if (ldc > 0)
1738  ldb = ldb / ldc;
1739 
1740  else
1741  ldb = 0;
1742  }
1743  ldb = ldb * working_space[ssize + j - i];//y[j]*h[j-i]/suma(h[j][k]x[k])
1744  }
1745  lda += ldb;
1746  }
1747  lda = lda * working_space[i];
1748  }
1749  working_space[3 * ssize + i] = lda;
1750  }
1751  for (i = 0; i < ssize; i++)
1752  working_space[i] = working_space[3 * ssize + i];
1753  }
1754  }
1755 
1756 //shift resulting spectrum
1757  for (i = 0; i < ssize; i++) {
1758  lda = working_space[i];
1759  j = i + posit;
1760  j = j % ssize;
1761  working_space[ssize + j] = lda;
1762  }
1763 
1764 //write back resulting spectrum
1765  for (i = 0; i < ssize; i++)
1766  source[i] = working_space[ssize + i];
1767  delete[]working_space;
1768  return 0;
1769 }
1770 
1771 
1772 ////////////////////////////////////////////////////////////////////////////////
1773 /// One-dimensional unfolding function
1774 ///
1775 /// This function unfolds source spectrum according to response matrix columns.
1776 /// The result is placed in the vector pointed by source pointer.
1777 /// The coefficients of the resulting vector represent contents of the columns
1778 /// (weights) in the input vector. On successful completion it returns 0. On
1779 /// error it returns pointer to the string describing error. If desired after
1780 /// every numberIterations one can apply boosting operation (exponential
1781 /// function with exponent given by boost coefficient) and repeat it
1782 /// numberRepetitions times. For details we refer to [1].
1783 ///
1784 /// #### Parameters:
1785 ///
1786 /// - source: pointer to the vector of source spectrum
1787 /// - respMatrix: pointer to the matrix of response spectra
1788 /// - ssizex: length of source spectrum and # of columns of the response
1789 /// matrix. ssizex must be >= ssizey.
1790 /// - ssizey: length of destination spectrum and # of rows of the response
1791 /// matrix.
1792 /// - numberIterations: number of iterations
1793 /// - numberRepetitions: number of repetitions for boosted deconvolution.
1794 /// It must be greater or equal to one.
1795 /// - boost: boosting coefficient, applies only if numberRepetitions is
1796 /// greater than one.
1797 ///
1798 /// ### Unfolding:
1799 ///
1800 /// The goal is the decomposition of spectrum to a given set of component spectra.
1801 ///
1802 /// The mathematical formulation of the discrete linear system is:
1803 ///
1804 /// \f[
1805 /// y(i) = \sum_{k=0}^{N_y-1} h(i,k)x(k), i = 0,1,2,...,N_x-1
1806 /// \f]
1807 /**
1808  \f[
1809  \begin{bmatrix}
1810  y(0) \\
1811  y(1) \\
1812  \dots \\
1813  y(N_x-1)
1814  \end{bmatrix}
1815  =
1816  \begin{bmatrix}
1817  h(0,0) & h(0,1) & \dots & h(0,N_y-1) \\
1818  h(1,0) & h(1,1) & \dots & h(1,N_y-1) \\
1819  \dots \\
1820  h(N_x-1,0) & h(N_x-1,1) & \dots & h(N_x-1,N_y-1)
1821  \end{bmatrix}
1822  \begin{bmatrix}
1823  x(0) \\
1824  x(1) \\
1825  \dots \\
1826  x(N_x-1)
1827  \end{bmatrix}
1828  \f]
1829 */
1830 /// #### References:
1831 ///
1832 /// 1. Jandel M., Morhac; M., Kliman J., Krupa L., Matouoek
1833 /// V., Hamilton J. H., Ramaya A. V.:
1834 /// Decomposition of continuum gamma-ray spectra using synthesised response matrix.
1835 /// NIM A 516 (2004), 172-183.
1836 ///
1837 /// ### Example of unfolding:
1838 ///
1839 /// ### Example 13 - script Unfolding.C:
1840 ///
1841 /// \image html TSpectrum_Unfolding3.gif Fig. 20 Response matrix composed of neutron spectra of pure chemical elements.
1842 /// \image html TSpectrum_Unfolding2.jpg Fig. 21 Source neutron spectrum to be decomposed
1843 /// \image html TSpectrum_Unfolding3.jpg Fig. 22 Spectrum after decomposition, contains 10 coefficients, which correspond to contents of chemical components (dominant 8-th and 10-th components, i.e. O, Si)
1844 ///
1845 /// #### Script:
1846 ///
1847 /// ~~~ {.cpp}
1848 /// // Example to illustrate unfolding function (class TSpectrum).
1849 /// // To execute this example, do
1850 /// // root > .x Unfolding.C
1851 ///
1852 /// #include <TSpectrum>
1853 ///
1854 /// void Unfolding() {
1855 /// Int_t i, j;
1856 /// Int_t nbinsx = 2048;
1857 /// Int_t nbinsy = 10;
1858 /// Double_t xmin = 0;
1859 /// Double_t xmax = nbinsx;
1860 /// Double_t ymin = 0;
1861 /// Double_t ymax = nbinsy;
1862 /// Double_t * source = new Double_t[nbinsx];
1863 /// Double_t ** response = new Double_t *[nbinsy];
1864 /// for (i=0;i<nbinsy;i++) response[i]=new Double_t[nbinsx];
1865 /// TH1F *h = new TH1F("h","",nbinsx,xmin,xmax);
1866 /// TH1F *d = new TH1F("d","Decomposition - unfolding",nbinsx,xmin,xmax);
1867 /// TH2F *decon_unf_resp = new TH2F("decon_unf_resp","Root File",nbinsy,ymin,ymax,nbinsx,xmin,xmax);
1868 /// TFile *f = new TFile("spectra/TSpectrum.root");
1869 /// h=(TH1F*) f->Get("decon_unf_in;1");
1870 /// TFile *fr = new TFile("spectra/TSpectrum.root");
1871 /// decon_unf_resp = (TH2F*) fr->Get("decon_unf_resp;1");
1872 /// for (i = 0; i < nbinsx; i++) source[i] = h->GetBinContent(i + 1);
1873 /// for (i = 0; i < nbinsy; i++){
1874 /// for (j = 0; j< nbinsx; j++){
1875 /// response[i][j] = decon_unf_resp->GetBinContent(i + 1, j + 1);
1876 /// }
1877 /// }
1878 /// TCanvas *Decon1 = gROOT->GetListOfCanvases()->FindObject("Decon1");
1879 /// if (!Decon1) Decon1 = new TCanvas("Decon1","Decon1",10,10,1000,700);
1880 /// h->Draw("L");
1881 /// TSpectrum *s = new TSpectrum();
1882 /// s->Unfolding(source,response,nbinsx,nbinsy,1000,1,1);
1883 /// for (i = 0; i < nbinsy; i++) d->SetBinContent(i + 1,source[i]);
1884 /// d->SetLineColor(kRed);
1885 /// d->SetAxisRange(0,nbinsy);
1886 /// d->Draw("");
1887 /// }
1888 /// ~~~
1889 
1890 const char *TSpectrum::Unfolding(Double_t *source,
1891  const Double_t **respMatrix,
1892  int ssizex, int ssizey,
1893  int numberIterations,
1894  int numberRepetitions, Double_t boost)
1895 {
1896  int i, j, k, lindex, lhx = 0, repet;
1897  Double_t lda, ldb, ldc, area;
1898  if (ssizex <= 0 || ssizey <= 0)
1899  return "Wrong Parameters";
1900  if (ssizex < ssizey)
1901  return "Sizex must be greater than sizey)";
1902  if (numberIterations <= 0)
1903  return "Number of iterations must be positive";
1904  Double_t *working_space =
1905  new Double_t[ssizex * ssizey + 2 * ssizey * ssizey + 4 * ssizex];
1906 
1907 /*read response matrix*/
1908  for (j = 0; j < ssizey && lhx != -1; j++) {
1909  area = 0;
1910  lhx = -1;
1911  for (i = 0; i < ssizex; i++) {
1912  lda = respMatrix[j][i];
1913  if (lda != 0) {
1914  lhx = i + 1;
1915  }
1916  working_space[j * ssizex + i] = lda;
1917  area = area + lda;
1918  }
1919  if (lhx != -1) {
1920  for (i = 0; i < ssizex; i++)
1921  working_space[j * ssizex + i] /= area;
1922  }
1923  }
1924  if (lhx == -1) {
1925  delete [] working_space;
1926  return ("ZERO COLUMN IN RESPONSE MATRIX");
1927  }
1928 
1929 /*read source vector*/
1930  for (i = 0; i < ssizex; i++)
1931  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
1932  source[i];
1933 
1934 /*create matrix at*a + at*y */
1935  for (i = 0; i < ssizey; i++) {
1936  for (j = 0; j < ssizey; j++) {
1937  lda = 0;
1938  for (k = 0; k < ssizex; k++) {
1939  ldb = working_space[ssizex * i + k];
1940  ldc = working_space[ssizex * j + k];
1941  lda = lda + ldb * ldc;
1942  }
1943  working_space[ssizex * ssizey + ssizey * i + j] = lda;
1944  }
1945  lda = 0;
1946  for (k = 0; k < ssizex; k++) {
1947  ldb = working_space[ssizex * i + k];
1948  ldc =
1949  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
1950  k];
1951  lda = lda + ldb * ldc;
1952  }
1953  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
1954  lda;
1955  }
1956 
1957 /*move vector at*y*/
1958  for (i = 0; i < ssizey; i++)
1959  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
1960  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
1961 
1962 /*create matrix at*a*at*a + vector at*a*at*y */
1963  for (i = 0; i < ssizey; i++) {
1964  for (j = 0; j < ssizey; j++) {
1965  lda = 0;
1966  for (k = 0; k < ssizey; k++) {
1967  ldb = working_space[ssizex * ssizey + ssizey * i + k];
1968  ldc = working_space[ssizex * ssizey + ssizey * j + k];
1969  lda = lda + ldb * ldc;
1970  }
1971  working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j] =
1972  lda;
1973  }
1974  lda = 0;
1975  for (k = 0; k < ssizey; k++) {
1976  ldb = working_space[ssizex * ssizey + ssizey * i + k];
1977  ldc =
1978  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
1979  k];
1980  lda = lda + ldb * ldc;
1981  }
1982  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
1983  lda;
1984  }
1985 
1986 /*move at*a*at*y*/
1987  for (i = 0; i < ssizey; i++)
1988  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
1989  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
1990 
1991 /*initialization in resulting vector */
1992  for (i = 0; i < ssizey; i++)
1993  working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = 1;
1994 
1995  /***START OF ITERATIONS***/
1996  for (repet = 0; repet < numberRepetitions; repet++) {
1997  if (repet != 0) {
1998  for (i = 0; i < ssizey; i++)
1999  working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = TMath::Power(working_space[ssizex * ssizey + 2 * ssizey * ssizey + i], boost);
2000  }
2001  for (lindex = 0; lindex < numberIterations; lindex++) {
2002  for (i = 0; i < ssizey; i++) {
2003  lda = 0;
2004  for (j = 0; j < ssizey; j++) {
2005  ldb =
2006  working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j];
2007  ldc = working_space[ssizex * ssizey + 2 * ssizey * ssizey + j];
2008  lda = lda + ldb * ldc;
2009  }
2010  ldb =
2011  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i];
2012  if (lda != 0) {
2013  lda = ldb / lda;
2014  }
2015 
2016  else
2017  lda = 0;
2018  ldb = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2019  lda = lda * ldb;
2020  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] = lda;
2021  }
2022  for (i = 0; i < ssizey; i++)
2023  working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
2024  working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2025  }
2026  }
2027 
2028 /*write back resulting spectrum*/
2029  for (i = 0; i < ssizex; i++) {
2030  if (i < ssizey)
2031  source[i] = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2032 
2033  else
2034  source[i] = 0;
2035  }
2036  delete[]working_space;
2037  return 0;
2038 }
2039 
2040 
2041 ////////////////////////////////////////////////////////////////////////////////
2042 /// One-dimensional high-resolution peak search function
2043 ///
2044 /// This function searches for peaks in source spectrum. It is based on
2045 /// deconvolution method. First the background is removed (if desired), then
2046 /// Markov smoothed spectrum is calculated (if desired), then the response
2047 /// function is generated according to given sigma and deconvolution is
2048 /// carried out. The order of peaks is arranged according to their heights in
2049 /// the spectrum after background elimination. The highest peak is the first in
2050 /// the list. On success it returns number of found peaks.
2051 ///
2052 /// #### Parameters:
2053 ///
2054 /// - source: pointer to the vector of source spectrum.
2055 /// - destVector: pointer to the vector of resulting deconvolved spectrum.
2056 /// - ssize: length of source spectrum.
2057 /// - sigma: sigma of searched peaks, for details we refer to manual.
2058 /// - threshold: threshold value in % for selected peaks, peaks with
2059 /// amplitude less than threshold*highest_peak/100
2060 /// are ignored, see manual.
2061 /// - backgroundRemove: logical variable, set if the removal of
2062 /// background before deconvolution is desired.
2063 /// - deconIterations-number of iterations in deconvolution operation.
2064 /// - markov: logical variable, if it is true, first the source spectrum
2065 /// is replaced by new spectrum calculated using Markov
2066 /// chains method.
2067 /// - averWindow: averaging window of searched peaks, for details
2068 /// we refer to manual (applies only for Markov method).
2069 ///
2070 /// ### Peaks searching:
2071 ///
2072 /// The goal of this function is to identify automatically the peaks in spectrum
2073 /// with the presence of the continuous background and statistical
2074 /// fluctuations - noise.
2075 ///
2076 /// The common problems connected with correct peak identification are:
2077 ///
2078 /// - non-sensitivity to noise, i.e., only statistically
2079 /// relevant peaks should be identified.
2080 /// - non-sensitivity of the algorithm to continuous
2081 /// background.
2082 /// - ability to identify peaks close to the edges of the
2083 /// spectrum region. Usually peak finders fail to detect them.
2084 /// - resolution, decomposition of Double_tts and multiplets.
2085 /// The algorithm should be able to recognise close positioned peaks.
2086 /// - ability to identify peaks with different sigma.
2087 ///
2088 /// \image html TSpectrum_Searching1.jpg Fig. 27 An example of one-dimensional synthetic spectrum with found peaks denoted by markers.
2089 ///
2090 /// #### References:
2091 ///
2092 /// 1. M.A. Mariscotti: A method for identification of peaks in the presence of
2093 /// background and its application to spectrum analysis. NIM 50 (1967),
2094 /// 309-320.
2095 /// 2. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky,
2096 /// I. Turzo.:Identification of peaks in
2097 /// multidimensional coincidence gamma-ray spectra. NIM, A443 (2000) 108-125.
2098 /// 3. Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM
2099 /// A 376 (1996), 451.
2100 ///
2101 /// Examples of peak searching method:
2102 ///
2103 /// The SearchHighRes function provides users with the possibility to vary the
2104 /// input parameters and with the access to the output deconvolved data in the
2105 /// destination spectrum. Based on the output data one can tune the parameters.
2106 ///
2107 /// ### Example 15 - script SearchHR1.C:
2108 ///
2109 /// One-dimensional spectrum with found peaks denoted by markers, 3 iterations
2110 /// steps in the deconvolution.
2111 ///
2112 /// #### Script:
2113 ///
2114 /// Begin_Macro(source)
2115 /// ../../../tutorials/spectrum/SearchHR1.C
2116 /// End_Macro
2117 ///
2118 /// ### Example 16 - script SearchHR3.C:
2119 ///
2120 /// Influence of number of iterations (3-red, 10-blue, 100- green, 1000-magenta),
2121 /// sigma=8, smoothing width=3.
2122 ///
2123 /// Begin_Macro(source)
2124 /// ../../../tutorials/spectrum/SearchHR3.C
2125 /// End_Macro
2126 
2127 Int_t TSpectrum::SearchHighRes(Double_t *source,Double_t *destVector, int ssize,
2128  Double_t sigma, Double_t threshold,
2129  bool backgroundRemove,int deconIterations,
2130  bool markov, int averWindow)
2131 {
2132  int i, j, numberIterations = (Int_t)(7 * sigma + 0.5);
2133  Double_t a, b, c;
2134  int k, lindex, posit, imin, imax, jmin, jmax, lh_gold, priz;
2135  Double_t lda, ldb, ldc, area, maximum, maximum_decon;
2136  int xmin, xmax, l, peak_index = 0, size_ext = ssize + 2 * numberIterations, shift = numberIterations, bw = 2, w;
2137  Double_t maxch;
2138  Double_t nom, nip, nim, sp, sm, plocha = 0;
2139  Double_t m0low=0,m1low=0,m2low=0,l0low=0,l1low=0,detlow,av,men;
2140  if (sigma < 1) {
2141  Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
2142  return 0;
2143  }
2144 
2145  if(threshold<=0 || threshold>=100){
2146  Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
2147  return 0;
2148  }
2149 
2150  j = (Int_t) (5.0 * sigma + 0.5);
2151  if (j >= PEAK_WINDOW / 2) {
2152  Error("SearchHighRes", "Too large sigma");
2153  return 0;
2154  }
2155 
2156  if (markov == true) {
2157  if (averWindow <= 0) {
2158  Error("SearchHighRes", "Averaging window must be positive");
2159  return 0;
2160  }
2161  }
2162 
2163  if(backgroundRemove == true){
2164  if(ssize < 2 * numberIterations + 1){
2165  Error("SearchHighRes", "Too large clipping window");
2166  return 0;
2167  }
2168  }
2169 
2170  k = int(2 * sigma+0.5);
2171  if(k >= 2){
2172  for(i = 0;i < k;i++){
2173  a = i,b = source[i];
2174  m0low += 1,m1low += a,m2low += a * a,l0low += b,l1low += a * b;
2175  }
2176  detlow = m0low * m2low - m1low * m1low;
2177  if(detlow != 0)
2178  l1low = (-l0low * m1low + l1low * m0low) / detlow;
2179 
2180  else
2181  l1low = 0;
2182  if(l1low > 0)
2183  l1low=0;
2184  }
2185 
2186  else{
2187  l1low = 0;
2188  }
2189 
2190  i = (Int_t)(7 * sigma + 0.5);
2191  i = 2 * i;
2192  Double_t *working_space = new Double_t [7 * (ssize + i)];
2193  for (j=0;j<7 * (ssize + i);j++) working_space[j] = 0;
2194  for(i = 0; i < size_ext; i++){
2195  if(i < shift){
2196  a = i - shift;
2197  working_space[i + size_ext] = source[0] + l1low * a;
2198  if(working_space[i + size_ext] < 0)
2199  working_space[i + size_ext]=0;
2200  }
2201 
2202  else if(i >= ssize + shift){
2203  a = i - (ssize - 1 + shift);
2204  working_space[i + size_ext] = source[ssize - 1];
2205  if(working_space[i + size_ext] < 0)
2206  working_space[i + size_ext]=0;
2207  }
2208 
2209  else
2210  working_space[i + size_ext] = source[i - shift];
2211  }
2212 
2213  if(backgroundRemove == true){
2214  for(i = 1; i <= numberIterations; i++){
2215  for(j = i; j < size_ext - i; j++){
2216  if(markov == false){
2217  a = working_space[size_ext + j];
2218  b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2219  if(b < a)
2220  a = b;
2221 
2222  working_space[j]=a;
2223  }
2224 
2225  else{
2226  a = working_space[size_ext + j];
2227  av = 0;
2228  men = 0;
2229  for (w = j - bw; w <= j + bw; w++){
2230  if ( w >= 0 && w < size_ext){
2231  av += working_space[size_ext + w];
2232  men +=1;
2233  }
2234  }
2235  av = av / men;
2236  b = 0;
2237  men = 0;
2238  for (w = j - i - bw; w <= j - i + bw; w++){
2239  if ( w >= 0 && w < size_ext){
2240  b += working_space[size_ext + w];
2241  men +=1;
2242  }
2243  }
2244  b = b / men;
2245  c = 0;
2246  men = 0;
2247  for (w = j + i - bw; w <= j + i + bw; w++){
2248  if ( w >= 0 && w < size_ext){
2249  c += working_space[size_ext + w];
2250  men +=1;
2251  }
2252  }
2253  c = c / men;
2254  b = (b + c) / 2;
2255  if (b < a)
2256  av = b;
2257  working_space[j]=av;
2258  }
2259  }
2260  for(j = i; j < size_ext - i; j++)
2261  working_space[size_ext + j] = working_space[j];
2262  }
2263  for(j = 0;j < size_ext; j++){
2264  if(j < shift){
2265  a = j - shift;
2266  b = source[0] + l1low * a;
2267  if (b < 0) b = 0;
2268  working_space[size_ext + j] = b - working_space[size_ext + j];
2269  }
2270 
2271  else if(j >= ssize + shift){
2272  a = j - (ssize - 1 + shift);
2273  b = source[ssize - 1];
2274  if (b < 0) b = 0;
2275  working_space[size_ext + j] = b - working_space[size_ext + j];
2276  }
2277 
2278  else{
2279  working_space[size_ext + j] = source[j - shift] - working_space[size_ext + j];
2280  }
2281  }
2282  for(j = 0;j < size_ext; j++){
2283  if(working_space[size_ext + j] < 0) working_space[size_ext + j] = 0;
2284  }
2285  }
2286 
2287  for(i = 0; i < size_ext; i++){
2288  working_space[i + 6*size_ext] = working_space[i + size_ext];
2289  }
2290 
2291  if(markov == true){
2292  for(j = 0; j < size_ext; j++)
2293  working_space[2 * size_ext + j] = working_space[size_ext + j];
2294  xmin = 0,xmax = size_ext - 1;
2295  for(i = 0, maxch = 0; i < size_ext; i++){
2296  working_space[i] = 0;
2297  if(maxch < working_space[2 * size_ext + i])
2298  maxch = working_space[2 * size_ext + i];
2299  plocha += working_space[2 * size_ext + i];
2300  }
2301  if(maxch == 0) {
2302  delete [] working_space;
2303  return 0;
2304  }
2305 
2306  nom = 1;
2307  working_space[xmin] = 1;
2308  for(i = xmin; i < xmax; i++){
2309  nip = working_space[2 * size_ext + i] / maxch;
2310  nim = working_space[2 * size_ext + i + 1] / maxch;
2311  sp = 0,sm = 0;
2312  for(l = 1; l <= averWindow; l++){
2313  if((i + l) > xmax)
2314  a = working_space[2 * size_ext + xmax] / maxch;
2315 
2316  else
2317  a = working_space[2 * size_ext + i + l] / maxch;
2318 
2319  b = a - nip;
2320  if(a + nip <= 0)
2321  a=1;
2322 
2323  else
2324  a = TMath::Sqrt(a + nip);
2325 
2326  b = b / a;
2327  b = TMath::Exp(b);
2328  sp = sp + b;
2329  if((i - l + 1) < xmin)
2330  a = working_space[2 * size_ext + xmin] / maxch;
2331 
2332  else
2333  a = working_space[2 * size_ext + i - l + 1] / maxch;
2334 
2335  b = a - nim;
2336  if(a + nim <= 0)
2337  a = 1;
2338 
2339  else
2340  a = TMath::Sqrt(a + nim);
2341 
2342  b = b / a;
2343  b = TMath::Exp(b);
2344  sm = sm + b;
2345  }
2346  a = sp / sm;
2347  a = working_space[i + 1] = working_space[i] * a;
2348  nom = nom + a;
2349  }
2350  for(i = xmin; i <= xmax; i++){
2351  working_space[i] = working_space[i] / nom;
2352  }
2353  for(j = 0; j < size_ext; j++)
2354  working_space[size_ext + j] = working_space[j] * plocha;
2355  for(j = 0; j < size_ext; j++){
2356  working_space[2 * size_ext + j] = working_space[size_ext + j];
2357  }
2358  if(backgroundRemove == true){
2359  for(i = 1; i <= numberIterations; i++){
2360  for(j = i; j < size_ext - i; j++){
2361  a = working_space[size_ext + j];
2362  b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2363  if(b < a)
2364  a = b;
2365  working_space[j] = a;
2366  }
2367  for(j = i; j < size_ext - i; j++)
2368  working_space[size_ext + j] = working_space[j];
2369  }
2370  for(j = 0; j < size_ext; j++){
2371  working_space[size_ext + j] = working_space[2 * size_ext + j] - working_space[size_ext + j];
2372  }
2373  }
2374  }
2375 //deconvolution starts
2376  area = 0;
2377  lh_gold = -1;
2378  posit = 0;
2379  maximum = 0;
2380 //generate response vector
2381  for(i = 0; i < size_ext; i++){
2382  lda = (Double_t)i - 3 * sigma;
2383  lda = lda * lda / (2 * sigma * sigma);
2384  j = (Int_t)(1000 * TMath::Exp(-lda));
2385  lda = j;
2386  if(lda != 0)
2387  lh_gold = i + 1;
2388 
2389  working_space[i] = lda;
2390  area = area + lda;
2391  if(lda > maximum){
2392  maximum = lda;
2393  posit = i;
2394  }
2395  }
2396 //read source vector
2397  for(i = 0; i < size_ext; i++)
2398  working_space[2 * size_ext + i] = TMath::Abs(working_space[size_ext + i]);
2399 //create matrix at*a(vector b)
2400  i = lh_gold - 1;
2401  if(i > size_ext)
2402  i = size_ext;
2403 
2404  imin = -i,imax = i;
2405  for(i = imin; i <= imax; i++){
2406  lda = 0;
2407  jmin = 0;
2408  if(i < 0)
2409  jmin = -i;
2410  jmax = lh_gold - 1 - i;
2411  if(jmax > (lh_gold - 1))
2412  jmax = lh_gold - 1;
2413 
2414  for(j = jmin;j <= jmax; j++){
2415  ldb = working_space[j];
2416  ldc = working_space[i + j];
2417  lda = lda + ldb * ldc;
2418  }
2419  working_space[size_ext + i - imin] = lda;
2420  }
2421 //create vector p
2422  i = lh_gold - 1;
2423  imin = -i,imax = size_ext + i - 1;
2424  for(i = imin; i <= imax; i++){
2425  lda = 0;
2426  for(j = 0; j <= (lh_gold - 1); j++){
2427  ldb = working_space[j];
2428  k = i + j;
2429  if(k >= 0 && k < size_ext){
2430  ldc = working_space[2 * size_ext + k];
2431  lda = lda + ldb * ldc;
2432  }
2433 
2434  }
2435  working_space[4 * size_ext + i - imin] = lda;
2436  }
2437 //move vector p
2438  for(i = imin; i <= imax; i++)
2439  working_space[2 * size_ext + i - imin] = working_space[4 * size_ext + i - imin];
2440 //initialization of resulting vector
2441  for(i = 0; i < size_ext; i++)
2442  working_space[i] = 1;
2443 //START OF ITERATIONS
2444  for(lindex = 0; lindex < deconIterations; lindex++){
2445  for(i = 0; i < size_ext; i++){
2446  if(TMath::Abs(working_space[2 * size_ext + i]) > 0.00001 && TMath::Abs(working_space[i]) > 0.00001){
2447  lda=0;
2448  jmin = lh_gold - 1;
2449  if(jmin > i)
2450  jmin = i;
2451 
2452  jmin = -jmin;
2453  jmax = lh_gold - 1;
2454  if(jmax > (size_ext - 1 - i))
2455  jmax=size_ext-1-i;
2456 
2457  for(j = jmin; j <= jmax; j++){
2458  ldb = working_space[j + lh_gold - 1 + size_ext];
2459  ldc = working_space[i + j];
2460  lda = lda + ldb * ldc;
2461  }
2462  ldb = working_space[2 * size_ext + i];
2463  if(lda != 0)
2464  lda = ldb / lda;
2465 
2466  else
2467  lda = 0;
2468 
2469  ldb = working_space[i];
2470  lda = lda * ldb;
2471  working_space[3 * size_ext + i] = lda;
2472  }
2473  }
2474  for(i = 0; i < size_ext; i++){
2475  working_space[i] = working_space[3 * size_ext + i];
2476  }
2477  }
2478 //shift resulting spectrum
2479  for(i=0;i<size_ext;i++){
2480  lda = working_space[i];
2481  j = i + posit;
2482  j = j % size_ext;
2483  working_space[size_ext + j] = lda;
2484  }
2485 //write back resulting spectrum
2486  maximum = 0, maximum_decon = 0;
2487  j = lh_gold - 1;
2488  for(i = 0; i < size_ext - j; i++){
2489  if(i >= shift && i < ssize + shift){
2490  working_space[i] = area * working_space[size_ext + i + j];
2491  if(maximum_decon < working_space[i])
2492  maximum_decon = working_space[i];
2493  if(maximum < working_space[6 * size_ext + i])
2494  maximum = working_space[6 * size_ext + i];
2495  }
2496 
2497  else
2498  working_space[i] = 0;
2499  }
2500  lda=1;
2501  if(lda>threshold)
2502  lda=threshold;
2503  lda=lda/100;
2504 
2505 //searching for peaks in deconvolved spectrum
2506  for(i = 1; i < size_ext - 1; i++){
2507  if(working_space[i] > working_space[i - 1] && working_space[i] > working_space[i + 1]){
2508  if(i >= shift && i < ssize + shift){
2509  if(working_space[i] > lda*maximum_decon && working_space[6 * size_ext + i] > threshold * maximum / 100.0){
2510  for(j = i - 1, a = 0, b = 0; j <= i + 1; j++){
2511  a += (Double_t)(j - shift) * working_space[j];
2512  b += working_space[j];
2513  }
2514  a = a / b;
2515  if(a < 0)
2516  a = 0;
2517 
2518  if(a >= ssize)
2519  a = ssize - 1;
2520  if(peak_index == 0){
2521  fPositionX[0] = a;
2522  peak_index = 1;
2523  }
2524 
2525  else{
2526  for(j = 0, priz = 0; j < peak_index && priz == 0; j++){
2527  if(working_space[6 * size_ext + shift + (Int_t)a] > working_space[6 * size_ext + shift + (Int_t)fPositionX[j]])
2528  priz = 1;
2529  }
2530  if(priz == 0){
2531  if(j < fMaxPeaks){
2532  fPositionX[j] = a;
2533  }
2534  }
2535 
2536  else{
2537  for(k = peak_index; k >= j; k--){
2538  if(k < fMaxPeaks){
2539  fPositionX[k] = fPositionX[k - 1];
2540  }
2541  }
2542  fPositionX[j - 1] = a;
2543  }
2544  if(peak_index < fMaxPeaks)
2545  peak_index += 1;
2546  }
2547  }
2548  }
2549  }
2550  }
2551 
2552  for(i = 0; i < ssize; i++) destVector[i] = working_space[i + shift];
2553  delete [] working_space;
2554  fNPeaks = peak_index;
2555  if(peak_index == fMaxPeaks)
2556  Warning("SearchHighRes", "Peak buffer full");
2557  return fNPeaks;
2558 }
2559 
2560 ////////////////////////////////////////////////////////////////////////////////
2561 /// Old name of SearcHighRes introduced for back compatibility.
2562 /// This function will be removed after the June 2006 release
2563 
2564 Int_t TSpectrum::Search1HighRes(Double_t *source,Double_t *destVector, int ssize,
2565  Double_t sigma, Double_t threshold,
2566  bool backgroundRemove,int deconIterations,
2567  bool markov, int averWindow)
2568 {
2569 
2570 
2571  return SearchHighRes(source,destVector,ssize,sigma,threshold,backgroundRemove,
2572  deconIterations,markov,averWindow);
2573 }
2574 
2575 ////////////////////////////////////////////////////////////////////////////////
2576 /// Static function, interface to TSpectrum::Search.
2577 
2578 Int_t TSpectrum::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
2580 
2581  TSpectrum s;
2582  return s.Search(hist,sigma,option,threshold);
2583 }
2584 
2585 ////////////////////////////////////////////////////////////////////////////////
2586 /// Static function, interface to TSpectrum::Background.
2587 
2588 TH1 *TSpectrum::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
2589 {
2590  TSpectrum s;
2591  return s.Background(hist,niter,option);
2592 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
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)
Old name of SearcHighRes introduced for back compatibility.
Definition: TSpectrum.cxx:2579
float xmin
Definition: THbookFile.cxx:93
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8253
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:444
const char Option_t
Definition: RtypesCore.h:62
Definition: Rtypes.h:56
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:640
TH1 * h
Definition: legend2.C:5
static TH1 * StaticBackground(const TH1 *hist, Int_t niter=20, Option_t *option="")
Static function, interface to TSpectrum::Background.
Definition: TSpectrum.cxx:2603
const char * DeconvolutionRL(Double_t *source, const Double_t *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional deconvolution function.
Definition: TSpectrum.cxx:1678
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4639
Basic string class.
Definition: TString.h:129
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1099
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:628
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:501
virtual Int_t GetDimension() const
Definition: TH1.h:263
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
Double_t * fPositionX
[fNPeaks] X position of peaks
Definition: TSpectrum.h:28
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)
One-dimensional high-resolution peak search function.
Definition: TSpectrum.cxx:2142
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
const Double_t sigma
virtual void Delete(Option_t *option="")
Delete this object.
Definition: TObject.cxx:176
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05)
One-dimensional peak search function.
Definition: TSpectrum.cxx:266
virtual ~TSpectrum()
Destructor.
Definition: TSpectrum.cxx:87
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:455
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:679
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:327
TH1 * fHistogram
resulting histogram
Definition: TSpectrum.h:31
Double_t * fPositionY
[fNPeaks] Y position of peaks
Definition: TSpectrum.h:29
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:873
const char * Deconvolution(Double_t *source, const Double_t *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional deconvolution function.
Definition: TSpectrum.cxx:1470
TLine * l
Definition: textangle.C:4
static void SetAverageWindow(Int_t w=3)
Static function: Set average window of searched peaks (see TSpectrum::SearchHighRes).
Definition: TSpectrum.cxx:99
static Int_t fgAverageWindow
Average window of searched peaks.
Definition: TSpectrum.h:32
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
float xmax
Definition: THbookFile.cxx:93
#define PEAK_WINDOW
Definition: TSpectrum.cxx:44
static void SetDeconIterations(Int_t n=3)
Static function: Set max number of decon iterations in deconvolution operation (see TSpectrum::Search...
Definition: TSpectrum.cxx:108
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:92
TSpectrum()
Constructor.
Definition: TSpectrum.cxx:50
Double_t Exp(Double_t x)
Definition: TMath.h:622
A PolyMarker is defined by an array on N points in a 2-D space.
Definition: TPolyMarker.h:31
#define ClassImp(name)
Definition: Rtypes.h:336
virtual void Print(Option_t *option="") const
Print the array of positions.
Definition: TSpectrum.cxx:219
double Double_t
Definition: RtypesCore.h:55
virtual TH1 * Background(const TH1 *hist, Int_t niter=20, Option_t *option="")
One-dimensional background estimation function.
Definition: TSpectrum.cxx:152
static Int_t fgIterations
Maximum number of decon iterations (default=3)
Definition: TSpectrum.h:33
Int_t fMaxPeaks
Maximum number of peaks to be found.
Definition: TSpectrum.h:25
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:572
The TH1 histogram class.
Definition: TH1.h:56
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
void SetResolution(Double_t resolution=1)
NOT USED resolution: determines resolution of the neighbouring peaks default value is 1 correspond to...
Definition: TSpectrum.cxx:351
Advanced Spectra Processing.
Definition: TSpectrum.h:18
#define hbname
Definition: THbookFile.cxx:123
Double_t * fPosition
[fNPeaks] array of current peak positions
Definition: TSpectrum.h:27
virtual void Add(TObject *obj)
Definition: TList.h:77
#define dest(otri, vertexptr)
Definition: triangle.c:1040
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
Static function, interface to TSpectrum::Search.
Definition: TSpectrum.cxx:2593
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2544
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
#define snprintf
Definition: civetweb.c:822
#define gPad
Definition: TVirtualPad.h:284
Int_t fNPeaks
number of peaks found
Definition: TSpectrum.h:26
Definition: first.py:1
Double_t Sqrt(Double_t x)
Definition: TMath.h:591
TList * GetListOfFunctions() const
Definition: TH1.h:224
const char * SmoothMarkov(Double_t *source, Int_t ssize, Int_t averWindow)
One-dimensional markov spectrum smoothing function.
Definition: TSpectrum.cxx:1203
const Bool_t kTRUE
Definition: RtypesCore.h:91
const Int_t n
Definition: legend1.C:16
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)
One-dimensional unfolding function.
Definition: TSpectrum.cxx:1903
Double_t fResolution
NOT USED resolution of the neighboring peaks
Definition: TSpectrum.h:30
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:859
TAxis * GetXaxis()
Definition: TH1.h:300
const char * Data() const
Definition: TString.h:347