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