Logo ROOT  

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