Logo ROOT   6.21/01
Reference Guide
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 
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 = nullptr;
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 = nullptr;
89  fNPeaks = 0;
90  SetResolution(resolution);
91 }
92 
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 /// Destructor.
96 
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 
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 
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 
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  for(i = 0;i < ssizex; i++){
890  for(j = 0;j < ssizey; j++)
891  delete[] working_space[i][j];
892  delete[] working_space[i];
893  }
894  delete [] working_space;
895  return 0;
896  }
897 
898  nom = 0;
899  working_space[xmin][ymin][zmin] = 1;
900  for(i = xmin;i < xmax; i++){
901  nip = source[i][ymin][zmin] / maxch;
902  nim = source[i + 1][ymin][zmin] / maxch;
903  sp = 0,sm = 0;
904  for(l = 1;l <= averWindow; l++){
905  if((i + l) > xmax)
906  a = source[xmax][ymin][zmin] / maxch;
907 
908  else
909  a = source[i + l][ymin][zmin] / maxch;
910 
911  b = a - nip;
912  if(a + nip <= 0)
913  a = 1;
914 
915  else
916  a = TMath::Sqrt(a + nip);
917 
918  b = b / a;
919  b = TMath::Exp(b);
920  sp = sp + b;
921  if(i - l + 1 < xmin)
922  a = source[xmin][ymin][zmin] / maxch;
923 
924  else
925  a = source[i - l + 1][ymin][zmin] / maxch;
926 
927  b = a - nim;
928  if(a + nim <= 0)
929  a = 1;
930 
931  else
932  a = TMath::Sqrt(a + nim);
933 
934  b = b / a;
935  b = TMath::Exp(b);
936  sm = sm + b;
937  }
938  a = sp / sm;
939  a = working_space[i + 1][ymin][zmin] = a * working_space[i][ymin][zmin];
940  nom = nom + a;
941  }
942  for(i = ymin;i < ymax; i++){
943  nip = source[xmin][i][zmin] / maxch;
944  nim = source[xmin][i + 1][zmin] / maxch;
945  sp = 0,sm = 0;
946  for(l = 1;l <= averWindow; l++){
947  if((i + l) > ymax)
948  a = source[xmin][ymax][zmin] / maxch;
949 
950  else
951  a = source[xmin][i + l][zmin] / maxch;
952 
953  b = a - nip;
954  if(a + nip <= 0)
955  a = 1;
956 
957  else
958  a = TMath::Sqrt(a + nip);
959 
960  b = b / a;
961  b = TMath::Exp(b);
962  sp = sp + b;
963  if(i - l + 1 < ymin)
964  a = source[xmin][ymin][zmin] / maxch;
965 
966  else
967  a = source[xmin][i - l + 1][zmin] / maxch;
968 
969  b = a - nim;
970  if(a + nim <= 0)
971  a = 1;
972 
973  else
974  a = TMath::Sqrt(a + nim);
975 
976  b = b / a;
977  b = TMath::Exp(b);
978  sm = sm + b;
979  }
980  a = sp / sm;
981  a = working_space[xmin][i + 1][zmin] = a * working_space[xmin][i][zmin];
982  nom = nom + a;
983  }
984  for(i = zmin;i < zmax; i++){
985  nip = source[xmin][ymin][i] / maxch;
986  nim = source[xmin][ymin][i + 1] / maxch;
987  sp = 0,sm = 0;
988  for(l = 1;l <= averWindow; l++){
989  if((i + l) > zmax)
990  a = source[xmin][ymin][zmax] / maxch;
991 
992  else
993  a = source[xmin][ymin][i + l] / maxch;
994 
995  b = a - nip;
996  if(a + nip <= 0)
997  a = 1;
998 
999  else
1000  a = TMath::Sqrt(a + nip);
1001 
1002  b = b / a;
1003  b = TMath::Exp(b);
1004  sp = sp + b;
1005  if(i - l + 1 < zmin)
1006  a = source[xmin][ymin][zmin] / maxch;
1007 
1008  else
1009  a = source[xmin][ymin][i - l + 1] / maxch;
1010 
1011  b = a - nim;
1012  if(a + nim <= 0)
1013  a = 1;
1014 
1015  else
1016  a = TMath::Sqrt(a + nim);
1017  b = b / a;
1018  b = TMath::Exp(b);
1019  sm = sm + b;
1020  }
1021  a = sp / sm;
1022  a = working_space[xmin][ymin][i + 1] = a * working_space[xmin][ymin][i];
1023  nom = nom + a;
1024  }
1025  for(i = xmin;i < xmax; i++){
1026  for(j = ymin;j < ymax; j++){
1027  nip = source[i][j + 1][zmin] / maxch;
1028  nim = source[i + 1][j + 1][zmin] / maxch;
1029  spx = 0,smx = 0;
1030  for(l = 1;l <= averWindow; l++){
1031  if(i + l > xmax)
1032  a = source[xmax][j][zmin] / maxch;
1033 
1034  else
1035  a = source[i + l][j][zmin] / maxch;
1036 
1037  b = a - nip;
1038  if(a + nip <= 0)
1039  a = 1;
1040 
1041  else
1042  a = TMath::Sqrt(a + nip);
1043 
1044  b = b / a;
1045  b = TMath::Exp(b);
1046  spx = spx + b;
1047  if(i - l + 1 < xmin)
1048  a = source[xmin][j][zmin] / maxch;
1049 
1050  else
1051  a = source[i - l + 1][j][zmin] / maxch;
1052 
1053  b = a - nim;
1054  if(a + nim <= 0)
1055  a = 1;
1056 
1057  else
1058  a = TMath::Sqrt(a + nim);
1059 
1060  b = b / a;
1061  b = TMath::Exp(b);
1062  smx = smx + b;
1063  }
1064  spy = 0,smy = 0;
1065  nip = source[i + 1][j][zmin] / maxch;
1066  nim = source[i + 1][j + 1][zmin] / maxch;
1067  for(l = 1;l <= averWindow; l++){
1068  if(j + l > ymax)
1069  a = source[i][ymax][zmin] / maxch;
1070 
1071  else
1072  a = source[i][j + l][zmin] / maxch;
1073 
1074  b = a - nip;
1075  if(a + nip <= 0)
1076  a = 1;
1077 
1078  else
1079  a = TMath::Sqrt(a + nip);
1080 
1081  b = b / a;
1082  b = TMath::Exp(b);
1083  spy = spy + b;
1084  if(j - l + 1 < ymin)
1085  a = source[i][ymin][zmin] / maxch;
1086 
1087  else
1088  a = source[i][j - l + 1][zmin] / maxch;
1089 
1090  b = a - nim;
1091  if(a + nim <= 0)
1092  a = 1;
1093 
1094  else
1095  a = TMath::Sqrt(a + nim);
1096 
1097  b = b / a;
1098  b = TMath::Exp(b);
1099  smy = smy + b;
1100  }
1101  a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
1102  working_space[i + 1][j + 1][zmin] = a;
1103  nom = nom + a;
1104  }
1105  }
1106  for(i = xmin;i < xmax; i++){
1107  for(j = zmin;j < zmax; j++){
1108  nip = source[i][ymin][j + 1] / maxch;
1109  nim = source[i + 1][ymin][j + 1] / maxch;
1110  spx = 0,smx = 0;
1111  for(l = 1;l <= averWindow; l++){
1112  if(i + l > xmax)
1113  a = source[xmax][ymin][j] / maxch;
1114 
1115  else
1116  a = source[i + l][ymin][j] / maxch;
1117 
1118  b = a - nip;
1119  if(a + nip <= 0)
1120  a = 1;
1121 
1122  else
1123  a = TMath::Sqrt(a + nip);
1124 
1125  b = b / a;
1126  b = TMath::Exp(b);
1127  spx = spx + b;
1128  if(i - l + 1 < xmin)
1129  a = source[xmin][ymin][j] / maxch;
1130 
1131  else
1132  a = source[i - l + 1][ymin][j] / maxch;
1133 
1134  b = a - nim;
1135  if(a + nim <= 0)
1136  a = 1;
1137 
1138  else
1139  a = TMath::Sqrt(a + nim);
1140 
1141  b = b / a;
1142  b = TMath::Exp(b);
1143  smx = smx + b;
1144  }
1145  spz = 0,smz = 0;
1146  nip = source[i + 1][ymin][j] / maxch;
1147  nim = source[i + 1][ymin][j + 1] / maxch;
1148  for(l = 1;l <= averWindow; l++){
1149  if(j + l > zmax)
1150  a = source[i][ymin][zmax] / maxch;
1151 
1152  else
1153  a = source[i][ymin][j + l] / maxch;
1154 
1155  b = a - nip;
1156  if(a + nip <= 0)
1157  a = 1;
1158 
1159  else
1160  a = TMath::Sqrt(a + nip);
1161 
1162  b = b / a;
1163  b = TMath::Exp(b);
1164  spz = spz + b;
1165  if(j - l + 1 < zmin)
1166  a = source[i][ymin][zmin] / maxch;
1167 
1168  else
1169  a = source[i][ymin][j - l + 1] / maxch;
1170 
1171  b = a - nim;
1172  if(a + nim <= 0)
1173  a = 1;
1174 
1175  else
1176  a = TMath::Sqrt(a + nim);
1177 
1178  b = b / a;
1179  b = TMath::Exp(b);
1180  smz = smz + b;
1181  }
1182  a = (spx * working_space[i][ymin][j + 1] + spz * working_space[i + 1][ymin][j]) / (smx + smz);
1183  working_space[i + 1][ymin][j + 1] = a;
1184  nom = nom + a;
1185  }
1186  }
1187  for(i = ymin;i < ymax; i++){
1188  for(j = zmin;j < zmax; j++){
1189  nip = source[xmin][i][j + 1] / maxch;
1190  nim = source[xmin][i + 1][j + 1] / maxch;
1191  spy = 0,smy = 0;
1192  for(l = 1;l <= averWindow; l++){
1193  if(i + l > ymax)
1194  a = source[xmin][ymax][j] / maxch;
1195 
1196  else
1197  a = source[xmin][i + l][j] / maxch;
1198 
1199  b = a - nip;
1200  if(a + nip <= 0)
1201  a = 1;
1202 
1203  else
1204  a = TMath::Sqrt(a + nip);
1205 
1206  b = b / a;
1207  b = TMath::Exp(b);
1208  spy = spy + b;
1209  if(i - l + 1 < ymin)
1210  a = source[xmin][ymin][j] / maxch;
1211 
1212  else
1213  a = source[xmin][i - l + 1][j] / maxch;
1214 
1215  b = a - nim;
1216  if(a + nim <= 0)
1217  a = 1;
1218 
1219  else
1220  a = TMath::Sqrt(a + nim);
1221 
1222  b = b / a;
1223  b = TMath::Exp(b);
1224  smy = smy + b;
1225  }
1226  spz = 0,smz = 0;
1227  nip = source[xmin][i + 1][j] / maxch;
1228  nim = source[xmin][i + 1][j + 1] / maxch;
1229  for(l = 1;l <= averWindow; l++){
1230  if(j + l > zmax)
1231  a = source[xmin][i][zmax] / maxch;
1232 
1233  else
1234  a = source[xmin][i][j + l] / maxch;
1235 
1236  b = a - nip;
1237  if(a + nip <= 0)
1238  a = 1;
1239 
1240  else
1241  a = TMath::Sqrt(a + nip);
1242 
1243  b = b / a;
1244  b = TMath::Exp(b);
1245  spz = spz + b;
1246  if(j - l + 1 < zmin)
1247  a = source[xmin][i][zmin] / maxch;
1248 
1249  else
1250  a = source[xmin][i][j - l + 1] / maxch;
1251 
1252  b = a - nim;
1253  if(a + nim <= 0)
1254  a = 1;
1255 
1256  else
1257  a = TMath::Sqrt(a + nim);
1258 
1259  b = b / a;
1260  b = TMath::Exp(b);
1261  smz = smz + b;
1262  }
1263  a = (spy * working_space[xmin][i][j + 1] + spz * working_space[xmin][i + 1][j]) / (smy + smz);
1264  working_space[xmin][i + 1][j + 1] = a;
1265  nom = nom + a;
1266  }
1267  }
1268  for(i = xmin;i < xmax; i++){
1269  for(j = ymin;j < ymax; j++){
1270  for(k = zmin;k < zmax; k++){
1271  nip = source[i][j + 1][k + 1] / maxch;
1272  nim = source[i + 1][j + 1][k + 1] / maxch;
1273  spx = 0,smx = 0;
1274  for(l = 1;l <= averWindow; l++){
1275  if(i + l > xmax)
1276  a = source[xmax][j][k] / maxch;
1277 
1278  else
1279  a = source[i + l][j][k] / maxch;
1280 
1281  b = a - nip;
1282  if(a + nip <= 0)
1283  a = 1;
1284 
1285  else
1286  a = TMath::Sqrt(a + nip);
1287 
1288  b = b / a;
1289  b = TMath::Exp(b);
1290  spx = spx + b;
1291  if(i - l + 1 < xmin)
1292  a = source[xmin][j][k] / maxch;
1293 
1294  else
1295  a = source[i - l + 1][j][k] / maxch;
1296 
1297  b = a - nim;
1298  if(a + nim <= 0)
1299  a = 1;
1300 
1301  else
1302  a = TMath::Sqrt(a + nim);
1303 
1304  b = b / a;
1305  b = TMath::Exp(b);
1306  smx = smx + b;
1307  }
1308  spy = 0,smy = 0;
1309  nip = source[i + 1][j][k + 1] / maxch;
1310  nim = source[i + 1][j + 1][k + 1] / maxch;
1311  for(l = 1;l <= averWindow; l++){
1312  if(j + l > ymax)
1313  a = source[i][ymax][k] / maxch;
1314 
1315  else
1316  a = source[i][j + l][k] / maxch;
1317 
1318  b = a - nip;
1319  if(a + nip <= 0)
1320  a = 1;
1321 
1322  else
1323  a = TMath::Sqrt(a + nip);
1324 
1325  b = b / a;
1326  b = TMath::Exp(b);
1327  spy = spy + b;
1328  if(j - l + 1 < ymin)
1329  a = source[i][ymin][k] / maxch;
1330 
1331  else
1332  a = source[i][j - l + 1][k] / maxch;
1333 
1334  b = a - nim;
1335  if(a + nim <= 0)
1336  a = 1;
1337 
1338  else
1339  a = TMath::Sqrt(a + nim);
1340 
1341  b = b / a;
1342  b = TMath::Exp(b);
1343  smy = smy + b;
1344  }
1345  spz = 0,smz = 0;
1346  nip = source[i + 1][j + 1][k] / maxch;
1347  nim = source[i + 1][j + 1][k + 1] / maxch;
1348  for(l = 1;l <= averWindow; l++){
1349  if(j + l > zmax)
1350  a = source[i][j][zmax] / maxch;
1351 
1352  else
1353  a = source[i][j][k + l] / maxch;
1354 
1355  b = a - nip;
1356  if(a + nip <= 0)
1357  a = 1;
1358 
1359  else
1360  a = TMath::Sqrt(a + nip);
1361 
1362  b = b / a;
1363  b = TMath::Exp(b);
1364  spz = spz + b;
1365  if(j - l + 1 < ymin)
1366  a = source[i][j][zmin] / maxch;
1367 
1368  else
1369  a = source[i][j][k - l + 1] / maxch;
1370 
1371  b = a - nim;
1372  if(a + nim <= 0)
1373  a = 1;
1374 
1375  else
1376  a = TMath::Sqrt(a + nim);
1377 
1378  b = b / a;
1379  b = TMath::Exp(b);
1380  smz = smz+b;
1381  }
1382  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);
1383  working_space[i + 1][j + 1][k + 1] = a;
1384  nom = nom + a;
1385  }
1386  }
1387  }
1388  for(i = xmin;i <= xmax; i++){
1389  for(j = ymin;j <= ymax; j++){
1390  for(k = zmin;k <= zmax; k++){
1391  working_space[i][j][k] = working_space[i][j][k] / nom;
1392  }
1393  }
1394  }
1395  for(i = 0;i < ssizex; i++){
1396  for(j = 0;j < ssizey; j++){
1397  for(k = 0;k < ssizez; k++){
1398  source[i][j][k] = plocha * working_space[i][j][k];
1399  }
1400  }
1401  }
1402  for(i = 0;i < ssizex; i++){
1403  for(j = 0;j < ssizey; j++)
1404  delete[] working_space[i][j];
1405  delete[] working_space[i];
1406  }
1407  delete[] working_space;
1408  return 0;
1409 }
1410 
1411 ////////////////////////////////////////////////////////////////////////////////
1412 /// This function calculates deconvolution from source spectrum
1413 /// according to response spectrum
1414 /// The result is placed in the cube pointed by source pointer.
1415 ///
1416 /// Function parameters:
1417 /// - source-pointer to the cube of source spectrum
1418 /// - resp-pointer to the cube of response spectrum
1419 /// - ssizex-x length of source and response spectra
1420 /// - ssizey-y length of source and response spectra
1421 /// - ssizey-y length of source and response spectra
1422 /// - numberIterations, for details we refer to manual
1423 /// - numberRepetitions, for details we refer to manual
1424 /// - boost, boosting factor, for details we refer to manual
1425 ///
1426 /// ### Deconvolution
1427 ///
1428 /// Goal: Improvement of the resolution in spectra, decomposition of multiplets
1429 ///
1430 /// Mathematical formulation of the 3-dimensional convolution system is
1431 ///
1432 /// \image html spectrum3_deconvolution_image001.gif
1433 ///
1434 /// where h(i,j,k) is the impulse response function, x, y are
1435 /// input and output fields, respectively, \f$ N_1, N_2, N3\f$, are the lengths of x and h fields
1436 ///
1437 /// - let us assume that we know the response and the output fields (spectra)
1438 /// of the above given system.
1439 ///
1440 /// - the deconvolution represents solution of the overdetermined system of
1441 /// linear equations, i.e., the calculation of the field -x.
1442 ///
1443 /// - from numerical stability point of view the operation of deconvolution is
1444 /// extremely critical (ill-posed problem) as well as time consuming operation.
1445 ///
1446 /// - the Gold deconvolution algorithm proves to work very well even for
1447 /// 2-dimensional systems. Generalization of the algorithm for 2-dimensional
1448 /// systems was presented in [1], and for multidimensional systems in [2].
1449 ///
1450 /// - for Gold deconvolution algorithm as well as for boosted deconvolution
1451 /// algorithm we refer also to TSpectrum and TSpectrum2
1452 ///
1453 /// #### References:
1454 ///
1455 /// [1] M.Morhac, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.: Efficient
1456 /// one- and two-dimensional Gold deconvolution and its application to gamma-ray
1457 /// spectra decomposition. NIM, A401 (1997) 385-408.
1458 ///
1459 /// [2] Morhac M., Matouoek V.,
1460 /// Kliman J., Efficient algorithm of multidimensional deconvolution and its
1461 /// application to nuclear data processing, Digital Signal Processing 13 (2003) 144.
1462 ///
1463 /// ### Example 1 - script Decon.c :
1464 ///
1465 /// response function (usually peak) should be shifted to the beginning of
1466 /// the coordinate system (see Fig. 1)
1467 ///
1468 /// \image html spectrum3_deconvolution_image003.jpg Fig. 1 Three-dimensional response spectrum
1469 /// \image html spectrum3_deconvolution_image004.jpg Fig. 2 Three-dimensional input spectrum (before deconvolution)
1470 /// \image html spectrum3_deconvolution_image005.jpg Fig. 3 Spectrum from Fig. 2 after deconvolution (100 iterations)
1471 ///
1472 /// #### Script:
1473 ///
1474 /// Example to illustrate the Gold deconvolution (class TSpectrum3).
1475 /// To execute this example, do:
1476 ///
1477 /// `root > .x Decon3.C`
1478 ///
1479 /// ~~~ {.cpp}
1480 /// #include <TSpectrum3>
1481 /// void Decon3() {
1482 /// Int_t i, j, k;
1483 /// Int_t nbinsx = 32;
1484 /// Int_t nbinsy = 32;
1485 /// Int_t nbinsz = 32;
1486 /// Int_t xmin = 0;
1487 /// Int_t xmax = nbinsx;
1488 /// Int_t ymin = 0;
1489 /// Int_t ymax = nbinsy;
1490 /// Int_t zmin = 0;
1491 /// Int_t zmax = nbinsz;
1492 /// Double_t*** source = newDouble_t**[nbinsx];
1493 /// Double_t*** resp = new Double_t**[nbinsx];
1494 /// for(i=0;i<nbinsx;i++){
1495 /// source[i]=new Double_t* [nbinsy];
1496 /// for(j=0;j<nbinsy;j++)
1497 /// source[i][j]=new Double_t[nbinsz];
1498 /// }
1499 /// for(i=0;i<nbinsx;i++){
1500 /// resp[i]=new Double_t*[nbinsy];
1501 /// for(j=0;j<nbinsy;j++)
1502 /// resp[i][j]=new Double_t[nbinsz];
1503 /// }
1504 /// TH3F *decon_in = new TH3F("decon_in","Deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1505 /// TH3F *decon_resp = new TH3F("decon_resp","Deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1506 /// TFile *f = new TFile("TSpectrum3.root");
1507 /// decon_in=(TH3F*) f->Get("decon_in;1");
1508 /// decon_resp=(TH3F*) f->Get("decon_resp;1");
1509 /// TCanvas *Deconvolution = new TCanvas("Deconvolution","Deconvolution of 3-dimensional spectra",10,10,1000,700);
1510 /// TSpectrum3 *s = new TSpectrum3();
1511 /// for (i = 0; i < nbinsx; i++){
1512 /// for (j = 0; j < nbinsy; j++){
1513 /// for (k = 0; k < nbinsz; k++){
1514 /// source[i][j][k] = decon_in->GetBinContent(i + 1,j + 1,k + 1);
1515 /// resp[i][j][k] = decon_resp->GetBinContent(i + 1,j + 1,k + 1);
1516 /// }
1517 /// }
1518 /// }
1519 /// s->Deconvolution(source,resp,nbinsx,nbinsy,nbinsz,100,1,1);
1520 /// for (i = 0; i < nbinsx; i++){
1521 /// for (j = 0; j < nbinsy; j++){
1522 /// for (k = 0; k < nbinsz; k++){
1523 /// decon_in->SetBinContent(i + 1,j + 1,k + 1, source[i][j][k]);
1524 /// }
1525 /// }
1526 /// }
1527 /// decon_in->Draw("");
1528 /// }
1529 /// ~~~
1530 ///
1531 /// ### Example 2 - script Decon_hr.c :
1532 ///
1533 /// This example illustrates repeated
1534 /// Gold deconvolution with boosting. After every 10 iterations we apply power
1535 /// function with exponent = 2 to the spectrum given in Fig. 2.
1536 ///
1537 /// \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.
1538 ///
1539 /// #### Script:
1540 ///
1541 /// Example to illustrate the Gold deconvolution (class TSpectrum3).
1542 /// To execute this example, do:
1543 ///
1544 /// `root > .x Decon3_hr.C`
1545 ///
1546 /// ~~~ {.cpp}
1547 /// void Decon3_hr() {
1548 /// Int_t i, j, k;
1549 /// Int_t nbinsx = 32;
1550 /// Int_t nbinsy = 32;
1551 /// Int_t nbinsz = 32;
1552 /// Int_t xmin = 0;
1553 /// Int_t xmax = nbinsx;
1554 /// Int_t ymin = 0;
1555 /// Int_t ymax = nbinsy;
1556 /// Int_t zmin = 0;
1557 /// Int_t zmax = nbinsz;
1558 /// Double_t*** source = new Double_t**[nbinsx];
1559 /// Double_t*** resp = new Double_t**[nbinsx];
1560 /// for(i=0;i<nbinsx;i++){
1561 /// source[i]=new Double_t*[nbinsy];
1562 /// for(j=0;j<nbinsy;j++)
1563 /// source[i][j]=new Double_t[nbinsz];
1564 /// }
1565 /// for(i=0;i<nbinsx;i++){
1566 /// resp[i]=new Double_t*[nbinsy];
1567 /// for(j=0;j<nbinsy;j++)
1568 /// resp[i][j]=new Double_t[nbinsz];
1569 /// }
1570 /// TH3F *decon_in = new TH3F("decon_in","Deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1571 /// TH3F *decon_resp = new TH3F("decon_resp","Deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1572 /// TFile *f = new TFile("TSpectrum3.root");
1573 /// decon_in=(TH3F*)f->Get("decon_in;1");
1574 /// decon_resp=(TH3F*)f->Get("decon_resp;1");
1575 /// TCanvas *Deconvolution = new TCanvas("Deconvolution","High resolution deconvolution of 3-dimensional spectra",10,10,1000,700);
1576 /// TSpectrum3 *s = new TSpectrum3();
1577 /// for (i = 0; i < nbinsx; i++){
1578 /// for (j = 0; j < nbinsy; j++){
1579 /// for (k = 0; k < nbinsz; k++){
1580 /// source[i][j][k] = decon_in->GetBinContent(i + 1,j + 1,k + 1);
1581 /// resp[i][j][k] = decon_resp->GetBinContent(i + 1,j + 1,k + 1);
1582 /// }
1583 /// }
1584 /// }
1585 /// s->Deconvolution(source,resp,nbinsx,nbinsy,nbinsz,10,10,2);
1586 /// for (i = 0; i < nbinsx; i++){
1587 /// for (j = 0; j < nbinsy; j++){
1588 /// for (k = 0; k < nbinsz; k++){
1589 /// decon_in->SetBinContent(i + 1,j + 1,k + 1, source[i][j][k]);
1590 /// }
1591 /// }
1592 /// }
1593 /// decon_in->Draw("");
1594 /// }
1595 /// ~~~
1596 
1597 const char *TSpectrum3::Deconvolution(Double_t***source, const Double_t***resp,
1598  Int_t ssizex, Int_t ssizey, Int_t ssizez,
1599  Int_t numberIterations,
1600  Int_t numberRepetitions,
1601  Double_t boost)
1602 {
1603  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;
1604  Double_t lda, ldb, ldc, area, maximum = 0;
1605  if (ssizex <= 0 || ssizey <= 0 || ssizez <= 0)
1606  return "Wrong parameters";
1607  if (numberIterations <= 0)
1608  return "Number of iterations must be positive";
1609  if (numberRepetitions <= 0)
1610  return "Number of repetitions must be positive";
1611  Double_t ***working_space=new Double_t** [ssizex];
1612  for(i=0;i<ssizex;i++){
1613  working_space[i]=new Double_t* [ssizey];
1614  for(j=0;j<ssizey;j++)
1615  working_space[i][j]=new Double_t [5*ssizez];
1616  }
1617  area = 0;
1618  lhx = -1, lhy = -1, lhz = -1;
1619  for (i = 0; i < ssizex; i++) {
1620  for (j = 0; j < ssizey; j++) {
1621  for (k = 0; k < ssizez; k++) {
1622  lda = resp[i][j][k];
1623  if (lda != 0) {
1624  if ((i + 1) > lhx)
1625  lhx = i + 1;
1626  if ((j + 1) > lhy)
1627  lhy = j + 1;
1628  if ((k + 1) > lhz)
1629  lhz = k + 1;
1630  }
1631  working_space[i][j][k] = lda;
1632  area = area + lda;
1633  if (lda > maximum) {
1634  maximum = lda;
1635  positx = i, posity = j, positz = k;
1636  }
1637  }
1638  }
1639  }
1640  if (lhx == -1 || lhy == -1 || lhz == -1) {
1641  for(i = 0;i < ssizex; i++){
1642  for(j = 0;j < ssizey; j++)
1643  delete[] working_space[i][j];
1644  delete[] working_space[i];
1645  }
1646  delete [] working_space;
1647  return ("Zero response data");
1648  }
1649 
1650 //calculate ht*y and write into p
1651  for (i3 = 0; i3 < ssizez; i3++) {
1652  for (i2 = 0; i2 < ssizey; i2++) {
1653  for (i1 = 0; i1 < ssizex; i1++) {
1654  ldc = 0;
1655  for (j3 = 0; j3 <= (lhz - 1); j3++) {
1656  for (j2 = 0; j2 <= (lhy - 1); j2++) {
1657  for (j1 = 0; j1 <= (lhx - 1); j1++) {
1658  k3 = i3 + j3, k2 = i2 + j2, k1 = i1 + j1;
1659  if (k3 >= 0 && k3 < ssizez && k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
1660  lda = working_space[j1][j2][j3];
1661  ldb = source[k1][k2][k3];
1662  ldc = ldc + lda * ldb;
1663  }
1664  }
1665  }
1666  }
1667  working_space[i1][i2][i3 + ssizez] = ldc;
1668  }
1669  }
1670  }
1671 
1672 //calculate matrix b=ht*h
1673  i1min = -(lhx - 1), i1max = lhx - 1;
1674  i2min = -(lhy - 1), i2max = lhy - 1;
1675  i3min = -(lhz - 1), i3max = lhz - 1;
1676  for (i3 = i3min; i3 <= i3max; i3++) {
1677  for (i2 = i2min; i2 <= i2max; i2++) {
1678  for (i1 = i1min; i1 <= i1max; i1++) {
1679  ldc = 0;
1680  j3min = -i3;
1681  if (j3min < 0)
1682  j3min = 0;
1683  j3max = lhz - 1 - i3;
1684  if (j3max > lhz - 1)
1685  j3max = lhz - 1;
1686  for (j3 = j3min; j3 <= j3max; j3++) {
1687  j2min = -i2;
1688  if (j2min < 0)
1689  j2min = 0;
1690  j2max = lhy - 1 - i2;
1691  if (j2max > lhy - 1)
1692  j2max = lhy - 1;
1693  for (j2 = j2min; j2 <= j2max; j2++) {
1694  j1min = -i1;
1695  if (j1min < 0)
1696  j1min = 0;
1697  j1max = lhx - 1 - i1;
1698  if (j1max > lhx - 1)
1699  j1max = lhx - 1;
1700  for (j1 = j1min; j1 <= j1max; j1++) {
1701  lda = working_space[j1][j2][j3];
1702  if (i1 + j1 < ssizex && i2 + j2 < ssizey)
1703  ldb = working_space[i1 + j1][i2 + j2][i3 + j3];
1704  else
1705  ldb = 0;
1706  ldc = ldc + lda * ldb;
1707  }
1708  }
1709  }
1710  working_space[i1 - i1min][i2 - i2min][i3 - i3min + 2 * ssizez ] = ldc;
1711  }
1712  }
1713  }
1714 
1715 //initialization in x1 matrix
1716  for (i3 = 0; i3 < ssizez; i3++) {
1717  for (i2 = 0; i2 < ssizey; i2++) {
1718  for (i1 = 0; i1 < ssizex; i1++) {
1719  working_space[i1][i2][i3 + 3 * ssizez] = 1;
1720  working_space[i1][i2][i3 + 4 * ssizez] = 0;
1721  }
1722  }
1723  }
1724 
1725  //START OF ITERATIONS
1726  for (repet = 0; repet < numberRepetitions; repet++) {
1727  if (repet != 0) {
1728  for (i = 0; i < ssizex; i++) {
1729  for (j = 0; j < ssizey; j++) {
1730  for (k = 0; k < ssizez; k++) {
1731  working_space[i][j][k + 3 * ssizez] = TMath::Power(working_space[i][j][k + 3 * ssizez],boost);
1732  }
1733  }
1734  }
1735  }
1736  for (lindex = 0; lindex < numberIterations; lindex++) {
1737  for (i3 = 0; i3 < ssizez; i3++) {
1738  for (i2 = 0; i2 < ssizey; i2++) {
1739  for (i1 = 0; i1 < ssizex; i1++) {
1740  ldb = 0;
1741  j3min = i3;
1742  if (j3min > lhz - 1)
1743  j3min = lhz - 1;
1744  j3min = -j3min;
1745  j3max = ssizez - i3 - 1;
1746  if (j3max > lhz - 1)
1747  j3max = lhz - 1;
1748  j2min = i2;
1749  if (j2min > lhy - 1)
1750  j2min = lhy - 1;
1751  j2min = -j2min;
1752  j2max = ssizey - i2 - 1;
1753  if (j2max > lhy - 1)
1754  j2max = lhy - 1;
1755  j1min = i1;
1756  if (j1min > lhx - 1)
1757  j1min = lhx - 1;
1758  j1min = -j1min;
1759  j1max = ssizex - i1 - 1;
1760  if (j1max > lhx - 1)
1761  j1max = lhx - 1;
1762  for (j3 = j3min; j3 <= j3max; j3++) {
1763  for (j2 = j2min; j2 <= j2max; j2++) {
1764  for (j1 = j1min; j1 <= j1max; j1++) {
1765  ldc = working_space[j1 - i1min][j2 - i2min][j3 - i3min + 2 * ssizez];
1766  lda = working_space[i1 + j1][i2 + j2][i3 + j3 + 3 * ssizez];
1767  ldb = ldb + lda * ldc;
1768  }
1769  }
1770  }
1771  lda = working_space[i1][i2][i3 + 3 * ssizez];
1772  ldc = working_space[i1][i2][i3 + 1 * ssizez];
1773  if (ldc * lda != 0 && ldb != 0) {
1774  lda = lda * ldc / ldb;
1775  }
1776 
1777  else
1778  lda = 0;
1779  working_space[i1][i2][i3 + 4 * ssizez] = lda;
1780  }
1781  }
1782  }
1783  for (i3 = 0; i3 < ssizez; i3++) {
1784  for (i2 = 0; i2 < ssizey; i2++) {
1785  for (i1 = 0; i1 < ssizex; i1++)
1786  working_space[i1][i2][i3 + 3 * ssizez] = working_space[i1][i2][i3 + 4 * ssizez];
1787  }
1788  }
1789  }
1790  }
1791  for (i = 0; i < ssizex; i++) {
1792  for (j = 0; j < ssizey; j++){
1793  for (k = 0; k < ssizez; k++)
1794  source[(i + positx) % ssizex][(j + posity) % ssizey][(k + positz) % ssizez] = area * working_space[i][j][k + 3 * ssizez];
1795  }
1796  }
1797  for(i = 0;i < ssizex; i++){
1798  for(j = 0;j < ssizey; j++)
1799  delete[] working_space[i][j];
1800  delete[] working_space[i];
1801  }
1802  delete [] working_space;
1803  return 0;
1804 }
1805 
1806 ////////////////////////////////////////////////////////////////////////////////
1807 /// This function searches for peaks in source spectrum
1808 /// It is based on deconvolution method. First the background is
1809 /// removed (if desired), then Markov spectrum is calculated
1810 /// (if desired), then the response function is generated
1811 /// according to given sigma and deconvolution is carried out.
1812 /// It returns number of found peaks.
1813 ///
1814 /// Function parameters:
1815 /// - source-pointer to the matrix of source spectrum
1816 /// - dest-pointer to the matrix of resulting deconvolved spectrum
1817 /// - ssizex-x length of source spectrum
1818 /// - ssizey-y length of source spectrum
1819 /// - ssizez-z length of source spectrum
1820 /// - sigma-sigma of searched peaks, for details we refer to manual
1821 /// - threshold-threshold value in % for selected peaks, peaks with
1822 /// amplitude less than threshold*highest_peak/100
1823 /// are ignored, see manual
1824 /// - backgroundRemove-logical variable, set if the removal of
1825 /// background before deconvolution is desired
1826 /// - deconIterations-number of iterations in deconvolution operation
1827 /// - markov-logical variable, if it is true, first the source spectrum
1828 /// is replaced by new spectrum calculated using Markov
1829 /// chains method.
1830 /// - averWindow-averaging window of searched peaks, for details
1831 /// we refer to manual (applies only for Markov method)
1832 ///
1833 /// ### Peaks searching
1834 ///
1835 /// Goal: to identify automatically the peaks in spectrum with the presence of
1836 /// the continuous background, one- and two-fold coincidences (ridges) and statistical
1837 /// fluctuations - noise.
1838 ///
1839 /// The common problems connected
1840 /// with correct peak identification in three-dimensional coincidence spectra are
1841 ///
1842 /// - non-sensitivity to noise, i.e.,
1843 /// only statistically relevant peaks should be identified
1844 /// - non-sensitivity of the
1845 /// algorithm to continuous background
1846 /// - non-sensitivity to one-fold coincidences
1847 /// (coincidences peak - peak - background in all dimensions) and their
1848 /// crossings
1849 /// - non-sensitivity to two-fold
1850 /// coincidences (coincidences peak - background - background in all
1851 /// dimensions) and their crossings
1852 /// - ability to identify peaks close
1853 /// to the edges of the spectrum region
1854 /// - resolution, decomposition of
1855 /// doublets and multiplets. The algorithm should be able to recognise close
1856 /// positioned peaks.
1857 ///
1858 /// #### References:
1859 ///
1860 /// [1] M.A. Mariscotti: A method for
1861 /// identification of peaks in the presence of background and its application to
1862 /// spectrum analysis. NIM 50 (1967), 309-320.
1863 ///
1864 /// [2] M.Morhac, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:Identification
1865 /// of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
1866 /// 108-125.
1867 ///
1868 /// [3] Z.K. Silagadze, A new algorithm for automatic photo-peak searches. NIM A 376 (1996), 451.
1869 ///
1870 /// ### Example of peak searching method
1871 ///
1872 /// SearchHighRes function provides users with the possibility
1873 /// to vary the input parameters and with the access to the output deconvolved data
1874 /// in the destination spectrum. Based on the output data one can tune the
1875 /// parameters.
1876 ///
1877 /// #### Example 1 - script Search3.c:
1878 ///
1879 /// \image html spectrum3_searching_image001.jpg Fig. 1 Three-dimensional spectrum with 5 peaks (sigma=2, threshold=5%, 3 iterations steps in the deconvolution)
1880 /// \image html spectrum3_searching_image003.jpg Fig. 2 Spectrum from Fig. 1 after background elimination and deconvolution
1881 ///
1882 /// #### Script:
1883 ///
1884 /// Example to illustrate high resolution peak searching function (class TSpectrum3).
1885 /// To execute this example, do:
1886 ///
1887 /// `root > .x Search3.C`
1888 ///
1889 /// ~~~ {.cpp}
1890 /// void Search3() {
1891 /// Int_t i, j, k, nfound;
1892 /// Int_t nbinsx = 32;
1893 /// Int_t nbinsy = 32;
1894 /// Int_t nbinsz = 32;
1895 /// Int_t xmin = 0;
1896 /// Int_t xmax = nbinsx;
1897 /// Int_t ymin = 0;
1898 /// Int_t ymax = nbinsy;
1899 /// Int_t zmin = 0;
1900 /// Int_t zmax = nbinsz;
1901 /// Double_t*** source = new Double_t**[nbinsx];
1902 /// Double_t*** dest = new Double_t**[nbinsx];
1903 /// for(i=0;i<nbinsx;i++){
1904 /// source[i]=new Double_t*[nbinsy];
1905 /// for(j=0;j<nbinsy;j++)
1906 /// source[i][j]=new Double_t[nbinsz];
1907 /// }
1908 /// for(i=0;i<nbinsx;i++){
1909 /// dest[i]=new Double_t*[nbinsy];
1910 /// for(j=0;j<nbinsy;j++)
1911 /// dest[i][j]=new Double_t [nbinsz];
1912 /// }
1913 /// TH3F *search = new TH3F("Search","Peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax,nbinsz,zmin,zmax);
1914 /// TFile *f = new TFile("TSpectrum3.root");
1915 /// search=(TH3F*)f->Get("search2;1");
1916 /// TCanvas *Search = new TCanvas("Search","Peak searching",10,10,1000,700);
1917 /// TSpectrum3 *s = new TSpectrum3();
1918 /// for (i = 0; i < nbinsx; i++){
1919 /// for (j = 0; j < nbinsy; j++){
1920 /// for (k = 0; k < nbinsz; k++){
1921 /// source[i][j][k] = search->GetBinContent(i + 1,j + 1,k + 1);
1922 /// }
1923 /// }
1924 /// }
1925 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, nbinsz, 2, 5, kTRUE, 3, kFALSE, 3);
1926 /// printf("Found %d candidate peaks\n",nfound);
1927 /// for (i = 0; i < nbinsx; i++){
1928 /// for (j = 0; j < nbinsy; j++){
1929 /// for (k = 0; k < nbinsz; k++){
1930 /// search->SetBinContent(i + 1,j + 1,k + 1, dest[i][j][k]);
1931 /// }
1932 /// }
1933 /// }
1934 /// Double_t *PosX = new Double_t[nfound];
1935 /// Double_t *PosY = new Double_t[nfound];
1936 /// Double_t *PosZ = new Double_t[nfound];
1937 /// PosX = s->GetPositionX();
1938 /// PosY = s->GetPositionY();
1939 /// PosZ = s->GetPositionZ();
1940 /// for(i=0;i<nfound;i++)
1941 /// 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));
1942 /// search->Draw("");
1943 /// }
1944 /// ~~~
1945 
1946 Int_t TSpectrum3::SearchHighRes(const Double_t***source,Double_t***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez,
1947  Double_t sigma, Double_t threshold,
1948  Bool_t backgroundRemove,Int_t deconIterations,
1949  Bool_t markov, Int_t averWindow)
1950 
1951 {
1952  Int_t number_of_iterations = (Int_t)(4 * sigma + 0.5);
1953  Int_t k,lindex;
1954  Double_t lda,ldb,ldc,area,maximum;
1955  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;
1956  Int_t ymin,ymax,zmin,zmax,i,j;
1957  Double_t a,b,maxch,plocha = 0,plocha_markov = 0;
1958  Double_t nom,nip,nim,sp,sm,spx,spy,smx,smy,spz,smz;
1959  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;
1960  Int_t x,y,z;
1961  Double_t pocet_sigma = 5;
1962  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;
1963  if(sigma < 1){
1964  Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
1965  return 0;
1966  }
1967 
1968  if(threshold<=0||threshold>=100){
1969  Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
1970  return 0;
1971  }
1972 
1973  j = (Int_t)(pocet_sigma*sigma+0.5);
1974  if (j >= PEAK_WINDOW / 2) {
1975  Error("SearchHighRes", "Too large sigma");
1976  return 0;
1977  }
1978 
1979  if (markov == true) {
1980  if (averWindow <= 0) {
1981  Error("SearchHighRes", "Averaging window must be positive");
1982  return 0;
1983  }
1984  }
1985 
1986  if(backgroundRemove == true){
1987  if(sizex_ext < 2 * number_of_iterations + 1 || sizey_ext < 2 * number_of_iterations + 1 || sizez_ext < 2 * number_of_iterations + 1){
1988  Error("SearchHighRes", "Too large clipping window");
1989  return 0;
1990  }
1991  }
1992 
1993  i = (Int_t)(4 * sigma + 0.5);
1994  i = 4 * i;
1995  Double_t ***working_space = new Double_t** [ssizex + i];
1996  for(j = 0;j < ssizex + i; j++){
1997  working_space[j] = new Double_t* [ssizey + i];
1998  for(k = 0;k < ssizey + i; k++)
1999  working_space[j][k] = new Double_t [5 * (ssizez + i)];
2000  }
2001  for(k = 0;k < sizez_ext; k++){
2002  for(j = 0;j < sizey_ext; j++){
2003  for(i = 0;i < sizex_ext; i++){
2004  if(i < shift){
2005  if(j < shift){
2006  if(k < shift)
2007  working_space[i][j][k + sizez_ext] = source[0][0][0];
2008 
2009  else if(k >= ssizez + shift)
2010  working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1];
2011 
2012  else
2013  working_space[i][j][k + sizez_ext] = source[0][0][k - shift];
2014  }
2015 
2016  else if(j >= ssizey + shift){
2017  if(k < shift)
2018  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0];
2019 
2020  else if(k >= ssizez + shift)
2021  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1];
2022 
2023  else
2024  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift];
2025  }
2026 
2027  else{
2028  if(k < shift)
2029  working_space[i][j][k + sizez_ext] = source[0][j - shift][0];
2030 
2031  else if(k >= ssizez + shift)
2032  working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1];
2033 
2034  else
2035  working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift];
2036  }
2037  }
2038 
2039  else if(i >= ssizex + shift){
2040  if(j < shift){
2041  if(k < shift)
2042  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0];
2043 
2044  else if(k >= ssizez + shift)
2045  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1];
2046 
2047  else
2048  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift];
2049  }
2050 
2051  else if(j >= ssizey + shift){
2052  if(k < shift)
2053  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0];
2054 
2055  else if(k >= ssizez + shift)
2056  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1];
2057 
2058  else
2059  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift];
2060  }
2061 
2062  else{
2063  if(k < shift)
2064  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0];
2065 
2066  else if(k >= ssizez + shift)
2067  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1];
2068 
2069  else
2070  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift];
2071  }
2072  }
2073 
2074  else{
2075  if(j < shift){
2076  if(k < shift)
2077  working_space[i][j][k + sizez_ext] = source[i - shift][0][0];
2078 
2079  else if(k >= ssizez + shift)
2080  working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1];
2081 
2082  else
2083  working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift];
2084  }
2085 
2086  else if(j >= ssizey + shift){
2087  if(k < shift)
2088  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0];
2089 
2090  else if(k >= ssizez + shift)
2091  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1];
2092 
2093  else
2094  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift];
2095  }
2096 
2097  else{
2098  if(k < shift)
2099  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0];
2100 
2101  else if(k >= ssizez + shift)
2102  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1];
2103 
2104  else
2105  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift];
2106  }
2107  }
2108  }
2109  }
2110  }
2111  if(backgroundRemove == true){
2112  for(i = 1;i <= number_of_iterations; i++){
2113  for(z = i;z < sizez_ext - i; z++){
2114  for(y = i;y < sizey_ext - i; y++){
2115  for(x = i;x < sizex_ext - i; x++){
2116  a = working_space[x][y][z + sizez_ext];
2117  p1 = working_space[x + i][y + i][z - i + sizez_ext];
2118  p2 = working_space[x - i][y + i][z - i + sizez_ext];
2119  p3 = working_space[x + i][y - i][z - i + sizez_ext];
2120  p4 = working_space[x - i][y - i][z - i + sizez_ext];
2121  p5 = working_space[x + i][y + i][z + i + sizez_ext];
2122  p6 = working_space[x - i][y + i][z + i + sizez_ext];
2123  p7 = working_space[x + i][y - i][z + i + sizez_ext];
2124  p8 = working_space[x - i][y - i][z + i + sizez_ext];
2125  s1 = working_space[x + i][y ][z - i + sizez_ext];
2126  s2 = working_space[x ][y + i][z - i + sizez_ext];
2127  s3 = working_space[x - i][y ][z - i + sizez_ext];
2128  s4 = working_space[x ][y - i][z - i + sizez_ext];
2129  s5 = working_space[x + i][y ][z + i + sizez_ext];
2130  s6 = working_space[x ][y + i][z + i + sizez_ext];
2131  s7 = working_space[x - i][y ][z + i + sizez_ext];
2132  s8 = working_space[x ][y - i][z + i + sizez_ext];
2133  s9 = working_space[x - i][y + i][z + sizez_ext];
2134  s10 = working_space[x - i][y - i][z +sizez_ext];
2135  s11 = working_space[x + i][y + i][z +sizez_ext];
2136  s12 = working_space[x + i][y - i][z +sizez_ext];
2137  r1 = working_space[x ][y ][z - i + sizez_ext];
2138  r2 = working_space[x ][y ][z + i + sizez_ext];
2139  r3 = working_space[x - i][y ][z + sizez_ext];
2140  r4 = working_space[x + i][y ][z + sizez_ext];
2141  r5 = working_space[x ][y + i][z + sizez_ext];
2142  r6 = working_space[x ][y - i][z + sizez_ext];
2143  b = (p1 + p3) / 2.0;
2144  if(b > s1)
2145  s1 = b;
2146 
2147  b = (p1 + p2) / 2.0;
2148  if(b > s2)
2149  s2 = b;
2150 
2151  b = (p2 + p4) / 2.0;
2152  if(b > s3)
2153  s3 = b;
2154 
2155  b = (p3 + p4) / 2.0;
2156  if(b > s4)
2157  s4 = b;
2158 
2159  b = (p5 + p7) / 2.0;
2160  if(b > s5)
2161  s5 = b;
2162 
2163  b = (p5 + p6) / 2.0;
2164  if(b > s6)
2165  s6 = b;
2166 
2167  b = (p6 + p8) / 2.0;
2168  if(b > s7)
2169  s7 = b;
2170 
2171  b = (p7 + p8) / 2.0;
2172  if(b > s8)
2173  s8 = b;
2174 
2175  b = (p2 + p6) / 2.0;
2176  if(b > s9)
2177  s9 = b;
2178 
2179  b = (p4 + p8) / 2.0;
2180  if(b > s10)
2181  s10 = b;
2182 
2183  b = (p1 + p5) / 2.0;
2184  if(b > s11)
2185  s11 = b;
2186 
2187  b = (p3 + p7) / 2.0;
2188  if(b > s12)
2189  s12 = b;
2190 
2191  s1 = s1 - (p1 + p3) / 2.0;
2192  s2 = s2 - (p1 + p2) / 2.0;
2193  s3 = s3 - (p2 + p4) / 2.0;
2194  s4 = s4 - (p3 + p4) / 2.0;
2195  s5 = s5 - (p5 + p7) / 2.0;
2196  s6 = s6 - (p5 + p6) / 2.0;
2197  s7 = s7 - (p6 + p8) / 2.0;
2198  s8 = s8 - (p7 + p8) / 2.0;
2199  s9 = s9 - (p2 + p6) / 2.0;
2200  s10 = s10 - (p4 + p8) / 2.0;
2201  s11 = s11 - (p1 + p5) / 2.0;
2202  s12 = s12 - (p3 + p7) / 2.0;
2203  b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
2204  if(b > r1)
2205  r1 = b;
2206 
2207  b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
2208  if(b > r2)
2209  r2 = b;
2210 
2211  b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
2212  if(b > r3)
2213  r3 = b;
2214 
2215  b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
2216  if(b > r4)
2217  r4 = b;
2218 
2219  b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
2220  if(b > r5)
2221  r5 = b;
2222 
2223  b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
2224  if(b > r6)
2225  r6 = b;
2226 
2227  r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
2228  r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
2229  r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
2230  r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
2231  r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
2232  r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
2233  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;
2234  if(b < a)
2235  a = b;
2236 
2237  working_space[x][y][z] = a;
2238  }
2239  }
2240  }
2241  for(z = i;z < sizez_ext - i; z++){
2242  for(y = i;y < sizey_ext - i; y++){
2243  for(x = i;x < sizex_ext - i; x++){
2244  working_space[x][y][z + sizez_ext] = working_space[x][y][z];
2245  }
2246  }
2247  }
2248  }
2249  for(k = 0;k < sizez_ext; k++){
2250  for(j = 0;j < sizey_ext; j++){
2251  for(i = 0;i < sizex_ext; i++){
2252  if(i < shift){
2253  if(j < shift){
2254  if(k < shift)
2255  working_space[i][j][k + sizez_ext] = source[0][0][0] - working_space[i][j][k + sizez_ext];
2256 
2257  else if(k >= ssizez + shift)
2258  working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
2259 
2260  else
2261  working_space[i][j][k + sizez_ext] = source[0][0][k - shift] - working_space[i][j][k + sizez_ext];
2262  }
2263 
2264  else if(j >= ssizey + shift){
2265  if(k < shift)
2266  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
2267 
2268  else if(k >= ssizez + shift)
2269  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
2270 
2271  else
2272  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
2273  }
2274 
2275  else{
2276  if(k < shift)
2277  working_space[i][j][k + sizez_ext] = source[0][j - shift][0] - working_space[i][j][k + sizez_ext];
2278 
2279  else if(k >= ssizez + shift)
2280  working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
2281 
2282  else
2283  working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
2284  }
2285  }
2286 
2287  else if(i >= ssizex + shift){
2288  if(j < shift){
2289  if(k < shift)
2290  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0] - working_space[i][j][k + sizez_ext];
2291 
2292  else if(k >= ssizez + shift)
2293  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
2294 
2295  else
2296  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift] - working_space[i][j][k + sizez_ext];
2297  }
2298 
2299  else if(j >= ssizey + shift){
2300  if(k < shift)
2301  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
2302 
2303  else if(k >= ssizez + shift)
2304  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
2305 
2306  else
2307  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
2308  }
2309 
2310  else{
2311  if(k < shift)
2312  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0] - working_space[i][j][k + sizez_ext];
2313 
2314  else if(k >= ssizez + shift)
2315  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
2316 
2317  else
2318  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
2319  }
2320  }
2321 
2322  else{
2323  if(j < shift){
2324  if(k < shift)
2325  working_space[i][j][k + sizez_ext] = source[i - shift][0][0] - working_space[i][j][k + sizez_ext];
2326 
2327  else if(k >= ssizez + shift)
2328  working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
2329 
2330  else
2331  working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift] - working_space[i][j][k + sizez_ext];
2332  }
2333 
2334  else if(j >= ssizey + shift){
2335  if(k < shift)
2336  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
2337 
2338  else if(k >= ssizez + shift)
2339  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
2340 
2341  else
2342  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
2343  }
2344 
2345  else{
2346  if(k < shift)
2347  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0] - working_space[i][j][k + sizez_ext];
2348 
2349  else if(k >= ssizez + shift)
2350  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
2351 
2352  else
2353  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
2354  }
2355  }
2356  }
2357  }
2358  }
2359  }
2360 
2361  if(markov == true){
2362  for(i = 0;i < sizex_ext; i++){
2363  for(j = 0;j < sizey_ext; j++){
2364  for(k = 0;k < sizez_ext; k++){
2365  working_space[i][j][k + 2 * sizez_ext] = working_space[i][j][sizez_ext + k];
2366  plocha_markov = plocha_markov + working_space[i][j][sizez_ext + k];
2367  }
2368  }
2369  }
2370  xmin = 0;
2371  xmax = sizex_ext - 1;
2372  ymin = 0;
2373  ymax = sizey_ext - 1;
2374  zmin = 0;
2375  zmax = sizez_ext - 1;
2376  for(i = 0,maxch = 0;i < sizex_ext; i++){
2377  for(j = 0;j < sizey_ext;j++){
2378  for(k = 0;k < sizez_ext;k++){
2379  working_space[i][j][k] = 0;
2380  if(maxch < working_space[i][j][k + 2 * sizez_ext])
2381  maxch = working_space[i][j][k + 2 * sizez_ext];
2382 
2383  plocha += working_space[i][j][k + 2 * sizez_ext];
2384  }
2385  }
2386  }
2387  if(maxch == 0) {
2388  k = (Int_t)(4 * sigma + 0.5);
2389  k = 4 * k;
2390  for(i = 0;i < ssizex + k; i++){
2391  for(j = 0;j < ssizey + k; j++)
2392  delete[] working_space[i][j];
2393  delete[] working_space[i];
2394  }
2395  delete [] working_space;
2396  return 0;
2397  }
2398  nom = 0;
2399  working_space[xmin][ymin][zmin] = 1;
2400  for(i = xmin;i < xmax; i++){
2401  nip = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
2402  nim = working_space[i + 1][ymin][zmin + 2 * sizez_ext] / maxch;
2403  sp = 0,sm = 0;
2404  for(l = 1;l <= averWindow; l++){
2405  if((i + l) > xmax)
2406  a = working_space[xmax][ymin][zmin + 2 * sizez_ext] / maxch;
2407 
2408  else
2409  a = working_space[i + l][ymin][zmin + 2 * sizez_ext] / maxch;
2410 
2411  b = a - nip;
2412  if(a + nip <= 0)
2413  a = 1;
2414 
2415  else
2416  a = TMath::Sqrt(a + nip);
2417 
2418  b = b / a;
2419  b = TMath::Exp(b);
2420  sp = sp + b;
2421  if(i - l + 1 < xmin)
2422  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
2423 
2424  else
2425  a = working_space[i - l + 1][ymin][zmin + 2 * sizez_ext] / maxch;
2426 
2427  b = a - nim;
2428  if(a + nim <= 0)
2429  a = 1;
2430 
2431  else
2432  a = TMath::Sqrt(a + nim);
2433 
2434  b = b / a;
2435  b = TMath::Exp(b);
2436  sm = sm + b;
2437  }
2438  a = sp / sm;
2439  a = working_space[i + 1][ymin][zmin] = a * working_space[i][ymin][zmin];
2440  nom = nom + a;
2441  }
2442  for(i = ymin;i < ymax; i++){
2443  nip = working_space[xmin][i][zmin + 2 * sizez_ext] / maxch;
2444  nim = working_space[xmin][i + 1][zmin + 2 * sizez_ext] / maxch;
2445  sp = 0,sm = 0;
2446  for(l = 1;l <= averWindow; l++){
2447  if((i + l) > ymax)
2448  a = working_space[xmin][ymax][zmin + 2 * sizez_ext] / maxch;
2449 
2450  else
2451  a = working_space[xmin][i + l][zmin + 2 * sizez_ext] / maxch;
2452 
2453  b = a - nip;
2454  if(a + nip <= 0)
2455  a = 1;
2456 
2457  else
2458  a = TMath::Sqrt(a + nip);
2459 
2460  b = b / a;
2461  b = TMath::Exp(b);
2462  sp = sp + b;
2463  if(i - l + 1 < ymin)
2464  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
2465 
2466  else
2467  a = working_space[xmin][i - l + 1][zmin + 2 * sizez_ext] / maxch;
2468 
2469  b = a - nim;
2470  if(a + nim <= 0)
2471  a = 1;
2472 
2473  else
2474  a = TMath::Sqrt(a + nim);
2475 
2476  b = b / a;
2477  b = TMath::Exp(b);
2478  sm = sm + b;
2479  }
2480  a = sp / sm;
2481  a = working_space[xmin][i + 1][zmin] = a * working_space[xmin][i][zmin];
2482  nom = nom + a;
2483  }
2484  for(i = zmin;i < zmax;i++){
2485  nip = working_space[xmin][ymin][i + 2 * sizez_ext] / maxch;
2486  nim = working_space[xmin][ymin][i + 1 + 2 * sizez_ext] / maxch;
2487  sp = 0,sm = 0;
2488  for(l = 1;l <= averWindow;l++){
2489  if((i + l) > zmax)
2490  a = working_space[xmin][ymin][zmax + 2 * sizez_ext] / maxch;
2491 
2492  else
2493  a = working_space[xmin][ymin][i + l + 2 * sizez_ext] / maxch;
2494 
2495  b = a - nip;
2496  if(a + nip <= 0)
2497  a = 1;
2498 
2499  else
2500  a = TMath::Sqrt(a + nip);
2501 
2502  b = b / a;
2503  b = TMath::Exp(b);
2504  sp = sp + b;
2505  if(i - l + 1 < zmin)
2506  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
2507 
2508  else
2509  a = working_space[xmin][ymin][i - l + 1 + 2 * sizez_ext] / maxch;
2510 
2511  b = a - nim;
2512  if(a + nim <= 0)
2513  a = 1;
2514 
2515  else
2516  a = TMath::Sqrt(a + nim);
2517 
2518  b = b / a;
2519  b = TMath::Exp(b);
2520  sm = sm + b;
2521  }
2522  a = sp / sm;
2523  a = working_space[xmin][ymin][i + 1] = a * working_space[xmin][ymin][i];
2524  nom = nom + a;
2525  }
2526  for(i = xmin;i < xmax; i++){
2527  for(j = ymin;j < ymax; j++){
2528  nip = working_space[i][j + 1][zmin + 2 * sizez_ext] / maxch;
2529  nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
2530  spx = 0,smx = 0;
2531  for(l = 1;l <= averWindow; l++){
2532  if(i + l > xmax)
2533  a = working_space[xmax][j][zmin + 2 * sizez_ext] / maxch;
2534 
2535  else
2536  a = working_space[i + l][j][zmin + 2 * sizez_ext] / maxch;
2537 
2538  b = a - nip;
2539  if(a + nip <= 0)
2540  a = 1;
2541 
2542  else
2543  a = TMath::Sqrt(a + nip);
2544 
2545  b = b / a;
2546  b = TMath::Exp(b);
2547  spx = spx + b;
2548  if(i - l + 1 < xmin)
2549  a = working_space[xmin][j][zmin + 2 * sizez_ext] / maxch;
2550 
2551  else
2552  a = working_space[i - l + 1][j][zmin + 2 * sizez_ext] / maxch;
2553 
2554  b = a - nim;
2555  if(a + nim <= 0)
2556  a = 1;
2557 
2558  else
2559  a = TMath::Sqrt(a + nim);
2560 
2561  b = b / a;
2562  b = TMath::Exp(b);
2563  smx = smx + b;
2564  }
2565  spy = 0,smy = 0;
2566  nip = working_space[i + 1][j][zmin + 2 * sizez_ext] / maxch;
2567  nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
2568  for(l = 1;l <= averWindow; l++){
2569  if(j + l > ymax)
2570  a = working_space[i][ymax][zmin + 2 * sizez_ext] / maxch;
2571 
2572  else
2573  a = working_space[i][j + l][zmin + 2 * sizez_ext] / maxch;
2574 
2575  b = a - nip;
2576  if(a + nip <= 0)
2577  a = 1;
2578 
2579  else
2580  a = TMath::Sqrt(a + nip);
2581 
2582  b = b / a;
2583  b = TMath::Exp(b);
2584  spy = spy + b;
2585  if(j - l + 1 < ymin)
2586  a = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
2587 
2588  else
2589  a = working_space[i][j - l + 1][zmin + 2 * sizez_ext] / maxch;
2590 
2591  b = a - nim;
2592  if(a + nim <= 0)
2593  a = 1;
2594 
2595  else
2596  a = TMath::Sqrt(a + nim);
2597 
2598  b = b / a;
2599  b = TMath::Exp(b);
2600  smy = smy + b;
2601  }
2602  a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
2603  working_space[i + 1][j + 1][zmin] = a;
2604  nom = nom + a;
2605  }
2606  }
2607  for(i = xmin;i < xmax;i++){
2608  for(j = zmin;j < zmax;j++){
2609  nip = working_space[i][ymin][j + 1 + 2 * sizez_ext] / maxch;
2610  nim = working_space[i + 1][ymin][j + 1 + 2 * sizez_ext] / maxch;
2611  spx = 0,smx = 0;
2612  for(l = 1;l <= averWindow; l++){
2613  if(i + l > xmax)
2614  a = working_space[xmax][ymin][j + 2 * sizez_ext] / maxch;
2615 
2616  else
2617  a = working_space[i + l][ymin][j + 2 * sizez_ext] / maxch;
2618 
2619  b = a - nip;
2620  if(a + nip <= 0)
2621  a = 1;
2622 
2623  else
2624  a = TMath::Sqrt(a + nip);
2625 
2626  b = b / a;
2627  b = TMath::Exp(b);
2628  spx = spx + b;
2629  if(i - l + 1 < xmin)
2630  a = working_space[xmin][ymin][j + 2 * sizez_ext] / maxch;
2631 
2632  else
2633  a = working_space[i - l + 1][ymin][j + 2 * sizez_ext] / maxch;
2634 
2635  b = a - nim;
2636  if(a + nim <= 0)
2637  a = 1;
2638 
2639  else
2640  a = TMath::Sqrt(a + nim);
2641 
2642  b = b / a;
2643  b = TMath::Exp(b);
2644  smx = smx + b;
2645  }
2646  spz = 0,smz = 0;
2647  nip = working_space[i + 1][ymin][j + 2 * sizez_ext] / maxch;
2648  nim = working_space[i + 1][ymin][j + 1 + 2 * sizez_ext] / maxch;
2649  for(l = 1;l <= averWindow; l++){
2650  if(j + l > zmax)
2651  a = working_space[i][ymin][zmax + 2 * sizez_ext] / maxch;
2652 
2653  else
2654  a = working_space[i][ymin][j + l + 2 * sizez_ext] / maxch;
2655 
2656  b = a - nip;
2657  if(a + nip <= 0)
2658  a = 1;
2659 
2660  else
2661  a = TMath::Sqrt(a + nip);
2662 
2663  b = b / a;
2664  b = TMath::Exp(b);
2665  spz = spz + b;
2666  if(j - l + 1 < zmin)
2667  a = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
2668 
2669  else
2670  a = working_space[i][ymin][j - l + 1 + 2 * sizez_ext] / maxch;
2671 
2672  b = a - nim;
2673  if(a + nim <= 0)
2674  a = 1;
2675 
2676  else
2677  a = TMath::Sqrt(a + nim);
2678 
2679  b = b / a;
2680  b = TMath::Exp(b);
2681  smz = smz + b;
2682  }
2683  a = (spx * working_space[i][ymin][j + 1] + spz * working_space[i + 1][ymin][j]) / (smx + smz);
2684  working_space[i + 1][ymin][j + 1] = a;
2685  nom = nom + a;
2686  }
2687  }
2688  for(i = ymin;i < ymax;i++){
2689  for(j = zmin;j < zmax;j++){
2690  nip = working_space[xmin][i][j + 1 + 2 * sizez_ext] / maxch;
2691  nim = working_space[xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
2692  spy = 0,smy = 0;
2693  for(l = 1;l <= averWindow; l++){
2694  if(i + l > ymax)
2695  a = working_space[xmin][ymax][j + 2 * sizez_ext] / maxch;
2696 
2697  else
2698  a = working_space[xmin][i + l][j + 2 * sizez_ext] / maxch;
2699 
2700  b = a - nip;
2701  if(a + nip <= 0)
2702  a = 1;
2703 
2704  else
2705  a = TMath::Sqrt(a + nip);
2706 
2707  b = b / a;
2708  b = TMath::Exp(b);
2709  spy = spy + b;
2710  if(i - l + 1 < ymin)
2711  a = working_space[xmin][ymin][j + 2 * sizez_ext] / maxch;
2712 
2713  else
2714  a = working_space[xmin][i - l + 1][j + 2 * sizez_ext] / maxch;
2715 
2716  b = a - nim;
2717  if(a + nim <= 0)
2718  a = 1;
2719 
2720  else
2721  a = TMath::Sqrt(a + nim);
2722 
2723  b = b / a;
2724  b = TMath::Exp(b);
2725  smy = smy + b;
2726  }
2727  spz = 0,smz = 0;
2728  nip = working_space[xmin][i + 1][j + 2 * sizez_ext] / maxch;
2729  nim = working_space[xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
2730  for(l = 1;l <= averWindow; l++){
2731  if(j + l > zmax)
2732  a = working_space[xmin][i][zmax + 2 * sizez_ext] / maxch;
2733 
2734  else
2735  a = working_space[xmin][i][j + l + 2 * sizez_ext] / maxch;
2736 
2737  b = a - nip;
2738  if(a + nip <= 0)
2739  a = 1;
2740 
2741  else
2742  a = TMath::Sqrt(a + nip);
2743 
2744  b = b / a;
2745  b = TMath::Exp(b);
2746  spz = spz + b;
2747  if(j - l + 1 < zmin)
2748  a = working_space[xmin][i][zmin + 2 * sizez_ext] / maxch;
2749 
2750  else
2751  a = working_space[xmin][i][j - l + 1 + 2 * sizez_ext] / maxch;
2752 
2753  b = a - nim;
2754  if(a + nim <= 0)
2755  a = 1;
2756 
2757  else
2758  a = TMath::Sqrt(a + nim);
2759 
2760  b = b / a;
2761  b = TMath::Exp(b);
2762  smz = smz + b;
2763  }
2764  a = (spy * working_space[xmin][i][j + 1] + spz * working_space[xmin][i + 1][j]) / (smy + smz);
2765  working_space[xmin][i + 1][j + 1] = a;
2766  nom = nom + a;
2767  }
2768  }
2769  for(i = xmin;i < xmax; i++){
2770  for(j = ymin;j < ymax; j++){
2771  for(k = zmin;k < zmax; k++){
2772  nip = working_space[i][j + 1][k + 1 + 2 * sizez_ext] / maxch;
2773  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
2774  spx = 0,smx = 0;
2775  for(l = 1;l <= averWindow; l++){
2776  if(i + l > xmax)
2777  a = working_space[xmax][j][k + 2 * sizez_ext] / maxch;
2778 
2779  else
2780  a = working_space[i + l][j][k + 2 * sizez_ext] / maxch;
2781 
2782  b = a - nip;
2783  if(a + nip <= 0)
2784  a = 1;
2785 
2786  else
2787  a = TMath::Sqrt(a + nip);
2788 
2789  b = b / a;
2790  b = TMath::Exp(b);
2791  spx = spx + b;
2792  if(i - l + 1 < xmin)
2793  a = working_space[xmin][j][k + 2 * sizez_ext] / maxch;
2794 
2795  else
2796  a = working_space[i - l + 1][j][k + 2 * sizez_ext] / maxch;
2797 
2798  b = a - nim;
2799  if(a + nim <= 0)
2800  a = 1;
2801 
2802  else
2803  a = TMath::Sqrt(a + nim);
2804 
2805  b = b / a;
2806  b = TMath::Exp(b);
2807  smx = smx + b;
2808  }
2809  spy = 0,smy = 0;
2810  nip = working_space[i + 1][j][k + 1 + 2 * sizez_ext] / maxch;
2811  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
2812  for(l = 1;l <= averWindow; l++){
2813  if(j + l > ymax)
2814  a = working_space[i][ymax][k + 2 * sizez_ext] / maxch;
2815 
2816  else
2817  a = working_space[i][j + l][k + 2 * sizez_ext] / maxch;
2818 
2819  b = a - nip;
2820  if(a + nip <= 0)
2821  a = 1;
2822 
2823  else
2824  a = TMath::Sqrt(a + nip);
2825 
2826  b = b / a;
2827  b = TMath::Exp(b);
2828  spy = spy + b;
2829  if(j - l + 1 < ymin)
2830  a = working_space[i][ymin][k + 2 * sizez_ext] / maxch;
2831 
2832  else
2833  a = working_space[i][j - l + 1][k + 2 * sizez_ext] / maxch;
2834 
2835  b = a - nim;
2836  if(a + nim <= 0)
2837  a = 1;
2838 
2839  else
2840  a = TMath::Sqrt(a + nim);
2841 
2842  b = b / a;
2843  b = TMath::Exp(b);
2844  smy = smy + b;
2845  }
2846  spz = 0,smz = 0;
2847  nip = working_space[i + 1][j + 1][k + 2 * sizez_ext] / maxch;
2848  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
2849  for(l = 1;l <= averWindow; l++ ){
2850  if(j + l > zmax)
2851  a = working_space[i][j][zmax + 2 * sizez_ext] / maxch;
2852 
2853  else
2854  a = working_space[i][j][k + l + 2 * sizez_ext] / maxch;
2855 
2856  b = a - nip;
2857  if(a + nip <= 0)
2858  a = 1;
2859 
2860  else
2861  a = TMath::Sqrt(a + nip);
2862 
2863  b = b / a;
2864  b = TMath::Exp(b);
2865  spz = spz + b;
2866  if(j - l + 1 < ymin)
2867  a = working_space[i][j][zmin + 2 * sizez_ext] / maxch;
2868 
2869  else
2870  a = working_space[i][j][k - l + 1 + 2 * sizez_ext] / maxch;
2871 
2872  b = a - nim;
2873  if(a + nim <= 0)
2874  a = 1;
2875 
2876  else
2877  a = TMath::Sqrt(a + nim);
2878 
2879  b = b / a;
2880  b = TMath::Exp(b);
2881  smz = smz + b;
2882  }
2883  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);
2884  working_space[i + 1][j + 1][k + 1] = a;
2885  nom = nom + a;
2886  }
2887  }
2888  }
2889  a = 0;
2890  for(i = xmin;i <= xmax; i++){
2891  for(j = ymin;j <= ymax; j++){
2892  for(k = zmin;k <= zmax; k++){
2893  working_space[i][j][k] = working_space[i][j][k] / nom;
2894  a+=working_space[i][j][k];
2895  }
2896  }
2897  }
2898  for(i = 0;i < sizex_ext; i++){
2899  for(j = 0;j < sizey_ext; j++){
2900  for(k = 0;k < sizez_ext; k++){
2901  working_space[i][j][k + sizez_ext] = working_space[i][j][k] * plocha_markov / a;
2902  }
2903  }
2904  }
2905  }
2906  //deconvolution starts
2907  area = 0;
2908  lhx = -1,lhy = -1,lhz = -1;
2909  positx = 0,posity = 0,positz = 0;
2910  maximum = 0;
2911  //generate response cube
2912  for(i = 0;i < sizex_ext; i++){
2913  for(j = 0;j < sizey_ext; j++){
2914  for(k = 0;k < sizez_ext; k++){
2915  lda = (Double_t)i - 3 * sigma;
2916  ldb = (Double_t)j - 3 * sigma;
2917  ldc = (Double_t)k - 3 * sigma;
2918  lda = (lda * lda + ldb * ldb + ldc * ldc) / (2 * sigma * sigma);
2919  l = (Int_t)(1000 * exp(-lda));
2920  lda = l;
2921  if(lda!=0){
2922  if((i + 1) > lhx)
2923  lhx = i + 1;
2924 
2925  if((j + 1) > lhy)
2926  lhy = j + 1;
2927 
2928  if((k + 1) > lhz)
2929  lhz = k + 1;
2930  }
2931  working_space[i][j][k] = lda;
2932  area = area + lda;
2933  if(lda > maximum){
2934  maximum = lda;
2935  positx = i,posity = j,positz = k;
2936  }
2937  }
2938  }
2939  }
2940  //read source cube
2941  for(i = 0;i < sizex_ext; i++){
2942  for(j = 0;j < sizey_ext; j++){
2943  for(k = 0;k < sizez_ext; k++){
2944  working_space[i][j][k + 2 * sizez_ext] = TMath::Abs(working_space[i][j][k + sizez_ext]);
2945  }
2946  }
2947  }
2948  //calculate ht*y and write into p
2949  for (i3 = 0; i3 < sizez_ext; i3++) {
2950  for (i2 = 0; i2 < sizey_ext; i2++) {
2951  for (i1 = 0; i1 < sizex_ext; i1++) {
2952  ldc = 0;
2953  for (j3 = 0; j3 <= (lhz - 1); j3++) {
2954  for (j2 = 0; j2 <= (lhy - 1); j2++) {
2955  for (j1 = 0; j1 <= (lhx - 1); j1++) {
2956  k3 = i3 + j3, k2 = i2 + j2, k1 = i1 + j1;
2957  if (k3 >= 0 && k3 < sizez_ext && k2 >= 0 && k2 < sizey_ext && k1 >= 0 && k1 < sizex_ext) {
2958  lda = working_space[j1][j2][j3];
2959  ldb = working_space[k1][k2][k3+2*sizez_ext];
2960  ldc = ldc + lda * ldb;
2961  }
2962  }
2963  }
2964  }
2965  working_space[i1][i2][i3 + sizez_ext] = ldc;
2966  }
2967  }
2968  }
2969 //calculate b=ht*h
2970  i1min = -(lhx - 1), i1max = lhx - 1;
2971  i2min = -(lhy - 1), i2max = lhy - 1;
2972  i3min = -(lhz - 1), i3max = lhz - 1;
2973  for (i3 = i3min; i3 <= i3max; i3++) {
2974  for (i2 = i2min; i2 <= i2max; i2++) {
2975  for (i1 = i1min; i1 <= i1max; i1++) {
2976  ldc = 0;
2977  j3min = -i3;
2978  if (j3min < 0)
2979  j3min = 0;
2980 
2981  j3max = lhz - 1 - i3;
2982  if (j3max > lhz - 1)
2983  j3max = lhz - 1;
2984 
2985  for (j3 = j3min; j3 <= j3max; j3++) {
2986  j2min = -i2;
2987  if (j2min < 0)
2988  j2min = 0;
2989 
2990  j2max = lhy - 1 - i2;
2991  if (j2max > lhy - 1)
2992  j2max = lhy - 1;
2993 
2994  for (j2 = j2min; j2 <= j2max; j2++) {
2995  j1min = -i1;
2996  if (j1min < 0)
2997  j1min = 0;
2998 
2999  j1max = lhx - 1 - i1;
3000  if (j1max > lhx - 1)
3001  j1max = lhx - 1;
3002 
3003  for (j1 = j1min; j1 <= j1max; j1++) {
3004  lda = working_space[j1][j2][j3];
3005  if (i1 + j1 < sizex_ext && i2 + j2 < sizey_ext)
3006  ldb = working_space[i1 + j1][i2 + j2][i3 + j3];
3007 
3008  else
3009  ldb = 0;
3010 
3011  ldc = ldc + lda * ldb;
3012  }
3013  }
3014  }
3015  working_space[i1 - i1min][i2 - i2min][i3 - i3min + 2 * sizez_ext ] = ldc;
3016  }
3017  }
3018  }
3019 //initialization in x1 cube
3020  for (i3 = 0; i3 < sizez_ext; i3++) {
3021  for (i2 = 0; i2 < sizey_ext; i2++) {
3022  for (i1 = 0; i1 < sizex_ext; i1++) {
3023  working_space[i1][i2][i3 + 3 * sizez_ext] = 1;
3024  working_space[i1][i2][i3 + 4 * sizez_ext] = 0;
3025  }
3026  }
3027  }
3028 
3029 //START OF ITERATIONS
3030  for (lindex=0;lindex<deconIterations;lindex++){
3031  for (i3 = 0; i3 < sizez_ext; i3++) {
3032  for (i2 = 0; i2 < sizey_ext; i2++) {
3033  for (i1 = 0; i1 < sizex_ext; i1++) {
3034  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){
3035  ldb = 0;
3036  j3min = i3;
3037  if (j3min > lhz - 1)
3038  j3min = lhz - 1;
3039 
3040  j3min = -j3min;
3041  j3max = sizez_ext - i3 - 1;
3042  if (j3max > lhz - 1)
3043  j3max = lhz - 1;
3044 
3045  j2min = i2;
3046  if (j2min > lhy - 1)
3047  j2min = lhy - 1;
3048 
3049  j2min = -j2min;
3050  j2max = sizey_ext - i2 - 1;
3051  if (j2max > lhy - 1)
3052  j2max = lhy - 1;
3053 
3054  j1min = i1;
3055  if (j1min > lhx - 1)
3056  j1min = lhx - 1;
3057 
3058  j1min = -j1min;
3059  j1max = sizex_ext - i1 - 1;
3060  if (j1max > lhx - 1)
3061  j1max = lhx - 1;
3062 
3063  for (j3 = j3min; j3 <= j3max; j3++) {
3064  for (j2 = j2min; j2 <= j2max; j2++) {
3065  for (j1 = j1min; j1 <= j1max; j1++) {
3066  ldc = working_space[j1 - i1min][j2 - i2min][j3 - i3min + 2 * sizez_ext];
3067  lda = working_space[i1 + j1][i2 + j2][i3 + j3 + 3 * sizez_ext];
3068  ldb = ldb + lda * ldc;
3069  }
3070  }
3071  }
3072  lda = working_space[i1][i2][i3 + 3 * sizez_ext];
3073  ldc = working_space[i1][i2][i3 + 1 * sizez_ext];
3074  if (ldc * lda != 0 && ldb != 0) {
3075  lda = lda * ldc / ldb;
3076  }
3077 
3078  else
3079  lda = 0;
3080  working_space[i1][i2][i3 + 4 * sizez_ext] = lda;
3081  }
3082  }
3083  }
3084  }
3085  for (i3 = 0; i3 < sizez_ext; i3++) {
3086  for (i2 = 0; i2 < sizey_ext; i2++) {
3087  for (i1 = 0; i1 < sizex_ext; i1++)
3088  working_space[i1][i2][i3 + 3 * sizez_ext] = working_space[i1][i2][i3 + 4 * sizez_ext];
3089  }
3090  }
3091  }
3092 //write back resulting spectrum
3093  maximum=0;
3094  for(i = 0;i < sizex_ext; i++){
3095  for(j = 0;j < sizey_ext; j++){
3096  for(k = 0;k < sizez_ext; k++){
3097  working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext] = area * working_space[i][j][k + 3 * sizez_ext];
3098  if(maximum < working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext])
3099  maximum = working_space[(i + positx) % sizex_ext][(j + posity) % sizey_ext][(k + positz) % sizez_ext];
3100  }
3101  }
3102  }
3103 //searching for peaks in deconvolved spectrum
3104  for(i = 1;i < sizex_ext - 1; i++){
3105  for(j = 1;j < sizey_ext - 1; j++){
3106  for(l = 1;l < sizez_ext - 1; l++){
3107  a = working_space[i][j][l];
3108  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]){
3109  if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift && l >= shift && l < ssizez + shift){
3110  if(working_space[i][j][l] > threshold * maximum / 100.0){
3111  if(peak_index < fMaxPeaks){
3112  for(k = i - 1,a = 0,b = 0;k <= i + 1; k++){
3113  a += (Double_t)(k - shift) * working_space[k][j][l];
3114  b += working_space[k][j][l];
3115  }
3116  fPositionX[peak_index] = a / b;
3117  for(k = j - 1,a = 0,b = 0;k <= j + 1; k++){
3118  a += (Double_t)(k - shift) * working_space[i][k][l];
3119  b += working_space[i][k][l];
3120  }
3121  fPositionY[peak_index] = a / b;
3122  for(k = l - 1,a = 0,b = 0;k <= l + 1; k++){
3123  a += (Double_t)(k - shift) * working_space[i][j][k];
3124  b += working_space[i][j][k];
3125  }
3126  fPositionZ[peak_index] = a / b;
3127  peak_index += 1;
3128  }
3129  }
3130  }
3131  }
3132  }
3133  }
3134  }
3135  for(i = 0;i < ssizex; i++){
3136  for(j = 0;j < ssizey; j++){
3137  for(k = 0;k < ssizez; k++){
3138  dest[i][j][k] = working_space[i + shift][j + shift][k + shift];
3139  }
3140  }
3141  }
3142  k = (Int_t)(4 * sigma + 0.5);
3143  k = 4 * k;
3144  for(i = 0;i < ssizex + k; i++){
3145  for(j = 0;j < ssizey + k; j++)
3146  delete[] working_space[i][j];
3147  delete[] working_space[i];
3148  }
3149  delete[] working_space;
3150  fNPeaks = peak_index;
3151  return fNPeaks;
3152 }
3153 
3154 ////////////////////////////////////////////////////////////////////////////////
3155 /// THREE-DIMENSIONAL CLASSICAL PEAK SEARCH FUNCTION
3156 /// This function searches for peaks in source spectrum using
3157 /// the algorithm based on smoothed second differences.
3158 ///
3159 /// Function parameters:
3160 /// - source-pointer to the matrix of source spectrum
3161 /// - ssizex-x length of source spectrum
3162 /// - ssizey-y length of source spectrum
3163 /// - ssizez-z length of source spectrum
3164 /// - sigma-sigma of searched peaks, for details we refer to manual
3165 /// - threshold-threshold value in % for selected peaks, peaks with
3166 /// amplitude less than threshold*highest_peak/100
3167 /// are ignored, see manual
3168 /// - markov-logical variable, if it is true, first the source spectrum
3169 /// is replaced by new spectrum calculated using Markov
3170 /// chains method.
3171 /// - averWindow-averaging window of searched peaks, for details
3172 /// we refer to manual (applies only for Markov method)
3173 
3174 Int_t TSpectrum3::SearchFast(const Double_t***source, Double_t***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez,
3175  Double_t sigma, Double_t threshold,
3176  Bool_t markov, Int_t averWindow)
3177 
3178 {
3179  Int_t i,j,k,l,li,lj,lk,lmin,lmax,xmin,xmax,ymin,ymax,zmin,zmax;
3180  Double_t maxch,plocha = 0,plocha_markov = 0;
3181  Double_t nom,nip,nim,sp,sm,spx,spy,smx,smy,spz,smz;
3182  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;
3183  Double_t a,b,s,f,maximum;
3184  Int_t x,y,z,peak_index=0;
3185  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;
3186  Double_t pocet_sigma = 5;
3187  Int_t number_of_iterations=(Int_t)(4 * sigma + 0.5);
3188  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;
3189  Double_t c[PEAK_WINDOW],s_f_ratio_peaks = 5;
3190  if(sigma < 1){
3191  Error("SearchFast", "Invalid sigma, must be greater than or equal to 1");
3192  return 0;
3193  }
3194 
3195  if(threshold<=0||threshold>=100){
3196  Error("SearchFast", "Invalid threshold, must be positive and less than 100");
3197  return 0;
3198  }
3199 
3200  j = (Int_t)(pocet_sigma*sigma+0.5);
3201  if (j >= PEAK_WINDOW / 2) {
3202  Error("SearchFast", "Too large sigma");
3203  return 0;
3204  }
3205 
3206  if (markov == true) {
3207  if (averWindow <= 0) {
3208  Error("SearchFast", "Averaging window must be positive");
3209  return 0;
3210  }
3211  }
3212 
3213  if(sizex_ext < 2 * number_of_iterations + 1 || sizey_ext < 2 * number_of_iterations + 1 || sizez_ext < 2 * number_of_iterations + 1){
3214  Error("SearchFast", "Too large clipping window");
3215  return 0;
3216  }
3217 
3218  i = (Int_t)(4 * sigma + 0.5);
3219  i = 4 * i;
3220  Double_t ***working_space = new Double_t** [ssizex + i];
3221  for(j = 0;j < ssizex + i; j++){
3222  working_space[j] = new Double_t* [ssizey + i];
3223  for(k = 0;k < ssizey + i; k++)
3224  working_space[j][k] = new Double_t [4 * (ssizez + i)];
3225  }
3226 
3227  for(k = 0;k < sizez_ext; k++){
3228  for(j = 0;j < sizey_ext; j++){
3229  for(i = 0;i < sizex_ext; i++){
3230  if(i < shift){
3231  if(j < shift){
3232  if(k < shift)
3233  working_space[i][j][k + sizez_ext] = source[0][0][0];
3234 
3235  else if(k >= ssizez + shift)
3236  working_space[i][j][k + sizez_ext] = source[0][0][ssizez - 1];
3237 
3238  else
3239  working_space[i][j][k + sizez_ext] = source[0][0][k - shift];
3240  }
3241 
3242  else if(j >= ssizey + shift){
3243  if(k < shift)
3244  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][0];
3245 
3246  else if(k >= ssizez + shift)
3247  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][ssizez - 1];
3248 
3249  else
3250  working_space[i][j][k + sizez_ext] = source[0][ssizey - 1][k - shift];
3251  }
3252 
3253  else{
3254  if(k < shift)
3255  working_space[i][j][k + sizez_ext] = source[0][j - shift][0];
3256 
3257  else if(k >= ssizez + shift)
3258  working_space[i][j][k + sizez_ext] = source[0][j - shift][ssizez - 1];
3259 
3260  else
3261  working_space[i][j][k + sizez_ext] = source[0][j - shift][k - shift];
3262  }
3263  }
3264 
3265  else if(i >= ssizex + shift){
3266  if(j < shift){
3267  if(k < shift)
3268  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][0];
3269 
3270  else if(k >= ssizez + shift)
3271  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][ssizez - 1];
3272 
3273  else
3274  working_space[i][j][k + sizez_ext] = source[ssizex - 1][0][k - shift];
3275  }
3276 
3277  else if(j >= ssizey + shift){
3278  if(k < shift)
3279  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][0];
3280 
3281  else if(k >= ssizez + shift)
3282  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1];
3283 
3284  else
3285  working_space[i][j][k + sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift];
3286  }
3287 
3288  else{
3289  if(k < shift)
3290  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][0];
3291 
3292  else if(k >= ssizez + shift)
3293  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1];
3294 
3295  else
3296  working_space[i][j][k + sizez_ext] = source[ssizex - 1][j - shift][k - shift];
3297  }
3298  }
3299 
3300  else{
3301  if(j < shift){
3302  if(k < shift)
3303  working_space[i][j][k + sizez_ext] = source[i - shift][0][0];
3304 
3305  else if(k >= ssizez + shift)
3306  working_space[i][j][k + sizez_ext] = source[i - shift][0][ssizez - 1];
3307 
3308  else
3309  working_space[i][j][k + sizez_ext] = source[i - shift][0][k - shift];
3310  }
3311 
3312  else if(j >= ssizey + shift){
3313  if(k < shift)
3314  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][0];
3315 
3316  else if(k >= ssizez + shift)
3317  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1];
3318 
3319  else
3320  working_space[i][j][k + sizez_ext] = source[i - shift][ssizey - 1][k - shift];
3321  }
3322 
3323  else{
3324  if(k < shift)
3325  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][0];
3326 
3327  else if(k >= ssizez + shift)
3328  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][ssizez - 1];
3329 
3330  else
3331  working_space[i][j][k + sizez_ext] = source[i - shift][j - shift][k - shift];
3332  }
3333  }
3334  }
3335  }
3336  }
3337  for(i = 1;i <= number_of_iterations; i++){
3338  for(z = i;z < sizez_ext - i; z++){
3339  for(y = i;y < sizey_ext - i; y++){
3340  for(x = i;x < sizex_ext - i; x++){
3341  a = working_space[x][y][z + sizez_ext];
3342  p1 = working_space[x + i][y + i][z - i + sizez_ext];
3343  p2 = working_space[x - i][y + i][z - i + sizez_ext];
3344  p3 = working_space[x + i][y - i][z - i + sizez_ext];
3345  p4 = working_space[x - i][y - i][z - i + sizez_ext];
3346  p5 = working_space[x + i][y + i][z + i + sizez_ext];
3347  p6 = working_space[x - i][y + i][z + i + sizez_ext];
3348  p7 = working_space[x + i][y - i][z + i + sizez_ext];
3349  p8 = working_space[x - i][y - i][z + i + sizez_ext];
3350  s1 = working_space[x + i][y ][z - i + sizez_ext];
3351  s2 = working_space[x ][y + i][z - i + sizez_ext];
3352  s3 = working_space[x - i][y ][z - i + sizez_ext];
3353  s4 = working_space[x ][y - i][z - i + sizez_ext];
3354  s5 = working_space[x + i][y ][z + i + sizez_ext];
3355  s6 = working_space[x ][y + i][z + i + sizez_ext];
3356  s7 = working_space[x - i][y ][z + i + sizez_ext];
3357  s8 = working_space[x ][y - i][z + i + sizez_ext];
3358  s9 = working_space[x - i][y + i][z + sizez_ext];
3359  s10 = working_space[x - i][y - i][z +sizez_ext];
3360  s11 = working_space[x + i][y + i][z +sizez_ext];
3361  s12 = working_space[x + i][y - i][z +sizez_ext];
3362  r1 = working_space[x ][y ][z - i + sizez_ext];
3363  r2 = working_space[x ][y ][z + i + sizez_ext];
3364  r3 = working_space[x - i][y ][z + sizez_ext];
3365  r4 = working_space[x + i][y ][z + sizez_ext];
3366  r5 = working_space[x ][y + i][z + sizez_ext];
3367  r6 = working_space[x ][y - i][z + sizez_ext];
3368  b = (p1 + p3) / 2.0;
3369  if(b > s1)
3370  s1 = b;
3371 
3372  b = (p1 + p2) / 2.0;
3373  if(b > s2)
3374  s2 = b;
3375 
3376  b = (p2 + p4) / 2.0;
3377  if(b > s3)
3378  s3 = b;
3379 
3380  b = (p3 + p4) / 2.0;
3381  if(b > s4)
3382  s4 = b;
3383 
3384  b = (p5 + p7) / 2.0;
3385  if(b > s5)
3386  s5 = b;
3387 
3388  b = (p5 + p6) / 2.0;
3389  if(b > s6)
3390  s6 = b;
3391 
3392  b = (p6 + p8) / 2.0;
3393  if(b > s7)
3394  s7 = b;
3395 
3396  b = (p7 + p8) / 2.0;
3397  if(b > s8)
3398  s8 = b;
3399 
3400  b = (p2 + p6) / 2.0;
3401  if(b > s9)
3402  s9 = b;
3403 
3404  b = (p4 + p8) / 2.0;
3405  if(b > s10)
3406  s10 = b;
3407 
3408  b = (p1 + p5) / 2.0;
3409  if(b > s11)
3410  s11 = b;
3411 
3412  b = (p3 + p7) / 2.0;
3413  if(b > s12)
3414  s12 = b;
3415 
3416  s1 = s1 - (p1 + p3) / 2.0;
3417  s2 = s2 - (p1 + p2) / 2.0;
3418  s3 = s3 - (p2 + p4) / 2.0;
3419  s4 = s4 - (p3 + p4) / 2.0;
3420  s5 = s5 - (p5 + p7) / 2.0;
3421  s6 = s6 - (p5 + p6) / 2.0;
3422  s7 = s7 - (p6 + p8) / 2.0;
3423  s8 = s8 - (p7 + p8) / 2.0;
3424  s9 = s9 - (p2 + p6) / 2.0;
3425  s10 = s10 - (p4 + p8) / 2.0;
3426  s11 = s11 - (p1 + p5) / 2.0;
3427  s12 = s12 - (p3 + p7) / 2.0;
3428  b = (s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
3429  if(b > r1)
3430  r1 = b;
3431 
3432  b = (s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0;
3433  if(b > r2)
3434  r2 = b;
3435 
3436  b = (s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0;
3437  if(b > r3)
3438  r3 = b;
3439 
3440  b = (s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0;
3441  if(b > r4)
3442  r4 = b;
3443 
3444  b = (s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0;
3445  if(b > r5)
3446  r5 = b;
3447 
3448  b = (s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0;
3449  if(b > r6)
3450  r6 = b;
3451 
3452  r1 = r1 - ((s1 + s3) / 2.0 + (s2 + s4) / 2.0 + (p1 + p2 + p3 + p4) / 4.0);
3453  r2 = r2 - ((s5 + s7) / 2.0 + (s6 + s8) / 2.0 + (p5 + p6 + p7 + p8) / 4.0);
3454  r3 = r3 - ((s3 + s7) / 2.0 + (s9 + s10) / 2.0 + (p2 + p4 + p6 + p8) / 4.0);
3455  r4 = r4 - ((s1 + s5) / 2.0 + (s11 + s12) / 2.0 + (p1 + p3 + p5 + p7) / 4.0);
3456  r5 = r5 - ((s9 + s11) / 2.0 + (s2 + s6) / 2.0 + (p1 + p2 + p5 + p6) / 4.0);
3457  r6 = r6 - ((s4 + s8) / 2.0 + (s10 + s12) / 2.0 + (p3 + p4 + p7 + p8) / 4.0);
3458  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;
3459  if(b < a)
3460  a = b;
3461 
3462  working_space[x][y][z] = a;
3463  }
3464  }
3465  }
3466  for(z = i;z < sizez_ext - i; z++){
3467  for(y = i;y < sizey_ext - i; y++){
3468  for(x = i;x < sizex_ext - i; x++){
3469  working_space[x][y][z + sizez_ext] = working_space[x][y][z];
3470  }
3471  }
3472  }
3473  }
3474  for(k = 0;k < sizez_ext; k++){
3475  for(j = 0;j < sizey_ext; j++){
3476  for(i = 0;i < sizex_ext; i++){
3477  if(i < shift){
3478  if(j < shift){
3479  if(k < shift)
3480  working_space[i][j][k + 3 * sizez_ext] = source[0][0][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][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3484 
3485  else
3486  working_space[i][j][k + 3 * sizez_ext] = source[0][0][k - shift] - working_space[i][j][k + sizez_ext];
3487  }
3488 
3489  else if(j >= ssizey + shift){
3490  if(k < shift)
3491  working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3492 
3493  else if(k >= ssizez + shift)
3494  working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3495 
3496  else
3497  working_space[i][j][k + 3 * sizez_ext] = source[0][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3498  }
3499 
3500  else{
3501  if(k < shift)
3502  working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][0] - working_space[i][j][k + sizez_ext];
3503 
3504  else if(k >= ssizez + shift)
3505  working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3506 
3507  else
3508  working_space[i][j][k + 3 * sizez_ext] = source[0][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3509  }
3510  }
3511 
3512  else if(i >= ssizex + shift){
3513  if(j < shift){
3514  if(k < shift)
3515  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][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][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3519 
3520  else
3521  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][0][k - shift] - working_space[i][j][k + sizez_ext];
3522  }
3523 
3524  else if(j >= ssizey + shift){
3525  if(k < shift)
3526  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3527 
3528  else if(k >= ssizez + shift)
3529  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3530 
3531  else
3532  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3533  }
3534 
3535  else{
3536  if(k < shift)
3537  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][0] - working_space[i][j][k + sizez_ext];
3538 
3539  else if(k >= ssizez + shift)
3540  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3541 
3542  else
3543  working_space[i][j][k + 3 * sizez_ext] = source[ssizex - 1][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3544  }
3545  }
3546 
3547  else{
3548  if(j < shift){
3549  if(k < shift)
3550  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][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][0][ssizez - 1] - working_space[i][j][k + sizez_ext];
3554 
3555  else
3556  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][0][k - shift] - working_space[i][j][k + sizez_ext];
3557  }
3558 
3559  else if(j >= ssizey + shift){
3560  if(k < shift)
3561  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][0] - working_space[i][j][k + sizez_ext];
3562 
3563  else if(k >= ssizez + shift)
3564  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][ssizez - 1] - working_space[i][j][k + sizez_ext];
3565 
3566  else
3567  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][ssizey - 1][k - shift] - working_space[i][j][k + sizez_ext];
3568  }
3569 
3570  else{
3571  if(k < shift)
3572  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][0] - working_space[i][j][k + sizez_ext];
3573 
3574  else if(k >= ssizez + shift)
3575  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][ssizez - 1] - working_space[i][j][k + sizez_ext];
3576 
3577  else
3578  working_space[i][j][k + 3 * sizez_ext] = source[i - shift][j - shift][k - shift] - working_space[i][j][k + sizez_ext];
3579  }
3580  }
3581  }
3582  }
3583  }
3584 
3585  for(i = 0;i < sizex_ext; i++){
3586  for(j = 0;j < sizey_ext; j++){
3587  for(k = 0;k < sizez_ext; k++){
3588  if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift && k >= shift && k < ssizez + shift){
3589  working_space[i][j][k + 2 * sizez_ext] = source[i - shift][j - shift][k - shift];
3590  plocha_markov = plocha_markov + source[i - shift][j - shift][k - shift];
3591  }
3592  else
3593  working_space[i][j][k + 2 * sizez_ext] = 0;
3594  }
3595  }
3596  }
3597 
3598  if(markov == true){
3599  xmin = 0;
3600  xmax = sizex_ext - 1;
3601  ymin = 0;
3602  ymax = sizey_ext - 1;
3603  zmin = 0;
3604  zmax = sizez_ext - 1;
3605  for(i = 0,maxch = 0;i < sizex_ext; i++){
3606  for(j = 0;j < sizey_ext;j++){
3607  for(k = 0;k < sizez_ext;k++){
3608  working_space[i][j][k] = 0;
3609  if(maxch < working_space[i][j][k + 2 * sizez_ext])
3610  maxch = working_space[i][j][k + 2 * sizez_ext];
3611 
3612  plocha += working_space[i][j][k + 2 * sizez_ext];
3613  }
3614  }
3615  }
3616  if(maxch == 0) {
3617  k = (Int_t)(4 * sigma + 0.5);
3618  k = 4 * k;
3619  for(i = 0;i < ssizex + k; i++){
3620  for(j = 0;j < ssizey + k; j++)
3621  delete[] working_space[i][j];
3622  delete[] working_space[i];
3623  }
3624  delete [] working_space;
3625  return 0;
3626  }
3627 
3628  nom = 0;
3629  working_space[xmin][ymin][zmin] = 1;
3630  for(i = xmin;i < xmax; i++){
3631  nip = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
3632  nim = working_space[i + 1][ymin][zmin + 2 * sizez_ext] / maxch;
3633  sp = 0,sm = 0;
3634  for(l = 1;l <= averWindow; l++){
3635  if((i + l) > xmax)
3636  a = working_space[xmax][ymin][zmin + 2 * sizez_ext] / maxch;
3637 
3638  else
3639  a = working_space[i + l][ymin][zmin + 2 * sizez_ext] / maxch;
3640 
3641  b = a - nip;
3642  if(a + nip <= 0)
3643  a = 1;
3644 
3645  else
3646  a = TMath::Sqrt(a + nip);
3647 
3648  b = b / a;
3649  b = TMath::Exp(b);
3650  sp = sp + b;
3651  if(i - l + 1 < xmin)
3652  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
3653 
3654  else
3655  a = working_space[i - l + 1][ymin][zmin + 2 * sizez_ext] / maxch;
3656 
3657  b = a - nim;
3658  if(a + nim <= 0)
3659  a = 1;
3660 
3661  else
3662  a = TMath::Sqrt(a + nim);
3663 
3664  b = b / a;
3665  b = TMath::Exp(b);
3666  sm = sm + b;
3667  }
3668  a = sp / sm;
3669  a = working_space[i + 1][ymin][zmin] = a * working_space[i][ymin][zmin];
3670  nom = nom + a;
3671  }
3672  for(i = ymin;i < ymax; i++){
3673  nip = working_space[xmin][i][zmin + 2 * sizez_ext] / maxch;
3674  nim = working_space[xmin][i + 1][zmin + 2 * sizez_ext] / maxch;
3675  sp = 0,sm = 0;
3676  for(l = 1;l <= averWindow; l++){
3677  if((i + l) > ymax)
3678  a = working_space[xmin][ymax][zmin + 2 * sizez_ext] / maxch;
3679 
3680  else
3681  a = working_space[xmin][i + l][zmin + 2 * sizez_ext] / maxch;
3682 
3683  b = a - nip;
3684  if(a + nip <= 0)
3685  a = 1;
3686 
3687  else
3688  a = TMath::Sqrt(a + nip);
3689 
3690  b = b / a;
3691  b = TMath::Exp(b);
3692  sp = sp + b;
3693  if(i - l + 1 < ymin)
3694  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
3695 
3696  else
3697  a = working_space[xmin][i - l + 1][zmin + 2 * sizez_ext] / maxch;
3698 
3699  b = a - nim;
3700  if(a + nim <= 0)
3701  a = 1;
3702 
3703  else
3704  a = TMath::Sqrt(a + nim);
3705 
3706  b = b / a;
3707  b = TMath::Exp(b);
3708  sm = sm + b;
3709  }
3710  a = sp / sm;
3711  a = working_space[xmin][i + 1][zmin] = a * working_space[xmin][i][zmin];
3712  nom = nom + a;
3713  }
3714  for(i = zmin;i < zmax;i++){
3715  nip = working_space[xmin][ymin][i + 2 * sizez_ext] / maxch;
3716  nim = working_space[xmin][ymin][i + 1 + 2 * sizez_ext] / maxch;
3717  sp = 0,sm = 0;
3718  for(l = 1;l <= averWindow;l++){
3719  if((i + l) > zmax)
3720  a = working_space[xmin][ymin][zmax + 2 * sizez_ext] / maxch;
3721 
3722  else
3723  a = working_space[xmin][ymin][i + l + 2 * sizez_ext] / maxch;
3724 
3725  b = a - nip;
3726  if(a + nip <= 0)
3727  a = 1;
3728 
3729  else
3730  a = TMath::Sqrt(a + nip);
3731 
3732  b = b / a;
3733  b = TMath::Exp(b);
3734  sp = sp + b;
3735  if(i - l + 1 < zmin)
3736  a = working_space[xmin][ymin][zmin + 2 * sizez_ext] / maxch;
3737 
3738  else
3739  a = working_space[xmin][ymin][i - l + 1 + 2 * sizez_ext] / maxch;
3740 
3741  b = a - nim;
3742  if(a + nim <= 0)
3743  a = 1;
3744 
3745  else
3746  a = TMath::Sqrt(a + nim);
3747 
3748  b = b / a;
3749  b = TMath::Exp(b);
3750  sm = sm + b;
3751  }
3752  a = sp / sm;
3753  a = working_space[xmin][ymin][i + 1] = a * working_space[xmin][ymin][i];
3754  nom = nom + a;
3755  }
3756  for(i = xmin;i < xmax; i++){
3757  for(j = ymin;j < ymax; j++){
3758  nip = working_space[i][j + 1][zmin + 2 * sizez_ext] / maxch;
3759  nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
3760  spx = 0,smx = 0;
3761  for(l = 1;l <= averWindow; l++){
3762  if(i + l > xmax)
3763  a = working_space[xmax][j][zmin + 2 * sizez_ext] / maxch;
3764 
3765  else
3766  a = working_space[i + l][j][zmin + 2 * sizez_ext] / maxch;
3767 
3768  b = a - nip;
3769  if(a + nip <= 0)
3770  a = 1;
3771 
3772  else
3773  a = TMath::Sqrt(a + nip);
3774 
3775  b = b / a;
3776  b = TMath::Exp(b);
3777  spx = spx + b;
3778  if(i - l + 1 < xmin)
3779  a = working_space[xmin][j][zmin + 2 * sizez_ext] / maxch;
3780 
3781  else
3782  a = working_space[i - l + 1][j][zmin + 2 * sizez_ext] / maxch;
3783 
3784  b = a - nim;
3785  if(a + nim <= 0)
3786  a = 1;
3787 
3788  else
3789  a = TMath::Sqrt(a + nim);
3790 
3791  b = b / a;
3792  b = TMath::Exp(b);
3793  smx = smx + b;
3794  }
3795  spy = 0,smy = 0;
3796  nip = working_space[i + 1][j][zmin + 2 * sizez_ext] / maxch;
3797  nim = working_space[i + 1][j + 1][zmin + 2 * sizez_ext] / maxch;
3798  for(l = 1;l <= averWindow; l++){
3799  if(j + l > ymax)
3800  a = working_space[i][ymax][zmin + 2 * sizez_ext] / maxch;
3801 
3802  else
3803  a = working_space[i][j + l][zmin + 2 * sizez_ext] / maxch;
3804 
3805  b = a - nip;
3806  if(a + nip <= 0)
3807  a = 1;
3808 
3809  else
3810  a = TMath::Sqrt(a + nip);
3811 
3812  b = b / a;
3813  b = TMath::Exp(b);
3814  spy = spy + b;
3815  if(j - l + 1 < ymin)
3816  a = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
3817 
3818  else
3819  a = working_space[i][j - l + 1][zmin + 2 * sizez_ext] / maxch;
3820 
3821  b = a - nim;
3822  if(a + nim <= 0)
3823  a = 1;
3824 
3825  else
3826  a = TMath::Sqrt(a + nim);
3827 
3828  b = b / a;
3829  b = TMath::Exp(b);
3830  smy = smy + b;
3831  }
3832  a = (spx * working_space[i][j + 1][zmin] + spy * working_space[i + 1][j][zmin]) / (smx + smy);
3833  working_space[i + 1][j + 1][zmin] = a;
3834  nom = nom + a;
3835  }
3836  }
3837  for(i = xmin;i < xmax;i++){
3838  for(j = zmin;j < zmax;j++){
3839  nip = working_space[i][ymin][j + 1 + 2 * sizez_ext] / maxch;
3840  nim = working_space[i + 1][ymin][j + 1 + 2 * sizez_ext] / maxch;
3841  spx = 0,smx = 0;
3842  for(l = 1;l <= averWindow; l++){
3843  if(i + l > xmax)
3844  a = working_space[xmax][ymin][j + 2 * sizez_ext] / maxch;
3845 
3846  else
3847  a = working_space[i + l][ymin][j + 2 * sizez_ext] / maxch;
3848 
3849  b = a - nip;
3850  if(a + nip <= 0)
3851  a = 1;
3852 
3853  else
3854  a = TMath::Sqrt(a + nip);
3855 
3856  b = b / a;
3857  b = TMath::Exp(b);
3858  spx = spx + b;
3859  if(i - l + 1 < xmin)
3860  a = working_space[xmin][ymin][j + 2 * sizez_ext] / maxch;
3861 
3862  else
3863  a = working_space[i - l + 1][ymin][j + 2 * sizez_ext] / maxch;
3864 
3865  b = a - nim;
3866  if(a + nim <= 0)
3867  a = 1;
3868 
3869  else
3870  a = TMath::Sqrt(a + nim);
3871 
3872  b = b / a;
3873  b = TMath::Exp(b);
3874  smx = smx + b;
3875  }
3876  spz = 0,smz = 0;
3877  nip = working_space[i + 1][ymin][j + 2 * sizez_ext] / maxch;
3878  nim = working_space[i + 1][ymin][j + 1 + 2 * sizez_ext] / maxch;
3879  for(l = 1;l <= averWindow; l++){
3880  if(j + l > zmax)
3881  a = working_space[i][ymin][zmax + 2 * sizez_ext] / maxch;
3882 
3883  else
3884  a = working_space[i][ymin][j + l + 2 * sizez_ext] / maxch;
3885 
3886  b = a - nip;
3887  if(a + nip <= 0)
3888  a = 1;
3889 
3890  else
3891  a = TMath::Sqrt(a + nip);
3892 
3893  b = b / a;
3894  b = TMath::Exp(b);
3895  spz = spz + b;
3896  if(j - l + 1 < zmin)
3897  a = working_space[i][ymin][zmin + 2 * sizez_ext] / maxch;
3898 
3899  else
3900  a = working_space[i][ymin][j - l + 1 + 2 * sizez_ext] / maxch;
3901 
3902  b = a - nim;
3903  if(a + nim <= 0)
3904  a = 1;
3905 
3906  else
3907  a = TMath::Sqrt(a + nim);
3908 
3909  b = b / a;
3910  b = TMath::Exp(b);
3911  smz = smz + b;
3912  }
3913  a = (spx * working_space[i][ymin][j + 1] + spz * working_space[i + 1][ymin][j]) / (smx + smz);
3914  working_space[i + 1][ymin][j + 1] = a;
3915  nom = nom + a;
3916  }
3917  }
3918  for(i = ymin;i < ymax;i++){
3919  for(j = zmin;j < zmax;j++){
3920  nip = working_space[xmin][i][j + 1 + 2 * sizez_ext] / maxch;
3921  nim = working_space[xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
3922  spy = 0,smy = 0;
3923  for(l = 1;l <= averWindow; l++){
3924  if(i + l > ymax)
3925  a = working_space[xmin][ymax][j + 2 * sizez_ext] / maxch;
3926 
3927  else
3928  a = working_space[xmin][i + l][j + 2 * sizez_ext] / maxch;
3929 
3930  b = a - nip;
3931  if(a + nip <= 0)
3932  a = 1;
3933 
3934  else
3935  a = TMath::Sqrt(a + nip);
3936 
3937  b = b / a;
3938  b = TMath::Exp(b);
3939  spy = spy + b;
3940  if(i - l + 1 < ymin)
3941  a = working_space[xmin][ymin][j + 2 * sizez_ext] / maxch;
3942 
3943  else
3944  a = working_space[xmin][i - l + 1][j + 2 * sizez_ext] / maxch;
3945 
3946  b = a - nim;
3947  if(a + nim <= 0)
3948  a = 1;
3949 
3950  else
3951  a = TMath::Sqrt(a + nim);
3952 
3953  b = b / a;
3954  b = TMath::Exp(b);
3955  smy = smy + b;
3956  }
3957  spz = 0,smz = 0;
3958  nip = working_space[xmin][i + 1][j + 2 * sizez_ext] / maxch;
3959  nim = working_space[xmin][i + 1][j + 1 + 2 * sizez_ext] / maxch;
3960  for(l = 1;l <= averWindow; l++){
3961  if(j + l > zmax)
3962  a = working_space[xmin][i][zmax + 2 * sizez_ext] / maxch;
3963 
3964  else
3965  a = working_space[xmin][i][j + l + 2 * sizez_ext] / maxch;
3966 
3967  b = a - nip;
3968  if(a + nip <= 0)
3969  a = 1;
3970 
3971  else
3972  a = TMath::Sqrt(a + nip);
3973 
3974  b = b / a;
3975  b = TMath::Exp(b);
3976  spz = spz + b;
3977  if(j - l + 1 < zmin)
3978  a = working_space[xmin][i][zmin + 2 * sizez_ext] / maxch;
3979 
3980  else
3981  a = working_space[xmin][i][j - l + 1 + 2 * sizez_ext] / maxch;
3982 
3983  b = a - nim;
3984  if(a + nim <= 0)
3985  a = 1;
3986 
3987  else
3988  a = TMath::Sqrt(a + nim);
3989 
3990  b = b / a;
3991  b = TMath::Exp(b);
3992  smz = smz + b;
3993  }
3994  a = (spy * working_space[xmin][i][j + 1] + spz * working_space[xmin][i + 1][j]) / (smy + smz);
3995  working_space[xmin][i + 1][j + 1] = a;
3996  nom = nom + a;
3997  }
3998  }
3999  for(i = xmin;i < xmax; i++){
4000  for(j = ymin;j < ymax; j++){
4001  for(k = zmin;k < zmax; k++){
4002  nip = working_space[i][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4003  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4004  spx = 0,smx = 0;
4005  for(l = 1;l <= averWindow; l++){
4006  if(i + l > xmax)
4007  a = working_space[xmax][j][k + 2 * sizez_ext] / maxch;
4008 
4009  else
4010  a = working_space[i + l][j][k + 2 * sizez_ext] / maxch;
4011 
4012  b = a - nip;
4013  if(a + nip <= 0)
4014  a = 1;
4015 
4016  else
4017  a = TMath::Sqrt(a + nip);
4018 
4019  b = b / a;
4020  b = TMath::Exp(b);
4021  spx = spx + b;
4022  if(i - l + 1 < xmin)
4023  a = working_space[xmin][j][k + 2 * sizez_ext] / maxch;
4024 
4025  else
4026  a = working_space[i - l + 1][j][k + 2 * sizez_ext] / maxch;
4027 
4028  b = a - nim;
4029  if(a + nim <= 0)
4030  a = 1;
4031 
4032  else
4033  a = TMath::Sqrt(a + nim);
4034 
4035  b = b / a;
4036  b = TMath::Exp(b);
4037  smx = smx + b;
4038  }
4039  spy = 0,smy = 0;
4040  nip = working_space[i + 1][j][k + 1 + 2 * sizez_ext] / maxch;
4041  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4042  for(l = 1;l <= averWindow; l++){
4043  if(j + l > ymax)
4044  a = working_space[i][ymax][k + 2 * sizez_ext] / maxch;
4045 
4046  else
4047  a = working_space[i][j + l][k + 2 * sizez_ext] / maxch;
4048 
4049  b = a - nip;
4050  if(a + nip <= 0)
4051  a = 1;
4052 
4053  else
4054  a = TMath::Sqrt(a + nip);
4055 
4056  b = b / a;
4057  b = TMath::Exp(b);
4058  spy = spy + b;
4059  if(j - l + 1 < ymin)
4060  a = working_space[i][ymin][k + 2 * sizez_ext] / maxch;
4061 
4062  else
4063  a = working_space[i][j - l + 1][k + 2 * sizez_ext] / maxch;
4064 
4065  b = a - nim;
4066  if(a + nim <= 0)
4067  a = 1;
4068 
4069  else
4070  a = TMath::Sqrt(a + nim);
4071 
4072  b = b / a;
4073  b = TMath::Exp(b);
4074  smy = smy + b;
4075  }
4076  spz = 0,smz = 0;
4077  nip = working_space[i + 1][j + 1][k + 2 * sizez_ext] / maxch;
4078  nim = working_space[i + 1][j + 1][k + 1 + 2 * sizez_ext] / maxch;
4079  for(l = 1;l <= averWindow; l++ ){
4080  if(j + l > zmax)
4081  a = working_space[i][j][zmax + 2 * sizez_ext] / maxch;
4082 
4083  else
4084  a = working_space[i][j][k + l + 2 * sizez_ext] / maxch;
4085 
4086  b = a - nip;
4087  if(a + nip <= 0)
4088  a = 1;
4089 
4090  else
4091  a = TMath::Sqrt(a + nip);
4092 
4093  b = b / a;
4094  b = TMath::Exp(b);
4095  spz = spz + b;
4096  if(j - l + 1 < ymin)
4097  a = working_space[i][j][zmin + 2 * sizez_ext] / maxch;
4098 
4099  else
4100  a = working_space[i][j][k - l + 1 + 2 * sizez_ext] / maxch;
4101 
4102  b = a - nim;
4103  if(a + nim <= 0)
4104  a = 1;
4105 
4106  else
4107  a = TMath::Sqrt(a + nim);
4108 
4109  b = b / a;
4110  b = TMath::Exp(b);
4111  smz = smz + b;
4112  }
4113  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);
4114  working_space[i + 1][j + 1][k + 1] = a;
4115  nom = nom + a;
4116  }
4117  }
4118  }
4119  a = 0;
4120  for(i = xmin;i <= xmax; i++){
4121  for(j = ymin;j <= ymax; j++){
4122  for(k = zmin;k <= zmax; k++){
4123  working_space[i][j][k] = working_space[i][j][k] / nom;
4124  a+=working_space[i][j][k];
4125  }
4126  }
4127  }
4128  for(i = 0;i < sizex_ext; i++){
4129  for(j = 0;j < sizey_ext; j++){
4130  for(k = 0;k < sizez_ext; k++){
4131  working_space[i][j][k + 2 * sizez_ext] = working_space[i][j][k] * plocha_markov / a;
4132  }
4133  }
4134  }
4135  }
4136 
4137  maximum = 0;
4138  for(k = 0;k < ssizez; k++){
4139  for(j = 0;j < ssizey; j++){
4140  for(i = 0;i < ssizex; i++){
4141  working_space[i][j][k] = 0;
4142  working_space[i][j][sizez_ext + k] = 0;
4143  if(working_space[i][j][k + 3 * sizez_ext] > maximum)
4144  maximum=working_space[i][j][k+3*sizez_ext];
4145  }
4146  }
4147  }
4148  for(i = 0;i < PEAK_WINDOW; i++){
4149  c[i] = 0;
4150  }
4151  j = (Int_t)(pocet_sigma * sigma + 0.5);
4152  for(i = -j;i <= j; i++){
4153  a=i;
4154  a = -a * a;
4155  b = 2.0 * sigma * sigma;
4156  a = a / b;
4157  a = TMath::Exp(a);
4158  s = i;
4159  s = s * s;
4160  s = s - sigma * sigma;
4161  s = s / (sigma * sigma * sigma * sigma);
4162  s = s * a;
4163  c[PEAK_WINDOW / 2 + i] = s;
4164  }
4165  norma = 0;
4166  for(i = 0;i < PEAK_WINDOW; i++){
4167  norma = norma + TMath::Abs(c[i]);
4168  }
4169  for(i = 0;i < PEAK_WINDOW; i++){
4170  c[i] = c[i] / norma;
4171  }
4172  a = pocet_sigma * sigma + 0.5;
4173  i = (Int_t)a;
4174  zmin = i;
4175  zmax = sizez_ext - i - 1;
4176  ymin = i;
4177  ymax = sizey_ext - i - 1;
4178  xmin = i;
4179  xmax = sizex_ext - i - 1;
4180  lmin = PEAK_WINDOW / 2 - i;
4181  lmax = PEAK_WINDOW / 2 + i;
4182  for(i = xmin;i <= xmax; i++){
4183  for(j = ymin;j <= ymax; j++){
4184  for(k = zmin;k <= zmax; k++){
4185  s = 0,f = 0;
4186  for(li = lmin;li <= lmax; li++){
4187  for(lj = lmin;lj <= lmax; lj++){
4188  for(lk = lmin;lk <= lmax; lk++){
4189  a = working_space[i + li - PEAK_WINDOW / 2][j + lj - PEAK_WINDOW / 2][k + lk - PEAK_WINDOW / 2 + 2 * sizez_ext];
4190  b = c[li] * c[lj] * c[lk];
4191  s += a * b;
4192  f += a * b * b;
4193  }
4194  }
4195  }
4196  working_space[i][j][k] = s;
4197  working_space[i][j][k + sizez_ext] = TMath::Sqrt(f);
4198  }
4199  }
4200  }
4201  for(x = xmin;x <= xmax; x++){
4202  for(y = ymin + 1;y < ymax; y++){
4203  for(z = zmin + 1;z < zmax; z++){
4204  val = working_space[x][y][z];
4205  val1 = working_space[x - 1][y - 1][z - 1];
4206  val2 = working_space[x ][y - 1][z - 1];
4207  val3 = working_space[x + 1][y - 1][z - 1];
4208  val4 = working_space[x - 1][y ][z - 1];
4209  val5 = working_space[x ][y ][z - 1];
4210  val6 = working_space[x + 1][y ][z - 1];
4211  val7 = working_space[x - 1][y + 1][z - 1];
4212  val8 = working_space[x ][y + 1][z - 1];
4213  val9 = working_space[x + 1][y + 1][z - 1];
4214  val10 = working_space[x - 1][y - 1][z ];
4215  val11 = working_space[x ][y - 1][z ];
4216  val12 = working_space[x + 1][y - 1][z ];
4217  val13 = working_space[x - 1][y ][z ];
4218  val14 = working_space[x + 1][y ][z ];
4219  val15 = working_space[x - 1][y + 1][z ];
4220  val16 = working_space[x ][y + 1][z ];
4221  val17 = working_space[x + 1][y + 1][z ];
4222  val18 = working_space[x - 1][y - 1][z + 1];
4223  val19 = working_space[x ][y - 1][z + 1];
4224  val20 = working_space[x + 1][y - 1][z + 1];
4225  val21 = working_space[x - 1][y ][z + 1];
4226  val22 = working_space[x ][y ][z + 1];
4227  val23 = working_space[x + 1][y ][z + 1];
4228  val24 = working_space[x - 1][y + 1][z + 1];
4229  val25 = working_space[x ][y + 1][z + 1];
4230  val26 = working_space[x + 1][y + 1][z + 1];
4231  f = -s_f_ratio_peaks * working_space[x][y][z + sizez_ext];
4232  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){
4233  s=0,f=0;
4234  for(li = lmin;li <= lmax; li++){
4235  a = working_space[x + li - PEAK_WINDOW / 2][y][z + 2 * sizez_ext];
4236  s += a * c[li];
4237  f += a * c[li] * c[li];
4238  }
4239  f = -s_f_ratio_peaks * sqrt(f);
4240  if(s < f){
4241  s = 0,f = 0;
4242  for(li = lmin;li <= lmax; li++){
4243  a = working_space[x][y + li - PEAK_WINDOW / 2][z + 2 * sizez_ext];
4244  s += a * c[li];
4245  f += a * c[li] * c[li];
4246  }
4247  f = -s_f_ratio_peaks * sqrt(f);
4248  if(s < f){
4249  s = 0,f = 0;
4250  for(li = lmin;li <= lmax; li++){
4251  a = working_space[x][y][z + li - PEAK_WINDOW / 2 + 2 * sizez_ext];
4252  s += a * c[li];
4253  f += a * c[li] * c[li];
4254  }
4255  f = -s_f_ratio_peaks * sqrt(f);
4256  if(s < f){
4257  s = 0,f = 0;
4258  for(li = lmin;li <= lmax; li++){
4259  for(lj = lmin;lj <= lmax; lj++){
4260  a = working_space[x + li - PEAK_WINDOW / 2][y + lj - PEAK_WINDOW / 2][z + 2 * sizez_ext];
4261  s += a * c[li] * c[lj];
4262  f += a * c[li] * c[li] * c[lj] * c[lj];
4263  }
4264  }
4265  f = s_f_ratio_peaks * sqrt(f);
4266  if(s > f){
4267  s = 0,f = 0;
4268  for(li = lmin;li <= lmax; li++){
4269  for(lj = lmin;lj <= lmax; lj++){
4270  a = working_space[x + li - PEAK_WINDOW / 2][y][z + lj - PEAK_WINDOW / 2 + 2 * sizez_ext];
4271  s += a * c[li] * c[lj];
4272  f += a * c[li] * c[li] * c[lj] * c[lj];
4273  }
4274  }
4275  f = s_f_ratio_peaks * sqrt(f);
4276  if(s > f){
4277  s = 0,f = 0;
4278  for(li = lmin;li <= lmax; li++){
4279  for(lj=lmin;lj<=lmax;lj++){
4280  a = working_space[x][y + li - PEAK_WINDOW / 2][z + lj - PEAK_WINDOW / 2 + 2 * sizez_ext];
4281  s += a * c[li] * c[lj];
4282  f += a * c[li] * c[li] * c[lj] * c[lj];
4283  }
4284  }
4285  f = s_f_ratio_peaks * sqrt(f);
4286  if(s > f){
4287  if(x >= shift && x < ssizex + shift && y >= shift && y < ssizey + shift && z >= shift && z < ssizez + shift){
4288  if(working_space[x][y][z + 3 * sizez_ext] > threshold * maximum / 100.0){
4289  if(peak_index<fMaxPeaks){
4290  for(k = x - 1,a = 0,b = 0;k <= x + 1; k++){
4291  a += (Double_t)(k - shift) * working_space[k][y][z];
4292  b += working_space[k][y][z];
4293  }
4294  fPositionX[peak_index] = a / b;
4295  for(k = y - 1,a = 0,b = 0;k <= y + 1; k++){
4296  a += (Double_t)(k - shift) * working_space[x][k][z];
4297  b += working_space[x][k][z];
4298  }
4299  fPositionY[peak_index] = a / b;
4300  for(k = z - 1,a = 0,b = 0;k <= z + 1; k++){
4301  a += (Double_t)(k - shift) * working_space[x][y][k];
4302  b += working_space[x][y][k];
4303  }
4304  fPositionZ[peak_index] = a / b;
4305  peak_index += 1;
4306  }
4307  }
4308  }
4309  }
4310  }
4311  }
4312  }
4313  }
4314  }
4315  }
4316  }
4317  }
4318  }
4319  for(i = 0;i < ssizex; i++){
4320  for(j = 0;j < ssizey; j++){
4321  for(k = 0;k < ssizez; k++){
4322  val = -working_space[i + shift][j + shift][k + shift];
4323  if( val < 0)
4324  val = 0;
4325  dest[i][j][k] = val;
4326  }
4327  }
4328  }
4329  k = (Int_t)(4 * sigma + 0.5);
4330  k = 4 * k;
4331  for(i = 0;i < ssizex + k; i++){
4332  for(j = 0;j < ssizey + k; j++)
4333  delete[] working_space[i][j];
4334  delete[] working_space[i];
4335  }
4336  delete[] working_space;
4337  fNPeaks = peak_index;
4338  return fNPeaks;
4339 }
void SetResolution(Double_t resolution=1)
NOT USED resolution: determines resolution of the neighbouring peaks default value is 1 correspond to...
Definition: TSpectrum3.cxx:227
Double_t * fPosition
[fNPeaks] array of current peak positions
Definition: TSpectrum3.h:22
float xmin
Definition: THbookFile.cxx:93
Int_t fNPeaks
number of peaks found
Definition: TSpectrum3.h:21
const char Option_t
Definition: RtypesCore.h:62
float ymin
Definition: THbookFile.cxx:93
Int_t SearchFast(const Double_t ***source, Double_t ***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez, Double_t sigma, Double_t threshold, Bool_t markov, Int_t averWindow)
THREE-DIMENSIONAL CLASSICAL PEAK SEARCH FUNCTION This function searches for peaks in source spectrum ...
const char * SmoothMarkov(Double_t ***source, Int_t ssizex, Int_t ssizey, Int_t ssizez, Int_t averWindow)
This function calculates smoothed spectrum from source spectrum based on Markov chain method...
Definition: TSpectrum3.cxx:859
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4882
const char * Deconvolution(Double_t ***source, const Double_t ***resp, Int_t ssizex, Int_t ssizey, Int_t ssizez, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
This function calculates deconvolution from source spectrum according to response spectrum The result...
#define f(i)
Definition: RSha256.hxx:104
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
virtual Int_t GetDimension() const
Definition: TH1.h:278
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:476
Double_t * fPositionZ
[fNPeaks] Z positions of peaks
Definition: TSpectrum3.h:25
const Double_t sigma
static constexpr double s
virtual const char * Background(const TH1 *hist, Int_t niter, Option_t *option="goff")
This function calculates background spectrum from source in h.
Definition: TSpectrum3.cxx:115
Double_t * fPositionY
[fNPeaks] Y positions of peaks
Definition: TSpectrum3.h:24
float ymax
Definition: THbookFile.cxx:93
Double_t fResolution
NOT USED resolution of the neighboring peaks
Definition: TSpectrum3.h:26
auto * a
Definition: textangle.C:12
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only r...
Definition: VectorUtil.h:327
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
#define s1(x)
Definition: RSha256.hxx:91
TAxis * GetYaxis()
Definition: TH1.h:317
Advanced 3-dimensional spectra processing functions.
Definition: TSpectrum3.h:18
float xmax
Definition: THbookFile.cxx:93
#define h(i)
Definition: RSha256.hxx:106
const Bool_t kFALSE
Definition: RtypesCore.h:88
#define d(i)
Definition: RSha256.hxx:102
Int_t fMaxPeaks
Maximum number of peaks to be found.
Definition: TSpectrum3.h:20
Double_t Exp(Double_t x)
Definition: TMath.h:717
Int_t SearchHighRes(const Double_t ***source, Double_t ***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez, Double_t sigma, Double_t threshold, Bool_t backgroundRemove, Int_t deconIterations, Bool_t markov, Int_t averWindow)
This function searches for peaks in source spectrum It is based on deconvolution method.
#define ClassImp(name)
Definition: Rtypes.h:365
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
The TH1 histogram class.
Definition: TH1.h:56
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
TH1 * fHistogram
resulting histogram
Definition: TSpectrum3.h:27
#define PEAK_WINDOW
Definition: TSpectrum3.cxx:51
Double_t * fPositionX
[fNPeaks] X positions of peaks
Definition: TSpectrum3.h:23
TAxis * GetZaxis()
Definition: TH1.h:318
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
auto * l
Definition: textangle.C:4
#define dest(otri, vertexptr)
Definition: triangle.c:1040
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
#define c(i)
Definition: RSha256.hxx:101
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Int_t GetNbins() const
Definition: TAxis.h:121
double exp(double)
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
This function searches for peaks in source spectrum in hin The number of found peaks and their positi...
Definition: TSpectrum3.cxx:160
const Bool_t kTRUE
Definition: RtypesCore.h:87
const Int_t n
Definition: legend1.C:16
virtual void Print(Option_t *option="") const
Print the array of positions.
Definition: TSpectrum3.cxx:126
TSpectrum3()
Constructor.
Definition: TSpectrum3.cxx:58
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
virtual ~TSpectrum3()
Destructor.
Definition: TSpectrum3.cxx:97