Logo ROOT  

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