Logo ROOT  

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