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