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