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