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