Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
49 Int_t TSpectrum2::fgIterations = 3;
50 Int_t TSpectrum2::fgAverageWindow = 3;
51 
52 ClassImp(TSpectrum2);
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 
92 TSpectrum2::~TSpectrum2()
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 
104 void TSpectrum2::SetAverageWindow(Int_t w)
105 {
106  fgAverageWindow = w;
107 }
108 
109 ////////////////////////////////////////////////////////////////////////////////
110 /// static function: Set max number of decon iterations in deconvolution operation
111 /// see TSpectrum2::SearchHighRes
112 
113 void TSpectrum2::SetDeconIterations(Int_t n)
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 
166 void TSpectrum2::Print(Option_t *) const
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 
206 Int_t TSpectrum2::Search(const TH1 * hin, Double_t sigma,
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();
219  Bool_t background = kTRUE;
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 
287 void TSpectrum2::SetResolution(Double_t resolution)
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 - Background_gamma64.C
350 ///
351 /// Begin_Macro(source)
352 /// ../../../tutorials/spectrum/Background_gamma64.C
353 /// End_Macro
354 ///
355 /// ### Example 2- Background_gamma256.C
356 ///
357 /// Begin_Macro(source)
358 /// ../../../tutorials/spectrum/Background_gamma256.C
359 /// End_Macro
360 ///
361 /// ### Example 3- Background_synt256.C
362 ///
363 /// Begin_Macro(source)
364 /// ../../../tutorials/spectrum/Background_synt256.C
365 /// End_Macro
366 
367 const char *TSpectrum2::Background(Double_t **spectrum,
368  Int_t ssizex, Int_t ssizey,
369  Int_t numberIterationsX,
370  Int_t numberIterationsY,
371  Int_t direction,
372  Int_t filterType)
373 {
374  Int_t i, x, y, sampling, r1, r2;
375  Double_t a, b, p1, p2, p3, p4, s1, s2, s3, s4;
376  if (ssizex <= 0 || ssizey <= 0)
377  return "Wrong parameters";
378  if (numberIterationsX < 1 || numberIterationsY < 1)
379  return "Width of Clipping Window Must Be Positive";
380  if (ssizex < 2 * numberIterationsX + 1
381  || ssizey < 2 * numberIterationsY + 1)
382  return ("Too Large Clipping Window");
383  Double_t **working_space = new Double_t*[ssizex];
384  for (i = 0; i < ssizex; i++)
385  working_space[i] = new Double_t[ssizey];
386  sampling =
387  (Int_t) TMath::Max(numberIterationsX, numberIterationsY);
388  if (direction == kBackIncreasingWindow) {
389  if (filterType == kBackSuccessiveFiltering) {
390  for (i = 1; i <= sampling; i++) {
391  r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
392  (Int_t) TMath::Min(i, numberIterationsY);
393  for (y = r2; y < ssizey - r2; y++) {
394  for (x = r1; x < ssizex - r1; x++) {
395  a = spectrum[x][y];
396  p1 = spectrum[x - r1][y - r2];
397  p2 = spectrum[x - r1][y + r2];
398  p3 = spectrum[x + r1][y - r2];
399  p4 = spectrum[x + r1][y + r2];
400  s1 = spectrum[x][y - r2];
401  s2 = spectrum[x - r1][y];
402  s3 = spectrum[x + r1][y];
403  s4 = spectrum[x][y + r2];
404  b = (p1 + p2) / 2.0;
405  if (b > s2)
406  s2 = b;
407  b = (p1 + p3) / 2.0;
408  if (b > s1)
409  s1 = b;
410  b = (p2 + p4) / 2.0;
411  if (b > s4)
412  s4 = b;
413  b = (p3 + p4) / 2.0;
414  if (b > s3)
415  s3 = b;
416  s1 = s1 - (p1 + p3) / 2.0;
417  s2 = s2 - (p1 + p2) / 2.0;
418  s3 = s3 - (p3 + p4) / 2.0;
419  s4 = s4 - (p2 + p4) / 2.0;
420  b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
421  p3 +
422  p4) / 4.0;
423  if (b < a && b > 0)
424  a = b;
425  working_space[x][y] = a;
426  }
427  }
428  for (y = r2; y < ssizey - r2; y++) {
429  for (x = r1; x < ssizex - r1; x++) {
430  spectrum[x][y] = working_space[x][y];
431  }
432  }
433  }
434  }
435 
436  else if (filterType == kBackOneStepFiltering) {
437  for (i = 1; i <= sampling; i++) {
438  r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
439  (Int_t) TMath::Min(i, numberIterationsY);
440  for (y = r2; y < ssizey - r2; y++) {
441  for (x = r1; x < ssizex - r1; x++) {
442  a = spectrum[x][y];
443  b = -(spectrum[x - r1][y - r2] +
444  spectrum[x - r1][y + r2] + spectrum[x + r1][y -
445  r2]
446  + spectrum[x + r1][y + r2]) / 4 +
447  (spectrum[x][y - r2] + spectrum[x - r1][y] +
448  spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
449  if (b < a && b > 0)
450  a = b;
451  working_space[x][y] = a;
452  }
453  }
454  for (y = i; y < ssizey - i; y++) {
455  for (x = i; x < ssizex - i; x++) {
456  spectrum[x][y] = working_space[x][y];
457  }
458  }
459  }
460  }
461  }
462 
463  else if (direction == kBackDecreasingWindow) {
464  if (filterType == kBackSuccessiveFiltering) {
465  for (i = sampling; i >= 1; i--) {
466  r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
467  (Int_t) TMath::Min(i, numberIterationsY);
468  for (y = r2; y < ssizey - r2; y++) {
469  for (x = r1; x < ssizex - r1; x++) {
470  a = spectrum[x][y];
471  p1 = spectrum[x - r1][y - r2];
472  p2 = spectrum[x - r1][y + r2];
473  p3 = spectrum[x + r1][y - r2];
474  p4 = spectrum[x + r1][y + r2];
475  s1 = spectrum[x][y - r2];
476  s2 = spectrum[x - r1][y];
477  s3 = spectrum[x + r1][y];
478  s4 = spectrum[x][y + r2];
479  b = (p1 + p2) / 2.0;
480  if (b > s2)
481  s2 = b;
482  b = (p1 + p3) / 2.0;
483  if (b > s1)
484  s1 = b;
485  b = (p2 + p4) / 2.0;
486  if (b > s4)
487  s4 = b;
488  b = (p3 + p4) / 2.0;
489  if (b > s3)
490  s3 = b;
491  s1 = s1 - (p1 + p3) / 2.0;
492  s2 = s2 - (p1 + p2) / 2.0;
493  s3 = s3 - (p3 + p4) / 2.0;
494  s4 = s4 - (p2 + p4) / 2.0;
495  b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
496  p3 +
497  p4) / 4.0;
498  if (b < a && b > 0)
499  a = b;
500  working_space[x][y] = a;
501  }
502  }
503  for (y = r2; y < ssizey - r2; y++) {
504  for (x = r1; x < ssizex - r1; x++) {
505  spectrum[x][y] = working_space[x][y];
506  }
507  }
508  }
509  }
510 
511  else if (filterType == kBackOneStepFiltering) {
512  for (i = sampling; i >= 1; i--) {
513  r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
514  (Int_t) TMath::Min(i, numberIterationsY);
515  for (y = r2; y < ssizey - r2; y++) {
516  for (x = r1; x < ssizex - r1; x++) {
517  a = spectrum[x][y];
518  b = -(spectrum[x - r1][y - r2] +
519  spectrum[x - r1][y + r2] + spectrum[x + r1][y -
520  r2]
521  + spectrum[x + r1][y + r2]) / 4 +
522  (spectrum[x][y - r2] + spectrum[x - r1][y] +
523  spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
524  if (b < a && b > 0)
525  a = b;
526  working_space[x][y] = a;
527  }
528  }
529  for (y = i; y < ssizey - i; y++) {
530  for (x = i; x < ssizex - i; x++) {
531  spectrum[x][y] = working_space[x][y];
532  }
533  }
534  }
535  }
536  }
537  for (i = 0; i < ssizex; i++)
538  delete[]working_space[i];
539  delete[]working_space;
540  return 0;
541 }
542 
543 ////////////////////////////////////////////////////////////////////////////////
544 /// This function calculates smoothed spectrum from source spectrum
545 /// based on Markov chain method.
546 /// The result is placed in the array pointed by source pointer.
547 ///
548 /// Function parameters:
549 /// - source-pointer to the array of source spectrum
550 /// - ssizex-x length of source
551 /// - ssizey-y length of source
552 /// - averWindow-width of averaging smoothing window
553 ///
554 /// ### Smoothing
555 ///
556 /// Goal: Suppression of statistical fluctuations the algorithm is based on discrete
557 /// Markov chain, which has very simple invariant distribution
558 /// \f[
559 /// 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
560 /// \f]
561 /// \f$U_1\f$ being defined from the normalization condition \f$ \sum_{i=1}^{n} U_i = 1\f$
562 /// n is the length of the smoothed spectrum and
563 /// \f[
564 /// 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]
565 /// \f]
566 /// is the probability of the change of the peak position from channel i to the channel i+1.
567 /// \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
568 /// of smoothing window. We have extended this algorithm to two dimensions.
569 ///
570 /// #### Reference:
571 ///
572 /// [1] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 (1996), 451.
573 ///
574 /// ### Example 4 - Smooth.C
575 ///
576 /// Begin_Macro(source)
577 /// ../../../tutorials/spectrum/Smooth.C
578 /// End_Macro
579 
580 const char* TSpectrum2::SmoothMarkov(Double_t **source, Int_t ssizex, Int_t ssizey, Int_t averWindow)
581 {
582 
583  Int_t xmin, xmax, ymin, ymax, i, j, l;
584  Double_t a, b, maxch;
585  Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy, plocha = 0;
586  if(averWindow <= 0)
587  return "Averaging Window must be positive";
588  Double_t **working_space = new Double_t*[ssizex];
589  for(i = 0; i < ssizex; i++)
590  working_space[i] = new Double_t[ssizey];
591  xmin = 0;
592  xmax = ssizex - 1;
593  ymin = 0;
594  ymax = ssizey - 1;
595  for(i = 0, maxch = 0; i < ssizex; i++){
596  for(j = 0; j < ssizey; j++){
597  working_space[i][j] = 0;
598  if(maxch < source[i][j])
599  maxch = source[i][j];
600 
601  plocha += source[i][j];
602  }
603  }
604  if(maxch == 0) {
605  delete [] working_space;
606  return 0;
607  }
608 
609  nom = 0;
610  working_space[xmin][ymin] = 1;
611  for(i = xmin; i < xmax; i++){
612  nip = source[i][ymin] / maxch;
613  nim = source[i + 1][ymin] / maxch;
614  sp = 0,sm = 0;
615  for(l = 1; l <= averWindow; l++){
616  if((i + l) > xmax)
617  a = source[xmax][ymin] / maxch;
618 
619  else
620  a = source[i + l][ymin] / maxch;
621  b = a - nip;
622  if(a + nip <= 0)
623  a = 1;
624 
625  else
626  a = TMath::Sqrt(a + nip);
627  b = b / a;
628  b = TMath::Exp(b);
629  sp = sp + b;
630  if(i - l + 1 < xmin)
631  a = source[xmin][ymin] / maxch;
632 
633  else
634  a = source[i - l + 1][ymin] / maxch;
635  b = a - nim;
636  if(a + nim <= 0)
637  a = 1;
638 
639  else
640  a = TMath::Sqrt(a + nim);
641  b = b / a;
642  b = TMath::Exp(b);
643  sm = sm + b;
644  }
645  a = sp / sm;
646  a = working_space[i + 1][ymin] = a * working_space[i][ymin];
647  nom = nom + a;
648  }
649  for(i = ymin; i < ymax; i++){
650  nip = source[xmin][i] / maxch;
651  nim = source[xmin][i + 1] / maxch;
652  sp = 0,sm = 0;
653  for(l = 1; l <= averWindow; l++){
654  if((i + l) > ymax)
655  a = source[xmin][ymax] / maxch;
656 
657  else
658  a = source[xmin][i + l] / maxch;
659  b = a - nip;
660  if(a + nip <= 0)
661  a = 1;
662 
663  else
664  a = TMath::Sqrt(a + nip);
665  b = b / a;
666  b = TMath::Exp(b);
667  sp = sp + b;
668  if(i - l + 1 < ymin)
669  a = source[xmin][ymin] / maxch;
670 
671  else
672  a = source[xmin][i - l + 1] / maxch;
673  b = a - nim;
674  if(a + nim <= 0)
675  a = 1;
676 
677  else
678  a = TMath::Sqrt(a + nim);
679  b = b / a;
680  b = TMath::Exp(b);
681  sm = sm + b;
682  }
683  a = sp / sm;
684  a = working_space[xmin][i + 1] = a * working_space[xmin][i];
685  nom = nom + a;
686  }
687  for(i = xmin; i < xmax; i++){
688  for(j = ymin; j < ymax; j++){
689  nip = source[i][j + 1] / maxch;
690  nim = source[i + 1][j + 1] / maxch;
691  spx = 0,smx = 0;
692  for(l = 1; l <= averWindow; l++){
693  if(i + l > xmax)
694  a = source[xmax][j] / maxch;
695 
696  else
697  a = source[i + l][j] / maxch;
698  b = a - nip;
699  if(a + nip <= 0)
700  a = 1;
701 
702  else
703  a = TMath::Sqrt(a + nip);
704  b = b / a;
705  b = TMath::Exp(b);
706  spx = spx + b;
707  if(i - l + 1 < xmin)
708  a = source[xmin][j] / maxch;
709 
710  else
711  a = source[i - l + 1][j] / maxch;
712  b = a - nim;
713  if(a + nim <= 0)
714  a = 1;
715 
716  else
717  a = TMath::Sqrt(a + nim);
718  b = b / a;
719  b = TMath::Exp(b);
720  smx = smx + b;
721  }
722  spy = 0,smy = 0;
723  nip = source[i + 1][j] / maxch;
724  nim = source[i + 1][j + 1] / maxch;
725  for (l = 1; l <= averWindow; l++) {
726  if (j + l > ymax) a = source[i][ymax]/maxch;
727  else a = source[i][j + l] / maxch;
728  b = a - nip;
729  if (a + nip <= 0) a = 1;
730  else a = TMath::Sqrt(a + nip);
731  b = b / a;
732  b = TMath::Exp(b);
733  spy = spy + b;
734  if (j - l + 1 < ymin) a = source[i][ymin] / maxch;
735  else a = source[i][j - l + 1] / maxch;
736  b = a - nim;
737  if (a + nim <= 0) a = 1;
738  else a = TMath::Sqrt(a + nim);
739  b = b / a;
740  b = TMath::Exp(b);
741  smy = smy + b;
742  }
743  a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx +smy);
744  working_space[i + 1][j + 1] = a;
745  nom = nom + a;
746  }
747  }
748  for(i = xmin; i <= xmax; i++){
749  for(j = ymin; j <= ymax; j++){
750  working_space[i][j] = working_space[i][j] / nom;
751  }
752  }
753  for(i = 0;i < ssizex; i++){
754  for(j = 0; j < ssizey; j++){
755  source[i][j] = plocha * working_space[i][j];
756  }
757  }
758  for (i = 0; i < ssizex; i++)
759  delete[]working_space[i];
760  delete[]working_space;
761  return 0;
762 }
763 
764 ////////////////////////////////////////////////////////////////////////////////
765 /// This function calculates deconvolution from source spectrum
766 /// according to response spectrum
767 /// The result is placed in the matrix pointed by source pointer.
768 ///
769 /// Function parameters:
770 /// - source-pointer to the matrix of source spectrum
771 /// - resp-pointer to the matrix of response spectrum
772 /// - ssizex-x length of source and response spectra
773 /// - ssizey-y length of source and response spectra
774 /// - numberIterations, for details we refer to manual
775 /// - numberRepetitions, for details we refer to manual
776 /// - boost, boosting factor, for details we refer to manual
777 ///
778 /// ### Deconvolution
779 ///
780 /// Goal: Improvement of the resolution in spectra, decomposition of multiplets
781 ///
782 /// Mathematical formulation of the 2-dimensional convolution system is
783 ///
784 ///\f[
785 /// 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)
786 ///\f]
787 ///\f[
788 /// i_1 = 0,1,2, ... ,N_1-1, i_2 = 0,1,2, ... ,N_2-1
789 ///\f]
790 ///
791 /// where h(i,j) is the impulse response function, x, y are input and output
792 /// matrices, respectively, \f$ N_1, N_2\f$ are the lengths of x and h matrices
793 ///
794 /// - let us assume that we know the response and the output matrices (spectra) of the above given system.
795 /// - the deconvolution represents solution of the overdetermined system of linear equations, i.e., the
796 /// calculation of the matrix x.
797 /// - from numerical stability point of view the operation of deconvolution is
798 /// extremely critical (ill-posed problem) as well as time consuming operation.
799 /// - the Gold deconvolution algorithm proves to work very well even for 2-dimensional
800 /// systems. Generalisation of the algorithm for 2-dimensional systems was presented in [1], [2].
801 /// - for Gold deconvolution algorithm as well as for boosted deconvolution algorithm we refer also to TSpectrum
802 ///
803 /// #### References:
804 ///
805 /// [1] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:
806 /// Efficient one- and two-dimensional Gold deconvolution and its application to
807 /// gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
808 ///
809 /// [2] Morhac; M., Matouoek V., Kliman J., Efficient algorithm of multidimensional
810 /// deconvolution and its application to nuclear data processing, Digital Signal
811 /// Processing 13 (2003) 144.
812 ///
813 /// ### Example 5 - Deconvolution2_1.c
814 ///
815 /// Begin_Macro(source)
816 /// ../../../tutorials/spectrum/Deconvolution2_1.C
817 /// End_Macro
818 ///
819 /// ### Example 6 - Deconvolution2_2.C
820 ///
821 /// Begin_Macro(source)
822 /// ../../../tutorials/spectrum/Deconvolution2_2.C
823 /// End_Macro
824 ///
825 /// ### Example 7 - Deconvolution2_HR.C
826 ///
827 /// Begin_Macro(source)
828 /// ../../../tutorials/spectrum/Deconvolution2_HR.C
829 /// End_Macro
830 
831 const char *TSpectrum2::Deconvolution(Double_t **source, Double_t **resp,
832  Int_t ssizex, Int_t ssizey,
833  Int_t numberIterations,
834  Int_t numberRepetitions,
835  Double_t boost)
836 {
837  Int_t i, j, lhx, lhy, i1, i2, j1, j2, k1, k2, lindex, i1min, i1max,
838  i2min, i2max, j1min, j1max, j2min, j2max, positx = 0, posity = 0, repet;
839  Double_t lda, ldb, ldc, area, maximum = 0;
840  if (ssizex <= 0 || ssizey <= 0)
841  return "Wrong parameters";
842  if (numberIterations <= 0)
843  return "Number of iterations must be positive";
844  if (numberRepetitions <= 0)
845  return "Number of repetitions must be positive";
846  Double_t **working_space = new Double_t *[ssizex];
847  for (i = 0; i < ssizex; i++)
848  working_space[i] = new Double_t[5 * ssizey];
849  area = 0;
850  lhx = -1, lhy = -1;
851  for (i = 0; i < ssizex; i++) {
852  for (j = 0; j < ssizey; j++) {
853  lda = resp[i][j];
854  if (lda != 0) {
855  if ((i + 1) > lhx)
856  lhx = i + 1;
857  if ((j + 1) > lhy)
858  lhy = j + 1;
859  }
860  working_space[i][j] = lda;
861  area = area + lda;
862  if (lda > maximum) {
863  maximum = lda;
864  positx = i, posity = j;
865  }
866  }
867  }
868  if (lhx == -1 || lhy == -1) {
869  delete [] working_space;
870  return ("Zero response data");
871  }
872 
873 //calculate ht*y and write into p
874  for (i2 = 0; i2 < ssizey; i2++) {
875  for (i1 = 0; i1 < ssizex; i1++) {
876  ldc = 0;
877  for (j2 = 0; j2 <= (lhy - 1); j2++) {
878  for (j1 = 0; j1 <= (lhx - 1); j1++) {
879  k2 = i2 + j2, k1 = i1 + j1;
880  if (k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
881  lda = working_space[j1][j2];
882  ldb = source[k1][k2];
883  ldc = ldc + lda * ldb;
884  }
885  }
886  }
887  working_space[i1][i2 + ssizey] = ldc;
888  }
889  }
890 
891 //calculate matrix b=ht*h
892  i1min = -(lhx - 1), i1max = lhx - 1;
893  i2min = -(lhy - 1), i2max = lhy - 1;
894  for (i2 = i2min; i2 <= i2max; i2++) {
895  for (i1 = i1min; i1 <= i1max; i1++) {
896  ldc = 0;
897  j2min = -i2;
898  if (j2min < 0)
899  j2min = 0;
900  j2max = lhy - 1 - i2;
901  if (j2max > lhy - 1)
902  j2max = lhy - 1;
903  for (j2 = j2min; j2 <= j2max; j2++) {
904  j1min = -i1;
905  if (j1min < 0)
906  j1min = 0;
907  j1max = lhx - 1 - i1;
908  if (j1max > lhx - 1)
909  j1max = lhx - 1;
910  for (j1 = j1min; j1 <= j1max; j1++) {
911  lda = working_space[j1][j2];
912  if (i1 + j1 < ssizex && i2 + j2 < ssizey)
913  ldb = working_space[i1 + j1][i2 + j2];
914  else
915  ldb = 0;
916  ldc = ldc + lda * ldb;
917  }
918  }
919  working_space[i1 - i1min][i2 - i2min + 2 * ssizey ] = ldc;
920  }
921  }
922 
923 //initialization in x1 matrix
924  for (i2 = 0; i2 < ssizey; i2++) {
925  for (i1 = 0; i1 < ssizex; i1++) {
926  working_space[i1][i2 + 3 * ssizey] = 1;
927  working_space[i1][i2 + 4 * ssizey] = 0;
928  }
929  }
930 
931  //START OF ITERATIONS
932  for (repet = 0; repet < numberRepetitions; repet++) {
933  if (repet != 0) {
934  for (i = 0; i < ssizex; i++) {
935  for (j = 0; j < ssizey; j++) {
936  working_space[i][j + 3 * ssizey] =
937  TMath::Power(working_space[i][j + 3 * ssizey], boost);
938  }
939  }
940  }
941  for (lindex = 0; lindex < numberIterations; lindex++) {
942  for (i2 = 0; i2 < ssizey; i2++) {
943  for (i1 = 0; i1 < ssizex; i1++) {
944  ldb = 0;
945  j2min = i2;
946  if (j2min > lhy - 1)
947  j2min = lhy - 1;
948  j2min = -j2min;
949  j2max = ssizey - i2 - 1;
950  if (j2max > lhy - 1)
951  j2max = lhy - 1;
952  j1min = i1;
953  if (j1min > lhx - 1)
954  j1min = lhx - 1;
955  j1min = -j1min;
956  j1max = ssizex - i1 - 1;
957  if (j1max > lhx - 1)
958  j1max = lhx - 1;
959  for (j2 = j2min; j2 <= j2max; j2++) {
960  for (j1 = j1min; j1 <= j1max; j1++) {
961  ldc = working_space[j1 - i1min][j2 - i2min + 2 * ssizey];
962  lda = working_space[i1 + j1][i2 + j2 + 3 * ssizey];
963  ldb = ldb + lda * ldc;
964  }
965  }
966  lda = working_space[i1][i2 + 3 * ssizey];
967  ldc = working_space[i1][i2 + 1 * ssizey];
968  if (ldc * lda != 0 && ldb != 0) {
969  lda = lda * ldc / ldb;
970  }
971 
972  else
973  lda = 0;
974  working_space[i1][i2 + 4 * ssizey] = lda;
975  }
976  }
977  for (i2 = 0; i2 < ssizey; i2++) {
978  for (i1 = 0; i1 < ssizex; i1++)
979  working_space[i1][i2 + 3 * ssizey] =
980  working_space[i1][i2 + 4 * ssizey];
981  }
982  }
983  }
984  for (i = 0; i < ssizex; i++) {
985  for (j = 0; j < ssizey; j++)
986  source[(i + positx) % ssizex][(j + posity) % ssizey] =
987  area * working_space[i][j + 3 * ssizey];
988  }
989  for (i = 0; i < ssizex; i++)
990  delete[]working_space[i];
991  delete[]working_space;
992  return 0;
993 }
994 
995 ////////////////////////////////////////////////////////////////////////////////
996 /// This function searches for peaks in source spectrum
997 /// It is based on deconvolution method. First the background is
998 /// removed (if desired), then Markov spectrum is calculated
999 /// (if desired), then the response function is generated
1000 /// according to given sigma and deconvolution is carried out.
1001 ///
1002 /// Function parameters:
1003 /// - source-pointer to the matrix of source spectrum
1004 /// - dest-pointer to the matrix of resulting deconvolved spectrum
1005 /// - ssizex-x length of source spectrum
1006 /// - ssizey-y length of source spectrum
1007 /// - sigma-sigma of searched peaks, for details we refer to manual
1008 /// - threshold-threshold value in % for selected peaks, peaks with
1009 /// amplitude less than threshold*highest_peak/100
1010 /// are ignored, see manual
1011 /// - backgroundRemove-logical variable, set if the removal of
1012 /// background before deconvolution is desired
1013 /// - deconIterations-number of iterations in deconvolution operation
1014 /// - markov-logical variable, if it is true, first the source spectrum
1015 /// is replaced by new spectrum calculated using Markov
1016 /// chains method.
1017 /// - averWindow-averaging window of searched peaks, for details
1018 /// we refer to manual (applies only for Markov method)
1019 ///
1020 /// ### Peaks searching
1021 ///
1022 /// Goal: to identify automatically the peaks in spectrum with the presence of the
1023 /// continuous background, one-fold coincidences (ridges) and statistical
1024 /// fluctuations - noise.
1025 ///
1026 /// The common problems connected with correct peak identification in two-dimensional coincidence spectra are
1027 ///
1028 /// - non-sensitivity to noise, i.e., only statistically relevant peaks should be identified
1029 /// - non-sensitivity of the algorithm to continuous background
1030 /// - non-sensitivity to one-fold coincidences (coincidences peak - background in both dimensions) and their crossings
1031 /// - ability to identify peaks close to the edges of the spectrum region. Usually peak finders fail to detect them
1032 /// - resolution, decomposition of doublets and multiplets. The algorithm should be able to recognise close positioned peaks.
1033 /// - ability to identify peaks with different sigma
1034 ///
1035 /// #### References:
1036 ///
1037 /// [1] M.A. Mariscotti: A method for identification of peaks in the presence of
1038 /// background and its application to spectrum analysis. NIM 50 (1967), 309-320.
1039 ///
1040 /// [2] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:Identification
1041 /// of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
1042 /// 108-125.
1043 ///
1044 /// [3] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
1045 /// (1996), 451.
1046 ///
1047 /// ### Examples of peak searching method
1048 ///
1049 /// SearchHighRes function provides users with the possibility
1050 /// to vary the input parameters and with the access to the output deconvolved data
1051 /// in the destination spectrum. Based on the output data one can tune the
1052 /// parameters.
1053 ///
1054 /// ### Example 8 - Src.C
1055 ///
1056 /// Begin_Macro(source)
1057 /// ../../../tutorials/spectrum/Src.C
1058 /// End_Macro
1059 ///
1060 /// ### Example 9 - Src2.C
1061 ///
1062 /// Begin_Macro(source)
1063 /// ../../../tutorials/spectrum/Src2.C
1064 /// End_Macro
1065 ///
1066 /// ### Example 10 - Src3.C
1067 ///
1068 /// Begin_Macro(source)
1069 /// ../../../tutorials/spectrum/Src3.C
1070 /// End_Macro
1071 ///
1072 /// ### Example 11 - Src4.C
1073 ///
1074 /// Begin_Macro(source)
1075 /// ../../../tutorials/spectrum/Src4.C
1076 /// End_Macro
1077 ///
1078 /// ### Example 12 - Src5.C
1079 ///
1080 /// Begin_Macro(source)
1081 /// ../../../tutorials/spectrum/Src5.C
1082 /// End_Macro
1083 
1084 Int_t TSpectrum2::SearchHighRes(Double_t **source, Double_t **dest, Int_t ssizex, Int_t ssizey,
1085  Double_t sigma, Double_t threshold,
1086  Bool_t backgroundRemove,Int_t deconIterations,
1087  Bool_t markov, Int_t averWindow)
1088 
1089 {
1090  Int_t number_of_iterations = (Int_t)(4 * sigma + 0.5);
1091  Int_t k, lindex, priz;
1092  Double_t lda, ldb, ldc, area, maximum;
1093  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;
1094  Int_t ymin, ymax, i, j;
1095  Double_t a, b, ax, ay, maxch, plocha = 0;
1096  Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy;
1097  Double_t p1, p2, p3, p4, s1, s2, s3, s4;
1098  Int_t x, y;
1099  Int_t lhx, lhy, i1, i2, j1, j2, k1, k2, i1min, i1max, i2min, i2max, j1min, j1max, j2min, j2max, positx, posity;
1100  if (sigma < 1) {
1101  Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
1102  return 0;
1103  }
1104 
1105  if(threshold<=0||threshold>=100){
1106  Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
1107  return 0;
1108  }
1109 
1110  j = (Int_t) (5.0 * sigma + 0.5);
1111  if (j >= PEAK_WINDOW / 2) {
1112  Error("SearchHighRes", "Too large sigma");
1113  return 0;
1114  }
1115 
1116  if (markov == true) {
1117  if (averWindow <= 0) {
1118  Error("SearchHighRes", "Averaging window must be positive");
1119  return 0;
1120  }
1121  }
1122  if(backgroundRemove == true){
1123  if(ssizex_ext < 2 * number_of_iterations + 1 || ssizey_ext < 2 * number_of_iterations + 1){
1124  Error("SearchHighRes", "Too large clipping window");
1125  return 0;
1126  }
1127  }
1128  i = (Int_t)(4 * sigma + 0.5);
1129  i = 4 * i;
1130  Double_t **working_space = new Double_t *[ssizex + i];
1131  for (j = 0; j < ssizex + i; j++) {
1132  Double_t *wsk = working_space[j] = new Double_t[16 * (ssizey + i)];
1133  for (k=0;k<16 * (ssizey + i);k++) wsk[k] = 0;
1134  }
1135  for(j = 0; j < ssizey_ext; j++){
1136  for(i = 0; i < ssizex_ext; i++){
1137  if(i < shift){
1138  if(j < shift)
1139  working_space[i][j + ssizey_ext] = source[0][0];
1140 
1141  else if(j >= ssizey + shift)
1142  working_space[i][j + ssizey_ext] = source[0][ssizey - 1];
1143 
1144  else
1145  working_space[i][j + ssizey_ext] = source[0][j - shift];
1146  }
1147 
1148  else if(i >= ssizex + shift){
1149  if(j < shift)
1150  working_space[i][j + ssizey_ext] = source[ssizex - 1][0];
1151 
1152  else if(j >= ssizey + shift)
1153  working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1];
1154 
1155  else
1156  working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift];
1157  }
1158 
1159  else{
1160  if(j < shift)
1161  working_space[i][j + ssizey_ext] = source[i - shift][0];
1162 
1163  else if(j >= ssizey + shift)
1164  working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1];
1165 
1166  else
1167  working_space[i][j + ssizey_ext] = source[i - shift][j - shift];
1168  }
1169  }
1170  }
1171  if(backgroundRemove == true){
1172  for(i = 1; i <= number_of_iterations; i++){
1173  for(y = i; y < ssizey_ext - i; y++){
1174  for(x = i; x < ssizex_ext - i; x++){
1175  a = working_space[x][y + ssizey_ext];
1176  p1 = working_space[x - i][y + ssizey_ext - i];
1177  p2 = working_space[x - i][y + ssizey_ext + i];
1178  p3 = working_space[x + i][y + ssizey_ext - i];
1179  p4 = working_space[x + i][y + ssizey_ext + i];
1180  s1 = working_space[x][y + ssizey_ext - i];
1181  s2 = working_space[x - i][y + ssizey_ext];
1182  s3 = working_space[x + i][y + ssizey_ext];
1183  s4 = working_space[x][y + ssizey_ext + i];
1184  b = (p1 + p2) / 2.0;
1185  if(b > s2)
1186  s2 = b;
1187  b = (p1 + p3) / 2.0;
1188  if(b > s1)
1189  s1 = b;
1190  b = (p2 + p4) / 2.0;
1191  if(b > s4)
1192  s4 = b;
1193  b = (p3 + p4) / 2.0;
1194  if(b > s3)
1195  s3 = b;
1196  s1 = s1 - (p1 + p3) / 2.0;
1197  s2 = s2 - (p1 + p2) / 2.0;
1198  s3 = s3 - (p3 + p4) / 2.0;
1199  s4 = s4 - (p2 + p4) / 2.0;
1200  b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
1201  if(b < a)
1202  a = b;
1203  working_space[x][y] = a;
1204  }
1205  }
1206  for(y = i;y < ssizey_ext - i; y++){
1207  for(x = i; x < ssizex_ext - i; x++){
1208  working_space[x][y + ssizey_ext] = working_space[x][y];
1209  }
1210  }
1211  }
1212  for(j = 0;j < ssizey_ext; j++){
1213  for(i = 0; i < ssizex_ext; i++){
1214  if(i < shift){
1215  if(j < shift)
1216  working_space[i][j + ssizey_ext] = source[0][0] - working_space[i][j + ssizey_ext];
1217 
1218  else if(j >= ssizey + shift)
1219  working_space[i][j + ssizey_ext] = source[0][ssizey - 1] - working_space[i][j + ssizey_ext];
1220 
1221  else
1222  working_space[i][j + ssizey_ext] = source[0][j - shift] - working_space[i][j + ssizey_ext];
1223  }
1224 
1225  else if(i >= ssizex + shift){
1226  if(j < shift)
1227  working_space[i][j + ssizey_ext] = source[ssizex - 1][0] - working_space[i][j + ssizey_ext];
1228 
1229  else if(j >= ssizey + shift)
1230  working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1] - working_space[i][j + ssizey_ext];
1231 
1232  else
1233  working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift] - working_space[i][j + ssizey_ext];
1234  }
1235 
1236  else{
1237  if(j < shift)
1238  working_space[i][j + ssizey_ext] = source[i - shift][0] - working_space[i][j + ssizey_ext];
1239 
1240  else if(j >= ssizey + shift)
1241  working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1] - working_space[i][j + ssizey_ext];
1242 
1243  else
1244  working_space[i][j + ssizey_ext] = source[i - shift][j - shift] - working_space[i][j + ssizey_ext];
1245  }
1246  }
1247  }
1248  }
1249  for(j = 0; j < ssizey_ext; j++){
1250  for(i = 0; i < ssizex_ext; i++){
1251  working_space[i][j + 15*ssizey_ext] = working_space[i][j + ssizey_ext];
1252  }
1253  }
1254  if(markov == true){
1255  for(i = 0;i < ssizex_ext; i++){
1256  for(j = 0; j < ssizey_ext; j++)
1257  working_space[i][j + 2 * ssizex_ext] = working_space[i][ssizey_ext + j];
1258  }
1259  xmin = 0;
1260  xmax = ssizex_ext - 1;
1261  ymin = 0;
1262  ymax = ssizey_ext - 1;
1263  for(i = 0, maxch = 0; i < ssizex_ext; i++){
1264  for(j = 0; j < ssizey_ext; j++){
1265  working_space[i][j] = 0;
1266  if(maxch < working_space[i][j + 2 * ssizey_ext])
1267  maxch = working_space[i][j + 2 * ssizey_ext];
1268  plocha += working_space[i][j + 2 * ssizey_ext];
1269  }
1270  }
1271  if(maxch == 0) {
1272  delete [] working_space;
1273  return 0;
1274  }
1275 
1276  nom=0;
1277  working_space[xmin][ymin] = 1;
1278  for(i = xmin; i < xmax; i++){
1279  nip = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1280  nim = working_space[i + 1][ymin + 2 * ssizey_ext] / maxch;
1281  sp = 0,sm = 0;
1282  for(l = 1;l <= averWindow; l++){
1283  if((i + l) > xmax)
1284  a = working_space[xmax][ymin + 2 * ssizey_ext] / maxch;
1285 
1286  else
1287  a = working_space[i + l][ymin + 2 * ssizey_ext] / maxch;
1288 
1289  b = a - nip;
1290  if(a + nip <= 0)
1291  a = 1;
1292 
1293  else
1294  a=TMath::Sqrt(a + nip);
1295  b = b / a;
1296  b = TMath::Exp(b);
1297  sp = sp + b;
1298  if(i - l + 1 < xmin)
1299  a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1300 
1301  else
1302  a = working_space[i - l + 1][ymin + 2 * ssizey_ext] / maxch;
1303  b = a - nim;
1304  if(a + nim <= 0)
1305  a = 1;
1306 
1307  else
1308  a=TMath::Sqrt(a + nim);
1309  b = b / a;
1310  b = TMath::Exp(b);
1311  sm = sm + b;
1312  }
1313  a = sp / sm;
1314  a = working_space[i + 1][ymin] = a * working_space[i][ymin];
1315  nom = nom + a;
1316  }
1317  for(i = ymin; i < ymax; i++){
1318  nip = working_space[xmin][i + 2 * ssizey_ext] / maxch;
1319  nim = working_space[xmin][i + 1 + 2 * ssizey_ext] / maxch;
1320  sp = 0,sm = 0;
1321  for(l = 1; l <= averWindow; l++){
1322  if((i + l) > ymax)
1323  a = working_space[xmin][ymax + 2 * ssizey_ext] / maxch;
1324 
1325  else
1326  a = working_space[xmin][i + l + 2 * ssizey_ext] / maxch;
1327  b = a - nip;
1328  if(a + nip <= 0)
1329  a=1;
1330 
1331  else
1332  a=TMath::Sqrt(a + nip);
1333  b = b / a;
1334  b = TMath::Exp(b);
1335  sp = sp + b;
1336  if(i - l + 1 < ymin)
1337  a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1338 
1339  else
1340  a = working_space[xmin][i - l + 1 + 2 * ssizey_ext] / maxch;
1341  b = a - nim;
1342  if(a + nim <= 0)
1343  a = 1;
1344 
1345  else
1346  a=TMath::Sqrt(a + nim);
1347  b = b / a;
1348  b = TMath::Exp(b);
1349  sm = sm + b;
1350  }
1351  a = sp / sm;
1352  a = working_space[xmin][i + 1] = a * working_space[xmin][i];
1353  nom = nom + a;
1354  }
1355  for(i = xmin; i < xmax; i++){
1356  for(j = ymin; j < ymax; j++){
1357  nip = working_space[i][j + 1 + 2 * ssizey_ext] / maxch;
1358  nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1359  spx = 0,smx = 0;
1360  for(l = 1; l <= averWindow; l++){
1361  if(i + l > xmax)
1362  a = working_space[xmax][j + 2 * ssizey_ext] / maxch;
1363 
1364  else
1365  a = working_space[i + l][j + 2 * ssizey_ext] / maxch;
1366  b = a - nip;
1367  if(a + nip <= 0)
1368  a = 1;
1369 
1370  else
1371  a=TMath::Sqrt(a + nip);
1372  b = b / a;
1373  b = TMath::Exp(b);
1374  spx = spx + b;
1375  if(i - l + 1 < xmin)
1376  a = working_space[xmin][j + 2 * ssizey_ext] / maxch;
1377 
1378  else
1379  a = working_space[i - l + 1][j + 2 * ssizey_ext] / maxch;
1380  b = a - nim;
1381  if(a + nim <= 0)
1382  a=1;
1383 
1384  else
1385  a=TMath::Sqrt(a + nim);
1386  b = b / a;
1387  b = TMath::Exp(b);
1388  smx = smx + b;
1389  }
1390  spy = 0,smy = 0;
1391  nip = working_space[i + 1][j + 2 * ssizey_ext] / maxch;
1392  nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1393  for(l = 1; l <= averWindow; l++){
1394  if(j + l > ymax)
1395  a = working_space[i][ymax + 2 * ssizey_ext] / maxch;
1396 
1397  else
1398  a = working_space[i][j + l + 2 * ssizey_ext] / maxch;
1399  b = a - nip;
1400  if(a + nip <= 0)
1401  a = 1;
1402 
1403  else
1404  a=TMath::Sqrt(a + nip);
1405  b = b / a;
1406  b = TMath::Exp(b);
1407  spy = spy + b;
1408  if(j - l + 1 < ymin)
1409  a = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1410 
1411  else
1412  a = working_space[i][j - l + 1 + 2 * ssizey_ext] / maxch;
1413  b=a-nim;
1414  if(a + nim <= 0)
1415  a = 1;
1416  else
1417  a=TMath::Sqrt(a + nim);
1418  b = b / a;
1419  b = TMath::Exp(b);
1420  smy = smy + b;
1421  }
1422  a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx + smy);
1423  working_space[i + 1][j + 1] = a;
1424  nom = nom + a;
1425  }
1426  }
1427  for(i = xmin; i <= xmax; i++){
1428  for(j = ymin; j <= ymax; j++){
1429  working_space[i][j] = working_space[i][j] / nom;
1430  }
1431  }
1432  for(i = 0; i < ssizex_ext; i++){
1433  for(j = 0; j < ssizey_ext; j++){
1434  working_space[i][j + ssizey_ext] = working_space[i][j] * plocha;
1435  working_space[i][2 * ssizey_ext + j] = working_space[i][ssizey_ext + j];
1436  }
1437  }
1438  }
1439  //deconvolution starts
1440  area = 0;
1441  lhx = -1,lhy = -1;
1442  positx = 0,posity = 0;
1443  maximum = 0;
1444  //generate response matrix
1445  for(i = 0; i < ssizex_ext; i++){
1446  for(j = 0; j < ssizey_ext; j++){
1447  lda = (Double_t)i - 3 * sigma;
1448  ldb = (Double_t)j - 3 * sigma;
1449  lda = (lda * lda + ldb * ldb) / (2 * sigma * sigma);
1450  k=(Int_t)(1000 * TMath::Exp(-lda));
1451  lda = k;
1452  if(lda != 0){
1453  if((i + 1) > lhx)
1454  lhx = i + 1;
1455 
1456  if((j + 1) > lhy)
1457  lhy = j + 1;
1458  }
1459  working_space[i][j] = lda;
1460  area = area + lda;
1461  if(lda > maximum){
1462  maximum = lda;
1463  positx = i,posity = j;
1464  }
1465  }
1466  }
1467  //read source matrix
1468  for(i = 0;i < ssizex_ext; i++){
1469  for(j = 0;j < ssizey_ext; j++){
1470  working_space[i][j + 14 * ssizey_ext] = TMath::Abs(working_space[i][j + ssizey_ext]);
1471  }
1472  }
1473  //calculate matrix b=ht*h
1474  i = lhx - 1;
1475  if(i > ssizex_ext)
1476  i = ssizex_ext;
1477 
1478  j = lhy - 1;
1479  if(j>ssizey_ext)
1480  j = ssizey_ext;
1481 
1482  i1min = -i,i1max = i;
1483  i2min = -j,i2max = j;
1484  for(i2 = i2min; i2 <= i2max; i2++){
1485  for(i1 = i1min; i1 <= i1max; i1++){
1486  ldc = 0;
1487  j2min = -i2;
1488  if(j2min < 0)
1489  j2min = 0;
1490 
1491  j2max = lhy - 1 - i2;
1492  if(j2max > lhy - 1)
1493  j2max = lhy - 1;
1494 
1495  for(j2 = j2min; j2 <= j2max; j2++){
1496  j1min = -i1;
1497  if(j1min < 0)
1498  j1min = 0;
1499 
1500  j1max = lhx - 1 - i1;
1501  if(j1max > lhx - 1)
1502  j1max = lhx - 1;
1503 
1504  for(j1 = j1min; j1 <= j1max; j1++){
1505  lda = working_space[j1][j2];
1506  ldb = working_space[i1 + j1][i2 + j2];
1507  ldc = ldc + lda * ldb;
1508  }
1509  }
1510  k = (i1 + ssizex_ext) / ssizex_ext;
1511  working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext] = ldc;
1512  }
1513  }
1514  //calculate at*y and write into p
1515  i = lhx - 1;
1516  if(i > ssizex_ext)
1517  i = ssizex_ext;
1518 
1519  j = lhy - 1;
1520  if(j > ssizey_ext)
1521  j = ssizey_ext;
1522 
1523  i2min = -j,i2max = ssizey_ext + j - 1;
1524  i1min = -i,i1max = ssizex_ext + i - 1;
1525  for(i2 = i2min; i2 <= i2max; i2++){
1526  for(i1=i1min;i1<=i1max;i1++){
1527  ldc=0;
1528  for(j2 = 0; j2 <= (lhy - 1); j2++){
1529  for(j1 = 0; j1 <= (lhx - 1); j1++){
1530  k2 = i2 + j2,k1 = i1 + j1;
1531  if(k2 >= 0 && k2 < ssizey_ext && k1 >= 0 && k1 < ssizex_ext){
1532  lda = working_space[j1][j2];
1533  ldb = working_space[k1][k2 + 14 * ssizey_ext];
1534  ldc = ldc + lda * ldb;
1535  }
1536  }
1537  }
1538  k = (i1 + ssizex_ext) / ssizex_ext;
1539  working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext] = ldc;
1540  }
1541  }
1542  //move matrix p
1543  for(i2 = 0; i2 < ssizey_ext; i2++){
1544  for(i1 = 0; i1 < ssizex_ext; i1++){
1545  k = (i1 + ssizex_ext) / ssizex_ext;
1546  ldb = working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext];
1547  working_space[i1][i2 + 14 * ssizey_ext] = ldb;
1548  }
1549  }
1550  //initialization in x1 matrix
1551  for(i2 = 0; i2 < ssizey_ext; i2++){
1552  for(i1 = 0; i1 < ssizex_ext; i1++){
1553  working_space[i1][i2 + ssizey_ext] = 1;
1554  working_space[i1][i2 + 2 * ssizey_ext] = 0;
1555  }
1556  }
1557  //START OF ITERATIONS
1558  for(lindex = 0; lindex < deconIterations; lindex++){
1559  for(i2 = 0; i2 < ssizey_ext; i2++){
1560  for(i1 = 0; i1 < ssizex_ext; i1++){
1561  lda = working_space[i1][i2 + ssizey_ext];
1562  ldc = working_space[i1][i2 + 14 * ssizey_ext];
1563  if(lda > 0.000001 && ldc > 0.000001){
1564  ldb=0;
1565  j2min=i2;
1566  if(j2min > lhy - 1)
1567  j2min = lhy - 1;
1568 
1569  j2min = -j2min;
1570  j2max = ssizey_ext - i2 - 1;
1571  if(j2max > lhy - 1)
1572  j2max = lhy - 1;
1573 
1574  j1min = i1;
1575  if(j1min > lhx - 1)
1576  j1min = lhx - 1;
1577 
1578  j1min = -j1min;
1579  j1max = ssizex_ext - i1 - 1;
1580  if(j1max > lhx - 1)
1581  j1max = lhx - 1;
1582 
1583  for(j2 = j2min; j2 <= j2max; j2++){
1584  for(j1 = j1min; j1 <= j1max; j1++){
1585  k = (j1 + ssizex_ext) / ssizex_ext;
1586  ldc = working_space[(j1 + ssizex_ext) % ssizex_ext][j2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext];
1587  lda = working_space[i1 + j1][i2 + j2 + ssizey_ext];
1588  ldb = ldb + lda * ldc;
1589  }
1590  }
1591  lda = working_space[i1][i2 + ssizey_ext];
1592  ldc = working_space[i1][i2 + 14 * ssizey_ext];
1593  if(ldc * lda != 0 && ldb != 0){
1594  lda =lda * ldc / ldb;
1595  }
1596 
1597  else
1598  lda=0;
1599  working_space[i1][i2 + 2 * ssizey_ext] = lda;
1600  }
1601  }
1602  }
1603  for(i2 = 0; i2 < ssizey_ext; i2++){
1604  for(i1 = 0; i1 < ssizex_ext; i1++)
1605  working_space[i1][i2 + ssizey_ext] = working_space[i1][i2 + 2 * ssizey_ext];
1606  }
1607  }
1608  //looking for maximum
1609  maximum=0;
1610  for(i = 0; i < ssizex_ext; i++){
1611  for(j = 0; j < ssizey_ext; j++){
1612  working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext] = area * working_space[i][j + ssizey_ext];
1613  if(maximum < working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext])
1614  maximum = working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext];
1615 
1616  }
1617  }
1618  //searching for peaks in deconvolved spectrum
1619  for(i = 1; i < ssizex_ext - 1; i++){
1620  for(j = 1; j < ssizey_ext - 1; j++){
1621  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]){
1622  if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift){
1623  if(working_space[i][j] > threshold * maximum / 100.0){
1624  if(peak_index < fMaxPeaks){
1625  for(k = i - 1,a = 0,b = 0; k <= i + 1; k++){
1626  a += (Double_t)(k - shift) * working_space[k][j];
1627  b += working_space[k][j];
1628  }
1629  ax=a/b;
1630  if(ax < 0)
1631  ax = 0;
1632 
1633  if(ax >= ssizex)
1634  ax = ssizex - 1;
1635 
1636  for(k = j - 1,a = 0,b = 0; k <= j + 1; k++){
1637  a += (Double_t)(k - shift) * working_space[i][k];
1638  b += working_space[i][k];
1639  }
1640  ay=a/b;
1641  if(ay < 0)
1642  ay = 0;
1643 
1644  if(ay >= ssizey)
1645  ay = ssizey - 1;
1646 
1647  if(peak_index == 0){
1648  fPositionX[0] = ax;
1649  fPositionY[0] = ay;
1650  peak_index = 1;
1651  }
1652 
1653  else{
1654  for(k = 0, priz = 0; k < peak_index && priz == 0; k++){
1655  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)])
1656  priz = 1;
1657  }
1658  if(priz == 0){
1659  if(k < fMaxPeaks){
1660  fPositionX[k] = ax;
1661  fPositionY[k] = ay;
1662  }
1663  }
1664 
1665  else{
1666  for(l = peak_index; l >= k; l--){
1667  if(l < fMaxPeaks){
1668  fPositionX[l] = fPositionX[l - 1];
1669  fPositionY[l] = fPositionY[l - 1];
1670  }
1671  }
1672  fPositionX[k - 1] = ax;
1673  fPositionY[k - 1] = ay;
1674  }
1675  if(peak_index < fMaxPeaks)
1676  peak_index += 1;
1677  }
1678  }
1679  }
1680  }
1681  }
1682  }
1683  }
1684  //writing back deconvolved spectrum
1685  for(i = 0; i < ssizex; i++){
1686  for(j = 0; j < ssizey; j++){
1687  dest[i][j] = working_space[i + shift][j + shift];
1688  }
1689  }
1690  i = (Int_t)(4 * sigma + 0.5);
1691  i = 4 * i;
1692  for (j = 0; j < ssizex + i; j++)
1693  delete[]working_space[j];
1694  delete[]working_space;
1695  fNPeaks = peak_index;
1696  return fNPeaks;
1697 }
1698 
1699 ////////////////////////////////////////////////////////////////////////////////
1700 /// static function (called by TH1), interface to TSpectrum2::Search
1701 
1702 Int_t TSpectrum2::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
1703 {
1704  TSpectrum2 s;
1705  return s.Search(hist,sigma,option,threshold);
1706 }
1707 
1708 ////////////////////////////////////////////////////////////////////////////////
1709 /// static function (called by TH1), interface to TSpectrum2::Background
1710 
1711 TH1 *TSpectrum2::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
1712 {
1713  TSpectrum2 s;
1714  return s.Background(hist,niter,option);
1715 }