Logo ROOT   6.14/05
Reference Guide
TSpectrum2.cxx
Go to the documentation of this file.
1 // @(#)root/spectrum:$Id$
2 // Author: Miroslav Morhac 17/01/2006
3 
4 /** \class TSpectrum2
5  \ingroup Spectrum
6  \brief Advanced 2-dimensional spectra processing
7  \author Miroslav Morhac
8 
9  This class contains advanced spectra processing functions.
10 
11  - One-dimensional background estimation functions
12  - Two-dimensional background estimation functions
13  - One-dimensional smoothing functions
14  - Two-dimensional smoothing functions
15  - One-dimensional deconvolution functions
16  - Two-dimensional deconvolution functions
17  - One-dimensional peak search functions
18  - Two-dimensional peak search functions
19 
20  The algorithms in this class have been published in the following references:
21 
22  1. M.Morhac et al.: Background elimination methods for multidimensional coincidence gamma-ray spectra. Nuclear Instruments and Methods in Physics Research A 401 (1997) 113-132.
23  2. M.Morhac et al.: Efficient one- and two-dimensional Gold deconvolution and its application to gamma-ray spectra decomposition. Nuclear Instruments and Methods in Physics Research A 401 (1997) 385-408.
24  3. M.Morhac et al.: Identification of peaks in multidimensional coincidence gamma-ray spectra. Nuclear Instruments and Methods in Research Physics A 443(2000), 108-125.
25 
26  These NIM papers are also available as doc or ps files from:
27 
28  - [SpectrumDec.ps.gz](ftp://root.cern.ch/root/SpectrumDec.ps.gz)
29  - [SpectrumSrc.ps.gz](ftp://root.cern.ch/root/SpectrumSrc.ps.gz)
30  - [SpectrumBck.ps.gz](ftp://root.cern.ch/root/SpectrumBck.ps.gz)
31 
32  See also the
33  [online documentation](https://root.cern.ch/guides/tspectrum-manual) and
34  [tutorials](https://root.cern.ch/doc/master/group__tutorial__spectrum.html).
35 
36  All the figures in this page were prepared using the DaqProVis
37  system, Data Acquisition, Processing and Visualization system,
38  developed at the Institute of Physics, Slovak Academy of Sciences, Bratislava,
39  Slovakia.
40 */
41 
42 #include "TSpectrum2.h"
43 #include "TPolyMarker.h"
44 #include "TList.h"
45 #include "TH1.h"
46 #include "TMath.h"
47 #define PEAK_WINDOW 1024
48 
51 
53 
54 ////////////////////////////////////////////////////////////////////////////////
55 /// Constructor.
56 
57 TSpectrum2::TSpectrum2() :TNamed("Spectrum", "Miroslav Morhac peak finder")
58 {
59  Int_t n = 100;
60  fMaxPeaks = n;
61  fPosition = new Double_t[n];
62  fPositionX = new Double_t[n];
63  fPositionY = new Double_t[n];
64  fResolution = 1;
65  fHistogram = 0;
66  fNPeaks = 0;
67 }
68 
69 ////////////////////////////////////////////////////////////////////////////////
70 /// - maxpositions: maximum number of peaks
71 /// - resolution: *NOT USED* determines resolution of the neighbouring peaks
72 /// default value is 1 correspond to 3 sigma distance
73 /// between peaks. Higher values allow higher resolution
74 /// (smaller distance between peaks.
75 /// May be set later through SetResolution.
76 
77 TSpectrum2::TSpectrum2(Int_t maxpositions, Double_t resolution) :TNamed("Spectrum", "Miroslav Morhac peak finder")
78 {
79  Int_t n = maxpositions;
80  fMaxPeaks = n;
81  fPosition = new Double_t[n];
82  fPositionX = new Double_t[n];
83  fPositionY = new Double_t[n];
84  fHistogram = 0;
85  fNPeaks = 0;
86  SetResolution(resolution);
87 }
88 
89 ////////////////////////////////////////////////////////////////////////////////
90 /// Destructor.
91 
93 {
94  delete [] fPosition;
95  delete [] fPositionX;
96  delete [] fPositionY;
97  delete fHistogram;
98 }
99 
100 ////////////////////////////////////////////////////////////////////////////////
101 /// static function: Set average window of searched peaks
102 /// see TSpectrum2::SearchHighRes
103 
105 {
106  fgAverageWindow = w;
107 }
108 
109 ////////////////////////////////////////////////////////////////////////////////
110 /// static function: Set max number of decon iterations in deconvolution operation
111 /// see TSpectrum2::SearchHighRes
112 
114 {
115  fgIterations = n;
116 }
117 
118 ////////////////////////////////////////////////////////////////////////////////
119 /// This function calculates the background spectrum in the input histogram h.
120 /// The background is returned as a histogram.
121 ///
122 /// Function parameters:
123 /// - h: input 2-d histogram
124 /// - numberIterations, (default value = 20)
125 /// Increasing numberIterations make the result smoother and lower.
126 /// - option: may contain one of the following options
127 /// - to set the direction parameter
128 /// "BackIncreasingWindow". By default the direction is BackDecreasingWindow
129 /// - filterOrder-order of clipping filter. Possible values:
130 /// - "BackOrder2" (default)
131 /// - "BackOrder4"
132 /// - "BackOrder6"
133 /// - "BackOrder8"
134 /// - "nosmoothing"- if selected, the background is not smoothed
135 /// By default the background is smoothed.
136 /// - smoothWindow-width of smoothing window. Possible values:
137 /// - "BackSmoothing3" (default)
138 /// - "BackSmoothing5"
139 /// - "BackSmoothing7"
140 /// - "BackSmoothing9"
141 /// - "BackSmoothing11"
142 /// - "BackSmoothing13"
143 /// - "BackSmoothing15"
144 /// - "Compton" if selected the estimation of Compton edge
145 /// will be included.
146 /// - "same" : if this option is specified, the resulting background
147 /// histogram is superimposed on the picture in the current pad.
148 ///
149 /// NOTE that the background is only evaluated in the current range of h.
150 /// ie, if h has a bin range (set via h->GetXaxis()->SetRange(binmin,binmax),
151 /// the returned histogram will be created with the same number of bins
152 /// as the input histogram h, but only bins from binmin to binmax will be filled
153 /// with the estimated background.
154 
155 TH1 *TSpectrum2::Background(const TH1 * h, Int_t number_of_iterations,
156  Option_t * option)
157 {
158  Error("Background","function not yet implemented: h=%s, iter=%d, option=%sn"
159  , h->GetName(), number_of_iterations, option);
160  return 0;
161 }
162 
163 ////////////////////////////////////////////////////////////////////////////////
164 /// Print the array of positions.
165 
167 {
168  printf("\nNumber of positions = %d\n",fNPeaks);
169  for (Int_t i=0;i<fNPeaks;i++) {
170  printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
171  }
172 }
173 
174 ////////////////////////////////////////////////////////////////////////////////
175 /// This function searches for peaks in source spectrum in hin
176 /// The number of found peaks and their positions are written into
177 /// the members fNpeaks and fPositionX.
178 /// The search is performed in the current histogram range.
179 ///
180 /// Function parameters:
181 /// - hin: pointer to the histogram of source spectrum
182 /// - sigma: sigma of searched peaks, for details we refer to manual
183 /// - threshold: (default=0.05) peaks with amplitude less than
184 /// threshold*highest_peak are discarded. 0<threshold<1
185 ///
186 /// By default, the background is removed before deconvolution.
187 /// Specify the option "nobackground" to not remove the background.
188 ///
189 /// By default the "Markov" chain algorithm is used.
190 /// Specify the option "noMarkov" to disable this algorithm
191 /// Note that by default the source spectrum is replaced by a new spectrum
192 ///
193 /// By default a polymarker object is created and added to the list of
194 /// functions of the histogram. The histogram is drawn with the specified
195 /// option and the polymarker object drawn on top of the histogram.
196 /// The polymarker coordinates correspond to the npeaks peaks found in
197 /// the histogram.
198 /// A pointer to the polymarker object can be retrieved later via:
199 /// ~~~ {.cpp}
200 /// TList *functions = hin->GetListOfFunctions();
201 /// TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker")
202 /// ~~~
203 /// Specify the option "goff" to disable the storage and drawing of the
204 /// polymarker.
205 
207  Option_t * option, Double_t threshold)
208 {
209  if (hin == 0)
210  return 0;
211  Int_t dimension = hin->GetDimension();
212  if (dimension != 2) {
213  Error("Search", "Must be a 2-d histogram");
214  return 0;
215  }
216 
217  TString opt = option;
218  opt.ToLower();
220  if (opt.Contains("nobackground")) {
221  background = kFALSE;
222  opt.ReplaceAll("nobackground","");
223  }
224  Bool_t markov = kTRUE;
225  if (opt.Contains("nomarkov")) {
226  markov = kFALSE;
227  opt.ReplaceAll("nomarkov","");
228  }
229 
230  Int_t sizex = hin->GetXaxis()->GetNbins();
231  Int_t sizey = hin->GetYaxis()->GetNbins();
232  Int_t i, j, binx,biny, npeaks;
233  Double_t ** source = new Double_t*[sizex];
234  Double_t ** dest = new Double_t*[sizex];
235  for (i = 0; i < sizex; i++) {
236  source[i] = new Double_t[sizey];
237  dest[i] = new Double_t[sizey];
238  for (j = 0; j < sizey; j++) {
239  source[i][j] = hin->GetBinContent(i + 1, j + 1);
240  }
241  }
242  //npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, kTRUE, 3, kTRUE, 10);
243  //the smoothing option is used for 1-d but not for 2-d histograms
244  npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, background, fgIterations, markov, fgAverageWindow);
245 
246  //The logic in the loop should be improved to use the fact
247  //that fPositionX,Y give a precise position inside a bin.
248  //The current algorithm takes the center of the bin only.
249  for (i = 0; i < npeaks; i++) {
250  binx = 1 + Int_t(fPositionX[i] + 0.5);
251  biny = 1 + Int_t(fPositionY[i] + 0.5);
252  fPositionX[i] = hin->GetXaxis()->GetBinCenter(binx);
253  fPositionY[i] = hin->GetYaxis()->GetBinCenter(biny);
254  }
255  for (i = 0; i < sizex; i++) {
256  delete [] source[i];
257  delete [] dest[i];
258  }
259  delete [] source;
260  delete [] dest;
261 
262  if (opt.Contains("goff"))
263  return npeaks;
264  if (!npeaks) return 0;
265  TPolyMarker * pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
266  if (pm) {
267  hin->GetListOfFunctions()->Remove(pm);
268  delete pm;
269  }
270  pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
271  hin->GetListOfFunctions()->Add(pm);
272  pm->SetMarkerStyle(23);
273  pm->SetMarkerColor(kRed);
274  pm->SetMarkerSize(1.3);
275  ((TH1*)hin)->Draw(option);
276  return npeaks;
277 }
278 
279 ////////////////////////////////////////////////////////////////////////////////
280 /// *NOT USED*
281 /// resolution: determines resolution of the neighboring peaks
282 /// default value is 1 correspond to 3 sigma distance
283 /// between peaks. Higher values allow higher resolution
284 /// (smaller distance between peaks.
285 /// May be set later through SetResolution.
286 
288 {
289  if (resolution > 1)
290  fResolution = resolution;
291  else
292  fResolution = 1;
293 }
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 /// This function calculates background spectrum from source spectrum.
297 /// The result is placed to the array pointed by spectrum pointer.
298 ///
299 /// Function parameters:
300 /// - spectrum-pointer to the array of source spectrum
301 /// - ssizex-x length of spectrum
302 /// - ssizey-y length of spectrum
303 /// - numberIterationsX-maximal x width of clipping window
304 /// - numberIterationsY-maximal y width of clipping window
305 /// for details we refer to manual
306 /// - direction- direction of change of clipping window
307 /// - possible values=kBackIncreasingWindow
308 /// kBackDecreasingWindow
309 /// - filterType-determines the algorithm of the filtering
310 /// - possible values=kBackSuccessiveFiltering
311 /// kBackOneStepFiltering
312 ///
313 /// ### Background estimation
314 ///
315 /// Goal: Separation of useful information (peaks) from useless information (background)
316 ///
317 /// - method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping algorithm [1]
318 ///
319 /// - there exist two algorithms for the estimation of new value in the channel \f$ i_1,i_2\f$
320 ///
321 /// Algorithm based on Successive Comparisons:
322 ///
323 /// It is an extension of one-dimensional SNIP algorithm to another dimension. For
324 /// details we refer to [2].
325 ///
326 /// Algorithm based on One Step Filtering:
327 ///
328 /// New value in the estimated channel is calculated as:
329 /// \f[ a = \nu_{p-1}(i_1,i_2)\f]
330 /// \f[
331 /// b = \frac{\nu_{p-1}(i_1-p,i_2-p) - 2\nu_{p-1}(i_1-p,i_2) + \nu_{p-1}(i_1-p,i_2+p) - 2\nu_{p-1}(i_1,i_2-p)}{4} +
332 /// \f]
333 /// \f[
334 /// \frac{-2\nu_{p-1}(i_1,i_2+p) + \nu_{p-1}(i_1+p,i_2-p) - 2\nu_{p-1}(i_1+p,i_2) + \nu_{p-1}(i_1+p,i_2+p)}{4}
335 /// \f]
336 /// \f[ \nu_{p-1}(i_1,i_2) = min(a,b)\f]
337 /// where p = 1, 2, ..., number_of_iterations.
338 ///
339 /// #### References:
340 ///
341 /// [1] C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the
342 /// quantitative analysis of PIXE spectra in geoscience applications. NIM, B34 (1988), 396-402.
343 ///
344 /// [2] M. Morhac;, J. Kliman, V.
345 /// Matouoek, M. Veselsky, I. Turzo.:
346 /// Background elimination methods for multidimensional gamma-ray spectra. NIM,
347 /// A401 (1997) 113-132.
348 ///
349 /// ### Example 1 - script Back_gamma64.c
350 ///
351 /// \image html TSpectrum2_Background1.jpg Fig. 1 Original two-dimensional gamma-gamma-ray spectrum
352 /// \image html TSpectrum2_Background2.jpg Fig. 2 Background estimated from data from Fig. 1 using decreasing clipping window with widths 4, 4 and algorithm based on successive comparisons. The estimate includes not only continuously changing background but also one-dimensional ridges.
353 /// \image html TSpectrum2_Background3.jpg Fig. 3 Resulting peaks after subtraction of the estimated background (Fig. 2) from original two-dimensional gamma-gamma-ray spectrum (Fig. 1).
354 ///
355 /// #### Script:
356 ///
357 /// Example to illustrate the background estimator (class TSpectrum). To execute this example, do:
358 ///
359 /// `root > .x Back_gamma64.C`
360 ///
361 /// ~~~ {.cpp}
362 /// #include <TSpectrum>
363 /// void Back_gamma64() {
364 /// Int_t i, j;
365 /// Double_t nbinsx = 64;
366 /// Double_t nbinsy = 64;
367 /// Double_t xmin = 0;
368 /// Double_t xmax = (Double_t)nbinsx;
369 /// Double_t ymin = 0;
370 /// Double_t ymax = (Double_t)nbinsy;
371 /// Double_t ** source = new Double_t*[nbinsx];
372 /// for (i=0;i<nbinsx;i++)
373 /// source[i]=new Double_t[nbinsy];
374 /// TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
375 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
376 /// back=(TH2F*) f->Get("back1;1");
377 /// TCanvas *Background = new TCanvas("Background","Estimation of background with increasing window",10,10,1000,700);
378 /// TSpectrum *s = new TSpectrum();
379 /// for (i = 0; i < nbinsx; i++){
380 /// for (j = 0; j < nbinsy; j++){
381 /// source[i][j] = back->GetBinContent(i + 1,j + 1);
382 /// }
383 /// }
384 /// s->Background(source,nbinsx,nbinsy,4,4,kBackDecreasingWindow,kBackSuccessiveFiltering);
385 /// for (i = 0; i < nbinsx; i++){
386 /// for (j = 0; j < nbinsy; j++)
387 /// back->SetBinContent(i + 1,j + 1, source[i][j]);
388 /// }
389 /// back->Draw("SURF");
390 /// }
391 /// ~~~
392 ///
393 /// ### Example 2- script Back_gamma256.c
394 ///
395 /// \image html TSpectrum2_Background4.jpg Fig. 4 Original two-dimensional gamma-gamma-ray spectrum 256x256 channels
396 /// \image html TSpectrum2_Background5.jpg Fig. 5 Peaks after subtraction of the estimated background (increasing clipping window with widths 8, 8 and algorithm based on successive comparisons) from original two-dimensional gamma-gamma-ray spectrum (Fig. 4).
397 ///
398 /// #### Script:
399 ///
400 // Example to illustrate the background estimator (class TSpectrum).
401 /// To execute this example, do
402 ///
403 /// `root > .x Back_gamma256.C`
404 ///
405 /// ~~~ {.cpp}
406 /// #include <TSpectrum>
407 /// void Back_gamma256() {
408 /// Int_t i, j;
409 /// Double_t nbinsx = 64;
410 /// Double_t nbinsy = 64;
411 /// Double_t xmin = 0;
412 /// Double_t xmax = (Double_t)nbinsx;
413 /// Double_t ymin = 0;
414 /// Double_t ymax = (Double_t)nbinsy;
415 /// Double_t** source = new Double_t*[nbinsx];
416 /// for (i=0;i<nbinsx;i++)
417 /// source[i]=new Double_t[nbinsy];
418 /// TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
419 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
420 /// back=(TH2F*) f->Get("back2;1");
421 /// TCanvas *Background = new TCanvas("Background","Estimation of background with increasing window",10,10,1000,700);
422 /// TSpectrum *s = new TSpectrum();
423 /// for (i = 0; i < nbinsx; i++){
424 /// for (j = 0; j < nbinsy; j++){
425 /// source[i][j] = back->GetBinContent(i + 1,j + 1);
426 /// }
427 /// }
428 /// s->Background(source,nbinsx,nbinsy,8,8,kBackIncreasingWindow,kBackSuccessiveFiltering);
429 /// for (i = 0; i < nbinsx; i++){
430 /// for (j = 0; j < nbinsy; j++)
431 /// back->SetBinContent(i + 1,j + 1, source[i][j]);
432 /// }
433 /// back->Draw("SURF");
434 /// }
435 /// ~~~
436 ///
437 /// Example 3- script Back_synt256.c
438 ///
439 /// \image html TSpectrum2_Background6.jpg Fig. 6 Original two-dimensional synthetic spectrum 256x256 channels
440 /// \image html TSpectrum2_Background7.jpg Fig. 7 Peaks after subtraction of the estimated background (increasing clipping window with widths 8, 8 and algorithm based on successive comparisons) from original two-dimensional gamma-gamma-ray spectrum (Fig. 6). One can observe artefacts (ridges) around the peaks due to the employed algorithm.
441 /// \image html TSpectrum2_Background8.jpg Fig. 8 Peaks after subtraction of the estimated background (increasing clipping window with widths 8, 8 and algorithm based on one step filtering) from original two-dimensional gamma-gamma-ray spectrum (Fig. 6). The artefacts from the above given Fig. 7 disappeared.
442 ///
443 /// #### Script:
444 ///
445 /// Example to illustrate the background estimator (class TSpectrum).
446 /// To execute this example, do
447 ///
448 /// `root > .x Back_synt256.C`
449 ///
450 /// ~~~ {.cpp}
451 /// #include <TSpectrum>
452 /// void Back_synt256() {
453 /// Int_t i, j;
454 /// Double_t nbinsx = 64;
455 /// Double_t nbinsy = 64;
456 /// Double_t xmin = 0;
457 /// Double_t xmax = (Double_t)nbinsx;
458 /// Double_t ymin = 0;
459 /// Double_t ymax = (Double_t)nbinsy;
460 /// Double_t** source = new Double_t*[nbinsx];
461 /// for (i=0;i<nbinsx;i++)
462 /// source[i]=new Double_t[nbinsy];
463 /// TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
464 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
465 /// back=(TH2F*) f->Get("back3;1");
466 /// TCanvas *Background = new TCanvas("Background","Estimation of background with increasing window",10,10,1000,700);
467 /// TSpectrum *s = new TSpectrum();
468 /// for (i = 0; i < nbinsx; i++){
469 /// for (j = 0; j < nbinsy; j++){
470 /// source[i][j] = back->GetBinContent(i + 1,j + 1);
471 /// }
472 /// }
473 /// s->Background(source,nbinsx,nbinsy,8,8,
474 /// kBackIncreasingWindow,kBackSuccessiveFiltering);//kBackOneStepFiltering
475 /// for (i = 0; i < nbinsx; i++){
476 /// for (j = 0; j < nbinsy; j++)
477 /// back->SetBinContent(i + 1,j + 1, source[i][j]);
478 /// }
479 /// back->Draw("SURF");
480 /// }
481 /// ~~~
482 
483 const char *TSpectrum2::Background(Double_t **spectrum,
484  Int_t ssizex, Int_t ssizey,
485  Int_t numberIterationsX,
486  Int_t numberIterationsY,
487  Int_t direction,
488  Int_t filterType)
489 {
490  Int_t i, x, y, sampling, r1, r2;
491  Double_t a, b, p1, p2, p3, p4, s1, s2, s3, s4;
492  if (ssizex <= 0 || ssizey <= 0)
493  return "Wrong parameters";
494  if (numberIterationsX < 1 || numberIterationsY < 1)
495  return "Width of Clipping Window Must Be Positive";
496  if (ssizex < 2 * numberIterationsX + 1
497  || ssizey < 2 * numberIterationsY + 1)
498  return ("Too Large Clipping Window");
499  Double_t **working_space = new Double_t*[ssizex];
500  for (i = 0; i < ssizex; i++)
501  working_space[i] = new Double_t[ssizey];
502  sampling =
503  (Int_t) TMath::Max(numberIterationsX, numberIterationsY);
504  if (direction == kBackIncreasingWindow) {
505  if (filterType == kBackSuccessiveFiltering) {
506  for (i = 1; i <= sampling; i++) {
507  r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
508  (Int_t) TMath::Min(i, numberIterationsY);
509  for (y = r2; y < ssizey - r2; y++) {
510  for (x = r1; x < ssizex - r1; x++) {
511  a = spectrum[x][y];
512  p1 = spectrum[x - r1][y - r2];
513  p2 = spectrum[x - r1][y + r2];
514  p3 = spectrum[x + r1][y - r2];
515  p4 = spectrum[x + r1][y + r2];
516  s1 = spectrum[x][y - r2];
517  s2 = spectrum[x - r1][y];
518  s3 = spectrum[x + r1][y];
519  s4 = spectrum[x][y + r2];
520  b = (p1 + p2) / 2.0;
521  if (b > s2)
522  s2 = b;
523  b = (p1 + p3) / 2.0;
524  if (b > s1)
525  s1 = b;
526  b = (p2 + p4) / 2.0;
527  if (b > s4)
528  s4 = b;
529  b = (p3 + p4) / 2.0;
530  if (b > s3)
531  s3 = b;
532  s1 = s1 - (p1 + p3) / 2.0;
533  s2 = s2 - (p1 + p2) / 2.0;
534  s3 = s3 - (p3 + p4) / 2.0;
535  s4 = s4 - (p2 + p4) / 2.0;
536  b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
537  p3 +
538  p4) / 4.0;
539  if (b < a && b > 0)
540  a = b;
541  working_space[x][y] = a;
542  }
543  }
544  for (y = r2; y < ssizey - r2; y++) {
545  for (x = r1; x < ssizex - r1; x++) {
546  spectrum[x][y] = working_space[x][y];
547  }
548  }
549  }
550  }
551 
552  else if (filterType == kBackOneStepFiltering) {
553  for (i = 1; i <= sampling; i++) {
554  r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
555  (Int_t) TMath::Min(i, numberIterationsY);
556  for (y = r2; y < ssizey - r2; y++) {
557  for (x = r1; x < ssizex - r1; x++) {
558  a = spectrum[x][y];
559  b = -(spectrum[x - r1][y - r2] +
560  spectrum[x - r1][y + r2] + spectrum[x + r1][y -
561  r2]
562  + spectrum[x + r1][y + r2]) / 4 +
563  (spectrum[x][y - r2] + spectrum[x - r1][y] +
564  spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
565  if (b < a && b > 0)
566  a = b;
567  working_space[x][y] = a;
568  }
569  }
570  for (y = i; y < ssizey - i; y++) {
571  for (x = i; x < ssizex - i; x++) {
572  spectrum[x][y] = working_space[x][y];
573  }
574  }
575  }
576  }
577  }
578 
579  else if (direction == kBackDecreasingWindow) {
580  if (filterType == kBackSuccessiveFiltering) {
581  for (i = sampling; i >= 1; i--) {
582  r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
583  (Int_t) TMath::Min(i, numberIterationsY);
584  for (y = r2; y < ssizey - r2; y++) {
585  for (x = r1; x < ssizex - r1; x++) {
586  a = spectrum[x][y];
587  p1 = spectrum[x - r1][y - r2];
588  p2 = spectrum[x - r1][y + r2];
589  p3 = spectrum[x + r1][y - r2];
590  p4 = spectrum[x + r1][y + r2];
591  s1 = spectrum[x][y - r2];
592  s2 = spectrum[x - r1][y];
593  s3 = spectrum[x + r1][y];
594  s4 = spectrum[x][y + r2];
595  b = (p1 + p2) / 2.0;
596  if (b > s2)
597  s2 = b;
598  b = (p1 + p3) / 2.0;
599  if (b > s1)
600  s1 = b;
601  b = (p2 + p4) / 2.0;
602  if (b > s4)
603  s4 = b;
604  b = (p3 + p4) / 2.0;
605  if (b > s3)
606  s3 = b;
607  s1 = s1 - (p1 + p3) / 2.0;
608  s2 = s2 - (p1 + p2) / 2.0;
609  s3 = s3 - (p3 + p4) / 2.0;
610  s4 = s4 - (p2 + p4) / 2.0;
611  b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
612  p3 +
613  p4) / 4.0;
614  if (b < a && b > 0)
615  a = b;
616  working_space[x][y] = a;
617  }
618  }
619  for (y = r2; y < ssizey - r2; y++) {
620  for (x = r1; x < ssizex - r1; x++) {
621  spectrum[x][y] = working_space[x][y];
622  }
623  }
624  }
625  }
626 
627  else if (filterType == kBackOneStepFiltering) {
628  for (i = sampling; i >= 1; i--) {
629  r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
630  (Int_t) TMath::Min(i, numberIterationsY);
631  for (y = r2; y < ssizey - r2; y++) {
632  for (x = r1; x < ssizex - r1; x++) {
633  a = spectrum[x][y];
634  b = -(spectrum[x - r1][y - r2] +
635  spectrum[x - r1][y + r2] + spectrum[x + r1][y -
636  r2]
637  + spectrum[x + r1][y + r2]) / 4 +
638  (spectrum[x][y - r2] + spectrum[x - r1][y] +
639  spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
640  if (b < a && b > 0)
641  a = b;
642  working_space[x][y] = a;
643  }
644  }
645  for (y = i; y < ssizey - i; y++) {
646  for (x = i; x < ssizex - i; x++) {
647  spectrum[x][y] = working_space[x][y];
648  }
649  }
650  }
651  }
652  }
653  for (i = 0; i < ssizex; i++)
654  delete[]working_space[i];
655  delete[]working_space;
656  return 0;
657 }
658 
659 ////////////////////////////////////////////////////////////////////////////////
660 /// This function calculates smoothed spectrum from source spectrum
661 /// based on Markov chain method.
662 /// The result is placed in the array pointed by source pointer.
663 ///
664 /// Function parameters:
665 /// - source-pointer to the array of source spectrum
666 /// - ssizex-x length of source
667 /// - ssizey-y length of source
668 /// - averWindow-width of averaging smoothing window
669 ///
670 /// ### Smoothing
671 ///
672 /// Goal: Suppression of statistical fluctuations the algorithm is based on discrete
673 /// Markov chain, which has very simple invariant distribution
674 /// \f[
675 /// U_2 = \frac{p_{1.2}}{p_{2,1}}U_1, U_3 = \frac{p_{2,3}}{p_{3,2}}U_2 U_1, ... , U_n = \frac{p_{n-1,n}}{p_{n,n-1}}U_{n-1} ... U_2 U_1
676 /// \f]
677 /// \f$U_1\f$ being defined from the normalization condition \f$ \sum_{i=1}^{n} U_i = 1\f$
678 /// n is the length of the smoothed spectrum and
679 /// \f[
680 /// p_{i,i\pm1} = A_i \sum_{k=1}^{m} exp\left[\frac{y(i\pm k)-y(i)}{y(i\pm k)+y(i)}\right]
681 /// \f]
682 /// is the probability of the change of the peak position from channel i to the channel i+1.
683 /// \f$A_i\f$ is the normalization constant so that\f$ p_{i,i-1}+p_{i,i+1}=1\f$ and m is a width
684 /// of smoothing window. We have extended this algorithm to two dimensions.
685 ///
686 /// #### Reference:
687 ///
688 /// [1] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 (1996), 451.
689 ///
690 /// ### Example 4 - script Smooth.c
691 ///
692 /// \image html TSpectrum2_Smoothing1.jpg Fig. 9 Original noisy spectrum.
693 /// \image html TSpectrum2_Smoothing2.jpg Fig. 10 Smoothed spectrum m=3 Peaks can hardly be observed. Peaks become apparent.
694 /// \image html TSpectrum2_Smoothing3.jpg Fig. 11 Smoothed spectrum m=5
695 /// \image html TSpectrum2_Smoothing4.jpg Fig.12 Smoothed spectrum m=7
696 ///
697 /// #### Script:
698 ///
699 /// Example to illustrate the Markov smoothing (class TSpectrum).
700 /// To execute this example, do
701 ///
702 /// `root > .x Smooth.C`
703 ///
704 /// ~~~ {.cpp}
705 /// #include <TSpectrum>
706 /// void Smooth() {
707 /// Int_t i, j;
708 /// Double_t nbinsx = 256;
709 /// Double_t nbinsy = 256;
710 /// Double_t xmin = 0;
711 /// Double_t xmax = (Double_t)nbinsx;
712 /// Double_t ymin = 0;
713 /// Double_t ymax = (Double_t)nbinsy;
714 /// Double_t** source = new Double_t*[nbinsx];
715 /// for (i=0;i<nbinsx;i++)
716 /// source[i] = new Double_t[nbinsy];
717 /// TH2F *smooth = new TH2F("smooth","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
718 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
719 /// smooth=(TH2F*) f->Get("smooth1;1");
720 /// TCanvas *Smoothing = new TCanvas("Smoothing","Markov smoothing",10,10,1000,700);
721 /// TSpectrum *s = new TSpectrum();
722 /// for (i = 0; i < nbinsx; i++){
723 /// for (j = 0; j < nbinsy; j++){
724 /// source[i][j] = smooth->GetBinContent(i + 1,j + 1);
725 /// }
726 /// }
727 /// s->SmoothMarkov(source,nbinsx,nbinsx,3); //5,7
728 /// for (i = 0; i < nbinsx; i++){
729 /// for (j = 0; j < nbinsy; j++)
730 /// smooth->SetBinContent(i + 1,j + 1, source[i][j]);
731 /// }
732 /// smooth->Draw("SURF");
733 /// }
734 /// ~~~
735 
736 const char* TSpectrum2::SmoothMarkov(Double_t **source, Int_t ssizex, Int_t ssizey, Int_t averWindow)
737 {
738 
739  Int_t xmin, xmax, ymin, ymax, i, j, l;
740  Double_t a, b, maxch;
741  Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy, plocha = 0;
742  if(averWindow <= 0)
743  return "Averaging Window must be positive";
744  Double_t **working_space = new Double_t*[ssizex];
745  for(i = 0; i < ssizex; i++)
746  working_space[i] = new Double_t[ssizey];
747  xmin = 0;
748  xmax = ssizex - 1;
749  ymin = 0;
750  ymax = ssizey - 1;
751  for(i = 0, maxch = 0; i < ssizex; i++){
752  for(j = 0; j < ssizey; j++){
753  working_space[i][j] = 0;
754  if(maxch < source[i][j])
755  maxch = source[i][j];
756 
757  plocha += source[i][j];
758  }
759  }
760  if(maxch == 0) {
761  delete [] working_space;
762  return 0;
763  }
764 
765  nom = 0;
766  working_space[xmin][ymin] = 1;
767  for(i = xmin; i < xmax; i++){
768  nip = source[i][ymin] / maxch;
769  nim = source[i + 1][ymin] / maxch;
770  sp = 0,sm = 0;
771  for(l = 1; l <= averWindow; l++){
772  if((i + l) > xmax)
773  a = source[xmax][ymin] / maxch;
774 
775  else
776  a = source[i + l][ymin] / maxch;
777  b = a - nip;
778  if(a + nip <= 0)
779  a = 1;
780 
781  else
782  a = TMath::Sqrt(a + nip);
783  b = b / a;
784  b = TMath::Exp(b);
785  sp = sp + b;
786  if(i - l + 1 < xmin)
787  a = source[xmin][ymin] / maxch;
788 
789  else
790  a = source[i - l + 1][ymin] / maxch;
791  b = a - nim;
792  if(a + nim <= 0)
793  a = 1;
794 
795  else
796  a = TMath::Sqrt(a + nim);
797  b = b / a;
798  b = TMath::Exp(b);
799  sm = sm + b;
800  }
801  a = sp / sm;
802  a = working_space[i + 1][ymin] = a * working_space[i][ymin];
803  nom = nom + a;
804  }
805  for(i = ymin; i < ymax; i++){
806  nip = source[xmin][i] / maxch;
807  nim = source[xmin][i + 1] / maxch;
808  sp = 0,sm = 0;
809  for(l = 1; l <= averWindow; l++){
810  if((i + l) > ymax)
811  a = source[xmin][ymax] / maxch;
812 
813  else
814  a = source[xmin][i + l] / maxch;
815  b = a - nip;
816  if(a + nip <= 0)
817  a = 1;
818 
819  else
820  a = TMath::Sqrt(a + nip);
821  b = b / a;
822  b = TMath::Exp(b);
823  sp = sp + b;
824  if(i - l + 1 < ymin)
825  a = source[xmin][ymin] / maxch;
826 
827  else
828  a = source[xmin][i - l + 1] / maxch;
829  b = a - nim;
830  if(a + nim <= 0)
831  a = 1;
832 
833  else
834  a = TMath::Sqrt(a + nim);
835  b = b / a;
836  b = TMath::Exp(b);
837  sm = sm + b;
838  }
839  a = sp / sm;
840  a = working_space[xmin][i + 1] = a * working_space[xmin][i];
841  nom = nom + a;
842  }
843  for(i = xmin; i < xmax; i++){
844  for(j = ymin; j < ymax; j++){
845  nip = source[i][j + 1] / maxch;
846  nim = source[i + 1][j + 1] / maxch;
847  spx = 0,smx = 0;
848  for(l = 1; l <= averWindow; l++){
849  if(i + l > xmax)
850  a = source[xmax][j] / maxch;
851 
852  else
853  a = source[i + l][j] / maxch;
854  b = a - nip;
855  if(a + nip <= 0)
856  a = 1;
857 
858  else
859  a = TMath::Sqrt(a + nip);
860  b = b / a;
861  b = TMath::Exp(b);
862  spx = spx + b;
863  if(i - l + 1 < xmin)
864  a = source[xmin][j] / maxch;
865 
866  else
867  a = source[i - l + 1][j] / maxch;
868  b = a - nim;
869  if(a + nim <= 0)
870  a = 1;
871 
872  else
873  a = TMath::Sqrt(a + nim);
874  b = b / a;
875  b = TMath::Exp(b);
876  smx = smx + b;
877  }
878  spy = 0,smy = 0;
879  nip = source[i + 1][j] / maxch;
880  nim = source[i + 1][j + 1] / maxch;
881  for (l = 1; l <= averWindow; l++) {
882  if (j + l > ymax) a = source[i][ymax]/maxch;
883  else a = source[i][j + l] / maxch;
884  b = a - nip;
885  if (a + nip <= 0) a = 1;
886  else a = TMath::Sqrt(a + nip);
887  b = b / a;
888  b = TMath::Exp(b);
889  spy = spy + b;
890  if (j - l + 1 < ymin) a = source[i][ymin] / maxch;
891  else a = source[i][j - l + 1] / maxch;
892  b = a - nim;
893  if (a + nim <= 0) a = 1;
894  else a = TMath::Sqrt(a + nim);
895  b = b / a;
896  b = TMath::Exp(b);
897  smy = smy + b;
898  }
899  a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx +smy);
900  working_space[i + 1][j + 1] = a;
901  nom = nom + a;
902  }
903  }
904  for(i = xmin; i <= xmax; i++){
905  for(j = ymin; j <= ymax; j++){
906  working_space[i][j] = working_space[i][j] / nom;
907  }
908  }
909  for(i = 0;i < ssizex; i++){
910  for(j = 0; j < ssizey; j++){
911  source[i][j] = plocha * working_space[i][j];
912  }
913  }
914  for (i = 0; i < ssizex; i++)
915  delete[]working_space[i];
916  delete[]working_space;
917  return 0;
918 }
919 
920 ////////////////////////////////////////////////////////////////////////////////
921 /// This function calculates deconvolution from source spectrum
922 /// according to response spectrum
923 /// The result is placed in the matrix pointed by source pointer.
924 ///
925 /// Function parameters:
926 /// - source-pointer to the matrix of source spectrum
927 /// - resp-pointer to the matrix of response spectrum
928 /// - ssizex-x length of source and response spectra
929 /// - ssizey-y length of source and response spectra
930 /// - numberIterations, for details we refer to manual
931 /// - numberRepetitions, for details we refer to manual
932 /// - boost, boosting factor, for details we refer to manual
933 ///
934 /// ### Deconvolution
935 ///
936 /// Goal: Improvement of the resolution in spectra, decomposition of multiplets
937 ///
938 /// Mathematical formulation of the 2-dimensional convolution system is
939 ///
940 ///\f[
941 /// y(i_1,i_2) = \sum_{k_1=0}^{N_1-1} \sum_{k_2=0}^{N_2-1} h(i_1-k_1,i_2-k_2)x(k_1,k_2)
942 ///\f]
943 ///\f[
944 /// i_1 = 0,1,2, ... ,N_1-1, i_2 = 0,1,2, ... ,N_2-1
945 ///\f]
946 ///
947 /// where h(i,j) is the impulse response function, x, y are input and output
948 /// matrices, respectively, \f$ N_1, N_2\f$ are the lengths of x and h matrices
949 ///
950 /// - let us assume that we know the response and the output matrices (spectra) of the above given system.
951 /// - the deconvolution represents solution of the overdetermined system of linear equations, i.e., the
952 /// calculation of the matrix x.
953 /// - from numerical stability point of view the operation of deconvolution is
954 /// extremely critical (ill-posed problem) as well as time consuming operation.
955 /// - the Gold deconvolution algorithm proves to work very well even for 2-dimensional
956 /// systems. Generalisation of the algorithm for 2-dimensional systems was presented in [1], [2].
957 /// - for Gold deconvolution algorithm as well as for boosted deconvolution algorithm we refer also to TSpectrum
958 ///
959 /// #### References:
960 ///
961 /// [1] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:
962 /// Efficient one- and two-dimensional Gold deconvolution and its application to
963 /// gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
964 ///
965 /// [2] Morhac; M., Matouoek V., Kliman J., Efficient algorithm of multidimensional
966 /// deconvolution and its application to nuclear data processing, Digital Signal
967 /// Processing 13 (2003) 144.
968 ///
969 /// ### Example 5 - script Decon.c
970 ///
971 /// response function (usually peak) should be shifted to the beginning of the coordinate system (see Fig. 13)
972 ///
973 /// \image html TSpectrum2_Deconvolution1.jpg Fig. 13 2-dimensional response spectrum
974 /// \image html TSpectrum2_Deconvolution2.jpg Fig. 14 2-dimensional gamma-gamma-ray input spectrum (before deconvolution)
975 /// \image html TSpectrum2_Deconvolution3.jpg Fig. 15 Spectrum from Fig. 14 after deconvolution (1000 iterations)
976 ///
977 /// #### Script:
978 ///
979 /// Example to illustrate the Gold deconvolution (class TSpectrum2).
980 /// To execute this example, do
981 ///
982 /// `root > .x Decon.C`
983 ///
984 /// ~~~ {.cpp}
985 /// #include <TSpectrum2>
986 /// void Decon() {
987 /// Int_t i, j;
988 /// Double_t nbinsx = 256;
989 /// Double_t nbinsy = 256;
990 /// Double_t xmin = 0;
991 /// Double_t xmax = (Double_t)nbinsx;
992 /// Double_t ymin = 0;
993 /// Double_t ymax = (Double_t)nbinsy;
994 /// Double_t** source = new Double_t*[nbinsx];
995 /// for (i=0;i<nbinsx;i++)
996 /// source[i]=new Double_t[nbinsy];
997 /// TH2F *decon = new TH2F("decon","Gold deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
998 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
999 /// decon=(TH2F*) f->Get("decon1;1");
1000 /// Double_t** response = new Double_t*[nbinsx];
1001 /// for (i=0;i<nbinsx;i++)
1002 /// response[i]=new Double_t[nbinsy];
1003 /// TH2F *resp = new TH2F("resp","Response matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1004 /// resp=(TH2F*) f->Get("resp1;1");
1005 /// TCanvas *Deconvol = new TCanvas("Deconvolution","Gold deconvolution",10,10,1000,700);
1006 /// TSpectrum *s = new TSpectrum();
1007 /// for (i = 0; i < nbinsx; i++){
1008 /// for (j = 0; j < nbinsy; j++){
1009 /// source[i][j] = decon->GetBinContent(i + 1,j + 1);
1010 /// }
1011 /// }
1012 /// for (i = 0; i < nbinsx; i++){
1013 /// for (j = 0; j < nbinsy; j++){
1014 /// response[i][j] = resp->GetBinContent(i + 1,j + 1);
1015 /// }
1016 /// }
1017 /// s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
1018 /// for (i = 0; i < nbinsx; i++){
1019 /// for (j = 0; j < nbinsy; j++)
1020 /// decon->SetBinContent(i + 1,j + 1, source[i][j]);
1021 /// }
1022 /// decon->Draw("SURF");
1023 /// }
1024 /// ~~~
1025 ///
1026 /// ### Example 6 - script Decon2.c
1027 ///
1028 /// \image html TSpectrum2_Deconvolution4.jpg Fig. 16 Response spectrum
1029 /// \image html TSpectrum2_Deconvolution5.jpg Fig. 17 Original synthetic input spectrum (before deconvolution). It is composed of 17 peaks. 5 peaks are overlapping in the outlined multiplet and two peaks in doublet.
1030 /// \image html TSpectrum2_Deconvolution6.jpg Fig. 18 Spectrum from Fig. 17 after deconvolution (1000 iterations). Resolution is improved but the peaks in multiplet remained unresolved.
1031 ///
1032 /// #### Script:
1033 ///
1034 /// Example to illustrate the Gold deconvolution (class TSpectrum2).
1035 /// To execute this example, do
1036 ///
1037 /// `root > .x Decon2.C`
1038 ///
1039 /// ~~~ {.cpp}
1040 /// #include <TSpectrum2>
1041 /// void Decon2() {
1042 /// Int_t i, j;
1043 /// Double_t nbinsx = 64;
1044 /// Double_t nbinsy = 64;
1045 /// Double_t xmin = 0;
1046 /// Double_t xmax = (Double_t)nbinsx;
1047 /// Double_t ymin = 0;
1048 /// Double_t ymax = (Double_t)nbinsy;
1049 /// Double_t** source = new Double_t*[nbinsx];
1050 /// for (i=0;i<nbinsx;i++)
1051 /// source[i]=new Double_t[nbinsy];
1052 /// TH2F *decon = new TH2F("decon","Gold deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1053 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
1054 /// decon=(TH2F*) f->Get("decon2;1");
1055 /// Double_t** response = new Double_t*[nbinsx];
1056 /// for (i=0;i<nbinsx;i++)
1057 /// response[i]=new Double_t[nbinsy];
1058 /// TH2F *resp = new TH2F("resp","Response matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1059 /// resp=(TH2F*) f->Get("resp2;1");
1060 /// TCanvas *Deconvol = new TCanvas("Deconvolution","Gold deconvolution",10,10,1000,700);
1061 /// TSpectrum *s = new TSpectrum();
1062 /// for (i = 0; i < nbinsx; i++){
1063 /// for (j = 0; j < nbinsy; j++){
1064 /// source[i][j] = decon->GetBinContent(i + 1,j + 1);
1065 /// }
1066 /// }
1067 /// for (i = 0; i < nbinsx; i++){
1068 /// for (j = 0; j < nbinsy; j++){
1069 /// response[i][j] = resp->GetBinContent(i + 1,j + 1);
1070 /// }
1071 /// }
1072 /// s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
1073 /// for (i = 0; i < nbinsx; i++){
1074 /// for (j = 0; j < nbinsy; j++)
1075 /// decon->SetBinContent(i + 1,j + 1, source[i][j]);
1076 /// }
1077 /// decon->Draw("SURF");
1078 /// }
1079 /// ~~~
1080 ///
1081 /// ### Example 7 - script Decon2HR.c
1082 ///
1083 /// \image html TSpectrum2_Deconvolution7.jpg Fig. 19 Spectrum from Fig. 17 after boosted deconvolution (50 iterations repeated 20 times, boosting coefficient was 1.2). All the peaks in multiplet as well as in doublet are completely decomposed.
1084 ///
1085 /// #### Script:
1086 ///
1087 /// Example to illustrate boosted Gold deconvolution (class TSpectrum2).
1088 /// To execute this example, do
1089 ///
1090 /// `root > .x Decon2HR.C`
1091 ///
1092 /// ~~~ {.cpp}
1093 /// #include <TSpectrum2>
1094 /// void Decon2HR() {
1095 /// Int_t i, j;
1096 /// Double_t nbinsx = 64;
1097 /// Double_t nbinsy = 64;
1098 /// Double_t xmin = 0;
1099 /// Double_t xmax = (Double_t)nbinsx;
1100 /// Double_t ymin = 0;
1101 /// Double_t ymax = (Double_t)nbinsy;
1102 /// Double_t** source = new Double_t*[nbinsx];
1103 /// for (i=0;i<nbinsx;i++)
1104 /// source[i]=new Double_t[nbinsy];
1105 /// TH2F *decon = new TH2F("decon","Boosted Gold deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1106 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
1107 /// decon=(TH2F*) f->Get("decon2;1");
1108 /// Double_t** response = new Double_t*[nbinsx];
1109 /// for (i=0;i<nbinsx;i++)
1110 /// response[i]=new Double_t[nbinsy];
1111 /// TH2F *resp = new TH2F("resp","Response matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1112 /// resp=(TH2F*) f->Get("resp2;1");
1113 /// TCanvas *Deconvol = new TCanvas("Deconvolution","Gold deconvolution",10,10,1000,700);
1114 /// TSpectrum *s = new TSpectrum();
1115 /// for (i = 0; i < nbinsx; i++){
1116 /// for (j = 0; j < nbinsy; j++){
1117 /// source[i][j] = decon->GetBinContent(i + 1,j + 1);
1118 /// }
1119 /// }
1120 /// for (i = 0; i < nbinsx; i++){
1121 /// for (j = 0; j < nbinsy; j++){
1122 /// response[i][j] = resp->GetBinContent(i + 1,j + 1);
1123 /// }
1124 /// }
1125 /// s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
1126 /// for (i = 0; i < nbinsx; i++){
1127 /// for (j = 0; j < nbinsy; j++)
1128 /// dec on->SetBinContent(i + 1,j + 1, source[i][j]);
1129 /// }
1130 /// decon->Draw("SURF");
1131 /// }
1132 /// ~~~
1133 
1134 const char *TSpectrum2::Deconvolution(Double_t **source, Double_t **resp,
1135  Int_t ssizex, Int_t ssizey,
1136  Int_t numberIterations,
1137  Int_t numberRepetitions,
1138  Double_t boost)
1139 {
1140  Int_t i, j, lhx, lhy, i1, i2, j1, j2, k1, k2, lindex, i1min, i1max,
1141  i2min, i2max, j1min, j1max, j2min, j2max, positx = 0, posity = 0, repet;
1142  Double_t lda, ldb, ldc, area, maximum = 0;
1143  if (ssizex <= 0 || ssizey <= 0)
1144  return "Wrong parameters";
1145  if (numberIterations <= 0)
1146  return "Number of iterations must be positive";
1147  if (numberRepetitions <= 0)
1148  return "Number of repetitions must be positive";
1149  Double_t **working_space = new Double_t *[ssizex];
1150  for (i = 0; i < ssizex; i++)
1151  working_space[i] = new Double_t[5 * ssizey];
1152  area = 0;
1153  lhx = -1, lhy = -1;
1154  for (i = 0; i < ssizex; i++) {
1155  for (j = 0; j < ssizey; j++) {
1156  lda = resp[i][j];
1157  if (lda != 0) {
1158  if ((i + 1) > lhx)
1159  lhx = i + 1;
1160  if ((j + 1) > lhy)
1161  lhy = j + 1;
1162  }
1163  working_space[i][j] = lda;
1164  area = area + lda;
1165  if (lda > maximum) {
1166  maximum = lda;
1167  positx = i, posity = j;
1168  }
1169  }
1170  }
1171  if (lhx == -1 || lhy == -1) {
1172  delete [] working_space;
1173  return ("Zero response data");
1174  }
1175 
1176 //calculate ht*y and write into p
1177  for (i2 = 0; i2 < ssizey; i2++) {
1178  for (i1 = 0; i1 < ssizex; i1++) {
1179  ldc = 0;
1180  for (j2 = 0; j2 <= (lhy - 1); j2++) {
1181  for (j1 = 0; j1 <= (lhx - 1); j1++) {
1182  k2 = i2 + j2, k1 = i1 + j1;
1183  if (k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
1184  lda = working_space[j1][j2];
1185  ldb = source[k1][k2];
1186  ldc = ldc + lda * ldb;
1187  }
1188  }
1189  }
1190  working_space[i1][i2 + ssizey] = ldc;
1191  }
1192  }
1193 
1194 //calculate matrix b=ht*h
1195  i1min = -(lhx - 1), i1max = lhx - 1;
1196  i2min = -(lhy - 1), i2max = lhy - 1;
1197  for (i2 = i2min; i2 <= i2max; i2++) {
1198  for (i1 = i1min; i1 <= i1max; i1++) {
1199  ldc = 0;
1200  j2min = -i2;
1201  if (j2min < 0)
1202  j2min = 0;
1203  j2max = lhy - 1 - i2;
1204  if (j2max > lhy - 1)
1205  j2max = lhy - 1;
1206  for (j2 = j2min; j2 <= j2max; j2++) {
1207  j1min = -i1;
1208  if (j1min < 0)
1209  j1min = 0;
1210  j1max = lhx - 1 - i1;
1211  if (j1max > lhx - 1)
1212  j1max = lhx - 1;
1213  for (j1 = j1min; j1 <= j1max; j1++) {
1214  lda = working_space[j1][j2];
1215  if (i1 + j1 < ssizex && i2 + j2 < ssizey)
1216  ldb = working_space[i1 + j1][i2 + j2];
1217  else
1218  ldb = 0;
1219  ldc = ldc + lda * ldb;
1220  }
1221  }
1222  working_space[i1 - i1min][i2 - i2min + 2 * ssizey ] = ldc;
1223  }
1224  }
1225 
1226 //initialization in x1 matrix
1227  for (i2 = 0; i2 < ssizey; i2++) {
1228  for (i1 = 0; i1 < ssizex; i1++) {
1229  working_space[i1][i2 + 3 * ssizey] = 1;
1230  working_space[i1][i2 + 4 * ssizey] = 0;
1231  }
1232  }
1233 
1234  //START OF ITERATIONS
1235  for (repet = 0; repet < numberRepetitions; repet++) {
1236  if (repet != 0) {
1237  for (i = 0; i < ssizex; i++) {
1238  for (j = 0; j < ssizey; j++) {
1239  working_space[i][j + 3 * ssizey] =
1240  TMath::Power(working_space[i][j + 3 * ssizey], boost);
1241  }
1242  }
1243  }
1244  for (lindex = 0; lindex < numberIterations; lindex++) {
1245  for (i2 = 0; i2 < ssizey; i2++) {
1246  for (i1 = 0; i1 < ssizex; i1++) {
1247  ldb = 0;
1248  j2min = i2;
1249  if (j2min > lhy - 1)
1250  j2min = lhy - 1;
1251  j2min = -j2min;
1252  j2max = ssizey - i2 - 1;
1253  if (j2max > lhy - 1)
1254  j2max = lhy - 1;
1255  j1min = i1;
1256  if (j1min > lhx - 1)
1257  j1min = lhx - 1;
1258  j1min = -j1min;
1259  j1max = ssizex - i1 - 1;
1260  if (j1max > lhx - 1)
1261  j1max = lhx - 1;
1262  for (j2 = j2min; j2 <= j2max; j2++) {
1263  for (j1 = j1min; j1 <= j1max; j1++) {
1264  ldc = working_space[j1 - i1min][j2 - i2min + 2 * ssizey];
1265  lda = working_space[i1 + j1][i2 + j2 + 3 * ssizey];
1266  ldb = ldb + lda * ldc;
1267  }
1268  }
1269  lda = working_space[i1][i2 + 3 * ssizey];
1270  ldc = working_space[i1][i2 + 1 * ssizey];
1271  if (ldc * lda != 0 && ldb != 0) {
1272  lda = lda * ldc / ldb;
1273  }
1274 
1275  else
1276  lda = 0;
1277  working_space[i1][i2 + 4 * ssizey] = lda;
1278  }
1279  }
1280  for (i2 = 0; i2 < ssizey; i2++) {
1281  for (i1 = 0; i1 < ssizex; i1++)
1282  working_space[i1][i2 + 3 * ssizey] =
1283  working_space[i1][i2 + 4 * ssizey];
1284  }
1285  }
1286  }
1287  for (i = 0; i < ssizex; i++) {
1288  for (j = 0; j < ssizey; j++)
1289  source[(i + positx) % ssizex][(j + posity) % ssizey] =
1290  area * working_space[i][j + 3 * ssizey];
1291  }
1292  for (i = 0; i < ssizex; i++)
1293  delete[]working_space[i];
1294  delete[]working_space;
1295  return 0;
1296 }
1297 
1298 ////////////////////////////////////////////////////////////////////////////////
1299 /// This function searches for peaks in source spectrum
1300 /// It is based on deconvolution method. First the background is
1301 /// removed (if desired), then Markov spectrum is calculated
1302 /// (if desired), then the response function is generated
1303 /// according to given sigma and deconvolution is carried out.
1304 ///
1305 /// Function parameters:
1306 /// - source-pointer to the matrix of source spectrum
1307 /// - dest-pointer to the matrix of resulting deconvolved spectrum
1308 /// - ssizex-x length of source spectrum
1309 /// - ssizey-y length of source spectrum
1310 /// - sigma-sigma of searched peaks, for details we refer to manual
1311 /// - threshold-threshold value in % for selected peaks, peaks with
1312 /// amplitude less than threshold*highest_peak/100
1313 /// are ignored, see manual
1314 /// - backgroundRemove-logical variable, set if the removal of
1315 /// background before deconvolution is desired
1316 /// - deconIterations-number of iterations in deconvolution operation
1317 /// - markov-logical variable, if it is true, first the source spectrum
1318 /// is replaced by new spectrum calculated using Markov
1319 /// chains method.
1320 /// - averWindow-averaging window of searched peaks, for details
1321 /// we refer to manual (applies only for Markov method)
1322 ///
1323 /// ### Peaks searching
1324 ///
1325 /// Goal: to identify automatically the peaks in spectrum with the presence of the
1326 /// continuous background, one-fold coincidences (ridges) and statistical
1327 /// fluctuations - noise.
1328 ///
1329 /// The common problems connected with correct peak identification in two-dimensional coincidence spectra are
1330 ///
1331 /// - non-sensitivity to noise, i.e., only statistically relevant peaks should be identified
1332 /// - non-sensitivity of the algorithm to continuous background
1333 /// - non-sensitivity to one-fold coincidences (coincidences peak - background in both dimensions) and their crossings
1334 /// - ability to identify peaks close to the edges of the spectrum region. Usually peak finders fail to detect them
1335 /// - resolution, decomposition of doublets and multiplets. The algorithm should be able to recognise close positioned peaks.
1336 /// - ability to identify peaks with different sigma
1337 ///
1338 /// #### References:
1339 ///
1340 /// [1] M.A. Mariscotti: A method for identification of peaks in the presence of
1341 /// background and its application to spectrum analysis. NIM 50 (1967), 309-320.
1342 ///
1343 /// [2] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:Identification
1344 /// of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
1345 /// 108-125.
1346 ///
1347 /// [3] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
1348 /// (1996), 451.
1349 ///
1350 /// ### Examples of peak searching method
1351 ///
1352 /// SearchHighRes function provides users with the possibility
1353 /// to vary the input parameters and with the access to the output deconvolved data
1354 /// in the destination spectrum. Based on the output data one can tune the
1355 /// parameters.
1356 ///
1357 /// ### Example 8 - script Src.c
1358 ///
1359 /// \image html TSpectrum2_Searching1.jpg Fig. 20 Two-dimensional spectrum with found peaks denoted by markers (sigma=2, threshold=5%, 3 iterations steps in the deconvolution)
1360 /// \image html TSpectrum2_Searching3.jpg Fig. 21 Spectrum from Fig. 20 after background elimination and deconvolution
1361 ///
1362 /// #### Script:
1363 ///
1364 /// Example to illustrate high resolution peak searching function (class TSpectrum).
1365 /// To execute this example, do
1366 ///
1367 /// `root > .x Src.C
1368 ///
1369 /// ~~~ {.cpp}
1370 /// #include <TSpectrum2>
1371 /// void Src() {
1372 /// Int_t i, j, nfound;
1373 /// Double_t nbinsx = 64;
1374 /// Double_t nbinsy = 64;
1375 /// Double_t xmin = 0;
1376 /// Double_t xmax = (Double_t)nbinsx;
1377 /// Double_t ymin = 0;
1378 /// Double_t ymax = (Double_t)nbinsy;
1379 /// Double_t** source = new Double_t*[nbinsx];
1380 /// for (i=0;i<nbinsx;i++)
1381 /// source[i]=new Double_t[nbinsy];
1382 /// Double_t** dest = new Double_t*[nbinsx];
1383 /// for (i=0;i<nbinsx;i++)
1384 /// dest[i]=new Double_t[nbinsy];
1385 /// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);/
1386 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
1387 /// search=(TH2F*) f->Get("search4;1");
1388 /// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1389 /// TSpectrum2 *s = new TSpectrum2();
1390 /// for (i = 0; i < nbinsx; i++){
1391 /// for (j = 0; j < nbinsy; j++){
1392 /// source[i][j] = search->GetBinContent(i + 1,j + 1);
1393 /// }
1394 /// }
1395 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);
1396 /// printf("Found %d candidate peaks\n",nfound);
1397 /// for(i=0;i<nfound;i++)
1398 /// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1399 /// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1400 /// }
1401 /// ~~~
1402 ///
1403 /// ### Example 9 - script Src2.c
1404 ///
1405 /// \image html TSpectrum2_Searching4.jpg Fig. 22 Two-dimensional noisy spectrum with found peaks denoted by markers (sigma=2, threshold=10%, 10 iterations steps in the deconvolution). One can observe that the algorithm is insensitive to the crossings of one-dimensional ridges. It identifies only two-coincidence peaks.
1406 /// \image html TSpectrum2_Searching5.jpg Fig. 23 Spectrum from Fig. 22 after background elimination and deconvolution
1407 ///
1408 /// #### Script:
1409 ///
1410 /// Example to illustrate high resolution peak searching function (class TSpectrum).
1411 /// To execute this example, do
1412 ///
1413 /// `root > .x Src2.C`
1414 ///
1415 /// ~~~ {.cpp}
1416 /// #include <TSpectrum2>
1417 /// void Src2() {
1418 /// Int_t i, j, nfound;
1419 /// Double_t nbinsx = 256;
1420 /// Double_t nbinsy = 256;
1421 /// Double_t xmin = 0;
1422 /// Double_t xmax = (Double_t)nbinsx;
1423 /// Double_t ymin = 0;
1424 /// Double_t ymax = (Double_t)nbinsy;
1425 /// Double_t** source = new Double_t*[nbinsx];
1426 /// for (i=0;i<nbinsx;i++)
1427 /// source[i]=new Double_t[nbinsy];
1428 /// Double_t** dest = new Double_t*[nbinsx];
1429 /// for (i=0;i<nbinsx;i++)
1430 /// dest[i]=new Double_t[nbinsy];
1431 /// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1432 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
1433 /// search=(TH2F*) f->Get("back3;1");
1434 /// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1435 /// TSpectrum2 *s = new TSpectrum2();
1436 /// for (i = 0; i < nbinsx; i++){
1437 /// for (j = 0; j < nbinsy; j++){
1438 /// source[i][j] = search->GetBinContent(i + 1,j + 1);
1439 /// }
1440 /// }
1441 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 10, kTRUE, 10, kFALSE, 3);
1442 /// printf("Found %d candidate peaks\n",nfound);
1443 /// for(i=0;i<nfound;i++)
1444 /// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1445 /// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1446 /// }
1447 /// ~~~
1448 ///
1449 /// ### Example 10 - script Src3.c
1450 ///
1451 /// \image html TSpectrum2_Searching6.jpg Fig. 24 Two-dimensional spectrum with 15 found peaks denoted by markers. Some peaks are positioned close to each other. It is necessary to increase number of iterations in the deconvolution. In next 3 Figs. we shall study the influence of this parameter.
1452 /// \image html TSpectrum2_Searching7.jpg Fig. 25 Spectrum from Fig. 24 after deconvolution (# of iterations = 3). Number of identified peaks = 13.
1453 /// \image html TSpectrum2_Searching8.jpg Fig. 26 Spectrum from Fig. 24 after deconvolution (# of iterations = 10). Number of identified peaks = 13.
1454 /// \image html TSpectrum2_Searching9.jpg Fig. 27 Spectrum from Fig. 24 after deconvolution (# of iterations = 100). Number of identified peaks = 15. Now the algorithm is able to decompose two doublets in the spectrum.
1455 ///
1456 /// #### Script:
1457 ///
1458 /// Example to illustrate high resolution peak searching function (class TSpectrum).
1459 /// To execute this example, do
1460 ///
1461 /// `root > .x Src3.C`
1462 ///
1463 /// ~~~ {.cpp}
1464 /// #include <TSpectrum2>
1465 /// void Src3() {
1466 /// Int_t i, j, nfound;
1467 /// Double_t nbinsx = 64;
1468 /// Double_t nbinsy = 64;
1469 /// Double_t xmin = 0;
1470 /// Double_t xmax = (Double_t)nbinsx;
1471 /// Double_t ymin = 0;
1472 /// Double_t ymax = (Double_t)nbinsy;
1473 /// Double_t** source = new Double_t*[nbinsx];
1474 /// for (i=0;i<nbinsx;i++)
1475 /// source[i]=new Double_t[nbinsy];
1476 /// Double_t** dest = new Double_t*[nbinsx];
1477 /// for (i=0;i<nbinsx;i++)
1478 /// dest[i]=new Double_t[nbinsy];
1479 /// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1480 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
1481 /// search=(TH2F*) f->Get("search1;1");
1482 /// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1483 /// TSpectrum2 *s = new TSpectrum2();
1484 /// for (i = 0; i < nbinsx; i++){
1485 /// for (j = 0; j < nbinsy; j++){
1486 /// source[i][j] = search->GetBinContent(i + 1,j + 1);
1487 /// }
1488 /// }
1489 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 2, kFALSE, 3, kFALSE, 1);//3, 10, 100
1490 /// printf("Found %d candidate peaks\n",nfound);
1491 /// for(i=0;i<nfound;i++)
1492 /// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1493 /// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1494 /// }
1495 /// ~~~
1496 ///
1497 /// ### Example 11 - script Src4.c
1498 ///
1499 /// \image html TSpectrum2_Searching10.jpg Fig. 28 Two-dimensional spectrum with peaks with different sigma denoted by markers (sigma=3, threshold=5%, 10 iterations steps in the deconvolution, Markov smoothing with window=3)
1500 /// \image html TSpectrum2_Searching12.jpg Fig. 29 Spectrum from Fig. 28 after smoothing and deconvolution.
1501 ///
1502 /// #### Script:
1503 ///
1504 /// Example to illustrate high resolution peak searching function (class TSpectrum).
1505 /// To execute this example, do
1506 ///
1507 /// `root > .x Src4.C`
1508 ///
1509 /// ~~~ {.cpp}
1510 /// #include <TSpectrum2>
1511 /// void Src4() {
1512 /// Int_t i, j, nfound;
1513 /// Double_t nbinsx = 64;
1514 /// Double_t nbinsy = 64;
1515 /// Double_t xmin = 0;
1516 /// Double_t xmax = (Double_t)nbinsx;
1517 /// Double_t ymin = 0;
1518 /// Double_t ymax = (Double_t)nbinsy;
1519 /// Double_t** source = new Double_t*[nbinsx];
1520 /// for (i=0;i<nbinsx;i++)
1521 /// source[i]=new Double_t[nbinsy];
1522 /// Double_t** dest = new Double_t*[nbinsx];
1523 /// for (i=0;i<nbinsx;i++)
1524 /// dest[i]=new Double_t[nbinsy];
1525 /// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1526 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
1527 /// search=(TH2F*) f->Get("search2;1");
1528 /// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1529 /// TSpectrum2 *s = new TSpectrum2();
1530 /// for (i = 0; i < nbinsx; i++){
1531 /// for (j = 0; j < nbinsy; j++){
1532 /// source[i][j] = search->GetBinContent(i + 1,j + 1);
1533 /// }
1534 /// }
1535 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 3, 5, kFALSE, 10, kTRUE, 3);
1536 /// printf("Found %d candidate peaks\n",nfound);
1537 /// for(i=0;i<nfound;i++)
1538 /// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1539 /// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1540 /// }
1541 /// ~~~
1542 ///
1543 /// ### Example 12 - script Src5.c`
1544 ///
1545 /// \image html TSpectrum2_Searching13.jpg Fig. 30 Two-dimensional spectrum with peaks positioned close to the edges denoted by markers (sigma=2, threshold=5%, 10 iterations steps in the deconvolution)
1546 /// \image html TSpectrum2_Searching14.jpg Fig. 31 Spectrum from Fig. 30 after deconvolution.
1547 ///
1548 /// #### Script:
1549 ///
1550 /// Example to illustrate high resolution peak searching function (class TSpectrum).
1551 /// To execute this example, do
1552 ///
1553 /// `root > .x Src5.C`
1554 ///
1555 /// ~~~ {.cpp}
1556 /// #include <TSpectrum2>
1557 /// void Src5() {
1558 /// Int_t i, j, nfound;
1559 /// Double_t nbinsx = 64;
1560 /// Double_t nbinsy = 64;
1561 /// Double_t xmin = 0;
1562 /// Double_t xmax = (Double_t)nbinsx;
1563 /// Double_t ymin = 0;
1564 /// Double_t ymax = (Double_t)nbinsy;
1565 /// Double_t** source = new Double_t*[nbinsx];
1566 /// for (i=0;i<nbinsx;i++)
1567 /// source[i]=new Double_t[nbinsy];
1568 /// Double_t** dest = new Double_t*[nbinsx];
1569 /// for (i=0;i<nbinsx;i++)
1570 /// dest[i]=new Double_t[nbinsy];
1571 /// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1572 /// TFile *f = new TFile("spectra2/TSpectrum2.root");
1573 /// search=(TH2F*) f->Get("search3;1");
1574 /// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1575 /// TSpectrum2 *s = new TSpectrum2();
1576 /// for (i = 0; i < nbinsx; i++){
1577 /// for (j = 0; j < nbinsy; j++){
1578 /// source[i][j] = search->GetBinContent(i + 1,j + 1);
1579 /// }
1580 /// }
1581 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kFALSE, 10, kFALSE, 1);
1582 /// printf("Found %d candidate peaks\n",nfound);
1583 /// for(i=0;i<nfound;i++)
1584 /// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1585 /// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1586 /// }
1587 /// ~~~
1588 
1590  Double_t sigma, Double_t threshold,
1591  Bool_t backgroundRemove,Int_t deconIterations,
1592  Bool_t markov, Int_t averWindow)
1593 
1594 {
1595  Int_t number_of_iterations = (Int_t)(4 * sigma + 0.5);
1596  Int_t k, lindex, priz;
1597  Double_t lda, ldb, ldc, area, maximum;
1598  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;
1599  Int_t ymin, ymax, i, j;
1600  Double_t a, b, ax, ay, maxch, plocha = 0;
1601  Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy;
1602  Double_t p1, p2, p3, p4, s1, s2, s3, s4;
1603  Int_t x, y;
1604  Int_t lhx, lhy, i1, i2, j1, j2, k1, k2, i1min, i1max, i2min, i2max, j1min, j1max, j2min, j2max, positx, posity;
1605  if (sigma < 1) {
1606  Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
1607  return 0;
1608  }
1609 
1610  if(threshold<=0||threshold>=100){
1611  Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
1612  return 0;
1613  }
1614 
1615  j = (Int_t) (5.0 * sigma + 0.5);
1616  if (j >= PEAK_WINDOW / 2) {
1617  Error("SearchHighRes", "Too large sigma");
1618  return 0;
1619  }
1620 
1621  if (markov == true) {
1622  if (averWindow <= 0) {
1623  Error("SearchHighRes", "Averaging window must be positive");
1624  return 0;
1625  }
1626  }
1627  if(backgroundRemove == true){
1628  if(ssizex_ext < 2 * number_of_iterations + 1 || ssizey_ext < 2 * number_of_iterations + 1){
1629  Error("SearchHighRes", "Too large clipping window");
1630  return 0;
1631  }
1632  }
1633  i = (Int_t)(4 * sigma + 0.5);
1634  i = 4 * i;
1635  Double_t **working_space = new Double_t *[ssizex + i];
1636  for (j = 0; j < ssizex + i; j++) {
1637  Double_t *wsk = working_space[j] = new Double_t[16 * (ssizey + i)];
1638  for (k=0;k<16 * (ssizey + i);k++) wsk[k] = 0;
1639  }
1640  for(j = 0; j < ssizey_ext; j++){
1641  for(i = 0; i < ssizex_ext; i++){
1642  if(i < shift){
1643  if(j < shift)
1644  working_space[i][j + ssizey_ext] = source[0][0];
1645 
1646  else if(j >= ssizey + shift)
1647  working_space[i][j + ssizey_ext] = source[0][ssizey - 1];
1648 
1649  else
1650  working_space[i][j + ssizey_ext] = source[0][j - shift];
1651  }
1652 
1653  else if(i >= ssizex + shift){
1654  if(j < shift)
1655  working_space[i][j + ssizey_ext] = source[ssizex - 1][0];
1656 
1657  else if(j >= ssizey + shift)
1658  working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1];
1659 
1660  else
1661  working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift];
1662  }
1663 
1664  else{
1665  if(j < shift)
1666  working_space[i][j + ssizey_ext] = source[i - shift][0];
1667 
1668  else if(j >= ssizey + shift)
1669  working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1];
1670 
1671  else
1672  working_space[i][j + ssizey_ext] = source[i - shift][j - shift];
1673  }
1674  }
1675  }
1676  if(backgroundRemove == true){
1677  for(i = 1; i <= number_of_iterations; i++){
1678  for(y = i; y < ssizey_ext - i; y++){
1679  for(x = i; x < ssizex_ext - i; x++){
1680  a = working_space[x][y + ssizey_ext];
1681  p1 = working_space[x - i][y + ssizey_ext - i];
1682  p2 = working_space[x - i][y + ssizey_ext + i];
1683  p3 = working_space[x + i][y + ssizey_ext - i];
1684  p4 = working_space[x + i][y + ssizey_ext + i];
1685  s1 = working_space[x][y + ssizey_ext - i];
1686  s2 = working_space[x - i][y + ssizey_ext];
1687  s3 = working_space[x + i][y + ssizey_ext];
1688  s4 = working_space[x][y + ssizey_ext + i];
1689  b = (p1 + p2) / 2.0;
1690  if(b > s2)
1691  s2 = b;
1692  b = (p1 + p3) / 2.0;
1693  if(b > s1)
1694  s1 = b;
1695  b = (p2 + p4) / 2.0;
1696  if(b > s4)
1697  s4 = b;
1698  b = (p3 + p4) / 2.0;
1699  if(b > s3)
1700  s3 = b;
1701  s1 = s1 - (p1 + p3) / 2.0;
1702  s2 = s2 - (p1 + p2) / 2.0;
1703  s3 = s3 - (p3 + p4) / 2.0;
1704  s4 = s4 - (p2 + p4) / 2.0;
1705  b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
1706  if(b < a)
1707  a = b;
1708  working_space[x][y] = a;
1709  }
1710  }
1711  for(y = i;y < ssizey_ext - i; y++){
1712  for(x = i; x < ssizex_ext - i; x++){
1713  working_space[x][y + ssizey_ext] = working_space[x][y];
1714  }
1715  }
1716  }
1717  for(j = 0;j < ssizey_ext; j++){
1718  for(i = 0; i < ssizex_ext; i++){
1719  if(i < shift){
1720  if(j < shift)
1721  working_space[i][j + ssizey_ext] = source[0][0] - working_space[i][j + ssizey_ext];
1722 
1723  else if(j >= ssizey + shift)
1724  working_space[i][j + ssizey_ext] = source[0][ssizey - 1] - working_space[i][j + ssizey_ext];
1725 
1726  else
1727  working_space[i][j + ssizey_ext] = source[0][j - shift] - working_space[i][j + ssizey_ext];
1728  }
1729 
1730  else if(i >= ssizex + shift){
1731  if(j < shift)
1732  working_space[i][j + ssizey_ext] = source[ssizex - 1][0] - working_space[i][j + ssizey_ext];
1733 
1734  else if(j >= ssizey + shift)
1735  working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1] - working_space[i][j + ssizey_ext];
1736 
1737  else
1738  working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift] - working_space[i][j + ssizey_ext];
1739  }
1740 
1741  else{
1742  if(j < shift)
1743  working_space[i][j + ssizey_ext] = source[i - shift][0] - working_space[i][j + ssizey_ext];
1744 
1745  else if(j >= ssizey + shift)
1746  working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1] - working_space[i][j + ssizey_ext];
1747 
1748  else
1749  working_space[i][j + ssizey_ext] = source[i - shift][j - shift] - working_space[i][j + ssizey_ext];
1750  }
1751  }
1752  }
1753  }
1754  for(j = 0; j < ssizey_ext; j++){
1755  for(i = 0; i < ssizex_ext; i++){
1756  working_space[i][j + 15*ssizey_ext] = working_space[i][j + ssizey_ext];
1757  }
1758  }
1759  if(markov == true){
1760  for(i = 0;i < ssizex_ext; i++){
1761  for(j = 0; j < ssizey_ext; j++)
1762  working_space[i][j + 2 * ssizex_ext] = working_space[i][ssizey_ext + j];
1763  }
1764  xmin = 0;
1765  xmax = ssizex_ext - 1;
1766  ymin = 0;
1767  ymax = ssizey_ext - 1;
1768  for(i = 0, maxch = 0; i < ssizex_ext; i++){
1769  for(j = 0; j < ssizey_ext; j++){
1770  working_space[i][j] = 0;
1771  if(maxch < working_space[i][j + 2 * ssizey_ext])
1772  maxch = working_space[i][j + 2 * ssizey_ext];
1773  plocha += working_space[i][j + 2 * ssizey_ext];
1774  }
1775  }
1776  if(maxch == 0) {
1777  delete [] working_space;
1778  return 0;
1779  }
1780 
1781  nom=0;
1782  working_space[xmin][ymin] = 1;
1783  for(i = xmin; i < xmax; i++){
1784  nip = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1785  nim = working_space[i + 1][ymin + 2 * ssizey_ext] / maxch;
1786  sp = 0,sm = 0;
1787  for(l = 1;l <= averWindow; l++){
1788  if((i + l) > xmax)
1789  a = working_space[xmax][ymin + 2 * ssizey_ext] / maxch;
1790 
1791  else
1792  a = working_space[i + l][ymin + 2 * ssizey_ext] / maxch;
1793 
1794  b = a - nip;
1795  if(a + nip <= 0)
1796  a = 1;
1797 
1798  else
1799  a=TMath::Sqrt(a + nip);
1800  b = b / a;
1801  b = TMath::Exp(b);
1802  sp = sp + b;
1803  if(i - l + 1 < xmin)
1804  a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1805 
1806  else
1807  a = working_space[i - l + 1][ymin + 2 * ssizey_ext] / maxch;
1808  b = a - nim;
1809  if(a + nim <= 0)
1810  a = 1;
1811 
1812  else
1813  a=TMath::Sqrt(a + nim);
1814  b = b / a;
1815  b = TMath::Exp(b);
1816  sm = sm + b;
1817  }
1818  a = sp / sm;
1819  a = working_space[i + 1][ymin] = a * working_space[i][ymin];
1820  nom = nom + a;
1821  }
1822  for(i = ymin; i < ymax; i++){
1823  nip = working_space[xmin][i + 2 * ssizey_ext] / maxch;
1824  nim = working_space[xmin][i + 1 + 2 * ssizey_ext] / maxch;
1825  sp = 0,sm = 0;
1826  for(l = 1; l <= averWindow; l++){
1827  if((i + l) > ymax)
1828  a = working_space[xmin][ymax + 2 * ssizey_ext] / maxch;
1829 
1830  else
1831  a = working_space[xmin][i + l + 2 * ssizey_ext] / maxch;
1832  b = a - nip;
1833  if(a + nip <= 0)
1834  a=1;
1835 
1836  else
1837  a=TMath::Sqrt(a + nip);
1838  b = b / a;
1839  b = TMath::Exp(b);
1840  sp = sp + b;
1841  if(i - l + 1 < ymin)
1842  a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1843 
1844  else
1845  a = working_space[xmin][i - l + 1 + 2 * ssizey_ext] / maxch;
1846  b = a - nim;
1847  if(a + nim <= 0)
1848  a = 1;
1849 
1850  else
1851  a=TMath::Sqrt(a + nim);
1852  b = b / a;
1853  b = TMath::Exp(b);
1854  sm = sm + b;
1855  }
1856  a = sp / sm;
1857  a = working_space[xmin][i + 1] = a * working_space[xmin][i];
1858  nom = nom + a;
1859  }
1860  for(i = xmin; i < xmax; i++){
1861  for(j = ymin; j < ymax; j++){
1862  nip = working_space[i][j + 1 + 2 * ssizey_ext] / maxch;
1863  nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1864  spx = 0,smx = 0;
1865  for(l = 1; l <= averWindow; l++){
1866  if(i + l > xmax)
1867  a = working_space[xmax][j + 2 * ssizey_ext] / maxch;
1868 
1869  else
1870  a = working_space[i + l][j + 2 * ssizey_ext] / maxch;
1871  b = a - nip;
1872  if(a + nip <= 0)
1873  a = 1;
1874 
1875  else
1876  a=TMath::Sqrt(a + nip);
1877  b = b / a;
1878  b = TMath::Exp(b);
1879  spx = spx + b;
1880  if(i - l + 1 < xmin)
1881  a = working_space[xmin][j + 2 * ssizey_ext] / maxch;
1882 
1883  else
1884  a = working_space[i - l + 1][j + 2 * ssizey_ext] / maxch;
1885  b = a - nim;
1886  if(a + nim <= 0)
1887  a=1;
1888 
1889  else
1890  a=TMath::Sqrt(a + nim);
1891  b = b / a;
1892  b = TMath::Exp(b);
1893  smx = smx + b;
1894  }
1895  spy = 0,smy = 0;
1896  nip = working_space[i + 1][j + 2 * ssizey_ext] / maxch;
1897  nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1898  for(l = 1; l <= averWindow; l++){
1899  if(j + l > ymax)
1900  a = working_space[i][ymax + 2 * ssizey_ext] / maxch;
1901 
1902  else
1903  a = working_space[i][j + l + 2 * ssizey_ext] / maxch;
1904  b = a - nip;
1905  if(a + nip <= 0)
1906  a = 1;
1907 
1908  else
1909  a=TMath::Sqrt(a + nip);
1910  b = b / a;
1911  b = TMath::Exp(b);
1912  spy = spy + b;
1913  if(j - l + 1 < ymin)
1914  a = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1915 
1916  else
1917  a = working_space[i][j - l + 1 + 2 * ssizey_ext] / maxch;
1918  b=a-nim;
1919  if(a + nim <= 0)
1920  a = 1;
1921  else
1922  a=TMath::Sqrt(a + nim);
1923  b = b / a;
1924  b = TMath::Exp(b);
1925  smy = smy + b;
1926  }
1927  a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx + smy);
1928  working_space[i + 1][j + 1] = a;
1929  nom = nom + a;
1930  }
1931  }
1932  for(i = xmin; i <= xmax; i++){
1933  for(j = ymin; j <= ymax; j++){
1934  working_space[i][j] = working_space[i][j] / nom;
1935  }
1936  }
1937  for(i = 0; i < ssizex_ext; i++){
1938  for(j = 0; j < ssizey_ext; j++){
1939  working_space[i][j + ssizey_ext] = working_space[i][j] * plocha;
1940  working_space[i][2 * ssizey_ext + j] = working_space[i][ssizey_ext + j];
1941  }
1942  }
1943  }
1944  //deconvolution starts
1945  area = 0;
1946  lhx = -1,lhy = -1;
1947  positx = 0,posity = 0;
1948  maximum = 0;
1949  //generate response matrix
1950  for(i = 0; i < ssizex_ext; i++){
1951  for(j = 0; j < ssizey_ext; j++){
1952  lda = (Double_t)i - 3 * sigma;
1953  ldb = (Double_t)j - 3 * sigma;
1954  lda = (lda * lda + ldb * ldb) / (2 * sigma * sigma);
1955  k=(Int_t)(1000 * TMath::Exp(-lda));
1956  lda = k;
1957  if(lda != 0){
1958  if((i + 1) > lhx)
1959  lhx = i + 1;
1960 
1961  if((j + 1) > lhy)
1962  lhy = j + 1;
1963  }
1964  working_space[i][j] = lda;
1965  area = area + lda;
1966  if(lda > maximum){
1967  maximum = lda;
1968  positx = i,posity = j;
1969  }
1970  }
1971  }
1972  //read source matrix
1973  for(i = 0;i < ssizex_ext; i++){
1974  for(j = 0;j < ssizey_ext; j++){
1975  working_space[i][j + 14 * ssizey_ext] = TMath::Abs(working_space[i][j + ssizey_ext]);
1976  }
1977  }
1978  //calculate matrix b=ht*h
1979  i = lhx - 1;
1980  if(i > ssizex_ext)
1981  i = ssizex_ext;
1982 
1983  j = lhy - 1;
1984  if(j>ssizey_ext)
1985  j = ssizey_ext;
1986 
1987  i1min = -i,i1max = i;
1988  i2min = -j,i2max = j;
1989  for(i2 = i2min; i2 <= i2max; i2++){
1990  for(i1 = i1min; i1 <= i1max; i1++){
1991  ldc = 0;
1992  j2min = -i2;
1993  if(j2min < 0)
1994  j2min = 0;
1995 
1996  j2max = lhy - 1 - i2;
1997  if(j2max > lhy - 1)
1998  j2max = lhy - 1;
1999 
2000  for(j2 = j2min; j2 <= j2max; j2++){
2001  j1min = -i1;
2002  if(j1min < 0)
2003  j1min = 0;
2004 
2005  j1max = lhx - 1 - i1;
2006  if(j1max > lhx - 1)
2007  j1max = lhx - 1;
2008 
2009  for(j1 = j1min; j1 <= j1max; j1++){
2010  lda = working_space[j1][j2];
2011  ldb = working_space[i1 + j1][i2 + j2];
2012  ldc = ldc + lda * ldb;
2013  }
2014  }
2015  k = (i1 + ssizex_ext) / ssizex_ext;
2016  working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext] = ldc;
2017  }
2018  }
2019  //calculate at*y and write into p
2020  i = lhx - 1;
2021  if(i > ssizex_ext)
2022  i = ssizex_ext;
2023 
2024  j = lhy - 1;
2025  if(j > ssizey_ext)
2026  j = ssizey_ext;
2027 
2028  i2min = -j,i2max = ssizey_ext + j - 1;
2029  i1min = -i,i1max = ssizex_ext + i - 1;
2030  for(i2 = i2min; i2 <= i2max; i2++){
2031  for(i1=i1min;i1<=i1max;i1++){
2032  ldc=0;
2033  for(j2 = 0; j2 <= (lhy - 1); j2++){
2034  for(j1 = 0; j1 <= (lhx - 1); j1++){
2035  k2 = i2 + j2,k1 = i1 + j1;
2036  if(k2 >= 0 && k2 < ssizey_ext && k1 >= 0 && k1 < ssizex_ext){
2037  lda = working_space[j1][j2];
2038  ldb = working_space[k1][k2 + 14 * ssizey_ext];
2039  ldc = ldc + lda * ldb;
2040  }
2041  }
2042  }
2043  k = (i1 + ssizex_ext) / ssizex_ext;
2044  working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext] = ldc;
2045  }
2046  }
2047  //move matrix p
2048  for(i2 = 0; i2 < ssizey_ext; i2++){
2049  for(i1 = 0; i1 < ssizex_ext; i1++){
2050  k = (i1 + ssizex_ext) / ssizex_ext;
2051  ldb = working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext];
2052  working_space[i1][i2 + 14 * ssizey_ext] = ldb;
2053  }
2054  }
2055  //initialization in x1 matrix
2056  for(i2 = 0; i2 < ssizey_ext; i2++){
2057  for(i1 = 0; i1 < ssizex_ext; i1++){
2058  working_space[i1][i2 + ssizey_ext] = 1;
2059  working_space[i1][i2 + 2 * ssizey_ext] = 0;
2060  }
2061  }
2062  //START OF ITERATIONS
2063  for(lindex = 0; lindex < deconIterations; lindex++){
2064  for(i2 = 0; i2 < ssizey_ext; i2++){
2065  for(i1 = 0; i1 < ssizex_ext; i1++){
2066  lda = working_space[i1][i2 + ssizey_ext];
2067  ldc = working_space[i1][i2 + 14 * ssizey_ext];
2068  if(lda > 0.000001 && ldc > 0.000001){
2069  ldb=0;
2070  j2min=i2;
2071  if(j2min > lhy - 1)
2072  j2min = lhy - 1;
2073 
2074  j2min = -j2min;
2075  j2max = ssizey_ext - i2 - 1;
2076  if(j2max > lhy - 1)
2077  j2max = lhy - 1;
2078 
2079  j1min = i1;
2080  if(j1min > lhx - 1)
2081  j1min = lhx - 1;
2082 
2083  j1min = -j1min;
2084  j1max = ssizex_ext - i1 - 1;
2085  if(j1max > lhx - 1)
2086  j1max = lhx - 1;
2087 
2088  for(j2 = j2min; j2 <= j2max; j2++){
2089  for(j1 = j1min; j1 <= j1max; j1++){
2090  k = (j1 + ssizex_ext) / ssizex_ext;
2091  ldc = working_space[(j1 + ssizex_ext) % ssizex_ext][j2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext];
2092  lda = working_space[i1 + j1][i2 + j2 + ssizey_ext];
2093  ldb = ldb + lda * ldc;
2094  }
2095  }
2096  lda = working_space[i1][i2 + ssizey_ext];
2097  ldc = working_space[i1][i2 + 14 * ssizey_ext];
2098  if(ldc * lda != 0 && ldb != 0){
2099  lda =lda * ldc / ldb;
2100  }
2101 
2102  else
2103  lda=0;
2104  working_space[i1][i2 + 2 * ssizey_ext] = lda;
2105  }
2106  }
2107  }
2108  for(i2 = 0; i2 < ssizey_ext; i2++){
2109  for(i1 = 0; i1 < ssizex_ext; i1++)
2110  working_space[i1][i2 + ssizey_ext] = working_space[i1][i2 + 2 * ssizey_ext];
2111  }
2112  }
2113  //looking for maximum
2114  maximum=0;
2115  for(i = 0; i < ssizex_ext; i++){
2116  for(j = 0; j < ssizey_ext; j++){
2117  working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext] = area * working_space[i][j + ssizey_ext];
2118  if(maximum < working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext])
2119  maximum = working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext];
2120 
2121  }
2122  }
2123  //searching for peaks in deconvolved spectrum
2124  for(i = 1; i < ssizex_ext - 1; i++){
2125  for(j = 1; j < ssizey_ext - 1; j++){
2126  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]){
2127  if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift){
2128  if(working_space[i][j] > threshold * maximum / 100.0){
2129  if(peak_index < fMaxPeaks){
2130  for(k = i - 1,a = 0,b = 0; k <= i + 1; k++){
2131  a += (Double_t)(k - shift) * working_space[k][j];
2132  b += working_space[k][j];
2133  }
2134  ax=a/b;
2135  if(ax < 0)
2136  ax = 0;
2137 
2138  if(ax >= ssizex)
2139  ax = ssizex - 1;
2140 
2141  for(k = j - 1,a = 0,b = 0; k <= j + 1; k++){
2142  a += (Double_t)(k - shift) * working_space[i][k];
2143  b += working_space[i][k];
2144  }
2145  ay=a/b;
2146  if(ay < 0)
2147  ay = 0;
2148 
2149  if(ay >= ssizey)
2150  ay = ssizey - 1;
2151 
2152  if(peak_index == 0){
2153  fPositionX[0] = ax;
2154  fPositionY[0] = ay;
2155  peak_index = 1;
2156  }
2157 
2158  else{
2159  for(k = 0, priz = 0; k < peak_index && priz == 0; k++){
2160  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)])
2161  priz = 1;
2162  }
2163  if(priz == 0){
2164  if(k < fMaxPeaks){
2165  fPositionX[k] = ax;
2166  fPositionY[k] = ay;
2167  }
2168  }
2169 
2170  else{
2171  for(l = peak_index; l >= k; l--){
2172  if(l < fMaxPeaks){
2173  fPositionX[l] = fPositionX[l - 1];
2174  fPositionY[l] = fPositionY[l - 1];
2175  }
2176  }
2177  fPositionX[k - 1] = ax;
2178  fPositionY[k - 1] = ay;
2179  }
2180  if(peak_index < fMaxPeaks)
2181  peak_index += 1;
2182  }
2183  }
2184  }
2185  }
2186  }
2187  }
2188  }
2189  //writing back deconvolved spectrum
2190  for(i = 0; i < ssizex; i++){
2191  for(j = 0; j < ssizey; j++){
2192  dest[i][j] = working_space[i + shift][j + shift];
2193  }
2194  }
2195  i = (Int_t)(4 * sigma + 0.5);
2196  i = 4 * i;
2197  for (j = 0; j < ssizex + i; j++)
2198  delete[]working_space[j];
2199  delete[]working_space;
2200  fNPeaks = peak_index;
2201  return fNPeaks;
2202 }
2203 
2204 ////////////////////////////////////////////////////////////////////////////////
2205 /// static function (called by TH1), interface to TSpectrum2::Search
2206 
2208 {
2209  TSpectrum2 s;
2210  return s.Search(hist,sigma,option,threshold);
2211 }
2212 
2213 ////////////////////////////////////////////////////////////////////////////////
2214 /// static function (called by TH1), interface to TSpectrum2::Background
2215 
2216 TH1 *TSpectrum2::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
2217 {
2218  TSpectrum2 s;
2219  return s.Background(hist,niter,option);
2220 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
const char * Deconvolution(Double_t **source, Double_t **resp, Int_t ssizex, Int_t ssizey, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
This function calculates deconvolution from source spectrum according to response spectrum The result...
float xmin
Definition: THbookFile.cxx:93
static double p3(double t, double a, double b, double c, double d)
Double_t * fPositionY
[fNPeaks] Y position of peaks
Definition: TSpectrum2.h:24
const char Option_t
Definition: RtypesCore.h:62
Definition: Rtypes.h:59
float ymin
Definition: THbookFile.cxx:93
Advanced 2-dimensional spectra processing.
Definition: TSpectrum2.h:18
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4770
static Int_t fgAverageWindow
Average window of searched peaks.
Definition: TSpectrum2.h:27
TH1 * fHistogram
resulting histogram
Definition: TSpectrum2.h:26
Basic string class.
Definition: TString.h:131
static void SetDeconIterations(Int_t n=3)
static function: Set max number of decon iterations in deconvolution operation see TSpectrum2::Search...
Definition: TSpectrum2.cxx:113
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1100
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Double_t fResolution
NOT USED resolution of the neighboring peaks
Definition: TSpectrum2.h:25
void SetResolution(Double_t resolution=1)
NOT USED resolution: determines resolution of the neighboring peaks default value is 1 correspond to ...
Definition: TSpectrum2.cxx:287
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:734
virtual TObject * FindObject(const char *name) const
Delete a TObjLink object.
Definition: TList.cxx:574
virtual Int_t GetDimension() const
Definition: TH1.h:277
Double_t x[n]
Definition: legend1.C:17
Int_t fNPeaks
number of peaks found
Definition: TSpectrum2.h:21
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
static double p2(double t, double a, double b, double c)
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
virtual void Print(Option_t *option="") const
Print the array of positions.
Definition: TSpectrum2.cxx:166
const Double_t sigma
Int_t fMaxPeaks
Maximum number of peaks to be found.
Definition: TSpectrum2.h:20
Double_t * fPositionX
[fNPeaks] X position of peaks
Definition: TSpectrum2.h:23
float ymax
Definition: THbookFile.cxx:93
auto * a
Definition: textangle.C:12
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:818
static void SetAverageWindow(Int_t w=3)
static function: Set average window of searched peaks see TSpectrum2::SearchHighRes ...
Definition: TSpectrum2.cxx:104
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only re...
Definition: VectorUtil.h:327
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
#define s1(x)
Definition: RSha256.hxx:91
TSpectrum2()
Constructor.
Definition: TSpectrum2.cxx:57
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
TAxis * GetYaxis()
Definition: TH1.h:316
float xmax
Definition: THbookFile.cxx:93
static double p1(double t, double a, double b)
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
#define h(i)
Definition: RSha256.hxx:106
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t Exp(Double_t x)
Definition: TMath.h:726
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
static function (called by TH1), interface to TSpectrum2::Search
A PolyMarker is defined by an array on N points in a 2-D space.
Definition: TPolyMarker.h:31
Double_t * fPosition
[fNPeaks] array of current peak positions
Definition: TSpectrum2.h:22
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
virtual ~TSpectrum2()
Destructor.
Definition: TSpectrum2.cxx:92
Double_t y[n]
Definition: legend1.C:17
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
The TH1 histogram class.
Definition: TH1.h:56
static constexpr double s
virtual void Add(TObject *obj)
Definition: TList.h:87
auto * l
Definition: textangle.C:4
const char * SmoothMarkov(Double_t **source, Int_t ssizex, Int_t ssizey, Int_t averWindow)
This function calculates smoothed spectrum from source spectrum based on Markov chain method...
Definition: TSpectrum2.cxx:736
#define dest(otri, vertexptr)
Definition: triangle.c:1040
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
static TH1 * StaticBackground(const TH1 *hist, Int_t niter=20, Option_t *option="")
static function (called by TH1), interface to TSpectrum2::Background
virtual TH1 * Background(const TH1 *hist, Int_t niter=20, Option_t *option="")
This function calculates the background spectrum in the input histogram h.
Definition: TSpectrum2.cxx:155
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
Int_t GetNbins() const
Definition: TAxis.h:121
TList * GetListOfFunctions() const
Definition: TH1.h:238
const Bool_t kTRUE
Definition: RtypesCore.h:87
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)
This function searches for peaks in source spectrum It is based on deconvolution method.
const Int_t n
Definition: legend1.C:16
#define PEAK_WINDOW
Definition: TSpectrum2.cxx:47
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:315
static Int_t fgIterations
Maximum number of decon iterations (default=3)
Definition: TSpectrum2.h:28
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05)
This function searches for peaks in source spectrum in hin The number of found peaks and their positi...
Definition: TSpectrum2.cxx:206