Logo ROOT  
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 #include "snprintf.h"
11 
12 /** \class TSpectrum
13  \ingroup Spectrum
14  \brief Advanced Spectra Processing
15  \author Miroslav Morhac
16 
17  This class contains advanced spectra processing functions for:
18 
19  - One-dimensional background estimation
20  - One-dimensional smoothing
21  - One-dimensional deconvolution
22  - One-dimensional peak search
23 
24  The algorithms in this class have been published in the following references:
25 
26  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.
27  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.
28  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.
29 
30  These NIM papers are also available as doc or ps files from:
31 
32  - [Spectrum.doc](https://root.cern.ch/download/Spectrum.doc)
33  - [SpectrumDec.ps.gz](https://root.cern.ch/download/SpectrumDec.ps.gz)
34  - [SpectrumSrc.ps.gz](https://root.cern.ch/download/SpectrumSrc.ps.gz)
35  - [SpectrumBck.ps.gz](https://root.cern.ch/download/SpectrumBck.ps.gz)
36 
37  See also the
38  [online documentation](https://root.cern.ch/guides/tspectrum-manual) and
39  [tutorials](https://root.cern.ch/doc/master/group__tutorial__spectrum.html).
40 */
41 
44 
45 #define PEAK_WINDOW 1024
47 
48 ////////////////////////////////////////////////////////////////////////////////
49 /// Constructor.
50 
51 TSpectrum::TSpectrum() :TNamed("Spectrum", "Miroslav Morhac peak finder")
52 {
53  Int_t n = 100;
54  fMaxPeaks = n;
55  fPosition = new Double_t[n];
56  fPositionX = new Double_t[n];
57  fPositionY = new Double_t[n];
58  fResolution = 1;
59  fHistogram = 0;
60  fNPeaks = 0;
61 }
62 
63 ////////////////////////////////////////////////////////////////////////////////
64 /// - maxpositions: maximum number of peaks
65 /// - resolution: *NOT USED* determines resolution of the neighbouring peaks
66 /// default value is 1 correspond to 3 sigma distance
67 /// between peaks. Higher values allow higher resolution
68 /// (smaller distance between peaks.
69 /// May be set later through SetResolution.
70 
71 TSpectrum::TSpectrum(Int_t maxpositions, Double_t resolution)
72  :TNamed("Spectrum", "Miroslav Morhac peak finder")
73 {
74  Int_t n = maxpositions;
75  if (n <= 0) n = 1;
76  fMaxPeaks = n;
77  fPosition = new Double_t[n];
78  fPositionX = new Double_t[n];
79  fPositionY = new Double_t[n];
80  fHistogram = 0;
81  fNPeaks = 0;
82  SetResolution(resolution);
83 }
84 
85 ////////////////////////////////////////////////////////////////////////////////
86 /// Destructor.
87 
89 {
90  delete [] fPosition;
91  delete [] fPositionX;
92  delete [] fPositionY;
93  delete fHistogram;
94 }
95 
96 ////////////////////////////////////////////////////////////////////////////////
97 /// Static function: Set average window of searched peaks
98 /// (see TSpectrum::SearchHighRes).
99 
101 {
102  fgAverageWindow = w;
103 }
104 
105 ////////////////////////////////////////////////////////////////////////////////
106 /// Static function: Set max number of decon iterations in deconvolution
107 /// operation (see TSpectrum::SearchHighRes).
108 
110 {
111  fgIterations = n;
112 }
113 
114 ////////////////////////////////////////////////////////////////////////////////
115 /// One-dimensional background estimation function.
116 ///
117 /// This function calculates the background spectrum in the input histogram h.
118 /// The background is returned as a histogram.
119 ///
120 /// #### Parameters:
121 ///
122 /// - h: input 1-d histogram
123 /// - numberIterations, (default value = 20)
124 /// Increasing numberIterations make the result smoother and lower.
125 /// - option: may contain one of the following options:
126 ///
127 /// - to set the direction parameter
128 /// "BackIncreasingWindow". By default the direction is BackDecreasingWindow
129 /// - filterOrder-order of clipping filter, (default "BackOrder2")
130 /// -possible values= "BackOrder4"
131 /// "BackOrder6"
132 /// "BackOrder8"
133 /// - "nosmoothing"- if selected, the background is not smoothed
134 /// By default the background is smoothed.
135 /// - smoothWindow-width of smoothing window, (default is "BackSmoothing3")
136 /// -possible values= "BackSmoothing5"
137 /// "BackSmoothing7"
138 /// "BackSmoothing9"
139 /// "BackSmoothing11"
140 /// "BackSmoothing13"
141 /// "BackSmoothing15"
142 /// - "Compton" if selected the estimation of Compton edge
143 /// will be included.
144 /// - "same" : if this option is specified, the resulting background
145 /// histogram is superimposed on the picture in the current pad.
146 ///
147 /// NOTE that the background is only evaluated in the current range of h.
148 /// ie, if h has a bin range (set via `h->GetXaxis()->SetRange(binmin,binmax)`,
149 /// the returned histogram will be created with the same number of bins
150 /// as the input histogram h, but only bins from `binmin` to `binmax` will be filled
151 /// with the estimated background.
152 
153 TH1 *TSpectrum::Background(const TH1 * h, int numberIterations,
154  Option_t * option)
155 {
156  if (h == 0) return 0;
157  Int_t dimension = h->GetDimension();
158  if (dimension > 1) {
159  Error("Search", "Only implemented for 1-d histograms");
160  return 0;
161  }
162  TString opt = option;
163  opt.ToLower();
164 
165  //set options
166  Int_t direction = kBackDecreasingWindow;
167  if (opt.Contains("backincreasingwindow")) direction = kBackIncreasingWindow;
168  Int_t filterOrder = kBackOrder2;
169  if (opt.Contains("backorder4")) filterOrder = kBackOrder4;
170  if (opt.Contains("backorder6")) filterOrder = kBackOrder6;
171  if (opt.Contains("backorder8")) filterOrder = kBackOrder8;
172  Bool_t smoothing = kTRUE;
173  if (opt.Contains("nosmoothing")) smoothing = kFALSE;
174  Int_t smoothWindow = kBackSmoothing3;
175  if (opt.Contains("backsmoothing5")) smoothWindow = kBackSmoothing5;
176  if (opt.Contains("backsmoothing7")) smoothWindow = kBackSmoothing7;
177  if (opt.Contains("backsmoothing9")) smoothWindow = kBackSmoothing9;
178  if (opt.Contains("backsmoothing11")) smoothWindow = kBackSmoothing11;
179  if (opt.Contains("backsmoothing13")) smoothWindow = kBackSmoothing13;
180  if (opt.Contains("backsmoothing15")) smoothWindow = kBackSmoothing15;
181  Bool_t compton = kFALSE;
182  if (opt.Contains("compton")) compton = kTRUE;
183 
184  Int_t first = h->GetXaxis()->GetFirst();
185  Int_t last = h->GetXaxis()->GetLast();
186  Int_t size = last-first+1;
187  Int_t i;
188  Double_t * source = new Double_t[size];
189  for (i = 0; i < size; i++) source[i] = h->GetBinContent(i + first);
190 
191  //find background (source is input and in output contains the background
192  Background(source,size,numberIterations, direction, filterOrder,smoothing,
193  smoothWindow,compton);
194 
195  //create output histogram containing background
196  //only bins in the range of the input histogram are filled
197  Int_t nch = strlen(h->GetName());
198  char *hbname = new char[nch+20];
199  snprintf(hbname,nch+20,"%s_background",h->GetName());
200  TH1 *hb = (TH1*)h->Clone(hbname);
201  hb->Reset();
202  hb->GetListOfFunctions()->Delete();
203  hb->SetLineColor(2);
204  for (i=0; i< size; i++) hb->SetBinContent(i+first,source[i]);
205  hb->SetEntries(size);
206 
207  //if option "same is specified, draw the result in the pad
208  if (opt.Contains("same")) {
209  if (gPad) delete gPad->GetPrimitive(hbname);
210  hb->Draw("same");
211  }
212  delete [] source;
213  delete [] hbname;
214  return hb;
215 }
216 
217 ////////////////////////////////////////////////////////////////////////////////
218 /// Print the array of positions.
219 
221 {
222  printf("\nNumber of positions = %d\n",fNPeaks);
223  for (Int_t i=0;i<fNPeaks;i++) {
224  printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
225  }
226 }
227 
228 ////////////////////////////////////////////////////////////////////////////////
229 /// One-dimensional peak search function
230 ///
231 /// This function searches for peaks in source spectrum in hin
232 /// The number of found peaks and their positions are written into
233 /// the members fNpeaks and fPositionX.
234 /// The search is performed in the current histogram range.
235 ///
236 /// #### Parameters:
237 ///
238 /// - hin: pointer to the histogram of source spectrum
239 /// - sigma: sigma of searched peaks, for details we refer to manual
240 /// - threshold: (default=0.05) peaks with amplitude less than
241 /// threshold*highest_peak are discarded. 0<threshold<1
242 ///
243 /// By default, the background is removed before deconvolution.
244 /// Specify the option "nobackground" to not remove the background.
245 ///
246 /// By default the "Markov" chain algorithm is used.
247 /// Specify the option "noMarkov" to disable this algorithm
248 /// Note that by default the source spectrum is replaced by a new spectrum
249 ///
250 /// By default a polymarker object is created and added to the list of
251 /// functions of the histogram. The histogram is drawn with the specified
252 /// option and the polymarker object drawn on top of the histogram.
253 /// The polymarker coordinates correspond to the npeaks peaks found in
254 /// the histogram.
255 ///
256 /// A pointer to the polymarker object can be retrieved later via:
257 /// ~~~ {.cpp}
258 /// TList *functions = hin->GetListOfFunctions();
259 /// TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
260 /// ~~~
261 /// Specify the option "goff" to disable the storage and drawing of the
262 /// polymarker.
263 ///
264 /// To disable the final drawing of the histogram with the search results (in case
265 /// you want to draw it yourself) specify "nodraw" in the options parameter.
266 
268  Double_t threshold)
269 {
270  if (hin == 0) return 0;
271  Int_t dimension = hin->GetDimension();
272  if (dimension > 2) {
273  Error("Search", "Only implemented for 1-d and 2-d histograms");
274  return 0;
275  }
276  if (threshold <=0 || threshold >= 1) {
277  Warning("Search","threshold must 0<threshold<1, threshold=0.05 assumed");
278  threshold = 0.05;
279  }
280  TString opt = option;
281  opt.ToLower();
282  Bool_t background = kTRUE;
283  if (opt.Contains("nobackground")) {
284  background = kFALSE;
285  opt.ReplaceAll("nobackground","");
286  }
287  Bool_t markov = kTRUE;
288  if (opt.Contains("nomarkov")) {
289  markov = kFALSE;
290  opt.ReplaceAll("nomarkov","");
291  }
292  Bool_t draw = kTRUE;
293  if (opt.Contains("nodraw")) {
294  draw = kFALSE;
295  opt.ReplaceAll("nodraw","");
296  }
297  if (dimension == 1) {
298  Int_t first = hin->GetXaxis()->GetFirst();
299  Int_t last = hin->GetXaxis()->GetLast();
300  Int_t size = last-first+1;
301  Int_t i, bin, npeaks;
302  Double_t * source = new Double_t[size];
303  Double_t * dest = new Double_t[size];
304  for (i = 0; i < size; i++) source[i] = hin->GetBinContent(i + first);
305  if (sigma < 1) {
306  sigma = size/fMaxPeaks;
307  if (sigma < 1) sigma = 1;
308  if (sigma > 8) sigma = 8;
309  }
310  npeaks = SearchHighRes(source, dest, size, sigma, 100*threshold,
311  background, fgIterations, markov, fgAverageWindow);
312 
313  for (i = 0; i < npeaks; i++) {
314  bin = first + Int_t(fPositionX[i] + 0.5);
315  fPositionX[i] = hin->GetBinCenter(bin);
316  fPositionY[i] = hin->GetBinContent(bin);
317  }
318  delete [] source;
319  delete [] dest;
320 
321  if (opt.Contains("goff"))
322  return npeaks;
323  if (!npeaks) return 0;
324  TPolyMarker * pm =
325  (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
326  if (pm) {
327  hin->GetListOfFunctions()->Remove(pm);
328  delete pm;
329  }
330  pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
331  hin->GetListOfFunctions()->Add(pm);
332  pm->SetMarkerStyle(23);
333  pm->SetMarkerColor(kRed);
334  pm->SetMarkerSize(1.3);
335  opt.ReplaceAll(" ","");
336  opt.ReplaceAll(",","");
337  if (draw)
338  ((TH1*)hin)->Draw(opt.Data());
339  return npeaks;
340  }
341  return 0;
342 }
343 
344 ////////////////////////////////////////////////////////////////////////////////
345 /// *NOT USED*
346 /// resolution: determines resolution of the neighbouring peaks
347 /// default value is 1 correspond to 3 sigma distance
348 /// between peaks. Higher values allow higher resolution
349 /// (smaller distance between peaks.
350 /// May be set later through SetResolution.
351 
353 {
354  if (resolution > 1)
355  fResolution = resolution;
356  else
357  fResolution = 1;
358 }
359 
360 ////////////////////////////////////////////////////////////////////////////////
361 /// This function calculates background spectrum from source spectrum.
362 /// The result is placed in the vector pointed by spe1945ctrum pointer.
363 /// The goal is to separate the useful information (peaks) from useless
364 /// information (background).
365 ///
366 /// - method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping
367 /// algorithm.
368 /// - new value in the channel "i" is calculated
369 ///
370 /// \f[
371 /// v_p(i) = min \left\{ v_{p-1}(i)^{\frac{\left[v_{p-1}(i+p)+v_{p-1}(i-p)\right]}{2}} \right\}
372 /// \f]
373 ///
374 /// where p = 1, 2, ..., numberIterations. In fact it represents second order
375 /// difference filter (-1,2,-1).
376 ///
377 /// One can also change the
378 /// direction of the change of the clipping window, the order of the clipping
379 /// filter, to include smoothing, to set width of smoothing window and to include
380 /// the estimation of Compton edges. On successful completion it returns 0. On
381 /// error it returns pointer to the string describing error.
382 ///
383 /// #### Parameters:
384 ///
385 /// - spectrum: pointer to the vector of source spectrum
386 /// - ssize: length of the spectrum vector
387 /// - numberIterations: maximal width of clipping window,
388 /// - direction: direction of change of clipping window.
389 /// Possible values: kBackIncreasingWindow, kBackDecreasingWindow
390 /// - filterOrder: order of clipping filter.
391 /// Possible values: kBackOrder2, kBackOrder4, kBackOrder6, kBackOrder8
392 /// - smoothing: logical variable whether the smoothing operation in the
393 /// estimation of background will be included.
394 /// Possible values: kFALSE, kTRUE
395 /// - smoothWindow: width of smoothing window.
396 /// Possible values: kBackSmoothing3, kBackSmoothing5, kBackSmoothing7,
397 /// kBackSmoothing9, kBackSmoothing11, kBackSmoothing13, kBackSmoothing15.
398 /// - compton: logical variable whether the estimation of Compton edge will be
399 /// included. Possible values: kFALSE, kTRUE.
400 ///
401 /// #### References:
402 ///
403 /// 1. C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the
404 /// quantitative analysis of PIXE spectra in geoscience applications. NIM, B34
405 /// (1988), 396-402.
406 ///
407 /// 2. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo:
408 /// Background elimination methods for multidimensional gamma-ray spectra. NIM,
409 /// A401 (1997) 113-132.
410 ///
411 /// 3. D. D. Burgess, R. J. Tervo: Background estimation for gamma-ray
412 /// spectroscopy. NIM 214 (1983), 431-434.
413 ///
414 /// ### Example 1 script Background_incr.C:
415 ///
416 /// Example of the estimation of background for number of iterations=6.
417 /// Original spectrum is shown in black color, estimated background in red color.
418 ///
419 /// Begin_Macro(source)
420 /// ../../../tutorials/spectrum/Background_incr.C
421 /// End_Macro
422 ///
423 /// ### Example 2 script Background_decr.C:
424 ///
425 /// In Example 1. one can notice that at the edges of the peaks the estimated
426 /// background goes under the peaks. An alternative approach is to decrease the
427 /// clipping window from a given value numberIterations to the value of one, which
428 /// is presented in this example.
429 ///
430 /// Example of the estimation of background for numberIterations=6 using
431 /// decreasing clipping window algorithm. Original spectrum is shown in black
432 /// color, estimated background in red color.
433 ///
434 /// Begin_Macro(source)
435 /// ../../../tutorials/spectrum/Background_decr.C
436 /// End_Macro
437 ///
438 /// ### Example 3 script Background_width.C:
439 ///
440 /// The question is how to choose the width of the clipping window, i.e.,
441 /// numberIterations parameter. The influence of this parameter on the estimated
442 /// background is illustrated in Example 3.
443 ///
444 /// Example of the influence of clipping window width on the estimated background
445 /// for numberIterations=4 (red line), 6 (orange line) 8 (green line) using decreasing
446 /// clipping window algorithm.
447 ///
448 /// in general one should set this parameter so that the value
449 /// 2*numberIterations+1 was greater than the widths of preserved objects (peaks).
450 ///
451 /// Begin_Macro(source)
452 /// ../../../tutorials/spectrum/Background_width.C
453 /// End_Macro
454 ///
455 /// ### Example 4 script Background_width2.C:
456 ///
457 /// another example for very complex spectrum is given here.
458 ///
459 /// Example of the influence of clipping window width on the estimated background
460 /// for numberIterations=10 (red line), 20 (blue line), 30 (green line) and
461 /// 40 (magenta line) using decreasing clipping window algorithm.
462 ///
463 /// Begin_Macro(source)
464 /// ../../../tutorials/spectrum/Background_width2.C
465 /// End_Macro
466 ///
467 /// ### Example 5 script Background_order.C:
468 ///
469 /// Second order difference filter removes linear (quasi-linear) background and
470 /// preserves symmetrical peaks. However if the shape of the background is more
471 /// complex one can employ higher-order clipping filters.
472 ///
473 /// Example of the influence of clipping filter difference order on the estimated
474 /// background for fNnumberIterations=40, 2-nd order red line, 4-th order blue line,
475 /// 6-th order green line and 8-th order magenta line, and using decreasing
476 /// clipping window algorithm.
477 ///
478 /// Begin_Macro(source)
479 /// ../../../tutorials/spectrum/Background_order.C
480 /// End_Macro
481 ///
482 /// ### Example 6 script Background_smooth.C:
483 ///
484 /// The estimate of the background can be influenced by noise present in the
485 /// spectrum. We proposed the algorithm of the background estimate with
486 /// simultaneous smoothing. In the original algorithm without smoothing, the
487 /// estimated background snatches the lower spikes in the noise. Consequently,
488 /// the areas of peaks are biased by this error.
489 ///
490 /// \image html TSpectrum_Background_smooth1.jpg Principle of background estimation algorithm with simultaneous smoothing.
491 ///
492 /// Begin_Macro(source)
493 /// ../../../tutorials/spectrum/Background_smooth.C
494 /// End_Macro
495 ///
496 /// ### Example 8 script Background_compton.C:
497 ///
498 /// Sometimes it is necessary to include also the Compton edges into the estimate of
499 /// the background. This example presents the synthetic spectrum
500 /// with Compton edges. The background was estimated using the 8-th order filter
501 /// with the estimation of the Compton edges using decreasing
502 /// clipping window algorithm (numberIterations=10) with smoothing
503 /// (smoothingWindow=5).
504 ///
505 /// Example of the estimate of the background with Compton edges (red line) for
506 /// numberIterations=10, 8-th order difference filter, using decreasing clipping
507 /// window algorithm and smoothing (smoothingWindow=5).
508 ///
509 /// Begin_Macro(source)
510 /// ../../../tutorials/spectrum/Background_compton.C
511 /// End_Macro
512 
513 const char *TSpectrum::Background(Double_t *spectrum, int ssize,
514  int numberIterations,
515  int direction, int filterOrder,
516  bool smoothing,int smoothWindow,
517  bool compton)
518 {
519  int i, j, w, bw, b1, b2, priz;
520  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;
521  if (ssize <= 0)
522  return "Wrong Parameters";
523  if (numberIterations < 1)
524  return "Width of Clipping Window Must Be Positive";
525  if (ssize < 2 * numberIterations + 1)
526  return "Too Large Clipping Window";
527  if (smoothing == kTRUE && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15)
528  return "Incorrect width of smoothing window";
529  Double_t *working_space = new Double_t[2 * ssize];
530  for (i = 0; i < ssize; i++){
531  working_space[i] = spectrum[i];
532  working_space[i + ssize] = spectrum[i];
533  }
534  bw=(smoothWindow-1)/2;
535  if (direction == kBackIncreasingWindow)
536  i = 1;
537  else if(direction == kBackDecreasingWindow)
538  i = numberIterations;
539  if (filterOrder == kBackOrder2) {
540  do{
541  for (j = i; j < ssize - i; j++) {
542  if (smoothing == kFALSE){
543  a = working_space[ssize + j];
544  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
545  if (b < a)
546  a = b;
547  working_space[j] = a;
548  }
549 
550  else if (smoothing == kTRUE){
551  a = working_space[ssize + j];
552  av = 0;
553  men = 0;
554  for (w = j - bw; w <= j + bw; w++){
555  if ( w >= 0 && w < ssize){
556  av += working_space[ssize + w];
557  men +=1;
558  }
559  }
560  av = av / men;
561  b = 0;
562  men = 0;
563  for (w = j - i - bw; w <= j - i + bw; w++){
564  if ( w >= 0 && w < ssize){
565  b += working_space[ssize + w];
566  men +=1;
567  }
568  }
569  b = b / men;
570  c = 0;
571  men = 0;
572  for (w = j + i - bw; w <= j + i + bw; w++){
573  if ( w >= 0 && w < ssize){
574  c += working_space[ssize + w];
575  men +=1;
576  }
577  }
578  c = c / men;
579  b = (b + c) / 2;
580  if (b < a)
581  av = b;
582  working_space[j]=av;
583  }
584  }
585  for (j = i; j < ssize - i; j++)
586  working_space[ssize + j] = working_space[j];
587  if (direction == kBackIncreasingWindow)
588  i+=1;
589  else if(direction == kBackDecreasingWindow)
590  i-=1;
591  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
592  }
593 
594  else if (filterOrder == kBackOrder4) {
595  do{
596  for (j = i; j < ssize - i; j++) {
597  if (smoothing == kFALSE){
598  a = working_space[ssize + j];
599  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
600  c = 0;
601  ai = i / 2;
602  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
603  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
604  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
605  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
606  if (b < c)
607  b = c;
608  if (b < a)
609  a = b;
610  working_space[j] = a;
611  }
612 
613  else if (smoothing == kTRUE){
614  a = working_space[ssize + j];
615  av = 0;
616  men = 0;
617  for (w = j - bw; w <= j + bw; w++){
618  if ( w >= 0 && w < ssize){
619  av += working_space[ssize + w];
620  men +=1;
621  }
622  }
623  av = av / men;
624  b = 0;
625  men = 0;
626  for (w = j - i - bw; w <= j - i + bw; w++){
627  if ( w >= 0 && w < ssize){
628  b += working_space[ssize + w];
629  men +=1;
630  }
631  }
632  b = b / men;
633  c = 0;
634  men = 0;
635  for (w = j + i - bw; w <= j + i + bw; w++){
636  if ( w >= 0 && w < ssize){
637  c += working_space[ssize + w];
638  men +=1;
639  }
640  }
641  c = c / men;
642  b = (b + c) / 2;
643  ai = i / 2;
644  b4 = 0, men = 0;
645  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
646  if (w >= 0 && w < ssize){
647  b4 += working_space[ssize + w];
648  men +=1;
649  }
650  }
651  b4 = b4 / men;
652  c4 = 0, men = 0;
653  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
654  if (w >= 0 && w < ssize){
655  c4 += working_space[ssize + w];
656  men +=1;
657  }
658  }
659  c4 = c4 / men;
660  d4 = 0, men = 0;
661  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
662  if (w >= 0 && w < ssize){
663  d4 += working_space[ssize + w];
664  men +=1;
665  }
666  }
667  d4 = d4 / men;
668  e4 = 0, men = 0;
669  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
670  if (w >= 0 && w < ssize){
671  e4 += working_space[ssize + w];
672  men +=1;
673  }
674  }
675  e4 = e4 / men;
676  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
677  if (b < b4)
678  b = b4;
679  if (b < a)
680  av = b;
681  working_space[j]=av;
682  }
683  }
684  for (j = i; j < ssize - i; j++)
685  working_space[ssize + j] = working_space[j];
686  if (direction == kBackIncreasingWindow)
687  i+=1;
688  else if(direction == kBackDecreasingWindow)
689  i-=1;
690  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
691  }
692 
693  else if (filterOrder == kBackOrder6) {
694  do{
695  for (j = i; j < ssize - i; j++) {
696  if (smoothing == kFALSE){
697  a = working_space[ssize + j];
698  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
699  c = 0;
700  ai = i / 2;
701  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
702  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
703  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
704  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
705  d = 0;
706  ai = i / 3;
707  d += working_space[ssize + j - (Int_t) (3 * ai)] / 20;
708  d -= 6 * working_space[ssize + j - (Int_t) (2 * ai)] / 20;
709  d += 15 * working_space[ssize + j - (Int_t) ai] / 20;
710  d += 15 * working_space[ssize + j + (Int_t) ai] / 20;
711  d -= 6 * working_space[ssize + j + (Int_t) (2 * ai)] / 20;
712  d += working_space[ssize + j + (Int_t) (3 * ai)] / 20;
713  if (b < d)
714  b = d;
715  if (b < c)
716  b = c;
717  if (b < a)
718  a = b;
719  working_space[j] = a;
720  }
721 
722  else if (smoothing == kTRUE){
723  a = working_space[ssize + j];
724  av = 0;
725  men = 0;
726  for (w = j - bw; w <= j + bw; w++){
727  if ( w >= 0 && w < ssize){
728  av += working_space[ssize + w];
729  men +=1;
730  }
731  }
732  av = av / men;
733  b = 0;
734  men = 0;
735  for (w = j - i - bw; w <= j - i + bw; w++){
736  if ( w >= 0 && w < ssize){
737  b += working_space[ssize + w];
738  men +=1;
739  }
740  }
741  b = b / men;
742  c = 0;
743  men = 0;
744  for (w = j + i - bw; w <= j + i + bw; w++){
745  if ( w >= 0 && w < ssize){
746  c += working_space[ssize + w];
747  men +=1;
748  }
749  }
750  c = c / men;
751  b = (b + c) / 2;
752  ai = i / 2;
753  b4 = 0, men = 0;
754  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
755  if (w >= 0 && w < ssize){
756  b4 += working_space[ssize + w];
757  men +=1;
758  }
759  }
760  b4 = b4 / men;
761  c4 = 0, men = 0;
762  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
763  if (w >= 0 && w < ssize){
764  c4 += working_space[ssize + w];
765  men +=1;
766  }
767  }
768  c4 = c4 / men;
769  d4 = 0, men = 0;
770  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
771  if (w >= 0 && w < ssize){
772  d4 += working_space[ssize + w];
773  men +=1;
774  }
775  }
776  d4 = d4 / men;
777  e4 = 0, men = 0;
778  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
779  if (w >= 0 && w < ssize){
780  e4 += working_space[ssize + w];
781  men +=1;
782  }
783  }
784  e4 = e4 / men;
785  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
786  ai = i / 3;
787  b6 = 0, men = 0;
788  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
789  if (w >= 0 && w < ssize){
790  b6 += working_space[ssize + w];
791  men +=1;
792  }
793  }
794  b6 = b6 / men;
795  c6 = 0, men = 0;
796  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
797  if (w >= 0 && w < ssize){
798  c6 += working_space[ssize + w];
799  men +=1;
800  }
801  }
802  c6 = c6 / men;
803  d6 = 0, men = 0;
804  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
805  if (w >= 0 && w < ssize){
806  d6 += working_space[ssize + w];
807  men +=1;
808  }
809  }
810  d6 = d6 / men;
811  e6 = 0, men = 0;
812  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
813  if (w >= 0 && w < ssize){
814  e6 += working_space[ssize + w];
815  men +=1;
816  }
817  }
818  e6 = e6 / men;
819  f6 = 0, men = 0;
820  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
821  if (w >= 0 && w < ssize){
822  f6 += working_space[ssize + w];
823  men +=1;
824  }
825  }
826  f6 = f6 / men;
827  g6 = 0, men = 0;
828  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
829  if (w >= 0 && w < ssize){
830  g6 += working_space[ssize + w];
831  men +=1;
832  }
833  }
834  g6 = g6 / men;
835  b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
836  if (b < b6)
837  b = b6;
838  if (b < b4)
839  b = b4;
840  if (b < a)
841  av = b;
842  working_space[j]=av;
843  }
844  }
845  for (j = i; j < ssize - i; j++)
846  working_space[ssize + j] = working_space[j];
847  if (direction == kBackIncreasingWindow)
848  i+=1;
849  else if(direction == kBackDecreasingWindow)
850  i-=1;
851  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
852  }
853 
854  else if (filterOrder == kBackOrder8) {
855  do{
856  for (j = i; j < ssize - i; j++) {
857  if (smoothing == kFALSE){
858  a = working_space[ssize + j];
859  b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
860  c = 0;
861  ai = i / 2;
862  c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
863  c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
864  c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
865  c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
866  d = 0;
867  ai = i / 3;
868  d += working_space[ssize + j - (Int_t) (3 * ai)] / 20;
869  d -= 6 * working_space[ssize + j - (Int_t) (2 * ai)] / 20;
870  d += 15 * working_space[ssize + j - (Int_t) ai] / 20;
871  d += 15 * working_space[ssize + j + (Int_t) ai] / 20;
872  d -= 6 * working_space[ssize + j + (Int_t) (2 * ai)] / 20;
873  d += working_space[ssize + j + (Int_t) (3 * ai)] / 20;
874  e = 0;
875  ai = i / 4;
876  e -= working_space[ssize + j - (Int_t) (4 * ai)] / 70;
877  e += 8 * working_space[ssize + j - (Int_t) (3 * ai)] / 70;
878  e -= 28 * working_space[ssize + j - (Int_t) (2 * ai)] / 70;
879  e += 56 * working_space[ssize + j - (Int_t) ai] / 70;
880  e += 56 * working_space[ssize + j + (Int_t) ai] / 70;
881  e -= 28 * working_space[ssize + j + (Int_t) (2 * ai)] / 70;
882  e += 8 * working_space[ssize + j + (Int_t) (3 * ai)] / 70;
883  e -= working_space[ssize + j + (Int_t) (4 * ai)] / 70;
884  if (b < e)
885  b = e;
886  if (b < d)
887  b = d;
888  if (b < c)
889  b = c;
890  if (b < a)
891  a = b;
892  working_space[j] = a;
893  }
894 
895  else if (smoothing == kTRUE){
896  a = working_space[ssize + j];
897  av = 0;
898  men = 0;
899  for (w = j - bw; w <= j + bw; w++){
900  if ( w >= 0 && w < ssize){
901  av += working_space[ssize + w];
902  men +=1;
903  }
904  }
905  av = av / men;
906  b = 0;
907  men = 0;
908  for (w = j - i - bw; w <= j - i + bw; w++){
909  if ( w >= 0 && w < ssize){
910  b += working_space[ssize + w];
911  men +=1;
912  }
913  }
914  b = b / men;
915  c = 0;
916  men = 0;
917  for (w = j + i - bw; w <= j + i + bw; w++){
918  if ( w >= 0 && w < ssize){
919  c += working_space[ssize + w];
920  men +=1;
921  }
922  }
923  c = c / men;
924  b = (b + c) / 2;
925  ai = i / 2;
926  b4 = 0, men = 0;
927  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
928  if (w >= 0 && w < ssize){
929  b4 += working_space[ssize + w];
930  men +=1;
931  }
932  }
933  b4 = b4 / men;
934  c4 = 0, men = 0;
935  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
936  if (w >= 0 && w < ssize){
937  c4 += working_space[ssize + w];
938  men +=1;
939  }
940  }
941  c4 = c4 / men;
942  d4 = 0, men = 0;
943  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
944  if (w >= 0 && w < ssize){
945  d4 += working_space[ssize + w];
946  men +=1;
947  }
948  }
949  d4 = d4 / men;
950  e4 = 0, men = 0;
951  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
952  if (w >= 0 && w < ssize){
953  e4 += working_space[ssize + w];
954  men +=1;
955  }
956  }
957  e4 = e4 / men;
958  b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
959  ai = i / 3;
960  b6 = 0, men = 0;
961  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
962  if (w >= 0 && w < ssize){
963  b6 += working_space[ssize + w];
964  men +=1;
965  }
966  }
967  b6 = b6 / men;
968  c6 = 0, men = 0;
969  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
970  if (w >= 0 && w < ssize){
971  c6 += working_space[ssize + w];
972  men +=1;
973  }
974  }
975  c6 = c6 / men;
976  d6 = 0, men = 0;
977  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
978  if (w >= 0 && w < ssize){
979  d6 += working_space[ssize + w];
980  men +=1;
981  }
982  }
983  d6 = d6 / men;
984  e6 = 0, men = 0;
985  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
986  if (w >= 0 && w < ssize){
987  e6 += working_space[ssize + w];
988  men +=1;
989  }
990  }
991  e6 = e6 / men;
992  f6 = 0, men = 0;
993  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
994  if (w >= 0 && w < ssize){
995  f6 += working_space[ssize + w];
996  men +=1;
997  }
998  }
999  f6 = f6 / men;
1000  g6 = 0, men = 0;
1001  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1002  if (w >= 0 && w < ssize){
1003  g6 += working_space[ssize + w];
1004  men +=1;
1005  }
1006  }
1007  g6 = g6 / men;
1008  b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1009  ai = i / 4;
1010  b8 = 0, men = 0;
1011  for (w = j - (Int_t)(4 * ai) - bw; w <= j - (Int_t)(4 * ai) + bw; w++){
1012  if (w >= 0 && w < ssize){
1013  b8 += working_space[ssize + w];
1014  men +=1;
1015  }
1016  }
1017  b8 = b8 / men;
1018  c8 = 0, men = 0;
1019  for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
1020  if (w >= 0 && w < ssize){
1021  c8 += working_space[ssize + w];
1022  men +=1;
1023  }
1024  }
1025  c8 = c8 / men;
1026  d8 = 0, men = 0;
1027  for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1028  if (w >= 0 && w < ssize){
1029  d8 += working_space[ssize + w];
1030  men +=1;
1031  }
1032  }
1033  d8 = d8 / men;
1034  e8 = 0, men = 0;
1035  for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1036  if (w >= 0 && w < ssize){
1037  e8 += working_space[ssize + w];
1038  men +=1;
1039  }
1040  }
1041  e8 = e8 / men;
1042  f8 = 0, men = 0;
1043  for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1044  if (w >= 0 && w < ssize){
1045  f8 += working_space[ssize + w];
1046  men +=1;
1047  }
1048  }
1049  f8 = f8 / men;
1050  g8 = 0, men = 0;
1051  for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1052  if (w >= 0 && w < ssize){
1053  g8 += working_space[ssize + w];
1054  men +=1;
1055  }
1056  }
1057  g8 = g8 / men;
1058  h8 = 0, men = 0;
1059  for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1060  if (w >= 0 && w < ssize){
1061  h8 += working_space[ssize + w];
1062  men +=1;
1063  }
1064  }
1065  h8 = h8 / men;
1066  i8 = 0, men = 0;
1067  for (w = j + (Int_t)(4 * ai) - bw; w <= j + (Int_t)(4 * ai) + bw; w++){
1068  if (w >= 0 && w < ssize){
1069  i8 += working_space[ssize + w];
1070  men +=1;
1071  }
1072  }
1073  i8 = i8 / men;
1074  b8 = ( -b8 + 8 * c8 - 28 * d8 + 56 * e8 - 56 * f8 - 28 * g8 + 8 * h8 - i8)/70;
1075  if (b < b8)
1076  b = b8;
1077  if (b < b6)
1078  b = b6;
1079  if (b < b4)
1080  b = b4;
1081  if (b < a)
1082  av = b;
1083  working_space[j]=av;
1084  }
1085  }
1086  for (j = i; j < ssize - i; j++)
1087  working_space[ssize + j] = working_space[j];
1088  if (direction == kBackIncreasingWindow)
1089  i += 1;
1090  else if(direction == kBackDecreasingWindow)
1091  i -= 1;
1092  }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1093  }
1094 
1095  if (compton == kTRUE) {
1096  for (i = 0, b2 = 0; i < ssize; i++){
1097  b1 = b2;
1098  a = working_space[i], b = spectrum[i];
1099  j = i;
1100  if (TMath::Abs(a - b) >= 1) {
1101  b1 = i - 1;
1102  if (b1 < 0)
1103  b1 = 0;
1104  yb1 = working_space[b1];
1105  for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){
1106  a = working_space[b2], b = spectrum[b2];
1107  c = c + b - yb1;
1108  if (TMath::Abs(a - b) < 1) {
1109  priz = 1;
1110  yb2 = b;
1111  }
1112  }
1113  if (b2 == ssize)
1114  b2 -= 1;
1115  yb2 = working_space[b2];
1116  if (yb1 <= yb2){
1117  for (j = b1, c = 0; j <= b2; j++){
1118  b = spectrum[j];
1119  c = c + b - yb1;
1120  }
1121  if (c > 1){
1122  c = (yb2 - yb1) / c;
1123  for (j = b1, d = 0; j <= b2 && j < ssize; j++){
1124  b = spectrum[j];
1125  d = d + b - yb1;
1126  a = c * d + yb1;
1127  working_space[ssize + j] = a;
1128  }
1129  }
1130  }
1131 
1132  else{
1133  for (j = b2, c = 0; j >= b1; j--){
1134  b = spectrum[j];
1135  c = c + b - yb2;
1136  }
1137  if (c > 1){
1138  c = (yb1 - yb2) / c;
1139  for (j = b2, d = 0;j >= b1 && j >= 0; j--){
1140  b = spectrum[j];
1141  d = d + b - yb2;
1142  a = c * d + yb2;
1143  working_space[ssize + j] = a;
1144  }
1145  }
1146  }
1147  i=b2;
1148  }
1149  }
1150  }
1151 
1152  for (j = 0; j < ssize; j++)
1153  spectrum[j] = working_space[ssize + j];
1154  delete[]working_space;
1155  return 0;
1156 }
1157 
1158 ////////////////////////////////////////////////////////////////////////////////
1159 /// One-dimensional markov spectrum smoothing function
1160 ///
1161 /// This function calculates smoothed spectrum from source spectrum based on
1162 /// Markov chain method. The result is placed in the array pointed by source
1163 /// pointer. On successful completion it returns 0. On error it returns pointer
1164 /// to the string describing error.
1165 ///
1166 /// #### Parameters:
1167 ///
1168 /// - source: pointer to the array of source spectrum
1169 /// - ssize: length of source array
1170 /// - averWindow: width of averaging smoothing window
1171 ///
1172 /// The goal of this function is the suppression of the statistical fluctuations.
1173 /// The algorithm is based on discrete Markov chain, which has very simple
1174 /// invariant distribution:
1175 ///
1176 /// \f[
1177 /// 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
1178 /// \f]
1179 /// \f$ U_1\f$ being defined from the normalization condition
1180 /// \f$ \sum_{i=1}^{n} U_i=1\f$. \f$ n \f$ is the length of the smoothed spectrum and
1181 /// \f[
1182 /// 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]
1183 /// \f]
1184 ///
1185 /// #### Reference:
1186 ///
1187 /// 1. Z.K. Silagadze, A new algorithm for automatic photopeak searches.
1188 /// NIM A 376 (1996), 451.
1189 ///
1190 /// ### Example 14 - script Smoothing.C
1191 ///
1192 /// Begin_Macro(source)
1193 /// ../../../tutorials/spectrum/Smoothing.C
1194 /// End_Macro
1195 
1196 const char* TSpectrum::SmoothMarkov(Double_t *source, int ssize, int averWindow)
1197 {
1198  int xmin, xmax, i, l;
1199  Double_t a, b, maxch;
1200  Double_t nom, nip, nim, sp, sm, area = 0;
1201  if(averWindow <= 0)
1202  return "Averaging Window must be positive";
1203  Double_t *working_space = new Double_t[ssize];
1204  xmin = 0,xmax = ssize - 1;
1205  for(i = 0, maxch = 0; i < ssize; i++){
1206  working_space[i]=0;
1207  if(maxch < source[i])
1208  maxch = source[i];
1209 
1210  area += source[i];
1211  }
1212  if(maxch == 0) {
1213  delete [] working_space;
1214  return 0 ;
1215  }
1216 
1217  nom = 1;
1218  working_space[xmin] = 1;
1219  for(i = xmin; i < xmax; i++){
1220  nip = source[i] / maxch;
1221  nim = source[i + 1] / maxch;
1222  sp = 0,sm = 0;
1223  for(l = 1; l <= averWindow; l++){
1224  if((i + l) > xmax)
1225  a = source[xmax] / maxch;
1226 
1227  else
1228  a = source[i + l] / maxch;
1229  b = a - nip;
1230  if(a + nip <= 0)
1231  a = 1;
1232 
1233  else
1234  a = TMath::Sqrt(a + nip);
1235  b = b / a;
1236  b = TMath::Exp(b);
1237  sp = sp + b;
1238  if((i - l + 1) < xmin)
1239  a = source[xmin] / maxch;
1240 
1241  else
1242  a = source[i - l + 1] / maxch;
1243  b = a - nim;
1244  if(a + nim <= 0)
1245  a = 1;
1246  else
1247  a = TMath::Sqrt(a + nim);
1248  b = b / a;
1249  b = TMath::Exp(b);
1250  sm = sm + b;
1251  }
1252  a = sp / sm;
1253  a = working_space[i + 1] = working_space[i] * a;
1254  nom = nom + a;
1255  }
1256  for(i = xmin; i <= xmax; i++){
1257  working_space[i] = working_space[i] / nom;
1258  }
1259  for(i = 0; i < ssize; i++)
1260  source[i] = working_space[i] * area;
1261  delete[]working_space;
1262  return 0;
1263 }
1264 
1265 ////////////////////////////////////////////////////////////////////////////////
1266 /// One-dimensional deconvolution function
1267 ///
1268 /// This function calculates deconvolution from source spectrum according to
1269 /// response spectrum using Gold deconvolution algorithm. The result is placed
1270 /// in the vector pointed by source pointer. On successful completion it
1271 /// returns 0. On error it returns pointer to the string describing error. If
1272 /// desired after every numberIterations one can apply boosting operation
1273 /// (exponential function with exponent given by boost coefficient) and repeat
1274 /// it numberRepetitions times.
1275 ///
1276 /// #### Parameters:
1277 ///
1278 /// - source: pointer to the vector of source spectrum
1279 /// - response: pointer to the vector of response spectrum
1280 /// - ssize: length of source and response spectra
1281 /// - numberIterations, for details we refer to the reference given below
1282 /// - numberRepetitions, for repeated boosted deconvolution
1283 /// - boost, boosting coefficient
1284 ///
1285 /// The goal of this function is the improvement of the resolution in spectra,
1286 /// decomposition of multiplets. The mathematical formulation of
1287 /// the convolution system is:
1288 ///
1289 /// \f[
1290 /// y(i) = \sum_{k=0}^{N-1} h(i-k)x(k), i=0,1,2,...,N-1
1291 /// \f]
1292 ///
1293 /// where h(i) is the impulse response function, x, y are input and output
1294 /// vectors, respectively, N is the length of x and h vectors. In matrix form
1295 /// we have:
1296 /**
1297  \f[
1298  \begin{bmatrix}
1299  y(0) \\
1300  y(1) \\
1301  \dots \\
1302  y(2N-2) \\
1303  y(2N-1)
1304  \end{bmatrix}
1305  =
1306  \begin{bmatrix}
1307  h(0) & 0 & 0 & \dots & 0 \\
1308  h(1) & h(0) & 0 & \dots & \dots \\
1309  \dots & h(1) & h(0) & \dots & \dots \\
1310  \dots & \dots & h(1) & \dots & \dots \\
1311  \dots & \dots & \dots & \dots & \dots \\
1312  h(N-1) & \dots & \dots &\dots & 0 \\
1313  0 & h(N-1) & \dots & \dots & h(0) \\
1314  0 & 0 & h(N-1) & \dots & h(1) \\
1315  \dots & \dots & \dots & \dots & \dots \\
1316  0 & 0 & 0 & \dots & h(N-1)
1317  \end{bmatrix}
1318  \begin{bmatrix}
1319  x(0) \\
1320  x(1) \\
1321  \dots \\
1322  x(N-1)
1323  \end{bmatrix}
1324  \f]
1325 */
1326 /// Let us assume that we know the response and the output vector (spectrum) of
1327 /// the above given system. The deconvolution represents solution of the
1328 /// overdetermined system of linear equations, i.e., the calculation of the
1329 /// vector x. From numerical stability point of view the operation of
1330 /// deconvolution is extremely critical (ill-posed problem) as well as time
1331 /// consuming operation. The Gold deconvolution algorithm proves to work very
1332 /// well, other methods (Fourier, VanCittert etc) oscillate. It is suitable to
1333 /// process positive definite data (e.g. histograms).
1334 ///
1335 /// #### Gold deconvolution algorithm:
1336 /**
1337  \f[
1338  y = Hx \\
1339  H^{T}y = H^{T}Hx \\
1340  y^{'} = H^{'}x \\
1341  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 \\
1342  k = 1,2,3,...,L
1343  x^{0} = [1,1, ..., 1]^T
1344  \f]
1345 */
1346 /// Where L is given number of iterations (numberIterations parameter).
1347 ///
1348 /// #### Boosted deconvolution:
1349 ///
1350 /// 1. Set the initial solution:
1351 /// \f$ x^{(0)} = [1,1,...,1]^{T} \f$
1352 /// 2. Set required number of repetitions R and iterations L.
1353 /// 3. Set r = 1.
1354 /// 4. Using Gold deconvolution algorithm for k=1,2,...,L find
1355 /// \f$ x^{(L)} \f$
1356 /// 5. If r = R stop calculation, else
1357 ///
1358 /// 1. Apply boosting operation, i.e., set
1359 /// \f$ x^{(0)}(i) = [x^{(L)}(i)]^{p} \f$
1360 /// i=0,1,...N-1 and p is boosting coefficient >0.
1361 /// 2. r = r + 1
1362 /// 3. continue in 4.
1363 ///
1364 /// #### References:
1365 ///
1366 /// 1. Gold R., ANL-6984, Argonne National Laboratories, Argonne Ill, 1964.
1367 /// 2. Coote G.E., Iterative smoothing and deconvolution of one- and two-dimensional
1368 /// elemental distribution data, NIM B 130 (1997) 118.
1369 /// 3. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky,
1370 /// I. Turzo: Efficient one- and two-dimensional Gold deconvolution and
1371 /// its application to gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
1372 /// 4. Morhac; M., Matouoek V., Kliman J., Efficient algorithm of multidimensional
1373 /// deconvolution and its application to nuclear data processing, Digital Signal
1374 /// Processing 13 (2003) 144.
1375 ///
1376 /// ### Example 8 - script Deconvolution.C :
1377 ///
1378 /// response function (usually peak) should be shifted left to the first
1379 /// non-zero channel (bin).
1380 ///
1381 /// \image html TSpectrum_Deconvolution2.jpg Principle how the response matrix is composed inside of the Deconvolution function.
1382 ///
1383 /// Begin_Macro(source)
1384 /// ../../../tutorials/spectrum/Deconvolution.C
1385 /// End_Macro
1386 ///
1387 /// ### Examples of Gold deconvolution method:
1388 ///
1389 /// First let us study the influence of the number of iterations on the
1390 /// deconvolved spectrum (Fig. 12).
1391 ///
1392 /// \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.
1393 ///
1394 /// For relatively narrow peaks in the above given example the Gold
1395 /// deconvolution method is able to decompose overlapping peaks practically to
1396 /// delta - functions. In the next example we have chosen a synthetic data
1397 /// (spectrum, 256 channels) consisting of 5 very closely positioned, relatively
1398 /// wide peaks (sigma =5), with added noise (Fig. 13). Thin lines represent
1399 /// pure Gaussians (see Table 1); thick line is a resulting spectrum with
1400 /// additive noise (10% of the amplitude of small peaks).
1401 ///
1402 /// \image html TSpectrum_Deconvolution_wide2.jpg Fig. 13 Testing example of synthetic spectrum composed of 5 Gaussians with added noise.
1403 ///
1404 /// | Peak # | Position | Height | Area |
1405 /// |----------|----------|--------|--------|
1406 /// | 1 | 50 | 500 | 10159 |
1407 /// | 2 | 70 | 3000 | 60957 |
1408 /// | 3 | 80 | 1000 | 20319 |
1409 /// | 4 | 100 | 5000 | 101596 |
1410 /// | 5 | 110 | 500 | 10159 |
1411 ///
1412 /// Table 1 Positions, heights and areas of peaks in the spectrum shown in Fig. 13.
1413 ///
1414 /// In ideal case, we should obtain the result given in Fig. 14. The areas of
1415 /// the Gaussian components of the spectrum are concentrated completely to
1416 /// delta-functions. When solving the overdetermined system of linear equations
1417 /// with data from Fig. 13 in the sense of minimum least squares criterion
1418 /// without any regularisation we obtain the result with large oscillations
1419 /// (Fig. 15). From mathematical point of view, it is the optimal solution in
1420 /// the unconstrained space of independent variables. From physical point of
1421 /// view we are interested only in a meaningful solution. Therefore, we have to
1422 /// employ regularisation techniques (e.g. Gold deconvolution) and/or to
1423 /// confine the space of allowed solutions to subspace of positive solutions.
1424 ///
1425 /// \image html TSpectrum_Deconvolution_wide3.jpg Fig. 14 The same spectrum like in Fig. 13, outlined bars show the contents of present components (peaks).
1426 /// \image html TSpectrum_Deconvolution_wide4.jpg Fig. 15 Least squares solution of the system of linear equations without regularisation.
1427 ///
1428 /// ### Example 9 - script Deconvolution_wide.C
1429 ///
1430 /// When we employ Gold deconvolution algorithm we obtain the result given in
1431 /// Fig. 16. One can observe that the resulting spectrum is smooth. On the
1432 /// other hand the method is not able to decompose completely the peaks in the
1433 /// spectrum.
1434 ///
1435 /// Example of Gold deconvolution for closely positioned wide peaks. The original
1436 /// source spectrum is drawn with black color, the spectrum after the deconvolution
1437 /// (10000 iterations) with red color.
1438 ///
1439 /// Begin_Macro(source)
1440 /// ../../../tutorials/spectrum/Deconvolution_wide.C
1441 /// End_Macro
1442 ///
1443 /// ### Example 10 - script Deconvolution_wide_boost.C :
1444 ///
1445 /// Further let us employ boosting operation into deconvolution (Fig. 17).
1446 ///
1447 /// The original source spectrum is drawn with black color, the spectrum after
1448 /// the deconvolution with red color. Number of iterations = 200, number of
1449 /// repetitions = 50 and boosting coefficient = 1.2.
1450 ///
1451 /// One can observe that peaks are decomposed practically to delta functions.
1452 /// Number of peaks is correct, positions of big peaks as well as their areas
1453 /// are relatively well estimated. However there is a considerable error in
1454 /// the estimation of the position of small right hand peak.
1455 ///
1456 /// Begin_Macro(source)
1457 /// ../../../tutorials/spectrum/Deconvolution_wide_boost.C
1458 /// End_Macro
1459 
1460 const char *TSpectrum::Deconvolution(Double_t *source, const Double_t *response,
1461  int ssize, int numberIterations,
1462  int numberRepetitions, Double_t boost )
1463 {
1464  if (ssize <= 0)
1465  return "Wrong Parameters";
1466 
1467  if (numberRepetitions <= 0)
1468  return "Wrong Parameters";
1469 
1470  // working_space-pointer to the working vector
1471  // (its size must be 4*ssize of source spectrum)
1472  Double_t *working_space = new Double_t[4 * ssize];
1473  int i, j, k, lindex, posit, lh_gold, l, repet;
1474  Double_t lda, ldb, ldc, area, maximum;
1475  area = 0;
1476  lh_gold = -1;
1477  posit = 0;
1478  maximum = 0;
1479 
1480 //read response vector
1481  for (i = 0; i < ssize; i++) {
1482  lda = response[i];
1483  if (lda != 0)
1484  lh_gold = i + 1;
1485  working_space[i] = lda;
1486  area += lda;
1487  if (lda > maximum) {
1488  maximum = lda;
1489  posit = i;
1490  }
1491  }
1492  if (lh_gold == -1) {
1493  delete [] working_space;
1494  return "ZERO RESPONSE VECTOR";
1495  }
1496 
1497 //read source vector
1498  for (i = 0; i < ssize; i++)
1499  working_space[2 * ssize + i] = source[i];
1500 
1501 // create matrix at*a and vector at*y
1502  for (i = 0; i < ssize; i++){
1503  lda = 0;
1504  for (j = 0; j < ssize; j++){
1505  ldb = working_space[j];
1506  k = i + j;
1507  if (k < ssize){
1508  ldc = working_space[k];
1509  lda = lda + ldb * ldc;
1510  }
1511  }
1512  working_space[ssize + i] = lda;
1513  lda = 0;
1514  for (k = 0; k < ssize; k++){
1515  l = k - i;
1516  if (l >= 0){
1517  ldb = working_space[l];
1518  ldc = working_space[2 * ssize + k];
1519  lda = lda + ldb * ldc;
1520  }
1521  }
1522  working_space[3 * ssize + i]=lda;
1523  }
1524 
1525 // move vector at*y
1526  for (i = 0; i < ssize; i++){
1527  working_space[2 * ssize + i] = working_space[3 * ssize + i];
1528  }
1529 
1530 //initialization of resulting vector
1531  for (i = 0; i < ssize; i++)
1532  working_space[i] = 1;
1533 
1534  //**START OF ITERATIONS**
1535  for (repet = 0; repet < numberRepetitions; repet++) {
1536  if (repet != 0) {
1537  for (i = 0; i < ssize; i++)
1538  working_space[i] = TMath::Power(working_space[i], boost);
1539  }
1540  for (lindex = 0; lindex < numberIterations; lindex++) {
1541  for (i = 0; i < ssize; i++) {
1542  if (working_space[2 * ssize + i] > 0.000001
1543  && working_space[i] > 0.000001) {
1544  lda = 0;
1545  for (j = 0; j < lh_gold; j++) {
1546  ldb = working_space[j + ssize];
1547  if (j != 0){
1548  k = i + j;
1549  ldc = 0;
1550  if (k < ssize)
1551  ldc = working_space[k];
1552  k = i - j;
1553  if (k >= 0)
1554  ldc += working_space[k];
1555  }
1556 
1557  else
1558  ldc = working_space[i];
1559  lda = lda + ldb * ldc;
1560  }
1561  ldb = working_space[2 * ssize + i];
1562  if (lda != 0)
1563  lda = ldb / lda;
1564 
1565  else
1566  lda = 0;
1567  ldb = working_space[i];
1568  lda = lda * ldb;
1569  working_space[3 * ssize + i] = lda;
1570  }
1571  }
1572  for (i = 0; i < ssize; i++)
1573  working_space[i] = working_space[3 * ssize + i];
1574  }
1575  }
1576 
1577 //shift resulting spectrum
1578  for (i = 0; i < ssize; i++) {
1579  lda = working_space[i];
1580  j = i + posit;
1581  j = j % ssize;
1582  working_space[ssize + j] = lda;
1583  }
1584 
1585 //write back resulting spectrum
1586  for (i = 0; i < ssize; i++)
1587  source[i] = area * working_space[ssize + i];
1588  delete[]working_space;
1589  return 0;
1590 }
1591 
1592 
1593 ////////////////////////////////////////////////////////////////////////////////
1594 /// One-dimensional deconvolution function.
1595 ///
1596 /// This function calculates deconvolution from source spectrum according to
1597 /// response spectrum using Richardson-Lucy deconvolution algorithm. The result
1598 /// is placed in the vector pointed by source pointer. On successful completion
1599 /// it returns 0. On error it returns pointer to the string describing error.
1600 /// If desired after every numberIterations one can apply boosting operation
1601 /// (exponential function with exponent given by boost coefficient) and repeat
1602 /// it numberRepetitions times (see Gold deconvolution).
1603 ///
1604 /// #### Parameters:
1605 ///
1606 /// - source: pointer to the vector of source spectrum
1607 /// - response: pointer to the vector of response spectrum
1608 /// - ssize: length of source and response spectra
1609 /// - numberIterations, for details we refer to the reference given above
1610 /// - numberRepetitions, for repeated boosted deconvolution
1611 /// - boost, boosting coefficient
1612 ///
1613 /// ### Richardson-Lucy deconvolution algorithm:
1614 ///
1615 /// For discrete systems it has the form:
1616 /**
1617  \f[
1618  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)} \\
1619  i \in \left<0,M-1\right>
1620  \f]
1621 */
1622 /// for positive input data and response matrix this iterative method forces
1623 /// the deconvoluted spectra to be non-negative. The Richardson-Lucy
1624 /// iteration converges to the maximum likelihood solution for Poisson statistics
1625 /// in the data.
1626 ///
1627 /// #### References:
1628 ///
1629 /// 1. Abreu M.C. et al., A four-dimensional deconvolution method to correct NA38
1630 /// experimental data, NIM A 405 (1998) 139.
1631 /// 2. Lucy L.B., A.J. 79 (1974) 745.
1632 /// 3. Richardson W.H., J. Opt. Soc. Am. 62 (1972) 55.
1633 ///
1634 /// ### Examples of Richardson-Lucy deconvolution method:
1635 ///
1636 /// ### Example 11 - script DeconvolutionRL_wide.C :
1637 ///
1638 /// When we employ Richardson-Lucy deconvolution algorithm to our data from
1639 /// Fig. 13 we obtain the following result. One can observe improvements
1640 /// as compared to the result achieved by Gold deconvolution. Nevertheless it is
1641 /// unable to decompose the multiplet.
1642 ///
1643 /// Example of Richardson-Lucy deconvolution for closely positioned wide peaks.
1644 /// The original source spectrum is drawn with black color, the spectrum after
1645 /// the deconvolution (10000 iterations) with red color.
1646 ///
1647 /// Begin_Macro(source)
1648 /// ../../../tutorials/spectrum/DeconvolutionRL_wide.C
1649 /// End_Macro
1650 ///
1651 /// ### Example 12 - script DeconvolutionRL_wide_boost.C :
1652 ///
1653 /// Further let us employ boosting operation into deconvolution.
1654 ///
1655 /// The original source spectrum is drawn with black color, the spectrum after
1656 /// the deconvolution with red color. Number of iterations = 200, number of
1657 /// repetitions = 50 and boosting coefficient = 1.2.
1658 ///
1659 /// One can observe improvements in the estimation of peak positions as compared
1660 /// to the results achieved by Gold deconvolution.
1661 ///
1662 /// Begin_Macro(source)
1663 /// ../../../tutorials/spectrum/DeconvolutionRL_wide_boost.C
1664 /// End_Macro
1665 
1666 const char *TSpectrum::DeconvolutionRL(Double_t *source, const Double_t *response,
1667  int ssize, int numberIterations,
1668  int numberRepetitions, Double_t boost )
1669 {
1670  if (ssize <= 0)
1671  return "Wrong Parameters";
1672 
1673  if (numberRepetitions <= 0)
1674  return "Wrong Parameters";
1675 
1676  // working_space-pointer to the working vector
1677  // (its size must be 4*ssize of source spectrum)
1678  Double_t *working_space = new Double_t[4 * ssize];
1679  int i, j, k, lindex, posit, lh_gold, repet, kmin, kmax;
1680  Double_t lda, ldb, ldc, maximum;
1681  lh_gold = -1;
1682  posit = 0;
1683  maximum = 0;
1684 
1685 //read response vector
1686  for (i = 0; i < ssize; i++) {
1687  lda = response[i];
1688  if (lda != 0)
1689  lh_gold = i + 1;
1690  working_space[ssize + i] = lda;
1691  if (lda > maximum) {
1692  maximum = lda;
1693  posit = i;
1694  }
1695  }
1696  if (lh_gold == -1) {
1697  delete [] working_space;
1698  return "ZERO RESPONSE VECTOR";
1699  }
1700 
1701 //read source vector
1702  for (i = 0; i < ssize; i++)
1703  working_space[2 * ssize + i] = source[i];
1704 
1705 //initialization of resulting vector
1706  for (i = 0; i < ssize; i++){
1707  if (i <= ssize - lh_gold)
1708  working_space[i] = 1;
1709 
1710  else
1711  working_space[i] = 0;
1712 
1713  }
1714  //**START OF ITERATIONS**
1715  for (repet = 0; repet < numberRepetitions; repet++) {
1716  if (repet != 0) {
1717  for (i = 0; i < ssize; i++)
1718  working_space[i] = TMath::Power(working_space[i], boost);
1719  }
1720  for (lindex = 0; lindex < numberIterations; lindex++) {
1721  for (i = 0; i <= ssize - lh_gold; i++){
1722  lda = 0;
1723  if (working_space[i] > 0){//x[i]
1724  for (j = i; j < i + lh_gold; j++){
1725  ldb = working_space[2 * ssize + j];//y[j]
1726  if (j < ssize){
1727  if (ldb > 0){//y[j]
1728  kmax = j;
1729  if (kmax > lh_gold - 1)
1730  kmax = lh_gold - 1;
1731  kmin = j + lh_gold - ssize;
1732  if (kmin < 0)
1733  kmin = 0;
1734  ldc = 0;
1735  for (k = kmax; k >= kmin; k--){
1736  ldc += working_space[ssize + k] * working_space[j - k];//h[k]*x[j-k]
1737  }
1738  if (ldc > 0)
1739  ldb = ldb / ldc;
1740 
1741  else
1742  ldb = 0;
1743  }
1744  ldb = ldb * working_space[ssize + j - i];//y[j]*h[j-i]/suma(h[j][k]x[k])
1745  }
1746  lda += ldb;
1747  }
1748  lda = lda * working_space[i];
1749  }
1750  working_space[3 * ssize + i] = lda;
1751  }
1752  for (i = 0; i < ssize; i++)
1753  working_space[i] = working_space[3 * ssize + i];
1754  }
1755  }
1756 
1757 //shift resulting spectrum
1758  for (i = 0; i < ssize; i++) {
1759  lda = working_space[i];
1760  j = i + posit;
1761  j = j % ssize;
1762  working_space[ssize + j] = lda;
1763  }
1764 
1765 //write back resulting spectrum
1766  for (i = 0; i < ssize; i++)
1767  source[i] = working_space[ssize + i];
1768  delete[]working_space;
1769  return 0;
1770 }
1771 
1772 
1773 ////////////////////////////////////////////////////////////////////////////////
1774 /// One-dimensional unfolding function
1775 ///
1776 /// This function unfolds source spectrum according to response matrix columns.
1777 /// The result is placed in the vector pointed by source pointer.
1778 /// The coefficients of the resulting vector represent contents of the columns
1779 /// (weights) in the input vector. On successful completion it returns 0. On
1780 /// error it returns pointer to the string describing error. If desired after
1781 /// every numberIterations one can apply boosting operation (exponential
1782 /// function with exponent given by boost coefficient) and repeat it
1783 /// numberRepetitions times. For details we refer to [1].
1784 ///
1785 /// #### Parameters:
1786 ///
1787 /// - source: pointer to the vector of source spectrum
1788 /// - respMatrix: pointer to the matrix of response spectra
1789 /// - ssizex: length of source spectrum and # of rows of the response
1790 /// matrix. ssizex must be >= ssizey.
1791 /// - ssizey: length of destination coefficients and # of columns of the response
1792 /// matrix.
1793 /// - numberIterations: number of iterations
1794 /// - numberRepetitions: number of repetitions for boosted deconvolution.
1795 /// It must be greater or equal to one.
1796 /// - boost: boosting coefficient, applies only if numberRepetitions is
1797 /// greater than one.
1798 ///
1799 /// ### Unfolding:
1800 ///
1801 /// The goal is the decomposition of spectrum to a given set of component spectra.
1802 ///
1803 /// The mathematical formulation of the discrete linear system is:
1804 ///
1805 /// \f[
1806 /// y(i) = \sum_{k=0}^{N_y-1} h(i,k)x(k), i = 0,1,2,...,N_x-1
1807 /// \f]
1808 /**
1809  \f[
1810  \begin{bmatrix}
1811  y(0) \\
1812  y(1) \\
1813  \dots \\
1814  y(N_x-1)
1815  \end{bmatrix}
1816  =
1817  \begin{bmatrix}
1818  h(0,0) & h(0,1) & \dots & h(0,N_y-1) \\
1819  h(1,0) & h(1,1) & \dots & h(1,N_y-1) \\
1820  \dots \\
1821  h(N_x-1,0) & h(N_x-1,1) & \dots & h(N_x-1,N_y-1)
1822  \end{bmatrix}
1823  \begin{bmatrix}
1824  x(0) \\
1825  x(1) \\
1826  \dots \\
1827  x(N_y-1)
1828  \end{bmatrix}
1829  \f]
1830 */
1831 ///
1832 /// #### References:
1833 ///
1834 /// 1. Jandel M., Morhac; M., Kliman J., Krupa L., Matouoek
1835 /// V., Hamilton J. H., Ramaya A. V.:
1836 /// Decomposition of continuum gamma-ray spectra using synthesised response matrix.
1837 /// NIM A 516 (2004), 172-183.
1838 ///
1839 /// ### Example of unfolding:
1840 ///
1841 /// ### Example 13 - script Unfolding.C:
1842 ///
1843 /// \image html TSpectrum_Unfolding3.gif Fig. 20 Response matrix composed of neutron spectra of pure chemical elements.
1844 /// \image html TSpectrum_Unfolding2.jpg Fig. 21 Source neutron spectrum to be decomposed
1845 /// \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)
1846 ///
1847 /// #### Script:
1848 ///
1849 /// ~~~ {.cpp}
1850 /// // Example to illustrate unfolding function (class TSpectrum).
1851 /// // To execute this example, do
1852 /// // root > .x Unfolding.C
1853 ///
1854 /// void Unfolding() {
1855 /// Int_t i, j;
1856 /// Int_t nbinsx = 2048;
1857 /// Int_t nbinsy = 10;
1858 /// double xmin = 0;
1859 /// double xmax = nbinsx;
1860 /// double ymin = 0;
1861 /// double ymax = nbinsy;
1862 /// double *source = new double[nbinsx];
1863 /// double **response = new double *[nbinsy];
1864 /// for (i=0;i<nbinsy;i++) response[i] = new double[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("TSpectrum.root");
1869 /// h = (TH1F*) f->Get("decon_unf_in;1");
1870 /// TFile *fr = new TFile("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 = (TCanvas *)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,(const double**)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 
2579 {
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 }
c
#define c(i)
Definition: RSha256.hxx:101
l
auto * l
Definition: textangle.C:4
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
n
const Int_t n
Definition: legend1.C:16
TSpectrum::fNPeaks
Int_t fNPeaks
number of peaks found
Definition: TSpectrum.h:26
first
Definition: first.py:1
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:100
e
#define e(i)
Definition: RSha256.hxx:103
TSpectrum::StaticBackground
static TH1 * StaticBackground(const TH1 *hist, Int_t niter=20, Option_t *option="")
Static function, interface to TSpectrum::Background.
Definition: TSpectrum.cxx:2588
snprintf
#define snprintf
Definition: civetweb.c:1540
Option_t
const char Option_t
Definition: RtypesCore.h:66
dest
#define dest(otri, vertexptr)
Definition: triangle.c:1041
TSpectrum::fResolution
Double_t fResolution
NOT USED resolution of the neighboring peaks
Definition: TSpectrum.h:30
TSpectrum::kBackSmoothing3
@ kBackSmoothing3
Definition: TSpectrum.h:43
TList::FindObject
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:578
TList::Delete
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition: TList.cxx:470
TString::Data
const char * Data() const
Definition: TString.h:369
TSpectrum::kBackOrder8
@ kBackOrder8
Definition: TSpectrum.h:40
TSpectrum::SetDeconIterations
static void SetDeconIterations(Int_t n=3)
Static function: Set max number of decon iterations in deconvolution operation (see TSpectrum::Search...
Definition: TSpectrum.cxx:109
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
xmax
float xmax
Definition: THbookFile.cxx:95
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
TSpectrum::Search
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:267
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:691
TAxis::GetFirst
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:458
TSpectrum::SmoothMarkov
const char * SmoothMarkov(Double_t *source, Int_t ssize, Int_t averWindow)
One-dimensional markov spectrum smoothing function.
Definition: TSpectrum.cxx:1196
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:727
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
Int_t
int Int_t
Definition: RtypesCore.h:45
hbname
#define hbname
Definition: THbookFile.cxx:125
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
TH1::SetBinContent
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:9056
TList.h
TH1::GetListOfFunctions
TList * GetListOfFunctions() const
Definition: TH1.h:243
PEAK_WINDOW
#define PEAK_WINDOW
Definition: TSpectrum.cxx:45
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
TString
Basic string class.
Definition: TString.h:136
TSpectrum::DeconvolutionRL
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:1666
TSpectrum::kBackSmoothing9
@ kBackSmoothing9
Definition: TSpectrum.h:46
b
#define b(i)
Definition: RSha256.hxx:100
bool
TString::ReplaceAll
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:692
TSpectrum.h
TSpectrum::fgIterations
static Int_t fgIterations
Maximum number of decon iterations (default=3)
Definition: TSpectrum.h:33
TSpectrum::StaticSearch
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:2578
TSpectrum::Deconvolution
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:1460
TH1::GetDimension
virtual Int_t GetDimension() const
Definition: TH1.h:282
TSpectrum::kBackIncreasingWindow
@ kBackIncreasingWindow
Definition: TSpectrum.h:41
TSpectrum::kBackDecreasingWindow
@ kBackDecreasingWindow
Definition: TSpectrum.h:42
TSpectrum::~TSpectrum
virtual ~TSpectrum()
Destructor.
Definition: TSpectrum.cxx:88
TSpectrum::kBackOrder2
@ kBackOrder2
Definition: TSpectrum.h:37
TSpectrum::fHistogram
TH1 * fHistogram
resulting histogram
Definition: TSpectrum.h:31
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4993
xmin
float xmin
Definition: THbookFile.cxx:95
h
#define h(i)
Definition: RSha256.hxx:106
TSpectrum::fPositionY
Double_t * fPositionY
[fNPeaks] Y position of peaks
Definition: TSpectrum.h:29
size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
TH1::SetEntries
virtual void SetEntries(Double_t n)
Definition: TH1.h:385
TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
a
auto * a
Definition: textangle.C:12
TSpectrum::fMaxPeaks
Int_t fMaxPeaks
Maximum number of peaks to be found.
Definition: TSpectrum.h:25
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
TSpectrum::TSpectrum
TSpectrum()
Constructor.
Definition: TSpectrum.cxx:51
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:101
TSpectrum::SearchHighRes
Int_t SearchHighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
One-dimensional high-resolution peak search function.
Definition: TSpectrum.cxx:2127
TH1::GetBinCenter
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8975
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:735
TAxis::GetLast
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:469
TVirtualPad.h
TSpectrum::fPositionX
Double_t * fPositionX
[fNPeaks] X position of peaks
Definition: TSpectrum.h:28
TSpectrum::Background
virtual TH1 * Background(const TH1 *hist, Int_t niter=20, Option_t *option="")
One-dimensional background estimation function.
Definition: TSpectrum.cxx:153
TSpectrum::Print
virtual void Print(Option_t *option="") const
Print the array of positions.
Definition: TSpectrum.cxx:220
kRed
@ kRed
Definition: Rtypes.h:66
TObject::Warning
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:879
TSpectrum::Search1HighRes
Int_t Search1HighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
Old name of SearcHighRes introduced for back compatibility.
Definition: TSpectrum.cxx:2564
ROOT::Math::VectorUtil::boost
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only re...
Definition: VectorUtil.h:365
TSpectrum::kBackOrder6
@ kBackOrder6
Definition: TSpectrum.h:39
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
TSpectrum::fPosition
Double_t * fPosition
[fNPeaks] array of current peak positions
Definition: TSpectrum.h:27
Double_t
double Double_t
Definition: RtypesCore.h:59
TSpectrum::kBackSmoothing5
@ kBackSmoothing5
Definition: TSpectrum.h:44
TSpectrum::SetAverageWindow
static void SetAverageWindow(Int_t w=3)
Static function: Set average window of searched peaks (see TSpectrum::SearchHighRes).
Definition: TSpectrum.cxx:100
TList::Remove
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:822
TSpectrum::SetResolution
void SetResolution(Double_t resolution=1)
NOT USED resolution: determines resolution of the neighbouring peaks default value is 1 correspond to...
Definition: TSpectrum.cxx:352
TList::Add
virtual void Add(TObject *obj)
Definition: TList.h:87
TSpectrum::kBackOrder4
@ kBackOrder4
Definition: TSpectrum.h:38
TH1::Reset
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH1.cxx:7070
TH1
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
TPolyMarker
A PolyMarker is defined by an array on N points in a 2-D space.
Definition: TPolyMarker.h:31
TSpectrum::kBackSmoothing13
@ kBackSmoothing13
Definition: TSpectrum.h:48
d
#define d(i)
Definition: RSha256.hxx:102
TSpectrum::kBackSmoothing7
@ kBackSmoothing7
Definition: TSpectrum.h:45
TPolyMarker.h
gPad
#define gPad
Definition: TVirtualPad.h:287
TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
TSpectrum::kBackSmoothing11
@ kBackSmoothing11
Definition: TSpectrum.h:47
TSpectrum::fgAverageWindow
static Int_t fgAverageWindow
Average window of searched peaks.
Definition: TSpectrum.h:32
TSpectrum::Unfolding
const char * Unfolding(Double_t *source, const Double_t **respMatrix, Int_t ssizex, Int_t ssizey, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional unfolding function.
Definition: TSpectrum.cxx:1890
TString::ToLower
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1150
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:320
TSpectrum::kBackSmoothing15
@ kBackSmoothing15
Definition: TSpectrum.h:49
TSpectrum
Advanced Spectra Processing.
Definition: TSpectrum.h:18
TH1.h
TMath.h
int
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3073