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