Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TSpectrum3.cxx
Go to the documentation of this file.
1 // @(#)root/spectrum:$Id$
2 // Author: Miroslav Morhac 25/09/2006
3 
4 /** \class TSpectrum3
5  \ingroup Spectrum
6  \brief Advanced 3-dimensional spectra processing functions
7  \author Miroslav Morhac
8 
9  This class contains advanced spectra processing functions.
10 
11  - Three-dimensional background estimation functions
12  - Three-dimensional smoothing functions
13  - Three-dimensional deconvolution functions
14  - Three-dimensional peak search functions
15 
16 
17  The algorithms in this class have been published in the following
18  references:
19 
20  [1] M.Morhac et al.: Background elimination methods for
21  multidimensional coincidence gamma-ray spectra. Nuclear
22  Instruments and Methods in Physics Research A 401 (1997) 113-132.
23 
24  [2] M.Morhac et al.: Efficient one- and two-dimensional Gold
25  deconvolution and its application to gamma-ray spectra
26  decomposition. Nuclear Instruments and Methods in Physics
27  Research A 401 (1997) 385-408.
28 
29  [3] M. Morhac et al.: Efficient algorithm of multidimensional
30  deconvolution and its application to nuclear data processing. Digital
31  Signal Processing, Vol. 13, No. 1, (2003), 144-171.
32 
33  [4] M.Morhac et al.: Identification of peaks in multidimensional
34  coincidence gamma-ray spectra. Nuclear Instruments and Methods in
35  Research Physics A 443(2000), 108-125.
36 
37  These NIM papers are also available as Postscript files from:
38 
39  - [SpectrumDec.ps.gz](ftp://root.cern.ch/root/SpectrumDec.ps.gz)
40  - [SpectrumSrc.ps.gz](ftp://root.cern.ch/root/SpectrumSrc.ps.gz)
41  - [SpectrumBck.ps.gz](ftp://root.cern.ch/root/SpectrumBck.ps.gz)
42 
43  See also the
44  [online documentation](https://root.cern.ch/guides/tspectrum-manual) and
45  [tutorials](https://root.cern.ch/doc/master/group__tutorial__spectrum.html).
46 */
47 
48 #include "TSpectrum3.h"
49 #include "TH1.h"
50 #include "TMath.h"
51 #define PEAK_WINDOW 1024
52 
53 ClassImp(TSpectrum3);
54 
55 ////////////////////////////////////////////////////////////////////////////////
56 /// Constructor.
57 
58 TSpectrum3::TSpectrum3() :TNamed("Spectrum", "Miroslav Morhac peak finder")
59 {
60  Int_t n = 100;
61  fMaxPeaks = n;
62  fPosition = new Double_t[n];
63  fPositionX = new Double_t[n];
64  fPositionY = new Double_t[n];
65  fPositionZ = new Double_t[n];
66  fResolution = 1;
67  fHistogram = 0;
68  fNPeaks = 0;
69 }
70 
71 
72 ////////////////////////////////////////////////////////////////////////////////
73 /// - maxpositions: maximum number of peaks
74 /// - resolution: *NOT USED* determines resolution of the neighbouring peaks
75 /// default value is 1 correspond to 3 sigma distance
76 /// between peaks. Higher values allow higher resolution
77 /// (smaller distance between peaks.
78 /// May be set later through SetResolution.
79 
80 TSpectrum3::TSpectrum3(Int_t maxpositions, Double_t resolution) :TNamed("Spectrum", "Miroslav Morhac peak finder")
81 {
82  Int_t n = TMath::Max(maxpositions, 100);
83  fMaxPeaks = n;
84  fPosition = new Double_t[n];
85  fPositionX = new Double_t[n];
86  fPositionY = new Double_t[n];
87  fPositionZ = new Double_t[n];
88  fHistogram = 0;
89  fNPeaks = 0;
90  SetResolution(resolution);
91 }
92 
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 /// Destructor.
96 
97 TSpectrum3::~TSpectrum3()
98 {
99  delete [] fPosition;
100  delete [] fPositionX;
101  delete [] fPositionY;
102  delete [] fPositionZ;
103  delete fHistogram;
104 }
105 
106 ////////////////////////////////////////////////////////////////////////////////
107 /// This function calculates background spectrum from source in h.
108 /// The result is placed in the vector pointed by spectrum pointer.
109 ///
110 /// Function parameters:
111 /// - spectrum: pointer to the vector of source spectrum
112 /// - size: length of spectrum and working space vectors
113 /// - number_of_iterations, for details we refer to manual
114 
115 const char *TSpectrum3::Background(const TH1 * h, Int_t number_of_iterations,
116  Option_t * option)
117 {
118  Error("Background","function not yet implemented: h=%s, iter=%d, option=%sn"
119  , h->GetName(), number_of_iterations, option);
120  return 0;
121 }
122 
123 ////////////////////////////////////////////////////////////////////////////////
124 /// Print the array of positions
125 
126 void TSpectrum3::Print(Option_t *) const
127 {
128  printf("\nNumber of positions = %d\n",fNPeaks);
129  for (Int_t i=0;i<fNPeaks;i++) {
130  printf(" x[%d] = %g, y[%d] = %g, z[%d] = %g\n",i,fPositionX[i],i,fPositionY[i],i,fPositionZ[i]);
131  }
132 }
133 
134 
135 
136 ////////////////////////////////////////////////////////////////////////////////
137 /// This function searches for peaks in source spectrum in hin
138 /// The number of found peaks and their positions are written into
139 /// the members fNpeaks and fPositionX.
140 ///
141 /// Function parameters:
142 /// - hin: pointer to the histogram of source spectrum
143 /// - sigma: sigma of searched peaks, for details we refer to manual
144 /// Note that sigma is in number of bins
145 /// - threshold: (default=0.05) peaks with amplitude less than
146 /// threshold*highest_peak are discarded.
147 ///
148 /// if option is not equal to "goff" (goff is the default), then
149 /// a polymarker object is created and added to the list of functions of
150 /// the histogram. The histogram is drawn with the specified option and
151 /// the polymarker object drawn on top of the histogram.
152 /// The polymarker coordinates correspond to the npeaks peaks found in
153 /// the histogram.
154 /// A pointer to the polymarker object can be retrieved later via:
155 /// ~~~ {.cpp}
156 /// TList *functions = hin->GetListOfFunctions();
157 /// TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker")
158 /// ~~~
159 
160 Int_t TSpectrum3::Search(const TH1 * hin, Double_t sigma,
161  Option_t * option, Double_t threshold)
162 {
163  if (hin == 0)
164  return 0;
165  Int_t dimension = hin->GetDimension();
166  if (dimension != 3) {
167  Error("Search", "Must be a 3-d histogram");
168  return 0;
169  }
170 
171  Int_t sizex = hin->GetXaxis()->GetNbins();
172  Int_t sizey = hin->GetYaxis()->GetNbins();
173  Int_t sizez = hin->GetZaxis()->GetNbins();
174  Int_t i, j, k, binx,biny,binz, npeaks;
175  Double_t*** source = new Double_t**[sizex];
176  Double_t*** dest = new Double_t**[sizex];
177  for (i = 0; i < sizex; i++) {
178  source[i] = new Double_t*[sizey];
179  dest[i] = new Double_t*[sizey];
180  for (j = 0; j < sizey; j++) {
181  source[i][j] = new Double_t[sizez];
182  dest[i][j] = new Double_t[sizez];
183  for (k = 0; k < sizez; k++)
184  source[i][j][k] = (Double_t) hin->GetBinContent(i + 1, j + 1, k + 1);
185  }
186  }
187  //the smoothing option is used for 1-d but not for 2-d histograms
188  npeaks = SearchHighRes((const Double_t***)source, dest, sizex, sizey, sizez, sigma, 100*threshold, kTRUE, 3, kFALSE, 3);
189 
190  //The logic in the loop should be improved to use the fact
191  //that fPositionX,Y give a precise position inside a bin.
192  //The current algorithm takes the center of the bin only.
193  for (i = 0; i < npeaks; i++) {
194  binx = 1 + Int_t(fPositionX[i] + 0.5);
195  biny = 1 + Int_t(fPositionY[i] + 0.5);
196  binz = 1 + Int_t(fPositionZ[i] + 0.5);
197  fPositionX[i] = hin->GetXaxis()->GetBinCenter(binx);
198  fPositionY[i] = hin->GetYaxis()->GetBinCenter(biny);
199  fPositionZ[i] = hin->GetZaxis()->GetBinCenter(binz);
200  }
201  for (i = 0; i < sizex; i++) {
202  for (j = 0; j < sizey; j++){
203  delete [] source[i][j];
204  delete [] dest[i][j];
205  }
206  delete [] source[i];
207  delete [] dest[i];
208  }
209  delete [] source;
210  delete [] dest;
211 
212  if (strstr(option, "goff"))
213  return npeaks;
214  if (!npeaks) return 0;
215  return npeaks;
216 }
217 
218 
219 ////////////////////////////////////////////////////////////////////////////////
220 /// *NOT USED*
221 /// resolution: determines resolution of the neighbouring peaks
222 /// default value is 1 correspond to 3 sigma distance
223 /// between peaks. Higher values allow higher resolution
224 /// (smaller distance between peaks.
225 /// May be set later through SetResolution.
226 
227 void TSpectrum3::SetResolution(Double_t resolution)
228 {
229  if (resolution > 1)
230  fResolution = resolution;
231  else
232  fResolution = 1;
233 }
234 
235 ////////////////////////////////////////////////////////////////////////////////
236 /// This function calculates background spectrum from source spectrum.
237 /// The result is placed to the array pointed by spectrum pointer.
238 ///
239 /// Function parameters:
240 /// - spectrum-pointer to the array of source spectrum
241 /// - ssizex-x length of spectrum
242 /// - ssizey-y length of spectrum
243 /// - ssizez-z length of spectrum
244 /// - numberIterationsX-maximal x width of clipping window
245 /// - numberIterationsY-maximal y width of clipping window
246 /// - numberIterationsZ-maximal z width of clipping window
247 /// for details we refer to manual
248 /// - direction- direction of change of clipping window
249 /// - possible values=kBackIncreasingWindow, kBackDecreasingWindow
250 /// - filterType-determines the algorithm of the filtering
251 /// -possible values=kBackSuccessiveFiltering, kBackOneStepFiltering
252 ///
253 /// ### Background estimation
254 ///
255 /// Goal: Separation of useful information (peaks) from useless information (background)
256 ///
257 /// - method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping
258 /// algorithm [1]
259 /// - there exist two algorithms for the estimation of new value in the
260 /// channel \f$i_1, i_2, i_3\f$
261 ///
262 /// #### Algorithm based on Successive Comparisons
263 ///
264 /// It is an extension of one-dimensional SNIP algorithm to another dimension.
265 /// For details we refer to [2].
266 ///
267 /// #### Algorithm based on One Step Filtering
268 ///
269 /// The algorithm is analogous to that for 2-dimensional data. For details we
270 /// refer to TSpectrum2. New value in the estimated channel is calculated as
271 /// \f$ a = \nu_{p-1}(i_1, i_2, i_3)\f$
272 ///
273 /// \image html spectrum3_background_image003.gif
274 /// \f[
275 /// \nu_p(i_1, i_2, i_3) = min (a,b)
276 /// \f]
277 ///
278 /// where p = 1, 2, ..., number_of_iterations.
279 ///
280 /// #### References:
281 ///
282 /// [1] C. G Ryan et al.: SNIP, a
283 /// statistics-sensitive background treatment for the quantitative analysis of PIXE
284 /// spectra in geoscience applications. NIM, B34 (1988), 396-402./
285 ///
286 /// [2] M.Morhac, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.: Background
287 /// elimination methods for multidimensional gamma-ray spectra. NIM, A401 (1997)
288 /// 113-132.
289 ///
290 /// Example 1- script Back3.c :
291 ///
292 /// \image html spectrum3_background_image005.jpg Fig. 1 Original three-dimensional gamma-gamma-gamma-ray spectrum
293 /// \image html spectrum3_background_image006.jpg Fig. 2 Background estimated from data from Fig. 1 using decreasing clipping window with widths 5, 5, 5 and algorithm based on successive comparisons. The estimate includes not only continuously changing background but also one- and two-dimensional ridges.
294 /// \image html spectrum3_background_image007.jpg Fig. 3 Resulting peaks after subtraction of the estimated background (Fig. 2) from original three-dimensional gamma-gamma-gamma-ray spectrum (Fig. 1).
295 ///
296 /// #### Script:
297 ///
298 /// Example to illustrate the background estimator (class TSpectrum3).
299 /// To execute this example, do:
300 ///
301 /// `root > .x Back3.C`
302 ///
303 /// ~~~ {.cpp}
304 /// void Back3() {
305 /// Int_t i, j, k;
306 /// Int_t nbinsx = 64;
307 /// Int_t nbinsy = 64;
308 /// Int_t nbinsz = 64;
309 /// Int_t xmin = 0;
310 /// Int_t xmax = nbinsx;
311 /// Int_t ymin = 0;
312 /// Int_t ymax = nbinsy;
313 /// Int_t zmin = 0;
314 /// Int_t zmax = nbinsz;
315 /// Double_t*** source = new Double_t**[nbinsx];
316 /// Double_t*** dest = new Double_t**[nbinsx];
317 /// for(i=0;i<nbinsx;i++){
318 /// source[i]=new Double_t*[nbinsy];
319 /// for(j=0;j<nbinsy;j++)
320 /// source[i][j]=new Double_t[nbinsz];
321 /// }
322 /// for(i=0;i<nbinsx;i++){
323 /// dest[i]=new Double_t*[nbinsy];
324 /// for(j=0;j<nbinsy;j++)
325 /// dest[i][j]=new Double_t[nbinsz];
326 /// }
327 /// TH3F *back = new TH3F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
328 /// TFile *f = new TFile("TSpectrum3.root");
329 /// back=(TH3F*)f->Get("back;1");
330 /// TCanvas *Background = new TCanvas("Background","Estimation of background with decreasing window",10,10,1000,700);
331 /// TSpectrum3 *s = new TSpectrum3();
332 /// for (i = 0; i < nbinsx; i++){
333 /// for (j = 0; j < nbinsy; j++){
334 /// for (k = 0; k < nbinsz; k++){
335 /// source[i][j][k] = back->GetBinContent(i + 1,j + 1,k + 1);
336 /// dest[i][j][k] = back->GetBinContent(i + 1,j + 1,k + 1);
337 /// }
338 /// }
339 /// }
340 /// s->Background(dest,nbinsx,nbinsy,nbinsz,5,5,5,s->kBackDecreasingWindow,s->kBackSuccessiveFiltering);
341 /// for (i = 0; i < nbinsx; i++){
342 /// for (j = 0; j < nbinsy; j++){
343 /// for (k = 0; k < nbinsz; k++){
344 /// back->SetBinContent(i + 1,j + 1,k + 1, dest[i][j][k]);
345 /// }
346 /// }
347 /// }
348 /// FILE *out;
349 /// char PATH[80];
350 /// strcpy(PATH,"spectra3/back_output_5ds.spe");
351 /// out=fopen(PATH,"wb");
352 /// for(i=0;i<nbinsx;i++){
353 /// for(j=0;j<nbinsy;j++){
354 /// fwrite(dest[i][j], sizeof(dest[0][0][0]),nbinsz,out);
355 /// }
356 /// }
357 /// fclose(out);
358 /// for (i = 0; i < nbinsx; i++){
359 /// for (j = 0; j <nbinsy; j++){
360 /// for (k = 0; k <nbinsz; k++){
361 /// source[i][j][k] = source[i][j][k] - dest[i][j][k];
362 /// }
363 /// }
364 /// }
365 /// for (i = 0; i < nbinsx; i++){
366 /// for (j = 0; j < nbinsy; j++){
367 /// for (k = 0; k < nbinsz; k++){
368 /// back->SetBinContent(i + 1,j + 1,k + 1, source[i][j][k]);
369 /// }
370 /// }
371 /// }
372 /// strcpy(PATH,"spectra3/back_peaks_5ds.spe");
373 /// out=fopen(PATH,"wb");
374 /// for(i=0;i<nbinsx;i++){
375 /// for(j=0;j<nbinsy;j++){
376 /// fwrite(source[i][j], sizeof(source[0][0][0]),nbinsz,out);
377 /// }
378 /// }
379 /// fclose(out);
380 /// back->Draw("");
381 /// }
382 /// ~~~
383 
384 const char *TSpectrum3::Background(Double_t***spectrum,
385  Int_t ssizex, Int_t ssizey, Int_t ssizez,
386  Int_t numberIterationsX,
387  Int_t numberIterationsY,
388  Int_t numberIterationsZ,
389  Int_t direction,
390  Int_t filterType)
391 {
392  Int_t i, j, x, y, z, sampling, q1, q2, q3;
393  Double_t a, b, c, d, p1, p2, p3, p4, p5, p6, p7, p8, s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12, r1, r2, r3, r4, r5, r6;
394  if (ssizex <= 0 || ssizey <= 0 || ssizez <= 0)
395  return "Wrong parameters";
396  if (numberIterationsX < 1 || numberIterationsY < 1 || numberIterationsZ < 1)
397  return "Width of Clipping Window Must Be Positive";
398  if (ssizex < 2 * numberIterationsX + 1 || ssizey < 2 * numberIterationsY + 1 || ssizey < 2 * numberIterationsZ + 1)
399  return ("Too Large Clipping Window");
400  Double_t*** working_space=new Double_t**[ssizex];
401  for(i=0;i<ssizex;i++){
402  working_space[i] =new Double_t*[ssizey];
403  for(j=0;j<ssizey;j++)
404  working_space[i][j]=new Double_t[ssizez];
405  }
406  sampling =(Int_t) TMath::Max(numberIterationsX, numberIterationsY);
407  sampling =(Int_t) TMath::Max(sampling, numberIterationsZ);
408  if (direction == kBackIncreasingWindow) {
409  if (filterType == kBackSuccessiveFiltering) {
410  for (i = 1; i <= sampling; i++) {
411  q1 = (Int_t) TMath::Min(i, numberIterationsX), q2 =(Int_t) TMath::Min(i, numberIterationsY), q3 =(Int_t) TMath::Min(i, numberIterationsZ);
412  for (z = q3; z < ssizez - q3; z++) {
413  for (y = q2; y < ssizey - q2; y++) {
414  for (x = q1; x < ssizex - q1; x++) {
415  a = spectrum[x][y][z];
416  p1 = spectrum[x + q1][y + q2][z - q3];
417  p2 = spectrum[x - q1][y + q2][z - q3];
418  p3 = spectrum[x + q1][y - q2][z - q3];
419  p4 = spectrum[x - q1][y - q2][z - q3];
420  p5 = spectrum[x + q1][y + q2][z + q3];
421  p6 = spectrum[x - q1][y + q2][z + q3];
422  p7 = spectrum[x + q1][y - q2][z + q3];
423  p8 = spectrum[x - q1][y - q2][z + q3];
424  s1 = spectrum[x + q1][y ][z - q3];
425  s2 = spectrum[x ][y + q2][z - q3];
426  s3 = spectrum[x - q1][y ][z - q3];
427  s4 = spectrum[x ][y - q2][z - q3];
428  s5 = spectrum[x + q1][y ][z + q3];
429  s6 = spectrum[x ][y + q2][z + q3];
430  s7 = spectrum[x - q1][y ][z + q3];
431  s8 = spectrum[x ][y - q2][z + q3];
432  s9 = spectrum[x - q1][y + q2][z ];
433  s10 = spectrum[x - q1][y - q2][z ];
434  s11 = spectrum[x + q1][y + q2][z ];
435  s12 = spectrum[x + q1][y - q2][z ];
436  r1 = spectrum[x ][y ][z - q3];
437  r2 = spectrum[x ][y ][z + q3];
438  r3 = spectrum[x - q1][y ][z ];
439  r4 = spectrum[x + q1][y ][z ];
440  r5 = spectrum[x ][y + q2][z ];
441  r6 = spectrum[x ][y - q2][z ];
442  b = (p1 + p3) / 2.0;
443  if(b > s1)
444  s1 = b;
445  b = (p1 + p2) / 2.0;
446  if(b > s2)
447  s2 = b;
448  b = (p2 + p4) / 2.0;
449  if(b > s3)
450  s3 = b;
451  b = (p3 + p4) / 2.0;
452  if(b > s4)
453  s4 = b;
454  b = (p5 + p7) / 2.0;
455  if(b > s5)
456  s5 = b;
457  b = (p5 + p6) / 2.0;
458  if(b > s6)
459  s6 = b;
460  b = (p6 + p8) / 2.0;
461  if(b > s7)
462  s7 = b;
463  b = (p7 + p8) / 2.0;
464  if(b > s8)
465  s8 = b;
466  b = (p2 + p6) / 2.0;
467  if(b > s9)
468  s9 = b;
469  b = (p4 + p8) / 2.0;
470  if(b > s10)
471  s10 = b;
472  b = (p1 + p5) / 2.0;
473  if(b > s11)
474  s11 = b;
475  b = (p3 + p7) / 2.0;
476  if(b > s12)
477  s12 = b;
478  s1 = s1 - (p1 + p3) / 2.0;
479  s2 = s2 - (p1 + p2) / 2.0;
480  s3 = s3 - (p2 + p4) / 2.0;
481  s4 = s4 - (p3 + p4) / 2.0;
482  s5 = s5 - (p5 + p7) / 2.0;
483  s6 = s6 - (p5 + p6) / 2.0;
484  s7 = s7 - (p6 + p8) / 2.0;
485  s8 = s8 - (p7 + p8) / 2.0;
486  s9 = s9 - (p2 + p6) / 2.0;
487  s10 = s10 - (p4 + p8) / 2.0;
488  s11 = s11 - (p1 + p5) / 2.0;
489  s12 = s12 - (p3 + p7) / 2.0;
490  b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
491  if(b > r1)
492  r1 = b;
493  b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
494  if(b > r2)
495  r2 = b;
496  b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
497  if(b > r3)
498  r3 = b;
499  b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
500  if(b > r4)
501  r4 = b;
502  b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
503  if(b > r5)
504  r5 = b;
505  b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
506  if(b > r6)
507  r6 = b;
508  r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
509  r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
510  r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
511  r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
512  r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
513  r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
514  b = (r1 + r2) / 2.0 + (r3 + r4) / 2.0 + (r5 + r6) / 2.0 + (s1 + s3 + s5 + s7) / 4.0 + (s2 + s4 + s6 + s8) / 4.0 + (s9 + s10 + s11 + s12) / 4.0 + (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.0;
515  if(b < a)
516  a = b;
517  working_space[x][y][z] = a;
518  }
519  }
520  }
521  for (z = q3; z < ssizez - q3; z++) {
522  for (y = q2; y < ssizey - q2; y++) {
523  for (x = q1; x < ssizex - q1; x++) {
524  spectrum[x][y][z] = working_space[x][y][z];
525  }
526  }
527  }
528  }
529  }
530 
531  else if (filterType == kBackOneStepFiltering) {
532  for (i = 1; i <= sampling; i++) {
533  q1 = (Int_t) TMath::Min(i, numberIterationsX), q2 =(Int_t) TMath::Min(i, numberIterationsY), q3 =(Int_t) TMath::Min(i, numberIterationsZ);
534  for (z = q3; z < ssizez - q3; z++) {
535  for (y = q2; y < ssizey - q2; y++) {
536  for (x = q1; x < ssizex - q1; x++) {
537  a = spectrum[x][y][z];
538  p1 = spectrum[x + q1][y + q2][z - q3];
539  p2 = spectrum[x - q1][y + q2][z - q3];
540  p3 = spectrum[x + q1][y - q2][z - q3];
541  p4 = spectrum[x - q1][y - q2][z - q3];
542  p5 = spectrum[x + q1][y + q2][z + q3];
543  p6 = spectrum[x - q1][y + q2][z + q3];
544  p7 = spectrum[x + q1][y - q2][z + q3];
545  p8 = spectrum[x - q1][y - q2][z + q3];
546  s1 = spectrum[x + q1][y ][z - q3];
547  s2 = spectrum[x ][y + q2][z - q3];
548  s3 = spectrum[x - q1][y ][z - q3];
549  s4 = spectrum[x ][y - q2][z - q3];
550  s5 = spectrum[x + q1][y ][z + q3];
551  s6 = spectrum[x ][y + q2][z + q3];
552  s7 = spectrum[x - q1][y ][z + q3];
553  s8 = spectrum[x ][y - q2][z + q3];
554  s9 = spectrum[x - q1][y + q2][z ];
555  s10 = spectrum[x - q1][y - q2][z ];
556  s11 = spectrum[x + q1][y + q2][z ];
557  s12 = spectrum[x + q1][y - q2][z ];
558  r1 = spectrum[x ][y ][z - q3];
559  r2 = spectrum[x ][y ][z + q3];
560  r3 = spectrum[x - q1][y ][z ];
561  r4 = spectrum[x + q1][y ][z ];
562  r5 = spectrum[x ][y + q2][z ];
563  r6 = spectrum[x ][y - q2][z ];
564  b=(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8 - (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4 + (r1 + r2 + r3 + r4 + r5 + r6) / 2;
565  c = -(s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4 + (r1 + r2 + r3 + r4 + r5 + r6) / 2;
566  d = -(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8 + (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 12;
567  if(b < a && b >= 0 && c >=0 && d >= 0)
568  a = b;
569  working_space[x][y][z] = a;
570  }
571  }
572  }
573  for (z = q3; z < ssizez - q3; z++) {
574  for (y = q2; y < ssizey - q2; y++) {
575  for (x = q1; x < ssizex - q1; x++) {
576  spectrum[x][y][z] = working_space[x][y][z];
577  }
578  }
579  }
580  }
581  }
582  }
583 
584  else if (direction == kBackDecreasingWindow) {
585  if (filterType == kBackSuccessiveFiltering) {
586  for (i = sampling; i >= 1; i--) {
587  q1 = (Int_t) TMath::Min(i, numberIterationsX), q2 =(Int_t) TMath::Min(i, numberIterationsY), q3 =(Int_t) TMath::Min(i, numberIterationsZ);
588  for (z = q3; z < ssizez - q3; z++) {
589  for (y = q2; y < ssizey - q2; y++) {
590  for (x = q1; x < ssizex - q1; x++) {
591  a = spectrum[x][y][z];
592  p1 = spectrum[x + q1][y + q2][z - q3];
593  p2 = spectrum[x - q1][y + q2][z - q3];
594  p3 = spectrum[x + q1][y - q2][z - q3];
595  p4 = spectrum[x - q1][y - q2][z - q3];
596  p5 = spectrum[x + q1][y + q2][z + q3];
597  p6 = spectrum[x - q1][y + q2][z + q3];
598  p7 = spectrum[x + q1][y - q2][z + q3];
599  p8 = spectrum[x - q1][y - q2][z + q3];
600  s1 = spectrum[x + q1][y ][z - q3];
601  s2 = spectrum[x ][y + q2][z - q3];
602  s3 = spectrum[x - q1][y ][z - q3];
603  s4 = spectrum[x ][y - q2][z - q3];
604  s5 = spectrum[x + q1][y ][z + q3];
605  s6 = spectrum[x ][y + q2][z + q3];
606  s7 = spectrum[x - q1][y ][z + q3];
607  s8 = spectrum[x ][y - q2][z + q3];
608  s9 = spectrum[x - q1][y + q2][z ];
609  s10 = spectrum[x - q1][y - q2][z ];
610  s11 = spectrum[x + q1][y + q2][z ];
611  s12 = spectrum[x + q1][y - q2][z ];
612  r1 = spectrum[x ][y ][z - q3];
613  r2 = spectrum[x ][y ][z + q3];
614  r3 = spectrum[x - q1][y ][z ];
615  r4 = spectrum[x + q1][y ][z ];
616  r5 = spectrum[x ][y + q2][z ];
617  r6 = spectrum[x ][y - q2][z ];
618  b = (p1 + p3) / 2.0;
619  if(b > s1)
620  s1 = b;
621  b = (p1 + p2) / 2.0;
622  if(b > s2)
623  s2 = b;
624  b = (p2 + p4) / 2.0;
625  if(b > s3)
626  s3 = b;
627  b = (p3 + p4) / 2.0;
628  if(b > s4)
629  s4 = b;
630  b = (p5 + p7) / 2.0;
631  if(b > s5)
632  s5 = b;
633  b = (p5 + p6) / 2.0;
634  if(b > s6)
635  s6 = b;
636  b = (p6 + p8) / 2.0;
637  if(b > s7)
638  s7 = b;
639  b = (p7 + p8) / 2.0;
640  if(b > s8)
641  s8 = b;
642  b = (p2 + p6) / 2.0;
643  if(b > s9)
644  s9 = b;
645  b = (p4 + p8) / 2.0;
646  if(b > s10)
647  s10 = b;
648  b = (p1 + p5) / 2.0;
649  if(b > s11)
650  s11 = b;
651  b = (p3 + p7) / 2.0;
652  if(b > s12)
653  s12 = b;
654  s1 = s1 - (p1 + p3) / 2.0;
655  s2 = s2 - (p1 + p2) / 2.0;
656  s3 = s3 - (p2 + p4) / 2.0;
657  s4 = s4 - (p3 + p4) / 2.0;
658  s5 = s5 - (p5 + p7) / 2.0;
659  s6 = s6 - (p5 + p6) / 2.0;
660  s7 = s7 - (p6 + p8) / 2.0;
661  s8 = s8 - (p7 + p8) / 2.0;
662  s9 = s9 - (p2 + p6) / 2.0;
663  s10 = s10 - (p4 + p8) / 2.0;
664  s11 = s11 - (p1 + p5) / 2.0;
665  s12 = s12 - (p3 + p7) / 2.0;
666  b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
667  if(b > r1)
668  r1 = b;
669  b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
670  if(b > r2)
671  r2 = b;
672  b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
673  if(b > r3)
674  r3 = b;
675  b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
676  if(b > r4)
677  r4 = b;
678  b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
679  if(b > r5)
680  r5 = b;
681  b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
682  if(b > r6)
683  r6 = b;
684  r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
685  r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
686  r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
687  r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
688  r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
689  r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
690  b = (r1 + r2) / 2.0 + (r3 + r4) / 2.0 + (r5 + r6) / 2.0 + (s1 + s3 + s5 + s7) / 4.0 + (s2 + s4 + s6 + s8) / 4.0 + (s9 + s10 + s11 + s12) / 4.0 + (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.0;
691  if(b < a)
692  a = b;
693  working_space[x][y][z] = a;
694  }
695  }
696  }
697  for (z = q3; z < ssizez - q3; z++) {
698  for (y = q2; y < ssizey - q2; y++) {
699  for (x = q1; x < ssizex - q1; x++) {
700  spectrum[x][y][z] = working_space[x][y][z];
701  }
702  }
703  }
704  }
705  }
706 
707  else if (filterType == kBackOneStepFiltering) {
708  for (i = sampling; i >= 1; i--) {
709  q1 = (Int_t) TMath::Min(i, numberIterationsX), q2 =(Int_t) TMath::Min(i, numberIterationsY), q3 =(Int_t) TMath::Min(i, numberIterationsZ);
710  for (z = q3; z < ssizez - q3; z++) {
711  for (y = q2; y < ssizey - q2; y++) {
712  for (x = q1; x < ssizex - q1; x++) {
713  a = spectrum[x][y][z];
714  p1 = spectrum[x + q1][y + q2][z - q3];
715  p2 = spectrum[x - q1][y + q2][z - q3];
716  p3 = spectrum[x + q1][y - q2][z - q3];
717  p4 = spectrum[x - q1][y - q2][z - q3];
718  p5 = spectrum[x + q1][y + q2][z + q3];
719  p6 = spectrum[x - q1][y + q2][z + q3];
720  p7 = spectrum[x + q1][y - q2][z + q3];
721  p8 = spectrum[x - q1][y - q2][z + q3];
722  s1 = spectrum[x + q1][y ][z - q3];
723  s2 = spectrum[x ][y + q2][z - q3];
724  s3 = spectrum[x - q1][y ][z - q3];
725  s4 = spectrum[x ][y - q2][z - q3];
726  s5 = spectrum[x + q1][y ][z + q3];
727  s6 = spectrum[x ][y + q2][z + q3];
728  s7 = spectrum[x - q1][y ][z + q3];
729  s8 = spectrum[x ][y - q2][z + q3];
730  s9 = spectrum[x - q1][y + q2][z ];
731  s10 = spectrum[x - q1][y - q2][z ];
732  s11 = spectrum[x + q1][y + q2][z ];
733  s12 = spectrum[x + q1][y - q2][z ];
734  r1 = spectrum[x ][y ][z - q3];
735  r2 = spectrum[x ][y ][z + q3];
736  r3 = spectrum[x - q1][y ][z ];
737  r4 = spectrum[x + q1][y ][z ];
738  r5 = spectrum[x ][y + q2][z ];
739  r6 = spectrum[x ][y - q2][z ];
740  b = (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8 - (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4 + (r1 + r2 + r3 + r4 + r5 + r6) / 2;
741  c = -(s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 4+(r1 + r2 + r3 + r4 + r5 + r6) / 2;
742  d = -(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8)/8 + (s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 + s9 + s10 + s11 + s12) / 12;
743  if(b < a && b >= 0 && c >=0 && d >= 0)
744  a = b;
745  working_space[x][y][z] = a;
746  }
747  }
748  }
749  for (z = q3; z < ssizez - q3; z++) {
750  for (y = q2; y < ssizey - q2; y++) {
751  for (x = q1; x < ssizex - q1; x++) {
752  spectrum[x][y][z] = working_space[x][y][z];
753  }
754  }
755  }
756  }
757  }
758  }
759  for(i = 0;i < ssizex; i++){
760  for(j = 0;j < ssizey; j++)
761  delete[] working_space[i][j];
762  delete[] working_space[i];
763  }
764  delete[] working_space;
765  return 0;
766 }
767 
768 ////////////////////////////////////////////////////////////////////////////////
769 /// This function calculates smoothed spectrum from source spectrum
770 /// based on Markov chain method.
771 /// The result is placed in the array pointed by spectrum pointer.
772 ///
773 /// Function parameters:
774 /// - source-pointer to the array of source spectrum
775 /// - working_space-pointer to the working array
776 /// - ssizex-x length of spectrum and working space arrays
777 /// - ssizey-y length of spectrum and working space arrays
778 /// - ssizey-z length of spectrum and working space arrays
779 /// - averWindow-width of averaging smoothing window
780 ///
781 /// ### Smoothing
782 ///
783 /// Goal: Suppression of statistical fluctuations
784 /// the algorithm is based on discrete Markov chain, which has very simple
785 /// invariant distribution
786 ///
787 /// \f[
788 /// 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
789 /// \f]
790 /// \f$U_1\f$ being defined from the normalization condition \f$ \sum_{i=1}^{n} U_i = 1\f$
791 /// n is the length of the smoothed spectrum and
792 /// \f[
793 /// 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]
794 /// \f]
795 ///
796 /// is the probability of the change of the peak position from channel i to the channel i+1.
797 /// \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
798 /// of smoothing window. We have extended this algorithm to three dimensions.
799 ///
800 /// #### Reference:
801 ///
802 /// [1] Z.K. Silagadze, A new
803 /// algorithm for automatic photo-peak searches. NIM A 376 (1996), 451-.
804 ///
805 /// ### Example 1 - script SmootMarkov3.c :
806 ///
807 /// \image html spectrum3_smoothing_image007.jpg Fig. 1 Original noisy spectrum.
808 /// \image html spectrum3_smoothing_image008.jpg Fig. 2 Smoothed spectrum with averaging window m=3.
809 ///
810 /// #### Script:
811 ///
812 /// Example to illustrate the Markov smoothing (class TSpectrum3).
813 /// To execute this example, do:
814 ///
815 /// `root > .x SmoothMarkov3.C`
816 ///
817 /// ~~~ {.cpp}
818 /// void SmoothMarkov3() {
819 /// Int_t i, j, k;
820 /// Int_t nbinsx = 64;
821 /// Int_t nbinsy = 64;
822 /// Int_t nbinsz = 64;
823 /// Int_t xmin = 0;
824 /// Int_t xmax = nbinsx;
825 /// Int_t ymin = 0;
826 /// Int_t ymax = nbinsy;
827 /// Int_t zmin = 0;
828 /// Int_t zmax = nbinsz;
829 /// Double_t*** source = new Double_t**[nbinsx];
830 /// for(i=0;i<nbinsx;i++){
831 /// source[i]=new Double_t*[nbinsy];
832 /// for(j=0;j<nbinsy;j++)
833 /// source[i][j]=new Double_t[nbinsz];
834 /// }
835 /// TH3F *sm = new TH3F("Smoothing","Markov smoothing",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
836 /// TFile *f = new TFile("TSpectrum3.root");
837 /// sm=(TH3F*)f->Get("back;1");
838 /// TCanvas *Background = new TCanvas("Smoothing","Markov smoothing",10,10,1000,700);
839 /// TSpectrum3 *s = new TSpectrum3();
840 /// for (i = 0; i < nbinsx; i++){
841 /// for (j = 0; j < nbinsy; j++){
842 /// for (k = 0; k < nbinsz; k++){
843 /// source[i][j][k] = sm->GetBinContent(i + 1,j + 1,k + 1);
844 /// }
845 /// }
846 /// }
847 /// s->SmoothMarkov(source,nbinsx,nbinsy,nbinsz,3);
848 /// for (i = 0; i < nbinsx; i++){
849 /// for (j = 0; j < nbinsy; j++){
850 /// for (k = 0; k < nbinsz; k++){
851 /// sm->SetBinContent(i + 1,j + 1,k + 1, source[i][j][k]);
852 /// }
853 /// }
854 /// }
855 /// sm->Draw("");
856 /// }
857 /// ~~~
858 
859 const char* TSpectrum3::SmoothMarkov(Double_t***source, Int_t ssizex, Int_t ssizey, Int_t ssizez, Int_t averWindow)
860 {
861  Int_t xmin,xmax,ymin,ymax,zmin,zmax,i,j,k,l;
862  Double_t a,b,maxch;
863  Double_t nom,nip,nim,sp,sm,spx,smx,spy,smy,spz,smz,plocha=0;
864  if(averWindow<=0)
865  return "Averaging Window must be positive";
866  Double_t***working_space = new Double_t**[ssizex];
867  for(i = 0;i < ssizex; i++){
868  working_space[i] = new Double_t*[ssizey];
869  for(j = 0;j < ssizey; j++)
870  working_space[i][j] = new Double_t[ssizez];
871  }
872  xmin = 0;
873  xmax = ssizex - 1;
874  ymin = 0;
875  ymax = ssizey - 1;
876  zmin = 0;
877  zmax = ssizez - 1;
878  for(i = 0,maxch = 0;i < ssizex; i++){
879  for(j = 0;j < ssizey; j++){
880  for(k = 0;k < ssizez; k++){
881  working_space[i][j][k] = 0;
882  if(maxch < source[i][j][k])
883  maxch = source[i][j][k];
884  plocha += source[i][j][k];
885  }
886  }
887  }
888  if(maxch == 0) {
889  delete [] working_space;
890  return 0;
891  }
892 
893  nom = 0;
894  working_space[xmin][ymin][zmin] = 1;
895  for(i = xmin;i < xmax; i++){
896  nip = source[i][ymin][zmin] / maxch;
897  nim = source[i + 1][ymin][zmin] / maxch;
898  sp = 0,sm = 0;
899  for(l = 1;l <= averWindow; l++){
900  if((i + l) > xmax)
901  a = source[xmax][ymin][zmin] / maxch;
902 
903  else
904  a = source[i + l][ymin][zmin] / maxch;
905 
906  b = a - nip;
907  if(a + nip <= 0)
908  a = 1;
909 
910  else
911  a = TMath::Sqrt(a + nip);
912 
913  b = b / a;
914  b = TMath::Exp(b);
915  sp = sp + b;
916  if(i - l + 1 < xmin)
917  a = source[xmin][ymin][zmin] / maxch;
918 
919  else
920  a = source[i - l + 1][ymin][zmin] / maxch;
921 
922  b = a - nim;
923  if(a + nim <= 0)
924  a = 1;
925 
926  else
927  a = TMath::Sqrt(a + nim);
928 
929  b = b / a;
930  b = TMath::Exp(b);
931  sm = sm + b;
932  }
933  a = sp / sm;
934  a = working_space[i + 1][ymin][zmin] = a * working_space[i][ymin][zmin];
935  nom = nom + a;
936  }
937  for(i = ymin;i < ymax; i++){
938  nip = source[xmin][i][zmin] / maxch;
939  nim = source[xmin][i + 1][zmin] / maxch;
940  sp = 0,sm = 0;
941  for(l = 1;l <= averWindow; l++){
942  if((i + l) > ymax)
943  a = source[xmin][ymax][zmin] / maxch;
944 
945  else
946  a = source[xmin][i + l][zmin] / maxch;
947 
948  b = a - nip;
949  if(a + nip <= 0)
950  a = 1;
951 
952  else
953  a = TMath::Sqrt(a + nip);
954 
955  b = b / a;
956  b = TMath::Exp(b);
957  sp = sp + b;
958  if(i - l + 1 < ymin)
959  a = source[xmin][ymin][zmin] / maxch;
960 
961  else
962  a = source[xmin][i - l + 1][zmin] / maxch;
963 
964  b = a - nim;
965  if(a + nim <= 0)
966  a = 1;
967 
968  else
969  a = TMath::Sqrt(a + nim);
970 
971  b = b / a;
972  b = TMath::Exp(b);
973  sm = sm + b;
974  }
975  a = sp / sm;
976  a = working_space[xmin][i + 1][zmin] = a * working_space[xmin][i][zmin];
977  nom = nom + a;
978  }
979  for(i = zmin;i < zmax; i++){
980  nip = source[xmin][ymin][i] / maxch;
981  nim = source[xmin][ymin][i + 1] / maxch;
982  sp = 0,sm = 0;
983  for(l = 1;l <= averWindow; l++){
984  if((i + l) > zmax)
985  a = source[xmin][ymin][zmax] / maxch;
986 
987  else
988  a = source[xmin][ymin][i + l] / maxch;
989 
990  b = a - nip;
991  if(a + nip <= 0)
992  a = 1;
993 
994  else
995  a = TMath::Sqrt(a + nip);
996 
997  b = b / a;
998  b = TMath::Exp(b);
999  sp = sp + b;
1000  if(i - l + 1 < zmin)
1001  a = source[xmin][ymin][zmin] / maxch;
1002 
1003  else
1004  a = source[xmin][ymin][i - l + 1] / maxch;
1005 
1006  b = a - nim;
1007  if(a + nim <= 0)
1008  a = 1;
1009 
1010  else
1011  a = TMath::Sqrt(a + nim);
1012  b = b / a;
1013  b = TMath::Exp(b);
1014  sm = sm + b;
1015  }
1016  a = sp / sm;
1017  a = working_space[xmin][ymin][i + 1] = a * working_space[xmin][ymin][i];
1018  nom = nom + a;
1019  }
1020  for(i = xmin;i < xmax; i++){
1021  for(j = ymin;j < ymax; j++){
1022  nip = source[i][j + 1][zmin] / maxch;
1023  nim = source[i + 1][j + 1][zmin] / maxch;
1024  spx = 0,smx = 0;
1025  for(l = 1;l <= averWindow; l++){
1026  if(i + l > xmax)
1027  a = source[xmax][j][zmin] / maxch;
1028 
1029  else
1030  a = source[i + l][j][zmin] / maxch;
1031 
1032  b = a - nip;
1033  if(a + nip <= 0)
1034  a = 1;
1035 
1036  else
1037  a = TMath::Sqrt(a + nip);
1038 
1039  b = b / a;
1040  b = TMath::Exp(b);
1041  spx = spx + b;
1042  if(i - l + 1 < xmin)
1043  a = source[xmin][j][zmin] / maxch;
1044 
1045  else
1046  a = source[i - l + 1][j][zmin] / maxch;
1047 
1048  b = a - nim;
1049  if(a + nim <= 0)
1050  a = 1;
1051 
1052  else
1053  a = TMath::Sqrt(a + nim);
1054 
1055  b = b / a;
1056  b = TMath::Exp(b);
1057  smx = smx + b;
1058  }
1059  spy = 0,smy = 0;
1060  nip = source[i + 1][j][zmin] / maxch;
1061  nim = source[i + 1][j + 1][zmin] / maxch;
1062  for(l = 1;l <= averWindow; l++){
1063  if(j + l > ymax)
1064  a = source[i][ymax][zmin] / maxch;
1065 
1066  else
1067  a = source[i][j + l][zmin] / maxch;
1068 
1069  b = a - nip;
1070  if(a + nip <= 0)
1071  a = 1;
1072 
1073  else
1074  a = TMath::Sqrt(a + nip);
1075 
1076  b = b / a;
1077  b = TMath::Exp(b);
1078  spy = spy + b;
1079  if(j - l + 1 < ymin)
1080  a = source[i][ymin][zmin] / maxch;
1081 
1082  else
1083  a = source[i][j - l + 1][zmin] / maxch;
1084 
1085  b = a - nim;
1086  if(a + nim <= 0)
1087  a = 1;
1088 
1089  else
1090  a = TMath::Sqrt(a + nim);
1091 
1092  b = b / a;
1093  b = TMath::Exp(b);
1094  smy = smy + b;
1095  }
1096  a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
1097  working_space[i + 1][j + 1][zmin] = a;
1098  nom = nom + a;
1099  }
1100  }
1101  for(i = xmin;i < xmax; i++){
1102  for(j = zmin;j < zmax; j++){
1103  nip = source[i][ymin][j + 1] / maxch;
1104  nim = source[i + 1][ymin][j + 1] / maxch;
1105  spx = 0,smx = 0;
1106  for(l = 1;l <= averWindow; l++){
1107  if(i + l > xmax)
1108  a = source[xmax][ymin][j] / maxch;
1109 
1110  else
1111  a = source[i + l][ymin][j] / maxch;
1112 
1113  b = a - nip;
1114  if(a + nip <= 0)
1115  a = 1;
1116 
1117  else
1118  a = TMath::Sqrt(a + nip);
1119 
1120  b = b / a;
1121  b = TMath::Exp(b);
1122  spx = spx + b;
1123  if(i - l + 1 < xmin)
1124  a = source[xmin][ymin][j] / maxch;
1125 
1126  else
1127  a = source[i - l + 1][ymin][j] / maxch;
1128 
1129  b = a - nim;
1130  if(a + nim <= 0)
1131  a = 1;
1132 
1133  else
1134  a = TMath::Sqrt(a + nim);
1135 
1136  b = b / a;
1137  b = TMath::Exp(b);
1138  smx = smx + b;
1139  }
1140  spz = 0,smz = 0;
1141  nip = source[i + 1][ymin][j] / maxch;
1142  nim = source[i + 1][ymin][j + 1] / maxch;
1143  for(l = 1;l <= averWindow; l++){
1144  if(j + l > zmax)
1145  a = source[i][ymin][zmax] / maxch;
1146 
1147  else
1148  a = source[i][ymin][j + l] / maxch;
1149 
1150  b = a - nip;
1151  if(a + nip <= 0)
1152  a = 1;
1153 
1154  else
1155  a = TMath::Sqrt(a + nip);
1156 
1157  b = b / a;
1158  b = TMath::Exp(b);
1159  spz = spz + b;
1160  if(j - l + 1 < zmin)
1161  a = source[i][ymin][zmin] / maxch;
1162 
1163  else
1164  a = source[i][ymin][j - l + 1] / maxch;
1165 
1166  b = a - nim;
1167  if(a + nim <= 0)
1168  a = 1;
1169 
1170  else
1171  a = TMath::Sqrt(a + nim);
1172 
1173  b = b / a;
1174  b = TMath::Exp(b);
1175  smz = smz + b;
1176  }
1177  a = (spx * working_space[i][ymin][j + 1] + spz * working_space[i + 1][ymin][j]) / (smx + smz);
1178  working_space[i + 1][ymin][j + 1] = a;
1179  nom = nom + a;
1180  }
1181  }
1182  for(i = ymin;i < ymax; i++){
1183  for(j = zmin;j < zmax; j++){
1184  nip = source[xmin][i][j + 1] / maxch;
1185  nim = source[xmin][i + 1][j + 1] / maxch;
1186  spy = 0,smy = 0;
1187  for(l = 1;l <= averWindow; l++){
1188  if(i + l > ymax)
1189  a = source[xmin][ymax][j] / maxch;
1190 
1191  else
1192  a = source[xmin][i + l][j] / maxch;
1193 
1194  b = a - nip;
1195  if(a + nip <= 0)
1196  a = 1;
1197 
1198  else
1199  a = TMath::Sqrt(a + nip);
1200 
1201  b = b / a;
1202  b = TMath::Exp(b);
1203  spy = spy + b;
1204  if(i - l + 1 < ymin)
1205  a = source[xmin][ymin][j] / maxch;
1206 
1207  else
1208  a = source[xmin][i - l + 1][j] / maxch;
1209 
1210  b = a - nim;
1211  if(a + nim <= 0)
1212  a = 1;
1213 
1214  else
1215  a = TMath::Sqrt(a + nim);
1216 
1217  b = b / a;
1218  b = TMath::Exp(b);
1219  smy = smy + b;
1220  }
1221  spz = 0,smz = 0;
1222  nip = source[xmin][i + 1][j] / maxch;
1223  nim = source[xmin][i + 1][j + 1] / maxch;
1224  for(l = 1;l <= averWindow; l++){
1225  if(j + l > zmax)
1226  a = source[xmin][i][zmax] / maxch;
1227 
1228  else
1229  a = source[xmin][i][j + l] / maxch;
1230 
1231  b = a - nip;
1232  if(a + nip <= 0)
1233  a = 1;
1234 
1235  else
1236  a = TMath::Sqrt(a + nip);
1237 
1238  b = b / a;
1239  b = TMath::Exp(b);
1240  spz = spz + b;
1241  if(j - l + 1 < zmin)
1242  a = source[xmin][i][zmin] / maxch;
1243 
1244  else
1245  a = source[xmin][i][j - l + 1] / maxch;
1246 
1247  b = a - nim;
1248  if(a + nim <= 0)
1249  a = 1;
1250 
1251  else
1252  a = TMath::Sqrt(a + nim);
1253 
1254  b = b / a;
1255  b = TMath::Exp(b);
1256  smz = smz + b;
1257  }
1258  a = (spy * working_space[xmin][i][j + 1] + spz * working_space[xmin][i + 1][j]) / (smy + smz);
1259  working_space[xmin][i + 1][j + 1] = a;
1260  nom = nom + a;
1261  }
1262  }
1263  for(i = xmin;i < xmax; i++){
1264  for(j = ymin;j < ymax; j++){
1265  for(k = zmin;k < zmax; k++){
1266  nip = source[i][j + 1][k + 1] / maxch;
1267  nim = source[i + 1][j + 1][k + 1] / maxch;
1268  spx = 0,smx = 0;
1269  for(l = 1;l <= averWindow; l++){
1270  if(i + l > xmax)
1271  a = source[xmax][j][k] / maxch;
1272 
1273  else
1274  a = source[i + l][j][k] / maxch;
1275 
1276  b = a - nip;
1277  if(a + nip <= 0)
1278  a = 1;
1279 
1280  else
1281  a = TMath::Sqrt(a + nip);
1282 
1283  b = b / a;
1284  b = TMath::Exp(b);
1285  spx = spx + b;
1286  if(i - l + 1 < xmin)
1287  a = source[xmin][j][k] / maxch;
1288 
1289  else
1290  a = source[i - l + 1][j][k] / maxch;
1291 
1292  b = a - nim;
1293  if(a + nim <= 0)
1294  a = 1;
1295 
1296  else
1297  a = TMath::Sqrt(a + nim);
1298 
1299  b = b / a;
1300  b = TMath::Exp(b);
1301  smx = smx + b;
1302  }
1303  spy = 0,smy = 0;
1304  nip = source[i + 1][j][k + 1] / maxch;
1305  nim = source[i + 1][j + 1][k + 1] / maxch;
1306  for(l = 1;l <= averWindow; l++){
1307  if(j + l > ymax)
1308  a = source[i][ymax][k] / maxch;
1309 
1310  else
1311  a = source[i][j + l][k] / maxch;
1312 
1313  b = a - nip;
1314  if(a + nip <= 0)
1315  a = 1;
1316 
1317  else
1318  a = TMath::Sqrt(a + nip);
1319 
1320  b = b / a;
1321  b = TMath::Exp(b);
1322  spy = spy + b;
1323  if(j - l + 1 < ymin)
1324  a = source[i][ymin][k] / maxch;
1325 
1326  else
1327  a = source[i][j - l + 1][k] / maxch;
1328 
1329  b = a - nim;
1330  if(a + nim <= 0)
1331  a = 1;
1332 
1333  else
1334  a = TMath::Sqrt(a + nim);
1335 
1336  b = b / a;
1337  b = TMath::Exp(b);
1338  smy = smy + b;
1339  }
1340  spz = 0,smz = 0;
1341  nip = source[i + 1][j + 1][k] / maxch;
1342  nim = source[i + 1][j + 1][k + 1] / maxch;
1343  for(l = 1;l <= averWindow; l++){
1344  if(j + l > zmax)
1345  a = source[i][j][zmax] / maxch;
1346 
1347  else
1348  a = source[i][j][k + l] / maxch;
1349 
1350  b = a - nip;
1351  if(a + nip <= 0)
1352  a = 1;
1353 
1354  else
1355  a = TMath::Sqrt(a + nip);
1356 
1357  b = b / a;
1358  b = TMath::Exp(b);
1359  spz = spz + b;
1360  if(j - l + 1 < ymin)
1361  a = source[i][j][zmin] / maxch;
1362 
1363  else
1364  a = source[i][j][k - l + 1] / maxch;
1365 
1366  b = a - nim;
1367  if(a + nim <= 0)
1368  a = 1;
1369 
1370  else
1371  a = TMath::Sqrt(a + nim);
1372 
1373  b = b / a;
1374  b = TMath::Exp(b);
1375  smz = smz+b;
1376  }
1377  a = (spx * working_space[i][j + 1][k + 1] + spy * working_space[i + 1][j][k + 1] + spz * working_space[i + 1][j + 1][k]) / (smx + smy + smz);
1378  working_space[i + 1][j + 1][k + 1] = a;
1379  nom = nom + a;
1380  }
1381  }
1382  }
1383  for(i = xmin;i <= xmax; i++){
1384  for(j = ymin;j <= ymax; j++){
1385  for(k = zmin;k <= zmax; k++){
1386  working_space[i][j][k] = working_space[i][j][k] / nom;
1387  }
1388  }
1389  }
1390  for(i = 0;i < ssizex; i++){
1391  for(j = 0;j < ssizey; j++){
1392  for(k = 0;k < ssizez; k++){
1393  source[i][j][k] = plocha * working_space[i][j][k];
1394  }
1395  }
1396  }
1397  for(i = 0;i < ssizex; i++){
1398  for(j = 0;j < ssizey; j++)
1399  delete[] working_space[i][j];
1400  delete[] working_space[i];
1401  }
1402  delete[] working_space;
1403  return 0;
1404 }
1405 
1406 ////////////////////////////////////////////////////////////////////////////////
1407 /// This function calculates deconvolution from source spectrum
1408 /// according to response spectrum
1409 /// The result is placed in the cube pointed by source pointer.
1410 ///
1411 /// Function parameters:
1412 /// - source-pointer to the cube of source spectrum
1413 /// - resp-pointer to the cube of response spectrum
1414 /// - ssizex-x length of source and response spectra
1415 /// - ssizey-y length of source and response spectra
1416 /// - ssizey-y length of source and response spectra
1417 /// - numberIterations, for details we refer to manual
1418 /// - numberRepetitions, for details we refer to manual
1419 /// - boost, boosting factor, for details we refer to manual
1420 ///
1421 /// ### Deconvolution
1422 ///
1423 /// Goal: Improvement of the resolution in spectra, decomposition of multiplets
1424 ///
1425 /// Mathematical formulation of the 3-dimensional convolution system is
1426 ///
1427 /// \image html spectrum3_deconvolution_image001.gif
1428 ///
1429 /// where h(i,j,k) is the impulse response function, x, y are
1430 /// input and output fields, respectively, \f$ N_1, N_2, N3\f$, are the lengths of x and h fields
1431 ///
1432 /// - let us assume that we know the response and the output fields (spectra)
1433 /// of the above given system.
1434 ///
1435 /// - the deconvolution represents solution of the overdetermined system of
1436 /// linear equations, i.e., the calculation of the field -x.
1437 ///
1438 /// - from numerical stability point of view the operation of deconvolution is
1439 /// extremely critical (ill-posed problem) as well as time consuming operation.
1440 ///
1441 /// - the Gold deconvolution algorithm proves to work very well even for
1442 /// 2-dimensional systems. Generalization of the algorithm for 2-dimensional
1443 /// systems was presented in [1], and for multidimensional systems in [2].
1444 ///
1445 /// - for Gold deconvolution algorithm as well as for boosted deconvolution
1446 /// algorithm we refer also to TSpectrum and TSpectrum2
1447 ///
1448 /// #### References:
1449 ///
1450 /// [1] M.Morhac, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.: Efficient
1451 /// one- and two-dimensional Gold deconvolution and its application to gamma-ray
1452 /// spectra decomposition. NIM, A401 (1997) 385-408.
1453 ///
1454 /// [2] Morhac M., Matouoek V.,
1455 /// Kliman J., Efficient algorithm of multidimensional deconvolution and its
1456 /// application to nuclear data processing, Digital Signal Processing 13 (2003) 144.
1457 ///
1458 /// ### Example 1 - script Decon.c :
1459 ///
1460 /// response function (usually peak) should be shifted to the beginning of
1461 /// the coordinate system (see Fig. 1)
1462 ///
1463 /// \image html spectrum3_deconvolution_image003.jpg Fig. 1 Three-dimensional response spectrum
1464 /// \image html spectrum3_deconvolution_image004.jpg Fig. 2 Three-dimensional input spectrum (before deconvolution)
1465 /// \image html spectrum3_deconvolution_image005.jpg Fig. 3 Spectrum from Fig. 2 after deconvolution (100 iterations)
1466 ///
1467 /// #### Script:
1468 ///
1469 /// Example to illustrate the Gold deconvolution (class TSpectrum3).
1470 /// To execute this example, do:
1471 ///
1472 /// `root > .x Decon3.C`
1473 ///
1474 /// ~~~ {.cpp}
1475 /// #include <TSpectrum3>
1476 /// void Decon3() {
1477 /// Int_t i, j, k;
1478 /// Int_t nbinsx = 32;
1479 /// Int_t nbinsy = 32;
1480 /// Int_t nbinsz = 32;
1481 /// Int_t xmin = 0;
1482 /// Int_t xmax = nbinsx;
1483 /// Int_t ymin = 0;
1484 /// Int_t ymax = nbinsy;
1485 /// Int_t zmin = 0;
1486 /// Int_t zmax = nbinsz;
1487 /// Double_t*** source = newDouble_t**[nbinsx];
1488 /// Double_t*** resp = new Double_t**[nbinsx];
1489 /// for(i=0;i<nbinsx;i++){
1490 /// source[i]=new Double_t* [nbinsy];
1491 /// for(j=0;j<nbinsy;j++)
1492 /// source[i][j]=new Double_t[nbinsz];
1493 /// }
1494 /// for(i=0;i<nbinsx;i++){
1495 /// resp[i]=new Double_t*[nbinsy];
1496 /// for(j=0;j<nbinsy;j++)
1497 /// resp[i][j]=new Double_t[nbinsz];
1498 /// }
1499 /// TH3F *decon_in = new TH3F("decon_in","Deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1500 /// TH3F *decon_resp = new TH3F("decon_resp","Deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1501 /// TFile *f = new TFile("TSpectrum3.root");
1502 /// decon_in=(TH3F*) f->Get("decon_in;1");
1503 /// decon_resp=(TH3F*) f->Get("decon_resp;1");
1504 /// TCanvas *Deconvolution = new TCanvas("Deconvolution","Deconvolution of 3-dimensional spectra",10,10,1000,700);
1505 /// TSpectrum3 *s = new TSpectrum3();
1506 /// for (i = 0; i < nbinsx; i++){
1507 /// for (j = 0; j < nbinsy; j++){
1508 /// for (k = 0; k < nbinsz; k++){
1509 /// source[i][j][k] = decon_in->GetBinContent(i + 1,j + 1,k + 1);
1510 /// resp[i][j][k] = decon_resp->GetBinContent(i + 1,j + 1,k + 1);
1511 /// }
1512 /// }
1513 /// }
1514 /// s->Deconvolution(source,resp,nbinsx,nbinsy,nbinsz,100,1,1);
1515 /// for (i = 0; i < nbinsx; i++){
1516 /// for (j = 0; j < nbinsy; j++){
1517 /// for (k = 0; k < nbinsz; k++){
1518 /// decon_in->SetBinContent(i + 1,j + 1,k + 1, source[i][j][k]);
1519 /// }
1520 /// }
1521 /// }
1522 /// decon_in->Draw("");
1523 /// }
1524 /// ~~~
1525 ///
1526 /// ### Example 2 - script Decon_hr.c :
1527 ///
1528 /// This example illustrates repeated
1529 /// Gold deconvolution with boosting. After every 10 iterations we apply power
1530 /// function with exponent = 2 to the spectrum given in Fig. 2.
1531 ///
1532 /// \image html spectrum3_deconvolution_image006.jpg Fig. 4 Spectrum from Fig. 2 after boosted deconvolution (10 iterations repeated 10 times). It decomposes completely cluster of peaks from Fig 2.
1533 ///
1534 /// #### Script:
1535 ///
1536 /// Example to illustrate the Gold deconvolution (class TSpectrum3).
1537 /// To execute this example, do:
1538 ///
1539 /// `root > .x Decon3_hr.C`
1540 ///
1541 /// ~~~ {.cpp}
1542 /// void Decon3_hr() {
1543 /// Int_t i, j, k;
1544 /// Int_t nbinsx = 32;
1545 /// Int_t nbinsy = 32;
1546 /// Int_t nbinsz = 32;
1547 /// Int_t xmin = 0;
1548 /// Int_t xmax = nbinsx;
1549 /// Int_t ymin = 0;
1550 /// Int_t ymax = nbinsy;
1551 /// Int_t zmin = 0;
1552 /// Int_t zmax = nbinsz;
1553 /// Double_t*** source = new Double_t**[nbinsx];
1554 /// Double_t*** resp = new Double_t**[nbinsx];
1555 /// for(i=0;i<nbinsx;i++){
1556 /// source[i]=new Double_t*[nbinsy];
1557 /// for(j=0;j<nbinsy;j++)
1558 /// source[i][j]=new Double_t[nbinsz];
1559 /// }
1560 /// for(i=0;i<nbinsx;i++){
1561 /// resp[i]=new Double_t*[nbinsy];
1562 /// for(j=0;j<nbinsy;j++)
1563 /// resp[i][j]=new Double_t[nbinsz];
1564 /// }
1565 /// TH3F *decon_in = new TH3F("decon_in","Deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1566 /// TH3F *decon_resp = new TH3F("decon_resp","Deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1567 /// TFile *f = new TFile("TSpectrum3.root");
1568 /// decon_in=(TH3F*)f->Get("decon_in;1");
1569 /// decon_resp=(TH3F*)f->Get("decon_resp;1");
1570 /// TCanvas *Deconvolution = new TCanvas("Deconvolution","High resolution deconvolution of 3-dimensional spectra",10,10,1000,700);
1571 /// TSpectrum3 *s = new TSpectrum3();
1572 /// for (i = 0; i < nbinsx; i++){
1573 /// for (j = 0; j < nbinsy; j++){
1574 /// for (k = 0; k < nbinsz; k++){
1575 /// source[i][j][k] = decon_in->GetBinContent(i + 1,j + 1,k + 1);
1576 /// resp[i][j][k] = decon_resp->GetBinContent(i + 1,j + 1,k + 1);
1577 /// }
1578 /// }
1579 /// }
1580 /// s->Deconvolution(source,resp,nbinsx,nbinsy,nbinsz,10,10,2);
1581 /// for (i = 0; i < nbinsx; i++){
1582 /// for (j = 0; j < nbinsy; j++){
1583 /// for (k = 0; k < nbinsz; k++){
1584 /// decon_in->SetBinContent(i + 1,j + 1,k + 1, source[i][j][k]);
1585 /// }
1586 /// }
1587 /// }
1588 /// decon_in->Draw("");
1589 /// }
1590 /// ~~~
1591 
1592 const char *TSpectrum3::Deconvolution(Double_t***source, const Double_t***resp,
1593  Int_t ssizex, Int_t ssizey, Int_t ssizez,
1594  Int_t numberIterations,
1595  Int_t numberRepetitions,
1596  Double_t boost)
1597 {
1598  Int_t i, j, k, lhx, lhy, lhz, i1, i2, i3, j1, j2, j3, k1, k2, k3, lindex, i1min, i1max, i2min, i2max, i3min, i3max, j1min, j1max, j2min, j2max, j3min, j3max, positx = 0, posity = 0, positz = 0, repet;
1599  Double_t lda, ldb, ldc, area, maximum = 0;
1600  if (ssizex <= 0 || ssizey <= 0 || ssizez <= 0)
1601  return "Wrong parameters";
1602  if (numberIterations <= 0)
1603  return "Number of iterations must be positive";
1604  if (numberRepetitions <= 0)
1605  return "Number of repetitions must be positive";
1606  Double_t ***working_space=new Double_t** [ssizex];
1607  for(i=0;i<ssizex;i++){
1608  working_space[i]=new Double_t* [ssizey];
1609  for(j=0;j<ssizey;j++)
1610  working_space[i][j]=new Double_t [5*ssizez];
1611  }
1612  area = 0;
1613  lhx = -1, lhy = -1, lhz = -1;
1614  for (i = 0; i < ssizex; i++) {
1615  for (j = 0; j < ssizey; j++) {
1616  for (k = 0; k < ssizez; k++) {
1617  lda = resp[i][j][k];
1618  if (lda != 0) {
1619  if ((i + 1) > lhx)
1620  lhx = i + 1;
1621  if ((j + 1) > lhy)
1622  lhy = j + 1;
1623  if ((k + 1) > lhz)
1624  lhz = k + 1;
1625  }
1626  working_space[i][j][k] = lda;
1627  area = area + lda;
1628  if (lda > maximum) {
1629  maximum = lda;
1630  positx = i, posity = j, positz = k;
1631  }
1632  }
1633  }
1634  }
1635  if (lhx == -1 || lhy == -1 || lhz == -1) {
1636  delete [] working_space;
1637  return ("Zero response data");
1638  }
1639 
1640 //calculate ht*y and write into p
1641  for (i3 = 0; i3 < ssizez; i3++) {
1642  for (i2 = 0; i2 < ssizey; i2++) {
1643  for (i1 = 0; i1 < ssizex; i1++) {
1644  ldc = 0;
1645  for (j3 = 0; j3 <= (lhz - 1); j3++) {
1646  for (j2 = 0; j2 <= (lhy - 1); j2++) {
1647  for (j1 = 0; j1 <= (lhx - 1); j1++) {
1648  k3 = i3 + j3, k2 = i2 + j2, k1 = i1 + j1;
1649  if (k3 >= 0 && k3 < ssizez && k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
1650  lda = working_space[j1][j2][j3];
1651  ldb = source[k1][k2][k3];
1652  ldc = ldc + lda * ldb;
1653  }
1654  }
1655  }
1656  }
1657  working_space[i1][i2][i3 + ssizez] = ldc;
1658  }
1659  }
1660  }
1661 
1662 //calculate matrix b=ht*h
1663  i1min = -(lhx - 1), i1max = lhx - 1;
1664  i2min = -(lhy - 1), i2max = lhy - 1;
1665  i3min = -(lhz - 1), i3max = lhz - 1;
1666  for (i3 = i3min; i3 <= i3max; i3++) {
1667  for (i2 = i2min; i2 <= i2max; i2++) {
1668  for (i1 = i1min; i1 <= i1max; i1++) {
1669  ldc = 0;
1670  j3min = -i3;
1671  if (j3min < 0)
1672  j3min = 0;
1673  j3max = lhz - 1 - i3;
1674  if (j3max > lhz - 1)
1675  j3max = lhz - 1;
1676  for (j3 = j3min; j3 <= j3max; j3++) {
1677  j2min = -i2;
1678  if (j2min < 0)
1679  j2min = 0;
1680  j2max = lhy - 1 - i2;
1681  if (j2max > lhy - 1)
1682  j2max = lhy - 1;
1683  for (j2 = j2min; j2 <= j2max; j2++) {
1684  j1min = -i1;
1685  if (j1min < 0)
1686  j1min = 0;
1687  j1max = lhx - 1 - i1;
1688  if (j1max > lhx - 1)
1689  j1max = lhx - 1;
1690  for (j1 = j1min; j1 <= j1max; j1++) {
1691  lda = working_space[j1][j2][j3];
1692  if (i1 + j1 < ssizex && i2 + j2 < ssizey)
1693  ldb = working_space[i1 + j1][i2 + j2][i3 + j3];
1694  else
1695  ldb = 0;
1696  ldc = ldc + lda * ldb;
1697  }
1698  }
1699  }
1700  working_space[i1 - i1min][i2 - i2min][i3 - i3min + 2 * ssizez ] = ldc;
1701  }
1702  }
1703  }
1704 
1705 //initialization in x1 matrix
1706  for (i3 = 0; i3 < ssizez; i3++) {
1707  for (i2 = 0; i2 < ssizey; i2++) {
1708  for (i1 = 0; i1 < ssizex; i1++) {
1709  working_space[i1][i2][i3 + 3 * ssizez] = 1;
1710  working_space[i1][i2][i3 + 4 * ssizez] = 0;
1711  }
1712  }
1713  }
1714 
1715  //START OF ITERATIONS
1716  for (repet = 0; repet < numberRepetitions; repet++) {
1717  if (repet != 0) {
1718  for (i = 0; i < ssizex; i++) {
1719  for (j = 0; j < ssizey; j++) {
1720  for (k = 0; k < ssizez; k++) {
1721  working_space[i][j][k + 3 * ssizez] = TMath::Power(working_space[i][j][k + 3 * ssizez],boost);
1722  }
1723  }
1724  }
1725  }
1726  for (lindex = 0; lindex < numberIterations; lindex++) {
1727  for (i3 = 0; i3 < ssizez; i3++) {
1728  for (i2 = 0; i2 < ssizey; i2++) {
1729  for (i1 = 0; i1 < ssizex; i1++) {
1730  ldb = 0;
1731  j3min = i3;
1732  if (j3min > lhz - 1)
1733  j3min = lhz - 1;
1734  j3min = -j3min;
1735  j3max = ssizez - i3 - 1;
1736  if (j3max > lhz - 1)
1737  j3max = lhz - 1;
1738  j2min = i2;
1739  if (j2min > lhy - 1)
1740  j2min = lhy - 1;
1741  j2min = -j2min;
1742  j2max = ssizey - i2 - 1;
1743  if (j2max > lhy - 1)
1744  j2max = lhy - 1;
1745  j1min = i1;
1746  if (j1min > lhx - 1)
1747  j1min = lhx - 1;
1748  j1min = -j1min;
1749  j1max = ssizex - i1 - 1;
1750  if (j1max > lhx - 1)
1751  j1max = lhx - 1;
1752  for (j3 = j3min; j3 <= j3max; j3++) {
1753  for (j2 = j2min; j2 <= j2max; j2++) {
1754  for (j1 = j1min; j1 <= j1max; j1++) {
1755  ldc = working_space[j1 - i1min][j2 - i2min][j3 - i3min + 2 * ssizez];
1756  lda = working_space[i1 + j1][i2 + j2][i3 + j3 + 3 * ssizez];
1757  ldb = ldb + lda * ldc;
1758  }
1759  }
1760  }
1761  lda = working_space[i1][i2][i3 + 3 * ssizez];
1762  ldc = working_space[i1][i2][i3 + 1 * ssizez];
1763  if (ldc * lda != 0 && ldb != 0) {
1764  lda = lda * ldc / ldb;
1765  }
1766 
1767  else
1768  lda = 0;
1769  working_space[i1][i2][i3 + 4 * ssizez] = lda;
1770  }
1771  }
1772  }
1773  for (i3 = 0; i3 < ssizez; i3++) {
1774  for (i2 = 0; i2 < ssizey; i2++) {
1775  for (i1 = 0; i1 < ssizex; i1++)
1776  working_space[i1][i2][i3 + 3 * ssizez] = working_space[i1][i2][i3 + 4 * ssizez];
1777  }
1778  }
1779  }
1780  }
1781  for (i = 0; i < ssizex; i++) {
1782  for (j = 0; j < ssizey; j++){
1783  for (k = 0; k < ssizez; k++)
1784  source[(i + positx) % ssizex][(j + posity) % ssizey][(k + positz) % ssizez] = area * working_space[i][j][k + 3 * ssizez];
1785  }
1786  }
1787  delete [] working_space;
1788  return 0;
1789 }
1790 
1791 ////////////////////////////////////////////////////////////////////////////////
1792 /// This function searches for peaks in source spectrum
1793 /// It is based on deconvolution method. First the background is
1794 /// removed (if desired), then Markov spectrum is calculated
1795 /// (if desired), then the response function is generated
1796 /// according to given sigma and deconvolution is carried out.
1797 /// It returns number of found peaks.
1798 ///
1799 /// Function parameters:
1800 /// - source-pointer to the matrix of source spectrum
1801 /// - dest-pointer to the matrix of resulting deconvolved spectrum
1802 /// - ssizex-x length of source spectrum
1803 /// - ssizey-y length of source spectrum
1804 /// - ssizez-z length of source spectrum
1805 /// - sigma-sigma of searched peaks, for details we refer to manual
1806 /// - threshold-threshold value in % for selected peaks, peaks with
1807 /// amplitude less than threshold*highest_peak/100
1808 /// are ignored, see manual
1809 /// - backgroundRemove-logical variable, set if the removal of
1810 /// background before deconvolution is desired
1811 /// - deconIterations-number of iterations in deconvolution operation
1812 /// - markov-logical variable, if it is true, first the source spectrum
1813 /// is replaced by new spectrum calculated using Markov
1814 /// chains method.
1815 /// - averWindow-averaging window of searched peaks, for details
1816 /// we refer to manual (applies only for Markov method)
1817 ///
1818 /// ### Peaks searching
1819 ///
1820 /// Goal: to identify automatically the peaks in spectrum with the presence of
1821 /// the continuous background, one- and two-fold coincidences (ridges) and statistical
1822 /// fluctuations - noise.
1823 ///
1824 /// The common problems connected
1825 /// with correct peak identification in three-dimensional coincidence spectra are
1826 ///
1827 /// - non-sensitivity to noise, i.e.,
1828 /// only statistically relevant peaks should be identified
1829 /// - non-sensitivity of the
1830 /// algorithm to continuous background
1831 /// - non-sensitivity to one-fold coincidences
1832 /// (coincidences peak - peak - background in all dimensions) and their
1833 /// crossings
1834 /// - non-sensitivity to two-fold
1835 /// coincidences (coincidences peak - background - background in all
1836 /// dimensions) and their crossings
1837 /// - ability to identify peaks close
1838 /// to the edges of the spectrum region
1839 /// - resolution, decomposition of
1840 /// doublets and multiplets. The algorithm should be able to recognise close
1841 /// positioned peaks.
1842 ///
1843 /// #### References:
1844 ///
1845 /// [1] M.A. Mariscotti: A method for
1846 /// identification of peaks in the presence of background and its application to
1847 /// spectrum analysis. NIM 50 (1967), 309-320.
1848 ///
1849 /// [2] M.Morhac, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:Identification
1850 /// of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
1851 /// 108-125.
1852 ///
1853 /// [3] Z.K. Silagadze, A new algorithm for automatic photo-peak searches. NIM A 376 (1996), 451.
1854 ///
1855 /// ### Example of peak searching method
1856 ///
1857 /// SearchHighRes function provides users with the possibility
1858 /// to vary the input parameters and with the access to the output deconvolved data
1859 /// in the destination spectrum. Based on the output data one can tune the
1860 /// parameters.
1861 ///
1862 /// #### Example 1 - script Search3.c:
1863 ///
1864 /// \image html spectrum3_searching_image001.jpg Fig. 1 Three-dimensional spectrum with 5 peaks (sigma=2, threshold=5%, 3 iterations steps in the deconvolution)
1865 /// \image html spectrum3_searching_image003.jpg Fig. 2 Spectrum from Fig. 1 after background elimination and deconvolution
1866 ///
1867 /// #### Script:
1868 ///
1869 /// Example to illustrate high resolution peak searching function (class TSpectrum3).
1870 /// To execute this example, do:
1871 ///
1872 /// `root > .x Search3.C`
1873 ///
1874 /// ~~~ {.cpp}
1875 /// void Search3() {
1876 /// Int_t i, j, k, nfound;
1877 /// Int_t nbinsx = 32;
1878 /// Int_t nbinsy = 32;
1879 /// Int_t nbinsz = 32;
1880 /// Int_t xmin = 0;
1881 /// Int_t xmax = nbinsx;
1882 /// Int_t ymin = 0;
1883 /// Int_t ymax = nbinsy;
1884 /// Int_t zmin = 0;
1885 /// Int_t zmax = nbinsz;
1886 /// Double_t*** source = new Double_t**[nbinsx];
1887 /// Double_t*** dest = new Double_t**[nbinsx];
1888 /// for(i=0;i<nbinsx;i++){
1889 /// source[i]=new Double_t*[nbinsy];
1890 /// for(j=0;j<nbinsy;j++)
1891 /// source[i][j]=new Double_t[nbinsz];
1892 /// }
1893 /// for(i=0;i<nbinsx;i++){
1894 /// dest[i]=new Double_t*[nbinsy];
1895 /// for(j=0;j<nbinsy;j++)
1896 /// dest[i][j]=new Double_t [nbinsz];
1897 /// }
1898 /// TH3F *search = new TH3F("Search","Peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1899 /// TFile *f = new TFile("TSpectrum3.root");
1900 /// search=(TH3F*)f->Get("search2;1");
1901 /// TCanvas *Search = new TCanvas("Search","Peak searching",10,10,1000,700);
1902 /// TSpectrum3 *s = new TSpectrum3();
1903 /// for (i = 0; i < nbinsx; i++){
1904 /// for (j = 0; j < nbinsy; j++){
1905 /// for (k = 0; k < nbinsz; k++){
1906 /// source[i][j][k] = search->GetBinContent(i + 1,j + 1,k + 1);
1907 /// }
1908 /// }
1909 /// }
1910 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, nbinsz, 2, 5, kTRUE, 3, kFALSE, 3);
1911 /// printf("Found %d candidate peaks\n",nfound);
1912 /// for (i = 0; i < nbinsx; i++){
1913 /// for (j = 0; j < nbinsy; j++){
1914 /// for (k = 0; k < nbinsz; k++){
1915 /// search->SetBinContent(i + 1,j + 1,k + 1, dest[i][j][k]);
1916 /// }
1917 /// }
1918 /// }
1919 /// Double_t *PosX = new Double_t[nfound];
1920 /// Double_t *PosY = new Double_t[nfound];
1921 /// Double_t *PosZ = new Double_t[nfound];
1922 /// PosX = s->GetPositionX();
1923 /// PosY = s->GetPositionY();
1924 /// PosZ = s->GetPositionZ();
1925 /// for(i=0;i<nfound;i++)
1926 /// printf("posx= %d, posy= %d, posz=%d\n",(Int_t)(PosX[i]+0.5), (Int_t)(PosY[i]+0.5),(Int_t)(PosZ[i]+0.5));
1927 /// search->Draw("");
1928 /// }
1929 /// ~~~
1930 
1931 Int_t TSpectrum3::SearchHighRes(const Double_t***source,Double_t***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez,
1932  Double_t sigma, Double_t threshold,
1933  Bool_t backgroundRemove,Int_t deconIterations,
1934  Bool_t markov, Int_t averWindow)
1935 
1936 {
1937  Int_t number_of_iterations = (Int_t)(4 * sigma + 0.5);
1938  Int_t k,lindex;
1939  Double_t lda,ldb,ldc,area,maximum;
1940  Int_t xmin,xmax,l,peak_index = 0,sizex_ext=ssizex + 4 * number_of_iterations,sizey_ext = ssizey + 4 * number_of_iterations,sizez_ext = ssizez + 4 * number_of_iterations,shift = 2 * number_of_iterations;
1941  Int_t ymin,ymax,zmin,zmax,i,j;
1942  Double_t a,b,maxch,plocha = 0,plocha_markov = 0;
1943  Double_t nom,nip,nim,sp,sm,spx,spy,smx,smy,spz,smz;
1944  Double_t p1,p2,p3,p4,p5,p6,p7,p8,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,r1,r2,r3,r4,r5,r6;
1945  Int_t x,y,z;
1946  Double_t pocet_sigma = 5;
1947  Int_t lhx,lhy,lhz,i1,i2,i3,j1,j2,j3,k1,k2,k3,i1min,i1max,i2min,i2max,i3min,i3max,j1min,j1max,j2min,j2max,j3min,j3max,positx,posity,positz;
1948  if(sigma < 1){
1949  Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
1950  return 0;
1951  }
1952 
1953  if(threshold<=0||threshold>=100){
1954  Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
1955  return 0;
1956  }
1957 
1958  j = (Int_t)(pocet_sigma*sigma+0.5);
1959  if (j >= PEAK_WINDOW / 2) {
1960  Error("SearchHighRes", "Too large sigma");
1961  return 0;
1962  }
1963 
1964  if (markov == true) {
1965  if (averWindow <= 0) {
1966  Error("SearchHighRes", "Averaging window must be positive");
1967  return 0;
1968  }
1969  }
1970 
1971  if(backgroundRemove == true){
1972  if(sizex_ext < 2 * number_of_iterations + 1 || sizey_ext < 2 * number_of_iterations + 1 || sizez_ext < 2 * number_of_iterations + 1){
1973  Error("SearchHighRes", "Too large clipping window");
1974  return 0;
1975  }
1976  }
1977 
1978  i = (Int_t)(4 * sigma + 0.5);
1979  i = 4 * i;
1980  Double_t ***working_space = new Double_t** [ssizex + i];
1981  for(j = 0;j < ssizex + i; j++){
1982  working_space[j] = new Double_t* [ssizey + i];
1983  for(k = 0;k < ssizey + i; k++)
1984  working_space[j][k] = new Double_t [5 * (ssizez + i)];
1985  }
1986  for(k = 0;k < sizez_ext; k++){
1987  for(j = 0;j < sizey_ext; j++){
1988  for(i = 0;i < sizex_ext; i++){
1989  if(i < shift){
1990  if(j < shift){
1991  if(k < shift)
1992  working_space[i][j][k + sizez_ext] = source[0][0][0];
1993 
1994  else if(k >= ssizez + shift)
1995  working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1];
1996 
1997  else
1998  working_space[i][j][k + sizez_ext] = source[0][0][k - shift];
1999  }
2000 
2001  else if(j >= ssizey + shift){
2002  if(k < shift)
2003  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0];
2004 
2005  else if(k >= ssizez + shift)
2006  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1];
2007 
2008  else
2009  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift];
2010  }
2011 
2012  else{
2013  if(k < shift)
2014  working_space[i][j][k + sizez_ext] = source[0][j - shift][0];
2015 
2016  else if(k >= ssizez + shift)
2017  working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1];
2018 
2019  else
2020  working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift];
2021  }
2022  }
2023 
2024  else if(i >= ssizex + shift){
2025  if(j < shift){
2026  if(k < shift)
2027  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0];
2028 
2029  else if(k >= ssizez + shift)
2030  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1];
2031 
2032  else
2033  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift];
2034  }
2035 
2036  else if(j >= ssizey + shift){
2037  if(k < shift)
2038  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0];
2039 
2040  else if(k >= ssizez + shift)
2041  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1];
2042 
2043  else
2044  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift];
2045  }
2046 
2047  else{
2048  if(k < shift)
2049  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0];
2050 
2051  else if(k >= ssizez + shift)
2052  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1];
2053 
2054  else
2055  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift];
2056  }
2057  }
2058 
2059  else{
2060  if(j < shift){
2061  if(k < shift)
2062  working_space[i][j][k + sizez_ext] = source[i - shift][0][0];
2063 
2064  else if(k >= ssizez + shift)
2065  working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1];
2066 
2067  else
2068  working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift];
2069  }
2070 
2071  else if(j >= ssizey + shift){
2072  if(k < shift)
2073  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0];
2074 
2075  else if(k >= ssizez + shift)
2076  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1];
2077 
2078  else
2079  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift];
2080  }
2081 
2082  else{
2083  if(k < shift)
2084  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0];
2085 
2086  else if(k >= ssizez + shift)
2087  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1];
2088 
2089  else
2090  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift];
2091  }
2092  }
2093  }
2094  }
2095  }
2096  if(backgroundRemove == true){
2097  for(i = 1;i <= number_of_iterations; i++){
2098  for(z = i;z < sizez_ext - i; z++){
2099  for(y = i;y < sizey_ext - i; y++){
2100  for(x = i;x < sizex_ext - i; x++){
2101  a = working_space[x][y][z + sizez_ext];
2102  p1 = working_space[x + i][y + i][z - i + sizez_ext];
2103  p2 = working_space[x - i][y + i][z - i + sizez_ext];
2104  p3 = working_space[x + i][y - i][z - i + sizez_ext];
2105  p4 = working_space[x - i][y - i][z - i + sizez_ext];
2106  p5 = working_space[x + i][y + i][z + i + sizez_ext];
2107  p6 = working_space[x - i][y + i][z + i + sizez_ext];
2108  p7 = working_space[x + i][y - i][z + i + sizez_ext];
2109  p8 = working_space[x - i][y - i][z + i + sizez_ext];
2110  s1 = working_space[x + i][y ][z - i + sizez_ext];
2111  s2 = working_space[x ][y + i][z - i + sizez_ext];
2112  s3 = working_space[x - i][y ][z - i + sizez_ext];
2113  s4 = working_space[x ][y - i][z - i + sizez_ext];
2114  s5 = working_space[x + i][y ][z + i + sizez_ext];
2115  s6 = working_space[x ][y + i][z + i + sizez_ext];
2116  s7 = working_space[x - i][y ][z + i + sizez_ext];
2117  s8 = working_space[x ][y - i][z + i + sizez_ext];
2118  s9 = working_space[x - i][y + i][z + sizez_ext];
2119  s10 = working_space[x - i][y - i][z +sizez_ext];
2120  s11 = working_space[x + i][y + i][z +sizez_ext];
2121  s12 = working_space[x + i][y - i][z +sizez_ext];
2122  r1 = working_space[x ][y ][z - i + sizez_ext];
2123  r2 = working_space[x ][y ][z + i + sizez_ext];
2124  r3 = working_space[x - i][y ][z + sizez_ext];
2125  r4 = working_space[x + i][y ][z + sizez_ext];
2126  r5 = working_space[x ][y + i][z + sizez_ext];
2127  r6 = working_space[x ][y - i][z + sizez_ext];
2128  b = (p1 + p3) / 2.0;
2129  if(b > s1)
2130  s1 = b;
2131 
2132  b = (p1 + p2) / 2.0;
2133  if(b > s2)
2134  s2 = b;
2135 
2136  b = (p2 + p4) / 2.0;
2137  if(b > s3)
2138  s3 = b;
2139 
2140  b = (p3 + p4) / 2.0;
2141  if(b > s4)
2142  s4 = b;
2143 
2144  b = (p5 + p7) / 2.0;
2145  if(b > s5)
2146  s5 = b;
2147 
2148  b = (p5 + p6) / 2.0;
2149  if(b > s6)
2150  s6 = b;
2151 
2152  b = (p6 + p8) / 2.0;
2153  if(b > s7)
2154  s7 = b;
2155 
2156  b = (p7 + p8) / 2.0;
2157  if(b > s8)
2158  s8 = b;
2159 
2160  b = (p2 + p6) / 2.0;
2161  if(b > s9)
2162  s9 = b;
2163 
2164  b = (p4 + p8) / 2.0;
2165  if(b > s10)
2166  s10 = b;
2167 
2168  b = (p1 + p5) / 2.0;
2169  if(b > s11)
2170  s11 = b;
2171 
2172  b = (p3 + p7) / 2.0;
2173  if(b > s12)
2174  s12 = b;
2175 
2176  s1 = s1 - (p1 + p3) / 2.0;
2177  s2 = s2 - (p1 + p2) / 2.0;
2178  s3 = s3 - (p2 + p4) / 2.0;
2179  s4 = s4 - (p3 + p4) / 2.0;
2180  s5 = s5 - (p5 + p7) / 2.0;
2181  s6 = s6 - (p5 + p6) / 2.0;
2182  s7 = s7 - (p6 + p8) / 2.0;
2183  s8 = s8 - (p7 + p8) / 2.0;
2184  s9 = s9 - (p2 + p6) / 2.0;
2185  s10 = s10 - (p4 + p8) / 2.0;
2186  s11 = s11 - (p1 + p5) / 2.0;
2187  s12 = s12 - (p3 + p7) / 2.0;
2188  b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
2189  if(b > r1)
2190  r1 = b;
2191 
2192  b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
2193  if(b > r2)
2194  r2 = b;
2195 
2196  b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
2197  if(b > r3)
2198  r3 = b;
2199 
2200  b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
2201  if(b > r4)
2202  r4 = b;
2203 
2204  b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
2205  if(b > r5)
2206  r5 = b;
2207 
2208  b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
2209  if(b > r6)
2210  r6 = b;
2211 
2212  r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
2213  r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
2214  r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
2215  r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
2216  r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
2217  r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
2218  b = (r1 + r2) / 2.0 + (r3 + r4) / 2.0 + (r5 + r6) / 2.0 + (s1 + s3 + s5 + s7) / 4.0 + (s2 + s4 + s6 + s8) / 4.0 + (s9 + s10 + s11 + s12) / 4.0 + (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.0;
2219  if(b < a)
2220  a = b;
2221 
2222  working_space[x][y][z] = a;
2223  }
2224  }
2225  }
2226  for(z = i;z < sizez_ext - i; z++){
2227  for(y = i;y < sizey_ext - i; y++){
2228  for(x = i;x < sizex_ext - i; x++){
2229  working_space[x][y][z + sizez_ext] = working_space[x][y][z];
2230  }
2231  }
2232  }
2233  }
2234  for(k = 0;k < sizez_ext; k++){
2235  for(j = 0;j < sizey_ext; j++){
2236  for(i = 0;i < sizex_ext; i++){
2237  if(i < shift){
2238  if(j < shift){
2239  if(k < shift)
2240  working_space[i][j][k + sizez_ext] = source[0][0][0] - working_space[i][j][k + sizez_ext];
2241 
2242  else if(k >= ssizez + shift)
2243  working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
2244 
2245  else
2246  working_space[i][j][k + sizez_ext] = source[0][0][k - shift] - working_space[i][j][k + sizez_ext];
2247  }
2248 
2249  else if(j >= ssizey + shift){
2250  if(k < shift)
2251  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
2252 
2253  else if(k >= ssizez + shift)
2254  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
2255 
2256  else
2257  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
2258  }
2259 
2260  else{
2261  if(k < shift)
2262  working_space[i][j][k + sizez_ext] = source[0][j - shift][0] - working_space[i][j][k + sizez_ext];
2263 
2264  else if(k >= ssizez + shift)
2265  working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
2266 
2267  else
2268  working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
2269  }
2270  }
2271 
2272  else if(i >= ssizex + shift){
2273  if(j < shift){
2274  if(k < shift)
2275  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0] - working_space[i][j][k + sizez_ext];
2276 
2277  else if(k >= ssizez + shift)
2278  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
2279 
2280  else
2281  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift] - working_space[i][j][k + sizez_ext];
2282  }
2283 
2284  else if(j >= ssizey + shift){
2285  if(k < shift)
2286  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
2287 
2288  else if(k >= ssizez + shift)
2289  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
2290 
2291  else
2292  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
2293  }
2294 
2295  else{
2296  if(k < shift)
2297  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0] - working_space[i][j][k + sizez_ext];
2298 
2299  else if(k >= ssizez + shift)
2300  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
2301 
2302  else
2303  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
2304  }
2305  }
2306 
2307  else{
2308  if(j < shift){
2309  if(k < shift)
2310  working_space[i][j][k + sizez_ext] = source[i - shift][0][0] - working_space[i][j][k + sizez_ext];
2311 
2312  else if(k >= ssizez + shift)
2313  working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
2314 
2315  else
2316  working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift] - working_space[i][j][k + sizez_ext];
2317  }
2318 
2319  else if(j >= ssizey + shift){
2320  if(k < shift)
2321  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
2322 
2323  else if(k >= ssizez + shift)
2324  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
2325 
2326  else
2327  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
2328  }
2329 
2330  else{
2331  if(k < shift)
2332  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0] - working_space[i][j][k + sizez_ext];
2333 
2334  else if(k >= ssizez + shift)
2335  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
2336 
2337  else
2338  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
2339  }
2340  }
2341  }
2342  }
2343  }
2344  }
2345 
2346  if(markov == true){
2347  for(i = 0;i < sizex_ext; i++){
2348  for(j = 0;j < sizey_ext; j++){
2349  for(k = 0;k < sizez_ext; k++){
2350  working_space[i][j][k + 2 * sizez_ext] = working_space[i][j][sizez_ext + k];
2351  plocha_markov = plocha_markov + working_space[i][j][sizez_ext + k];
2352  }
2353  }
2354  }
2355  xmin = 0;
2356  xmax = sizex_ext - 1;
2357  ymin = 0;
2358  ymax = sizey_ext - 1;
2359  zmin = 0;
2360  zmax = sizez_ext - 1;
2361  for(i = 0,maxch = 0;i < sizex_ext; i++){
2362  for(j = 0;j < sizey_ext;j++){
2363  for(k = 0;k < sizez_ext;k++){
2364  working_space[i][j][k] = 0;
2365  if(maxch < working_space[i][j][k + 2 * sizez_ext])
2366  maxch = working_space[i][j][k + 2 * sizez_ext];
2367 
2368  plocha += working_space[i][j][k + 2 * sizez_ext];
2369  }
2370  }
2371  }
2372  if(maxch == 0) {
2373  delete [] working_space;
2374  return 0;
2375  }
2376  nom = 0;
2377  working_space[xmin][ymin][zmin] = 1;
2378  for(i = xmin;i < xmax; i++){
2379  nip = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
2380  nim = working_space[i + 1][ymin][zmin + 2 * sizez_ext] / maxch;
2381  sp = 0,sm = 0;
2382  for(l = 1;l <= averWindow; l++){
2383  if((i + l) > xmax)
2384  a = working_space[xmax][ymin][zmin + 2 * sizez_ext] / maxch;
2385 
2386  else
2387  a = working_space[i + l][ymin][zmin + 2 * sizez_ext] / maxch;
2388 
2389  b = a - nip;
2390  if(a + nip <= 0)
2391  a = 1;
2392 
2393  else
2394  a = TMath::Sqrt(a + nip);
2395 
2396  b = b / a;
2397  b = TMath::Exp(b);
2398  sp = sp + b;
2399  if(i - l + 1 < xmin)
2400  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
2401 
2402  else
2403  a = working_space[i - l + 1][ymin][zmin + 2 * sizez_ext] / maxch;
2404 
2405  b = a - nim;
2406  if(a + nim <= 0)
2407  a = 1;
2408 
2409  else
2410  a = TMath::Sqrt(a + nim);
2411 
2412  b = b / a;
2413  b = TMath::Exp(b);
2414  sm = sm + b;
2415  }
2416  a = sp / sm;
2417  a = working_space[i + 1][ymin][zmin] = a * working_space[i][ymin][zmin];
2418  nom = nom + a;
2419  }
2420  for(i = ymin;i < ymax; i++){
2421  nip = working_space[xmin][i][zmin + 2 * sizez_ext] / maxch;
2422  nim = working_space[xmin][i + 1][zmin + 2 * sizez_ext] / maxch;
2423  sp = 0,sm = 0;
2424  for(l = 1;l <= averWindow; l++){
2425  if((i + l) > ymax)
2426  a = working_space[xmin][ymax][zmin + 2 * sizez_ext] / maxch;
2427 
2428  else
2429  a = working_space[xmin][i + l][zmin + 2 * sizez_ext] / maxch;
2430 
2431  b = a - nip;
2432  if(a + nip <= 0)
2433  a = 1;
2434 
2435  else
2436  a = TMath::Sqrt(a + nip);
2437 
2438  b = b / a;
2439  b = TMath::Exp(b);
2440  sp = sp + b;
2441  if(i - l + 1 < ymin)
2442  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
2443 
2444  else
2445  a = working_space[xmin][i - l + 1][zmin + 2 * sizez_ext] / maxch;
2446 
2447  b = a - nim;
2448  if(a + nim <= 0)
2449  a = 1;
2450 
2451  else
2452  a = TMath::Sqrt(a + nim);
2453 
2454  b = b / a;
2455  b = TMath::Exp(b);
2456  sm = sm + b;
2457  }
2458  a = sp / sm;
2459  a = working_space[xmin][i + 1][zmin] = a * working_space[xmin][i][zmin];
2460  nom = nom + a;
2461  }
2462  for(i = zmin;i < zmax;i++){
2463  nip = working_space[xmin][ymin][i + 2 * sizez_ext] / maxch;
2464  nim = working_space[xmin][ymin][i + 1 + 2 * sizez_ext] / maxch;
2465  sp = 0,sm = 0;
2466  for(l = 1;l <= averWindow;l++){
2467  if((i + l) > zmax)
2468  a = working_space[xmin][ymin][zmax + 2 * sizez_ext] / maxch;
2469 
2470  else
2471  a = working_space[xmin][ymin][i + l + 2 * sizez_ext] / maxch;
2472 
2473  b = a - nip;
2474  if(a + nip <= 0)
2475  a = 1;
2476 
2477  else
2478  a = TMath::Sqrt(a + nip);
2479 
2480  b = b / a;
2481  b = TMath::Exp(b);
2482  sp = sp + b;
2483  if(i - l + 1 < zmin)
2484  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
2485 
2486  else
2487  a = working_space[xmin][ymin][i - l + 1 + 2 * sizez_ext] / maxch;
2488 
2489  b = a - nim;
2490  if(a + nim <= 0)
2491  a = 1;
2492 
2493  else
2494  a = TMath::Sqrt(a + nim);
2495 
2496  b = b / a;
2497  b = TMath::Exp(b);
2498  sm = sm + b;
2499  }
2500  a = sp / sm;
2501  a = working_space[xmin][ymin][i + 1] = a * working_space[xmin][ymin][i];
2502  nom = nom + a;
2503  }
2504  for(i = xmin;i < xmax; i++){
2505  for(j = ymin;j < ymax; j++){
2506  nip = working_space[i][j + 1][zmin + 2 * sizez_ext] / maxch;
2507  nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
2508  spx = 0,smx = 0;
2509  for(l = 1;l <= averWindow; l++){
2510  if(i + l > xmax)
2511  a = working_space[xmax][j][zmin + 2 * sizez_ext] / maxch;
2512 
2513  else
2514  a = working_space[i + l][j][zmin + 2 * sizez_ext] / maxch;
2515 
2516  b = a - nip;
2517  if(a + nip <= 0)
2518  a = 1;
2519 
2520  else
2521  a = TMath::Sqrt(a + nip);
2522 
2523  b = b / a;
2524  b = TMath::Exp(b);
2525  spx = spx + b;
2526  if(i - l + 1 < xmin)
2527  a = working_space[xmin][j][zmin + 2 * sizez_ext] / maxch;
2528 
2529  else
2530  a = working_space[i - l + 1][j][zmin + 2 * sizez_ext] / maxch;
2531 
2532  b = a - nim;
2533  if(a + nim <= 0)
2534  a = 1;
2535 
2536  else
2537  a = TMath::Sqrt(a + nim);
2538 
2539  b = b / a;
2540  b = TMath::Exp(b);
2541  smx = smx + b;
2542  }
2543  spy = 0,smy = 0;
2544  nip = working_space[i + 1][j][zmin + 2 * sizez_ext] / maxch;
2545  nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
2546  for(l = 1;l <= averWindow; l++){
2547  if(j + l > ymax)
2548  a = working_space[i][ymax][zmin + 2 * sizez_ext] / maxch;
2549 
2550  else
2551  a = working_space[i][j + l][zmin + 2 * sizez_ext] / maxch;
2552 
2553  b = a - nip;
2554  if(a + nip <= 0)
2555  a = 1;
2556 
2557  else
2558  a = TMath::Sqrt(a + nip);
2559 
2560  b = b / a;
2561  b = TMath::Exp(b);
2562  spy = spy + b;
2563  if(j - l + 1 < ymin)
2564  a = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
2565 
2566  else
2567  a = working_space[i][j - l + 1][zmin + 2 * sizez_ext] / maxch;
2568 
2569  b = a - nim;
2570  if(a + nim <= 0)
2571  a = 1;
2572 
2573  else
2574  a = TMath::Sqrt(a + nim);
2575 
2576  b = b / a;
2577  b = TMath::Exp(b);
2578  smy = smy + b;
2579  }
2580  a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
2581  working_space[i + 1][j + 1][zmin] = a;
2582  nom = nom + a;
2583  }
2584  }
2585  for(i = xmin;i < xmax;i++){
2586  for(j = zmin;j < zmax;j++){
2587  nip = working_space[i][ymin][j + 1 + 2 * sizez_ext] / maxch;
2588  nim = working_space[i + 1][ymin][j + 1 + 2 * sizez_ext] / maxch;
2589  spx = 0,smx = 0;
2590  for(l = 1;l <= averWindow; l++){
2591  if(i + l > xmax)
2592  a = working_space[xmax][ymin][j + 2 * sizez_ext] / maxch;
2593 
2594  else
2595  a = working_space[i + l][ymin][j + 2 * sizez_ext] / maxch;
2596 
2597  b = a - nip;
2598  if(a + nip <= 0)
2599  a = 1;
2600 
2601  else
2602  a = TMath::Sqrt(a + nip);
2603 
2604  b = b / a;
2605  b = TMath::Exp(b);
2606  spx = spx + b;
2607  if(i - l + 1 < xmin)
2608  a = working_space[xmin][ymin][j + 2 * sizez_ext] / maxch;
2609 
2610  else
2611  a = working_space[i - l + 1][ymin][j + 2 * sizez_ext] / maxch;
2612 
2613  b = a - nim;
2614  if(a + nim <= 0)
2615  a = 1;
2616 
2617  else
2618  a = TMath::Sqrt(a + nim);
2619 
2620  b = b / a;
2621  b = TMath::Exp(b);
2622  smx = smx + b;
2623  }
2624  spz = 0,smz = 0;
2625  nip = working_space[i + 1][ymin][j + 2 * sizez_ext] / maxch;
2626  nim = working_space[i + 1][ymin][j + 1 + 2 * sizez_ext] / maxch;
2627  for(l = 1;l <= averWindow; l++){
2628  if(j + l > zmax)
2629  a = working_space[i][ymin][zmax + 2 * sizez_ext] / maxch;
2630 
2631  else
2632  a = working_space[i][ymin][j + l + 2 * sizez_ext] / maxch;
2633 
2634  b = a - nip;
2635  if(a + nip <= 0)
2636  a = 1;
2637 
2638  else
2639  a = TMath::Sqrt(a + nip);
2640 
2641  b = b / a;
2642  b = TMath::Exp(b);
2643  spz = spz + b;
2644  if(j - l + 1 < zmin)
2645  a = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
2646 
2647  else
2648  a = working_space[i][ymin][j - l + 1 + 2 * sizez_ext] / maxch;
2649 
2650  b = a - nim;
2651  if(a + nim <= 0)
2652  a = 1;
2653 
2654  else
2655  a = TMath::Sqrt(a + nim);
2656 
2657  b = b / a;
2658  b = TMath::Exp(b);
2659  smz = smz + b;
2660  }
2661  a = (spx * working_space[i][ymin][j + 1] + spz * working_space[i + 1][ymin][j]) / (smx + smz);
2662  working_space[i + 1][ymin][j + 1] = a;
2663  nom = nom + a;
2664  }
2665  }
2666  for(i = ymin;i < ymax;i++){
2667  for(j = zmin;j < zmax;j++){
2668  nip = working_space[xmin][i][j + 1 + 2 * sizez_ext] / maxch;
2669  nim = working_space[xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
2670  spy = 0,smy = 0;
2671  for(l = 1;l <= averWindow; l++){
2672  if(i + l > ymax)
2673  a = working_space[xmin][ymax][j + 2 * sizez_ext] / maxch;
2674 
2675  else
2676  a = working_space[xmin][i + l][j + 2 * sizez_ext] / maxch;
2677 
2678  b = a - nip;
2679  if(a + nip <= 0)
2680  a = 1;
2681 
2682  else
2683  a = TMath::Sqrt(a + nip);
2684 
2685  b = b / a;
2686  b = TMath::Exp(b);
2687  spy = spy + b;
2688  if(i - l + 1 < ymin)
2689  a = working_space[xmin][ymin][j + 2 * sizez_ext] / maxch;
2690 
2691  else
2692  a = working_space[xmin][i - l + 1][j + 2 * sizez_ext] / maxch;
2693 
2694  b = a - nim;
2695  if(a + nim <= 0)
2696  a = 1;
2697 
2698  else
2699  a = TMath::Sqrt(a + nim);
2700 
2701  b = b / a;
2702  b = TMath::Exp(b);
2703  smy = smy + b;
2704  }
2705  spz = 0,smz = 0;
2706  nip = working_space[xmin][i + 1][j + 2 * sizez_ext] / maxch;
2707  nim = working_space[xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
2708  for(l = 1;l <= averWindow; l++){
2709  if(j + l > zmax)
2710  a = working_space[xmin][i][zmax + 2 * sizez_ext] / maxch;
2711 
2712  else
2713  a = working_space[xmin][i][j + l + 2 * sizez_ext] / maxch;
2714 
2715  b = a - nip;
2716  if(a + nip <= 0)
2717  a = 1;
2718 
2719  else
2720  a = TMath::Sqrt(a + nip);
2721 
2722  b = b / a;
2723  b = TMath::Exp(b);
2724  spz = spz + b;
2725  if(j - l + 1 < zmin)
2726  a = working_space[xmin][i][zmin + 2 * sizez_ext] / maxch;
2727 
2728  else
2729  a = working_space[xmin][i][j - l + 1 + 2 * sizez_ext] / maxch;
2730 
2731  b = a - nim;
2732  if(a + nim <= 0)
2733  a = 1;
2734 
2735  else
2736  a = TMath::Sqrt(a + nim);
2737 
2738  b = b / a;
2739  b = TMath::Exp(b);
2740  smz = smz + b;
2741  }
2742  a = (spy * working_space[xmin][i][j + 1] + spz * working_space[xmin][i + 1][j]) / (smy + smz);
2743  working_space[xmin][i + 1][j + 1] = a;
2744  nom = nom + a;
2745  }
2746  }
2747  for(i = xmin;i < xmax; i++){
2748  for(j = ymin;j < ymax; j++){
2749  for(k = zmin;k < zmax; k++){
2750  nip = working_space[i][j + 1][k + 1 + 2 * sizez_ext] / maxch;
2751  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
2752  spx = 0,smx = 0;
2753  for(l = 1;l <= averWindow; l++){
2754  if(i + l > xmax)
2755  a = working_space[xmax][j][k + 2 * sizez_ext] / maxch;
2756 
2757  else
2758  a = working_space[i + l][j][k + 2 * sizez_ext] / maxch;
2759 
2760  b = a - nip;
2761  if(a + nip <= 0)
2762  a = 1;
2763 
2764  else
2765  a = TMath::Sqrt(a + nip);
2766 
2767  b = b / a;
2768  b = TMath::Exp(b);
2769  spx = spx + b;
2770  if(i - l + 1 < xmin)
2771  a = working_space[xmin][j][k + 2 * sizez_ext] / maxch;
2772 
2773  else
2774  a = working_space[i - l + 1][j][k + 2 * sizez_ext] / maxch;
2775 
2776  b = a - nim;
2777  if(a + nim <= 0)
2778  a = 1;
2779 
2780  else
2781  a = TMath::Sqrt(a + nim);
2782 
2783  b = b / a;
2784  b = TMath::Exp(b);
2785  smx = smx + b;
2786  }
2787  spy = 0,smy = 0;
2788  nip = working_space[i + 1][j][k + 1 + 2 * sizez_ext] / maxch;
2789  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
2790  for(l = 1;l <= averWindow; l++){
2791  if(j + l > ymax)
2792  a = working_space[i][ymax][k + 2 * sizez_ext] / maxch;
2793 
2794  else
2795  a = working_space[i][j + l][k + 2 * sizez_ext] / maxch;
2796 
2797  b = a - nip;
2798  if(a + nip <= 0)
2799  a = 1;
2800 
2801  else
2802  a = TMath::Sqrt(a + nip);
2803 
2804  b = b / a;
2805  b = TMath::Exp(b);
2806  spy = spy + b;
2807  if(j - l + 1 < ymin)
2808  a = working_space[i][ymin][k + 2 * sizez_ext] / maxch;
2809 
2810  else
2811  a = working_space[i][j - l + 1][k + 2 * sizez_ext] / maxch;
2812 
2813  b = a - nim;
2814  if(a + nim <= 0)
2815  a = 1;
2816 
2817  else
2818  a = TMath::Sqrt(a + nim);
2819 
2820  b = b / a;
2821  b = TMath::Exp(b);
2822  smy = smy + b;
2823  }
2824  spz = 0,smz = 0;
2825  nip = working_space[i + 1][j + 1][k + 2 * sizez_ext] / maxch;
2826  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
2827  for(l = 1;l <= averWindow; l++ ){
2828  if(j + l > zmax)
2829  a = working_space[i][j][zmax + 2 * sizez_ext] / maxch;
2830 
2831  else
2832  a = working_space[i][j][k + l + 2 * sizez_ext] / maxch;
2833 
2834  b = a - nip;
2835  if(a + nip <= 0)
2836  a = 1;
2837 
2838  else
2839  a = TMath::Sqrt(a + nip);
2840 
2841  b = b / a;
2842  b = TMath::Exp(b);
2843  spz = spz + b;
2844  if(j - l + 1 < ymin)
2845  a = working_space[i][j][zmin + 2 * sizez_ext] / maxch;
2846 
2847  else
2848  a = working_space[i][j][k - l + 1 + 2 * sizez_ext] / maxch;
2849 
2850  b = a - nim;
2851  if(a + nim <= 0)
2852  a = 1;
2853 
2854  else
2855  a = TMath::Sqrt(a + nim);
2856 
2857  b = b / a;
2858  b = TMath::Exp(b);
2859  smz = smz + b;
2860  }
2861  a = (spx * working_space[i][j + 1][k + 1] + spy * working_space[i + 1][j][k + 1] + spz * working_space[i + 1][j + 1][k]) / (smx + smy + smz);
2862  working_space[i + 1][j + 1][k + 1] = a;
2863  nom = nom + a;
2864  }
2865  }
2866  }
2867  a = 0;
2868  for(i = xmin;i <= xmax; i++){
2869  for(j = ymin;j <= ymax; j++){
2870  for(k = zmin;k <= zmax; k++){
2871  working_space[i][j][k] = working_space[i][j][k] / nom;
2872  a+=working_space[i][j][k];
2873  }
2874  }
2875  }
2876  for(i = 0;i < sizex_ext; i++){
2877  for(j = 0;j < sizey_ext; j++){
2878  for(k = 0;k < sizez_ext; k++){
2879  working_space[i][j][k + sizez_ext] = working_space[i][j][k] * plocha_markov / a;
2880  }
2881  }
2882  }
2883  }
2884  //deconvolution starts
2885  area = 0;
2886  lhx = -1,lhy = -1,lhz = -1;
2887  positx = 0,posity = 0,positz = 0;
2888  maximum = 0;
2889  //generate response cube
2890  for(i = 0;i < sizex_ext; i++){
2891  for(j = 0;j < sizey_ext; j++){
2892  for(k = 0;k < sizez_ext; k++){
2893  lda = (Double_t)i - 3 * sigma;
2894  ldb = (Double_t)j - 3 * sigma;
2895  ldc = (Double_t)k - 3 * sigma;
2896  lda = (lda * lda + ldb * ldb + ldc * ldc) / (2 * sigma * sigma);
2897  l = (Int_t)(1000 * exp(-lda));
2898  lda = l;
2899  if(lda!=0){
2900  if((i + 1) > lhx)
2901  lhx = i + 1;
2902 
2903  if((j + 1) > lhy)
2904  lhy = j + 1;
2905 
2906  if((k + 1) > lhz)
2907  lhz = k + 1;
2908  }
2909  working_space[i][j][k] = lda;
2910  area = area + lda;
2911  if(lda > maximum){
2912  maximum = lda;
2913  positx = i,posity = j,positz = k;
2914  }
2915  }
2916  }
2917  }
2918  //read source cube
2919  for(i = 0;i < sizex_ext; i++){
2920  for(j = 0;j < sizey_ext; j++){
2921  for(k = 0;k < sizez_ext; k++){
2922  working_space[i][j][k + 2 * sizez_ext] = TMath::Abs(working_space[i][j][k + sizez_ext]);
2923  }
2924  }
2925  }
2926  //calculate ht*y and write into p
2927  for (i3 = 0; i3 < sizez_ext; i3++) {
2928  for (i2 = 0; i2 < sizey_ext; i2++) {
2929  for (i1 = 0; i1 < sizex_ext; i1++) {
2930  ldc = 0;
2931  for (j3 = 0; j3 <= (lhz - 1); j3++) {
2932  for (j2 = 0; j2 <= (lhy - 1); j2++) {
2933  for (j1 = 0; j1 <= (lhx - 1); j1++) {
2934  k3 = i3 + j3, k2 = i2 + j2, k1 = i1 + j1;
2935  if (k3 >= 0 && k3 < sizez_ext && k2 >= 0 && k2 < sizey_ext && k1 >= 0 && k1 < sizex_ext) {
2936  lda = working_space[j1][j2][j3];
2937  ldb = working_space[k1][k2][k3+2*sizez_ext];
2938  ldc = ldc + lda * ldb;
2939  }
2940  }
2941  }
2942  }
2943  working_space[i1][i2][i3 + sizez_ext] = ldc;
2944  }
2945  }
2946  }
2947 //calculate b=ht*h
2948  i1min = -(lhx - 1), i1max = lhx - 1;
2949  i2min = -(lhy - 1), i2max = lhy - 1;
2950  i3min = -(lhz - 1), i3max = lhz - 1;
2951  for (i3 = i3min; i3 <= i3max; i3++) {
2952  for (i2 = i2min; i2 <= i2max; i2++) {
2953  for (i1 = i1min; i1 <= i1max; i1++) {
2954  ldc = 0;
2955  j3min = -i3;
2956  if (j3min < 0)
2957  j3min = 0;
2958 
2959  j3max = lhz - 1 - i3;
2960  if (j3max > lhz - 1)
2961  j3max = lhz - 1;
2962 
2963  for (j3 = j3min; j3 <= j3max; j3++) {
2964  j2min = -i2;
2965  if (j2min < 0)
2966  j2min = 0;
2967 
2968  j2max = lhy - 1 - i2;
2969  if (j2max > lhy - 1)
2970  j2max = lhy - 1;
2971 
2972  for (j2 = j2min; j2 <= j2max; j2++) {
2973  j1min = -i1;
2974  if (j1min < 0)
2975  j1min = 0;
2976 
2977  j1max = lhx - 1 - i1;
2978  if (j1max > lhx - 1)
2979  j1max = lhx - 1;
2980 
2981  for (j1 = j1min; j1 <= j1max; j1++) {
2982  lda = working_space[j1][j2][j3];
2983  if (i1 + j1 < sizex_ext && i2 + j2 < sizey_ext)
2984  ldb = working_space[i1 + j1][i2 + j2][i3 + j3];
2985 
2986  else
2987  ldb = 0;
2988 
2989  ldc = ldc + lda * ldb;
2990  }
2991  }
2992  }
2993  working_space[i1 - i1min][i2 - i2min][i3 - i3min + 2 * sizez_ext ] = ldc;
2994  }
2995  }
2996  }
2997 //initialization in x1 cube
2998  for (i3 = 0; i3 < sizez_ext; i3++) {
2999  for (i2 = 0; i2 < sizey_ext; i2++) {
3000  for (i1 = 0; i1 < sizex_ext; i1++) {
3001  working_space[i1][i2][i3 + 3 * sizez_ext] = 1;
3002  working_space[i1][i2][i3 + 4 * sizez_ext] = 0;
3003  }
3004  }
3005  }
3006 
3007 //START OF ITERATIONS
3008  for (lindex=0;lindex<deconIterations;lindex++){
3009  for (i3 = 0; i3 < sizez_ext; i3++) {
3010  for (i2 = 0; i2 < sizey_ext; i2++) {
3011  for (i1 = 0; i1 < sizex_ext; i1++) {
3012  if (TMath::Abs(working_space[i1][i2][i3 + 3 * sizez_ext])>1e-6 && TMath::Abs(working_space[i1][i2][i3 + 1 * sizez_ext])>1e-6){
3013  ldb = 0;
3014  j3min = i3;
3015  if (j3min > lhz - 1)
3016  j3min = lhz - 1;
3017 
3018  j3min = -j3min;
3019  j3max = sizez_ext - i3 - 1;
3020  if (j3max > lhz - 1)
3021  j3max = lhz - 1;
3022 
3023  j2min = i2;
3024  if (j2min > lhy - 1)
3025  j2min = lhy - 1;
3026 
3027  j2min = -j2min;
3028  j2max = sizey_ext - i2 - 1;
3029  if (j2max > lhy - 1)
3030  j2max = lhy - 1;
3031 
3032  j1min = i1;
3033  if (j1min > lhx - 1)
3034  j1min = lhx - 1;
3035 
3036  j1min = -j1min;
3037  j1max = sizex_ext - i1 - 1;
3038  if (j1max > lhx - 1)
3039  j1max = lhx - 1;
3040 
3041  for (j3 = j3min; j3 <= j3max; j3++) {
3042  for (j2 = j2min; j2 <= j2max; j2++) {
3043  for (j1 = j1min; j1 <= j1max; j1++) {
3044  ldc = working_space[j1 - i1min][j2 - i2min][j3 - i3min + 2 * sizez_ext];
3045  lda = working_space[i1 + j1][i2 + j2][i3 + j3 + 3 * sizez_ext];
3046  ldb = ldb + lda * ldc;
3047  }
3048  }
3049  }
3050  lda = working_space[i1][i2][i3 + 3 * sizez_ext];
3051  ldc = working_space[i1][i2][i3 + 1 * sizez_ext];
3052  if (ldc * lda != 0 && ldb != 0) {
3053  lda = lda * ldc / ldb;
3054  }
3055 
3056  else
3057  lda = 0;
3058  working_space[i1][i2][i3 + 4 * sizez_ext] = lda;
3059  }
3060  }
3061  }
3062  }
3063  for (i3 = 0; i3 < sizez_ext; i3++) {
3064  for (i2 = 0; i2 < sizey_ext; i2++) {
3065  for (i1 = 0; i1 < sizex_ext; i1++)
3066  working_space[i1][i2][i3 + 3 * sizez_ext] = working_space[i1][i2][i3 + 4 * sizez_ext];
3067  }
3068  }
3069  }
3070 //write back resulting spectrum
3071  maximum=0;
3072  for(i = 0;i < sizex_ext; i++){
3073  for(j = 0;j < sizey_ext; j++){
3074  for(k = 0;k < sizez_ext; k++){
3075  working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext] = area * working_space[i][j][k + 3 * sizez_ext];
3076  if(maximum < working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext])
3077  maximum = working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext];
3078  }
3079  }
3080  }
3081 //searching for peaks in deconvolved spectrum
3082  for(i = 1;i < sizex_ext - 1; i++){
3083  for(j = 1;j < sizey_ext - 1; j++){
3084  for(l = 1;l < sizez_ext - 1; l++){
3085  a = working_space[i][j][l];
3086  if(a > working_space[i][j][l - 1] && a > working_space[i - 1][j][l - 1] && a > working_space[i - 1][j - 1][l - 1] && a > working_space[i][j - 1][l - 1] && a > working_space[i + 1][j - 1][l - 1] && a > working_space[i + 1][j][l - 1] && a > working_space[i + 1][j + 1][l - 1] && a > working_space[i][j + 1][l - 1] && a > working_space[i - 1][j + 1][l - 1] && a > working_space[i - 1][j][l] && a > working_space[i - 1][j - 1][l] && a > working_space[i][j - 1][l] && a > working_space[i + 1][j - 1][l] && a > working_space[i + 1][j][l] && a > working_space[i + 1][j + 1][l] && a > working_space[i][j + 1][l] && a > working_space[i - 1][j + 1][l] && a > working_space[i][j][l + 1] && a > working_space[i - 1][j][l + 1] && a > working_space[i - 1][j - 1][l + 1] && a > working_space[i][j - 1][l + 1] && a > working_space[i + 1][j - 1][l + 1] && a > working_space[i + 1][j][l + 1] && a > working_space[i + 1][j + 1][l + 1] && a > working_space[i][j + 1][l + 1] && a > working_space[i - 1][j + 1][l + 1]){
3087  if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift && l >= shift && l < ssizez + shift){
3088  if(working_space[i][j][l] > threshold * maximum / 100.0){
3089  if(peak_index < fMaxPeaks){
3090  for(k = i - 1,a = 0,b = 0;k <= i + 1; k++){
3091  a += (Double_t)(k - shift) * working_space[k][j][l];
3092  b += working_space[k][j][l];
3093  }
3094  fPositionX[peak_index] = a / b;
3095  for(k = j - 1,a = 0,b = 0;k <= j + 1; k++){
3096  a += (Double_t)(k - shift) * working_space[i][k][l];
3097  b += working_space[i][k][l];
3098  }
3099  fPositionY[peak_index] = a / b;
3100  for(k = l - 1,a = 0,b = 0;k <= l + 1; k++){
3101  a += (Double_t)(k - shift) * working_space[i][j][k];
3102  b += working_space[i][j][k];
3103  }
3104  fPositionZ[peak_index] = a / b;
3105  peak_index += 1;
3106  }
3107  }
3108  }
3109  }
3110  }
3111  }
3112  }
3113  for(i = 0;i < ssizex; i++){
3114  for(j = 0;j < ssizey; j++){
3115  for(k = 0;k < ssizez; k++){
3116  dest[i][j][k] = working_space[i + shift][j + shift][k + shift];
3117  }
3118  }
3119  }
3120  k = (Int_t)(4 * sigma + 0.5);
3121  k = 4 * k;
3122  for(i = 0;i < ssizex + k; i++){
3123  for(j = 0;j < ssizey + k; j++)
3124  delete[] working_space[i][j];
3125  delete[] working_space[i];
3126  }
3127  delete[] working_space;
3128  fNPeaks = peak_index;
3129  return fNPeaks;
3130 }
3131 
3132 ////////////////////////////////////////////////////////////////////////////////
3133 /// THREE-DIMENSIONAL CLASSICAL PEAK SEARCH FUNCTION
3134 /// This function searches for peaks in source spectrum using
3135 /// the algorithm based on smoothed second differences.
3136 ///
3137 /// Function parameters:
3138 /// - source-pointer to the matrix of source spectrum
3139 /// - ssizex-x length of source spectrum
3140 /// - ssizey-y length of source spectrum
3141 /// - ssizez-z length of source spectrum
3142 /// - sigma-sigma of searched peaks, for details we refer to manual
3143 /// - threshold-threshold value in % for selected peaks, peaks with
3144 /// amplitude less than threshold*highest_peak/100
3145 /// are ignored, see manual
3146 /// - markov-logical variable, if it is true, first the source spectrum
3147 /// is replaced by new spectrum calculated using Markov
3148 /// chains method.
3149 /// - averWindow-averaging window of searched peaks, for details
3150 /// we refer to manual (applies only for Markov method)
3151 
3152 Int_t TSpectrum3::SearchFast(const Double_t***source, Double_t***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez,
3153  Double_t sigma, Double_t threshold,
3154  Bool_t markov, Int_t averWindow)
3155 
3156 {
3157  Int_t i,j,k,l,li,lj,lk,lmin,lmax,xmin,xmax,ymin,ymax,zmin,zmax;
3158  Double_t maxch,plocha = 0,plocha_markov = 0;
3159  Double_t nom,nip,nim,sp,sm,spx,spy,smx,smy,spz,smz;
3160  Double_t norma,val,val1,val2,val3,val4,val5,val6,val7,val8,val9,val10,val11,val12,val13,val14,val15,val16,val17,val18,val19,val20,val21,val22,val23,val24,val25,val26;
3161  Double_t a,b,s,f,maximum;
3162  Int_t x,y,z,peak_index=0;
3163  Double_t p1,p2,p3,p4,p5,p6,p7,p8,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,r1,r2,r3,r4,r5,r6;
3164  Double_t pocet_sigma = 5;
3165  Int_t number_of_iterations=(Int_t)(4 * sigma + 0.5);
3166  Int_t sizex_ext=ssizex + 4 * number_of_iterations,sizey_ext = ssizey + 4 * number_of_iterations,sizez_ext = ssizez + 4 * number_of_iterations,shift = 2 * number_of_iterations;
3167  Double_t c[PEAK_WINDOW],s_f_ratio_peaks = 5;
3168  if(sigma < 1){
3169  Error("SearchFast", "Invalid sigma, must be greater than or equal to 1");
3170  return 0;
3171  }
3172 
3173  if(threshold<=0||threshold>=100){
3174  Error("SearchFast", "Invalid threshold, must be positive and less than 100");
3175  return 0;
3176  }
3177 
3178  j = (Int_t)(pocet_sigma*sigma+0.5);
3179  if (j >= PEAK_WINDOW / 2) {
3180  Error("SearchFast", "Too large sigma");
3181  return 0;
3182  }
3183 
3184  if (markov == true) {
3185  if (averWindow <= 0) {
3186  Error("SearchFast", "Averaging window must be positive");
3187  return 0;
3188  }
3189  }
3190 
3191  if(sizex_ext < 2 * number_of_iterations + 1 || sizey_ext < 2 * number_of_iterations + 1 || sizez_ext < 2 * number_of_iterations + 1){
3192  Error("SearchFast", "Too large clipping window");
3193  return 0;
3194  }
3195 
3196  i = (Int_t)(4 * sigma + 0.5);
3197  i = 4 * i;
3198  Double_t ***working_space = new Double_t** [ssizex + i];
3199  for(j = 0;j < ssizex + i; j++){
3200  working_space[j] = new Double_t* [ssizey + i];
3201  for(k = 0;k < ssizey + i; k++)
3202  working_space[j][k] = new Double_t [4 * (ssizez + i)];
3203  }
3204 
3205  for(k = 0;k < sizez_ext; k++){
3206  for(j = 0;j < sizey_ext; j++){
3207  for(i = 0;i < sizex_ext; i++){
3208  if(i < shift){
3209  if(j < shift){
3210  if(k < shift)
3211  working_space[i][j][k + sizez_ext] = source[0][0][0];
3212 
3213  else if(k >= ssizez + shift)
3214  working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1];
3215 
3216  else
3217  working_space[i][j][k + sizez_ext] = source[0][0][k - shift];
3218  }
3219 
3220  else if(j >= ssizey + shift){
3221  if(k < shift)
3222  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0];
3223 
3224  else if(k >= ssizez + shift)
3225  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1];
3226 
3227  else
3228  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift];
3229  }
3230 
3231  else{
3232  if(k < shift)
3233  working_space[i][j][k + sizez_ext] = source[0][j - shift][0];
3234 
3235  else if(k >= ssizez + shift)
3236  working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1];
3237 
3238  else
3239  working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift];
3240  }
3241  }
3242 
3243  else if(i >= ssizex + shift){
3244  if(j < shift){
3245  if(k < shift)
3246  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0];
3247 
3248  else if(k >= ssizez + shift)
3249  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1];
3250 
3251  else
3252  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift];
3253  }
3254 
3255  else if(j >= ssizey + shift){
3256  if(k < shift)
3257  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0];
3258 
3259  else if(k >= ssizez + shift)
3260  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1];
3261 
3262  else
3263  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift];
3264  }
3265 
3266  else{
3267  if(k < shift)
3268  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0];
3269 
3270  else if(k >= ssizez + shift)
3271  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1];
3272 
3273  else
3274  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift];
3275  }
3276  }
3277 
3278  else{
3279  if(j < shift){
3280  if(k < shift)
3281  working_space[i][j][k + sizez_ext] = source[i - shift][0][0];
3282 
3283  else if(k >= ssizez + shift)
3284  working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1];
3285 
3286  else
3287  working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift];
3288  }
3289 
3290  else if(j >= ssizey + shift){
3291  if(k < shift)
3292  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0];
3293 
3294  else if(k >= ssizez + shift)
3295  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1];
3296 
3297  else
3298  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift];
3299  }
3300 
3301  else{
3302  if(k < shift)
3303  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0];
3304 
3305  else if(k >= ssizez + shift)
3306  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1];
3307 
3308  else
3309  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift];
3310  }
3311  }
3312  }
3313  }
3314  }
3315  for(i = 1;i <= number_of_iterations; i++){
3316  for(z = i;z < sizez_ext - i; z++){
3317  for(y = i;y < sizey_ext - i; y++){
3318  for(x = i;x < sizex_ext - i; x++){
3319  a = working_space[x][y][z + sizez_ext];
3320  p1 = working_space[x + i][y + i][z - i + sizez_ext];
3321  p2 = working_space[x - i][y + i][z - i + sizez_ext];
3322  p3 = working_space[x + i][y - i][z - i + sizez_ext];
3323  p4 = working_space[x - i][y - i][z - i + sizez_ext];
3324  p5 = working_space[x + i][y + i][z + i + sizez_ext];
3325  p6 = working_space[x - i][y + i][z + i + sizez_ext];
3326  p7 = working_space[x + i][y - i][z + i + sizez_ext];
3327  p8 = working_space[x - i][y - i][z + i + sizez_ext];
3328  s1 = working_space[x + i][y ][z - i + sizez_ext];
3329  s2 = working_space[x ][y + i][z - i + sizez_ext];
3330  s3 = working_space[x - i][y ][z - i + sizez_ext];
3331  s4 = working_space[x ][y - i][z - i + sizez_ext];
3332  s5 = working_space[x + i][y ][z + i + sizez_ext];
3333  s6 = working_space[x ][y + i][z + i + sizez_ext];
3334  s7 = working_space[x - i][y ][z + i + sizez_ext];
3335  s8 = working_space[x ][y - i][z + i + sizez_ext];
3336  s9 = working_space[x - i][y + i][z + sizez_ext];
3337  s10 = working_space[x - i][y - i][z +sizez_ext];
3338  s11 = working_space[x + i][y + i][z +sizez_ext];
3339  s12 = working_space[x + i][y - i][z +sizez_ext];
3340  r1 = working_space[x ][y ][z - i + sizez_ext];
3341  r2 = working_space[x ][y ][z + i + sizez_ext];
3342  r3 = working_space[x - i][y ][z + sizez_ext];
3343  r4 = working_space[x + i][y ][z + sizez_ext];
3344  r5 = working_space[x ][y + i][z + sizez_ext];
3345  r6 = working_space[x ][y - i][z + sizez_ext];
3346  b = (p1 + p3) / 2.0;
3347  if(b > s1)
3348  s1 = b;
3349 
3350  b = (p1 + p2) / 2.0;
3351  if(b > s2)
3352  s2 = b;
3353 
3354  b = (p2 + p4) / 2.0;
3355  if(b > s3)
3356  s3 = b;
3357 
3358  b = (p3 + p4) / 2.0;
3359  if(b > s4)
3360  s4 = b;
3361 
3362  b = (p5 + p7) / 2.0;
3363  if(b > s5)
3364  s5 = b;
3365 
3366  b = (p5 + p6) / 2.0;
3367  if(b > s6)
3368  s6 = b;
3369 
3370  b = (p6 + p8) / 2.0;
3371  if(b > s7)
3372  s7 = b;
3373 
3374  b = (p7 + p8) / 2.0;
3375  if(b > s8)
3376  s8 = b;
3377 
3378  b = (p2 + p6) / 2.0;
3379  if(b > s9)
3380  s9 = b;
3381 
3382  b = (p4 + p8) / 2.0;
3383  if(b > s10)
3384  s10 = b;
3385 
3386  b = (p1 + p5) / 2.0;
3387  if(b > s11)
3388  s11 = b;
3389 
3390  b = (p3 + p7) / 2.0;
3391  if(b > s12)
3392  s12 = b;
3393 
3394  s1 = s1 - (p1 + p3) / 2.0;
3395  s2 = s2 - (p1 + p2) / 2.0;
3396  s3 = s3 - (p2 + p4) / 2.0;
3397  s4 = s4 - (p3 + p4) / 2.0;
3398  s5 = s5 - (p5 + p7) / 2.0;
3399  s6 = s6 - (p5 + p6) / 2.0;
3400  s7 = s7 - (p6 + p8) / 2.0;
3401  s8 = s8 - (p7 + p8) / 2.0;
3402  s9 = s9 - (p2 + p6) / 2.0;
3403  s10 = s10 - (p4 + p8) / 2.0;
3404  s11 = s11 - (p1 + p5) / 2.0;
3405  s12 = s12 - (p3 + p7) / 2.0;
3406  b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
3407  if(b > r1)
3408  r1 = b;
3409 
3410  b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
3411  if(b > r2)
3412  r2 = b;
3413 
3414  b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
3415  if(b > r3)
3416  r3 = b;
3417 
3418  b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
3419  if(b > r4)
3420  r4 = b;
3421 
3422  b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
3423  if(b > r5)
3424  r5 = b;
3425 
3426  b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
3427  if(b > r6)
3428  r6 = b;
3429 
3430  r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
3431  r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
3432  r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
3433  r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
3434  r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
3435  r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
3436  b = (r1 + r2) / 2.0 + (r3 + r4) / 2.0 + (r5 + r6) / 2.0 + (s1 + s3 + s5 + s7) / 4.0 + (s2 + s4 + s6 + s8) / 4.0 + (s9 + s10 + s11 + s12) / 4.0 + (p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.0;
3437  if(b < a)
3438  a = b;
3439 
3440  working_space[x][y][z] = a;
3441  }
3442  }
3443  }
3444  for(z = i;z < sizez_ext - i; z++){
3445  for(y = i;y < sizey_ext - i; y++){
3446  for(x = i;x < sizex_ext - i; x++){
3447  working_space[x][y][z + sizez_ext] = working_space[x][y][z];
3448  }
3449  }
3450  }
3451  }
3452  for(k = 0;k < sizez_ext; k++){
3453  for(j = 0;j < sizey_ext; j++){
3454  for(i = 0;i < sizex_ext; i++){
3455  if(i < shift){
3456  if(j < shift){
3457  if(k < shift)
3458  working_space[i][j][k + 3 * sizez_ext] = source[0][0][0] - working_space[i][j][k + sizez_ext];
3459 
3460  else if(k >= ssizez + shift)
3461  working_space[i][j][k + 3 * sizez_ext] = source[0][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3462 
3463  else
3464  working_space[i][j][k + 3 * sizez_ext] = source[0][0][k - shift] - working_space[i][j][k + sizez_ext];
3465  }
3466 
3467  else if(j >= ssizey + shift){
3468  if(k < shift)
3469  working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3470 
3471  else if(k >= ssizez + shift)
3472  working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3473 
3474  else
3475  working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3476  }
3477 
3478  else{
3479  if(k < shift)
3480  working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][0] - working_space[i][j][k + sizez_ext];
3481 
3482  else if(k >= ssizez + shift)
3483  working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3484 
3485  else
3486  working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3487  }
3488  }
3489 
3490  else if(i >= ssizex + shift){
3491  if(j < shift){
3492  if(k < shift)
3493  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][0] - working_space[i][j][k + sizez_ext];
3494 
3495  else if(k >= ssizez + shift)
3496  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3497 
3498  else
3499  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][k - shift] - working_space[i][j][k + sizez_ext];
3500  }
3501 
3502  else if(j >= ssizey + shift){
3503  if(k < shift)
3504  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3505 
3506  else if(k >= ssizez + shift)
3507  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3508 
3509  else
3510  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3511  }
3512 
3513  else{
3514  if(k < shift)
3515  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][0] - working_space[i][j][k + sizez_ext];
3516 
3517  else if(k >= ssizez + shift)
3518  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3519 
3520  else
3521  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3522  }
3523  }
3524 
3525  else{
3526  if(j < shift){
3527  if(k < shift)
3528  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][0] - working_space[i][j][k + sizez_ext];
3529 
3530  else if(k >= ssizez + shift)
3531  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3532 
3533  else
3534  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][k - shift] - working_space[i][j][k + sizez_ext];
3535  }
3536 
3537  else if(j >= ssizey + shift){
3538  if(k < shift)
3539  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3540 
3541  else if(k >= ssizez + shift)
3542  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3543 
3544  else
3545  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3546  }
3547 
3548  else{
3549  if(k < shift)
3550  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][0] - working_space[i][j][k + sizez_ext];
3551 
3552  else if(k >= ssizez + shift)
3553  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3554 
3555  else
3556  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3557  }
3558  }
3559  }
3560  }
3561  }
3562 
3563  for(i = 0;i < sizex_ext; i++){
3564  for(j = 0;j < sizey_ext; j++){
3565  for(k = 0;k < sizez_ext; k++){
3566  if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift && k >= shift && k < ssizez + shift){
3567  working_space[i][j][k + 2 * sizez_ext] = source[i - shift][j - shift][k - shift];
3568  plocha_markov = plocha_markov + source[i - shift][j - shift][k - shift];
3569  }
3570  else
3571  working_space[i][j][k + 2 * sizez_ext] = 0;
3572  }
3573  }
3574  }
3575 
3576  if(markov == true){
3577  xmin = 0;
3578  xmax = sizex_ext - 1;
3579  ymin = 0;
3580  ymax = sizey_ext - 1;
3581  zmin = 0;
3582  zmax = sizez_ext - 1;
3583  for(i = 0,maxch = 0;i < sizex_ext; i++){
3584  for(j = 0;j < sizey_ext;j++){
3585  for(k = 0;k < sizez_ext;k++){
3586  working_space[i][j][k] = 0;
3587  if(maxch < working_space[i][j][k + 2 * sizez_ext])
3588  maxch = working_space[i][j][k + 2 * sizez_ext];
3589 
3590  plocha += working_space[i][j][k + 2 * sizez_ext];
3591  }
3592  }
3593  }
3594  if(maxch == 0) {
3595  delete [] working_space;
3596  return 0;
3597  }
3598 
3599  nom = 0;
3600  working_space[xmin][ymin][zmin] = 1;
3601  for(i = xmin;i < xmax; i++){
3602  nip = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
3603  nim = working_space[i + 1][ymin][zmin + 2 * sizez_ext] / maxch;
3604  sp = 0,sm = 0;
3605  for(l = 1;l <= averWindow; l++){
3606  if((i + l) > xmax)
3607  a = working_space[xmax][ymin][zmin + 2 * sizez_ext] / maxch;
3608 
3609  else
3610  a = working_space[i + l][ymin][zmin + 2 * sizez_ext] / maxch;
3611 
3612  b = a - nip;
3613  if(a + nip <= 0)
3614  a = 1;
3615 
3616  else
3617  a = TMath::Sqrt(a + nip);
3618 
3619  b = b / a;
3620  b = TMath::Exp(b);
3621  sp = sp + b;
3622  if(i - l + 1 < xmin)
3623  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
3624 
3625  else
3626  a = working_space[i - l + 1][ymin][zmin + 2 * sizez_ext] / maxch;
3627 
3628  b = a - nim;
3629  if(a + nim <= 0)
3630  a = 1;
3631 
3632  else
3633  a = TMath::Sqrt(a + nim);
3634 
3635  b = b / a;
3636  b = TMath::Exp(b);
3637  sm = sm + b;
3638  }
3639  a = sp / sm;
3640  a = working_space[i + 1][ymin][zmin] = a * working_space[i][ymin][zmin];
3641  nom = nom + a;
3642  }
3643  for(i = ymin;i < ymax; i++){
3644  nip = working_space[xmin][i][zmin + 2 * sizez_ext] / maxch;
3645  nim = working_space[xmin][i + 1][zmin + 2 * sizez_ext] / maxch;
3646  sp = 0,sm = 0;
3647  for(l = 1;l <= averWindow; l++){
3648  if((i + l) > ymax)
3649  a = working_space[xmin][ymax][zmin + 2 * sizez_ext] / maxch;
3650 
3651  else
3652  a = working_space[xmin][i + l][zmin + 2 * sizez_ext] / maxch;
3653 
3654  b = a - nip;
3655  if(a + nip <= 0)
3656  a = 1;
3657 
3658  else
3659  a = TMath::Sqrt(a + nip);
3660 
3661  b = b / a;
3662  b = TMath::Exp(b);
3663  sp = sp + b;
3664  if(i - l + 1 < ymin)
3665  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
3666 
3667  else
3668  a = working_space[xmin][i - l + 1][zmin + 2 * sizez_ext] / maxch;
3669 
3670  b = a - nim;
3671  if(a + nim <= 0)
3672  a = 1;
3673 
3674  else
3675  a = TMath::Sqrt(a + nim);
3676 
3677  b = b / a;
3678  b = TMath::Exp(b);
3679  sm = sm + b;
3680  }
3681  a = sp / sm;
3682  a = working_space[xmin][i + 1][zmin] = a * working_space[xmin][i][zmin];
3683  nom = nom + a;
3684  }
3685  for(i = zmin;i < zmax;i++){
3686  nip = working_space[xmin][ymin][i + 2 * sizez_ext] / maxch;
3687  nim = working_space[xmin][ymin][i + 1 + 2 * sizez_ext] / maxch;
3688  sp = 0,sm = 0;
3689  for(l = 1;l <= averWindow;l++){
3690  if((i + l) > zmax)
3691  a = working_space[xmin][ymin][zmax + 2 * sizez_ext] / maxch;
3692 
3693  else
3694  a = working_space[xmin][ymin][i + l + 2 * sizez_ext] / maxch;
3695 
3696  b = a - nip;
3697  if(a + nip <= 0)
3698  a = 1;
3699 
3700  else
3701  a = TMath::Sqrt(a + nip);
3702 
3703  b = b / a;
3704  b = TMath::Exp(b);
3705  sp = sp + b;
3706  if(i - l + 1 < zmin)
3707  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
3708 
3709  else
3710  a = working_space[xmin][ymin][i - l + 1 + 2 * sizez_ext] / maxch;
3711 
3712  b = a - nim;
3713  if(a + nim <= 0)
3714  a = 1;
3715 
3716  else
3717  a = TMath::Sqrt(a + nim);
3718 
3719  b = b / a;
3720  b = TMath::Exp(b);
3721  sm = sm + b;
3722  }
3723  a = sp / sm;
3724  a = working_space[xmin][ymin][i + 1] = a * working_space[xmin][ymin][i];
3725  nom = nom + a;
3726  }
3727  for(i = xmin;i < xmax; i++){
3728  for(j = ymin;j < ymax; j++){
3729  nip = working_space[i][j + 1][zmin + 2 * sizez_ext] / maxch;
3730  nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
3731  spx = 0,smx = 0;
3732  for(l = 1;l <= averWindow; l++){
3733  if(i + l > xmax)
3734  a = working_space[xmax][j][zmin + 2 * sizez_ext] / maxch;
3735 
3736  else
3737  a = working_space[i + l][j][zmin + 2 * sizez_ext] / maxch;
3738 
3739  b = a - nip;
3740  if(a + nip <= 0)
3741  a = 1;
3742 
3743  else
3744  a = TMath::Sqrt(a + nip);
3745 
3746  b = b / a;
3747  b = TMath::Exp(b);
3748  spx = spx + b;
3749  if(i - l + 1 < xmin)
3750  a = working_space[xmin][j][zmin + 2 * sizez_ext] / maxch;
3751 
3752  else
3753  a = working_space[i - l + 1][j][zmin + 2 * sizez_ext] / maxch;
3754 
3755  b = a - nim;
3756  if(a + nim <= 0)
3757  a = 1;
3758 
3759  else
3760  a = TMath::Sqrt(a + nim);
3761 
3762  b = b / a;
3763  b = TMath::Exp(b);
3764  smx = smx + b;
3765  }
3766  spy = 0,smy = 0;
3767  nip = working_space[i + 1][j][zmin + 2 * sizez_ext] / maxch;
3768  nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
3769  for(l = 1;l <= averWindow; l++){
3770  if(j + l > ymax)
3771  a = working_space[i][ymax][zmin + 2 * sizez_ext] / maxch;
3772 
3773  else
3774  a = working_space[i][j + l][zmin + 2 * sizez_ext] / maxch;
3775 
3776  b = a - nip;
3777  if(a + nip <= 0)
3778  a = 1;
3779 
3780  else
3781  a = TMath::Sqrt(a + nip);
3782 
3783  b = b / a;
3784  b = TMath::Exp(b);
3785  spy = spy + b;
3786  if(j - l + 1 < ymin)
3787  a = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
3788 
3789  else
3790  a = working_space[i][j - l + 1][zmin + 2 * sizez_ext] / maxch;
3791 
3792  b = a - nim;
3793  if(a + nim <= 0)
3794  a = 1;
3795 
3796  else
3797  a = TMath::Sqrt(a + nim);
3798 
3799  b = b / a;
3800  b = TMath::Exp(b);
3801  smy = smy + b;
3802  }
3803  a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
3804  working_space[i + 1][j + 1][zmin] = a;
3805  nom = nom + a;
3806  }
3807  }
3808  for(i = xmin;i < xmax;i++){
3809  for(j = zmin;j < zmax;j++){
3810  nip = working_space[i][ymin][j + 1 + 2 * sizez_ext] / maxch;
3811  nim = working_space[i + 1][ymin][j + 1 + 2 * sizez_ext] / maxch;
3812  spx = 0,smx = 0;
3813  for(l = 1;l <= averWindow; l++){
3814  if(i + l > xmax)
3815  a = working_space[xmax][ymin][j + 2 * sizez_ext] / maxch;
3816 
3817  else
3818  a = working_space[i + l][ymin][j + 2 * sizez_ext] / maxch;
3819 
3820  b = a - nip;
3821  if(a + nip <= 0)
3822  a = 1;
3823 
3824  else
3825  a = TMath::Sqrt(a + nip);
3826 
3827  b = b / a;
3828  b = TMath::Exp(b);
3829  spx = spx + b;
3830  if(i - l + 1 < xmin)
3831  a = working_space[xmin][ymin][j + 2 * sizez_ext] / maxch;
3832 
3833  else
3834  a = working_space[i - l + 1][ymin][j + 2 * sizez_ext] / maxch;
3835 
3836  b = a - nim;
3837  if(a + nim <= 0)
3838  a = 1;
3839 
3840  else
3841  a = TMath::Sqrt(a + nim);
3842 
3843  b = b / a;
3844  b = TMath::Exp(b);
3845  smx = smx + b;
3846  }
3847  spz = 0,smz = 0;
3848  nip = working_space[i + 1][ymin][j + 2 * sizez_ext] / maxch;
3849  nim = working_space[i + 1][ymin][j + 1 + 2 * sizez_ext] / maxch;
3850  for(l = 1;l <= averWindow; l++){
3851  if(j + l > zmax)
3852  a = working_space[i][ymin][zmax + 2 * sizez_ext] / maxch;
3853 
3854  else
3855  a = working_space[i][ymin][j + l + 2 * sizez_ext] / maxch;
3856 
3857  b = a - nip;
3858  if(a + nip <= 0)
3859  a = 1;
3860 
3861  else
3862  a = TMath::Sqrt(a + nip);
3863 
3864  b = b / a;
3865  b = TMath::Exp(b);
3866  spz = spz + b;
3867  if(j - l + 1 < zmin)
3868  a = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
3869 
3870  else
3871  a = working_space[i][ymin][j - l + 1 + 2 * sizez_ext] / maxch;
3872 
3873  b = a - nim;
3874  if(a + nim <= 0)
3875  a = 1;
3876 
3877  else
3878  a = TMath::Sqrt(a + nim);
3879 
3880  b = b / a;
3881  b = TMath::Exp(b);
3882  smz = smz + b;
3883  }
3884  a = (spx * working_space[i][ymin][j + 1] + spz * working_space[i + 1][ymin][j]) / (smx + smz);
3885  working_space[i + 1][ymin][j + 1] = a;
3886  nom = nom + a;
3887  }
3888  }
3889  for(i = ymin;i < ymax;i++){
3890  for(j = zmin;j < zmax;j++){
3891  nip = working_space[xmin][i][j + 1 + 2 * sizez_ext] / maxch;
3892  nim = working_space[xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
3893  spy = 0,smy = 0;
3894  for(l = 1;l <= averWindow; l++){
3895  if(i + l > ymax)
3896  a = working_space[xmin][ymax][j + 2 * sizez_ext] / maxch;
3897 
3898  else
3899  a = working_space[xmin][i + l][j + 2 * sizez_ext] / maxch;
3900 
3901  b = a - nip;
3902  if(a + nip <= 0)
3903  a = 1;
3904 
3905  else
3906  a = TMath::Sqrt(a + nip);
3907 
3908  b = b / a;
3909  b = TMath::Exp(b);
3910  spy = spy + b;
3911  if(i - l + 1 < ymin)
3912  a = working_space[xmin][ymin][j + 2 * sizez_ext] / maxch;
3913 
3914  else
3915  a = working_space[xmin][i - l + 1][j + 2 * sizez_ext] / maxch;
3916 
3917  b = a - nim;
3918  if(a + nim <= 0)
3919  a = 1;
3920 
3921  else
3922  a = TMath::Sqrt(a + nim);
3923 
3924  b = b / a;
3925  b = TMath::Exp(b);
3926  smy = smy + b;
3927  }
3928  spz = 0,smz = 0;
3929  nip = working_space[xmin][i + 1][j + 2 * sizez_ext] / maxch;
3930  nim = working_space[xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
3931  for(l = 1;l <= averWindow; l++){
3932  if(j + l > zmax)
3933  a = working_space[xmin][i][zmax + 2 * sizez_ext] / maxch;
3934 
3935  else
3936  a = working_space[xmin][i][j + l + 2 * sizez_ext] / maxch;
3937 
3938  b = a - nip;
3939  if(a + nip <= 0)
3940  a = 1;
3941 
3942  else
3943  a = TMath::Sqrt(a + nip);
3944 
3945  b = b / a;
3946  b = TMath::Exp(b);
3947  spz = spz + b;
3948  if(j - l + 1 < zmin)
3949  a = working_space[xmin][i][zmin + 2 * sizez_ext] / maxch;
3950 
3951  else
3952  a = working_space[xmin][i][j - l + 1 + 2 * sizez_ext] / maxch;
3953 
3954  b = a - nim;
3955  if(a + nim <= 0)
3956  a = 1;
3957 
3958  else
3959  a = TMath::Sqrt(a + nim);
3960 
3961  b = b / a;
3962  b = TMath::Exp(b);
3963  smz = smz + b;
3964  }
3965  a = (spy * working_space[xmin][i][j + 1] + spz * working_space[xmin][i + 1][j]) / (smy + smz);
3966  working_space[xmin][i + 1][j + 1] = a;
3967  nom = nom + a;
3968  }
3969  }
3970  for(i = xmin;i < xmax; i++){
3971  for(j = ymin;j < ymax; j++){
3972  for(k = zmin;k < zmax; k++){
3973  nip = working_space[i][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3974  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
3975  spx = 0,smx = 0;
3976  for(l = 1;l <= averWindow; l++){
3977  if(i + l > xmax)
3978  a = working_space[xmax][j][k + 2 * sizez_ext] / maxch;
3979 
3980  else
3981  a = working_space[i + l][j][k + 2 * sizez_ext] / maxch;
3982 
3983  b = a - nip;
3984  if(a + nip <= 0)
3985  a = 1;
3986 
3987  else
3988  a = TMath::Sqrt(a + nip);
3989 
3990  b = b / a;
3991  b = TMath::Exp(b);
3992  spx = spx + b;
3993  if(i - l + 1 < xmin)
3994  a = working_space[xmin][j][k + 2 * sizez_ext] / maxch;
3995 
3996  else
3997  a = working_space[i - l + 1][j][k + 2 * sizez_ext] / maxch;
3998 
3999  b = a - nim;
4000  if(a + nim <= 0)
4001  a = 1;
4002 
4003  else
4004  a = TMath::Sqrt(a + nim);
4005 
4006  b = b / a;
4007  b = TMath::Exp(b);
4008  smx = smx + b;
4009  }
4010  spy = 0,smy = 0;
4011  nip = working_space[i + 1][j][k + 1 + 2 * sizez_ext] / maxch;
4012  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4013  for(l = 1;l <= averWindow; l++){
4014  if(j + l > ymax)
4015  a = working_space[i][ymax][k + 2 * sizez_ext] / maxch;
4016 
4017  else
4018  a = working_space[i][j + l][k + 2 * sizez_ext] / maxch;
4019 
4020  b = a - nip;
4021  if(a + nip <= 0)
4022  a = 1;
4023 
4024  else
4025  a = TMath::Sqrt(a + nip);
4026 
4027  b = b / a;
4028  b = TMath::Exp(b);
4029  spy = spy + b;
4030  if(j - l + 1 < ymin)
4031  a = working_space[i][ymin][k + 2 * sizez_ext] / maxch;
4032 
4033  else
4034  a = working_space[i][j - l + 1][k + 2 * sizez_ext] / maxch;
4035 
4036  b = a - nim;
4037  if(a + nim <= 0)
4038  a = 1;
4039 
4040  else
4041  a = TMath::Sqrt(a + nim);
4042 
4043  b = b / a;
4044  b = TMath::Exp(b);
4045  smy = smy + b;
4046  }
4047  spz = 0,smz = 0;
4048  nip = working_space[i + 1][j + 1][k + 2 * sizez_ext] / maxch;
4049  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4050  for(l = 1;l <= averWindow; l++ ){
4051  if(j + l > zmax)
4052  a = working_space[i][j][zmax + 2 * sizez_ext] / maxch;
4053 
4054  else
4055  a = working_space[i][j][k + l + 2 * sizez_ext] / maxch;
4056 
4057  b = a - nip;
4058  if(a + nip <= 0)
4059  a = 1;
4060 
4061  else
4062  a = TMath::Sqrt(a + nip);
4063 
4064  b = b / a;
4065  b = TMath::Exp(b);
4066  spz = spz + b;
4067  if(j - l + 1 < ymin)
4068  a = working_space[i][j][zmin + 2 * sizez_ext] / maxch;
4069 
4070  else
4071  a = working_space[i][j][k - l + 1 + 2 * sizez_ext] / maxch;
4072 
4073  b = a - nim;
4074  if(a + nim <= 0)
4075  a = 1;
4076 
4077  else
4078  a = TMath::Sqrt(a + nim);
4079 
4080  b = b / a;
4081  b = TMath::Exp(b);
4082  smz = smz + b;
4083  }
4084  a = (spx * working_space[i][j + 1][k + 1] + spy * working_space[i + 1][j][k + 1] + spz * working_space[i + 1][j + 1][k]) / (smx + smy + smz);
4085  working_space[i + 1][j + 1][k + 1] = a;
4086  nom = nom + a;
4087  }
4088  }
4089  }
4090  a = 0;
4091  for(i = xmin;i <= xmax; i++){
4092  for(j = ymin;j <= ymax; j++){
4093  for(k = zmin;k <= zmax; k++){
4094  working_space[i][j][k] = working_space[i][j][k] / nom;
4095  a+=working_space[i][j][k];
4096  }
4097  }
4098  }
4099  for(i = 0;i < sizex_ext; i++){
4100  for(j = 0;j < sizey_ext; j++){
4101  for(k = 0;k < sizez_ext; k++){
4102  working_space[i][j][k + 2 * sizez_ext] = working_space[i][j][k] * plocha_markov / a;
4103  }
4104  }
4105  }
4106  }
4107 
4108  maximum = 0;
4109  for(k = 0;k < ssizez; k++){
4110  for(j = 0;j < ssizey; j++){
4111  for(i = 0;i < ssizex; i++){
4112  working_space[i][j][k] = 0;
4113  working_space[i][j][sizez_ext + k] = 0;
4114  if(working_space[i][j][k + 3 * sizez_ext] > maximum)
4115  maximum=working_space[i][j][k+3*sizez_ext];
4116  }
4117  }
4118  }
4119  for(i = 0;i < PEAK_WINDOW; i++){
4120  c[i] = 0;
4121  }
4122  j = (Int_t)(pocet_sigma * sigma + 0.5);
4123  for(i = -j;i <= j; i++){
4124  a=i;
4125  a = -a * a;
4126  b = 2.0 * sigma * sigma;
4127  a = a / b;
4128  a = TMath::Exp(a);
4129  s = i;
4130  s = s * s;
4131  s = s - sigma * sigma;
4132  s = s / (sigma * sigma * sigma * sigma);
4133  s = s * a;
4134  c[PEAK_WINDOW / 2 + i] = s;
4135  }
4136  norma = 0;
4137  for(i = 0;i < PEAK_WINDOW; i++){
4138  norma = norma + TMath::Abs(c[i]);
4139  }
4140  for(i = 0;i < PEAK_WINDOW; i++){
4141  c[i] = c[i] / norma;
4142  }
4143  a = pocet_sigma * sigma + 0.5;
4144  i = (Int_t)a;
4145  zmin = i;
4146  zmax = sizez_ext - i - 1;
4147  ymin = i;
4148  ymax = sizey_ext - i - 1;
4149  xmin = i;
4150  xmax = sizex_ext - i - 1;
4151  lmin = PEAK_WINDOW / 2 - i;
4152  lmax = PEAK_WINDOW / 2 + i;
4153  for(i = xmin;i <= xmax; i++){
4154  for(j = ymin;j <= ymax; j++){
4155  for(k = zmin;k <= zmax; k++){
4156  s = 0,f = 0;
4157  for(li = lmin;li <= lmax; li++){
4158  for(lj = lmin;lj <= lmax; lj++){
4159  for(lk = lmin;lk <= lmax; lk++){
4160  a = working_space[i + li - PEAK_WINDOW / 2][j + lj - PEAK_WINDOW / 2][k + lk - PEAK_WINDOW / 2 + 2 * sizez_ext];
4161  b = c[li] * c[lj] * c[lk];
4162  s += a * b;
4163  f += a * b * b;
4164  }
4165  }
4166  }
4167  working_space[i][j][k] = s;
4168  working_space[i][j][k + sizez_ext] = TMath::Sqrt(f);
4169  }
4170  }
4171  }
4172  for(x = xmin;x <= xmax; x++){
4173  for(y = ymin + 1;y < ymax; y++){
4174  for(z = zmin + 1;z < zmax; z++){
4175  val = working_space[x][y][z];
4176  val1 = working_space[x - 1][y - 1][z - 1];
4177  val2 = working_space[x ][y - 1][z - 1];
4178  val3 = working_space[x + 1][y - 1][z - 1];
4179  val4 = working_space[x - 1][y ][z - 1];
4180  val5 = working_space[x ][y ][z - 1];
4181  val6 = working_space[x + 1][y ][z - 1];
4182  val7 = working_space[x - 1][y + 1][z - 1];
4183  val8 = working_space[x ][y + 1][z - 1];
4184  val9 = working_space[x + 1][y + 1][z - 1];
4185  val10 = working_space[x - 1][y - 1][z ];
4186  val11 = working_space[x ][y - 1][z ];
4187  val12 = working_space[x + 1][y - 1][z ];
4188  val13 = working_space[x - 1][y ][z ];
4189  val14 = working_space[x + 1][y ][z ];
4190  val15 = working_space[x - 1][y + 1][z ];
4191  val16 = working_space[x ][y + 1][z ];
4192  val17 = working_space[x + 1][y + 1][z ];
4193  val18 = working_space[x - 1][y - 1][z + 1];
4194  val19 = working_space[x ][y - 1][z + 1];
4195  val20 = working_space[x + 1][y - 1][z + 1];
4196  val21 = working_space[x - 1][y ][z + 1];
4197  val22 = working_space[x ][y ][z + 1];
4198  val23 = working_space[x + 1][y ][z + 1];
4199  val24 = working_space[x - 1][y + 1][z + 1];
4200  val25 = working_space[x ][y + 1][z + 1];
4201  val26 = working_space[x + 1][y + 1][z + 1];
4202  f = -s_f_ratio_peaks * working_space[x][y][z + sizez_ext];
4203  if(val < f && val < val1 && val < val2 && val < val3 && val < val4 && val < val5 && val < val6 && val < val7 && val < val8 && val < val9 && val < val10 && val < val11 && val < val12 && val < val13 && val < val14 && val < val15 && val < val16 && val < val17 && val < val18 && val < val19 && val < val20 && val < val21 && val < val22 && val < val23 && val < val24 && val < val25 && val < val26){
4204  s=0,f=0;
4205  for(li = lmin;li <= lmax; li++){
4206  a = working_space[x + li - PEAK_WINDOW / 2][y][z + 2 * sizez_ext];
4207  s += a * c[li];
4208  f += a * c[li] * c[li];
4209  }
4210  f = -s_f_ratio_peaks * sqrt(f);
4211  if(s < f){
4212  s = 0,f = 0;
4213  for(li = lmin;li <= lmax; li++){
4214  a = working_space[x][y + li - PEAK_WINDOW / 2][z + 2 * sizez_ext];
4215  s += a * c[li];
4216  f += a * c[li] * c[li];
4217  }
4218  f = -s_f_ratio_peaks * sqrt(f);
4219  if(s < f){
4220  s = 0,f = 0;
4221  for(li = lmin;li <= lmax; li++){
4222  a = working_space[x][y][z + li - PEAK_WINDOW / 2 + 2 * sizez_ext];
4223  s += a * c[li];
4224  f += a * c[li] * c[li];
4225  }
4226  f = -s_f_ratio_peaks * sqrt(f);
4227  if(s < f){
4228  s = 0,f = 0;
4229  for(li = lmin;li <= lmax; li++){
4230  for(lj = lmin;lj <= lmax; lj++){
4231  a = working_space[x + li - PEAK_WINDOW / 2][y + lj - PEAK_WINDOW / 2][z + 2 * sizez_ext];
4232  s += a * c[li] * c[lj];
4233  f += a * c[li] * c[li] * c[lj] * c[lj];
4234  }
4235  }
4236  f = s_f_ratio_peaks * sqrt(f);
4237  if(s > f){
4238  s = 0,f = 0;
4239  for(li = lmin;li <= lmax; li++){
4240  for(lj = lmin;lj <= lmax; lj++){
4241  a = working_space[x + li - PEAK_WINDOW / 2][y][z + lj - PEAK_WINDOW / 2 + 2 * sizez_ext];
4242  s += a * c[li] * c[lj];
4243  f += a * c[li] * c[li] * c[lj] * c[lj];
4244  }
4245  }
4246  f = s_f_ratio_peaks * sqrt(f);
4247  if(s > f){
4248  s = 0,f = 0;
4249  for(li = lmin;li <= lmax; li++){
4250  for(lj=lmin;lj<=lmax;lj++){
4251  a = working_space[x][y + li - PEAK_WINDOW / 2][z + lj - PEAK_WINDOW / 2 + 2 * sizez_ext];
4252  s += a * c[li] * c[lj];
4253  f += a * c[li] * c[li] * c[lj] * c[lj];
4254  }
4255  }
4256  f = s_f_ratio_peaks * sqrt(f);
4257  if(s > f){
4258  if(x >= shift && x < ssizex + shift && y >= shift && y < ssizey + shift && z >= shift && z < ssizez + shift){
4259  if(working_space[x][y][z + 3 * sizez_ext] > threshold * maximum / 100.0){
4260  if(peak_index<fMaxPeaks){
4261  for(k = x - 1,a = 0,b = 0;k <= x + 1; k++){
4262  a += (Double_t)(k - shift) * working_space[k][y][z];
4263  b += working_space[k][y][z];
4264  }
4265  fPositionX[peak_index] = a / b;
4266  for(k = y - 1,a = 0,b = 0;k <= y + 1; k++){
4267  a += (Double_t)(k - shift) * working_space[x][k][z];
4268  b += working_space[x][k][z];
4269  }
4270  fPositionY[peak_index] = a / b;
4271  for(k = z - 1,a = 0,b = 0;k <= z + 1; k++){
4272  a += (Double_t)(k - shift) * working_space[x][y][k];
4273  b += working_space[x][y][k];
4274  }
4275  fPositionZ[peak_index] = a / b;
4276  peak_index += 1;
4277  }
4278  }
4279  }
4280  }
4281  }
4282  }
4283  }
4284  }
4285  }
4286  }
4287  }
4288  }
4289  }
4290  for(i = 0;i < ssizex; i++){
4291  for(j = 0;j < ssizey; j++){
4292  for(k = 0;k < ssizez; k++){
4293  val = -working_space[i + shift][j + shift][k + shift];
4294  if( val < 0)
4295  val = 0;
4296  dest[i][j][k] = val;
4297  }
4298  }
4299  }
4300  k = (Int_t)(4 * sigma + 0.5);
4301  k = 4 * k;
4302  for(i = 0;i < ssizex + k; i++){
4303  for(j = 0;j < ssizey + k; j++)
4304  delete[] working_space[i][j];
4305  delete[] working_space[i];
4306  }
4307  delete[] working_space;
4308  fNPeaks = peak_index;
4309  return fNPeaks;
4310 }