Logo ROOT   6.18/05
Reference Guide
TSpectrum2.cxx
Go to the documentation of this file.
1// @(#)root/spectrum:$Id$
2// Author: Miroslav Morhac 17/01/2006
3
4/** \class TSpectrum2
5 \ingroup Spectrum
6 \brief Advanced 2-dimensional spectra processing
7 \author Miroslav Morhac
8
9 This class contains advanced spectra processing functions.
10
11 - One-dimensional background estimation functions
12 - Two-dimensional background estimation functions
13 - One-dimensional smoothing functions
14 - Two-dimensional smoothing functions
15 - One-dimensional deconvolution functions
16 - Two-dimensional deconvolution functions
17 - One-dimensional peak search functions
18 - Two-dimensional peak search functions
19
20 The algorithms in this class have been published in the following references:
21
22 1. M.Morhac et al.: Background elimination methods for multidimensional coincidence gamma-ray spectra. Nuclear Instruments and Methods in Physics Research A 401 (1997) 113-132.
23 2. M.Morhac et al.: Efficient one- and two-dimensional Gold deconvolution and its application to gamma-ray spectra decomposition. Nuclear Instruments and Methods in Physics Research A 401 (1997) 385-408.
24 3. M.Morhac et al.: Identification of peaks in multidimensional coincidence gamma-ray spectra. Nuclear Instruments and Methods in Research Physics A 443(2000), 108-125.
25
26 These NIM papers are also available as doc or ps files from:
27
28 - [SpectrumDec.ps.gz](ftp://root.cern.ch/root/SpectrumDec.ps.gz)
29 - [SpectrumSrc.ps.gz](ftp://root.cern.ch/root/SpectrumSrc.ps.gz)
30 - [SpectrumBck.ps.gz](ftp://root.cern.ch/root/SpectrumBck.ps.gz)
31
32 See also the
33 [online documentation](https://root.cern.ch/guides/tspectrum-manual) and
34 [tutorials](https://root.cern.ch/doc/master/group__tutorial__spectrum.html).
35
36 All the figures in this page were prepared using the DaqProVis
37 system, Data Acquisition, Processing and Visualization system,
38 developed at the Institute of Physics, Slovak Academy of Sciences, Bratislava,
39 Slovakia.
40*/
41
42#include "TSpectrum2.h"
43#include "TPolyMarker.h"
44#include "TList.h"
45#include "TH1.h"
46#include "TMath.h"
47#define PEAK_WINDOW 1024
48
49Int_t TSpectrum2::fgIterations = 3;
50Int_t TSpectrum2::fgAverageWindow = 3;
51
52ClassImp(TSpectrum2);
53
54////////////////////////////////////////////////////////////////////////////////
55/// Constructor.
56
57TSpectrum2::TSpectrum2() :TNamed("Spectrum", "Miroslav Morhac peak finder")
58{
59 Int_t n = 100;
60 fMaxPeaks = n;
61 fPosition = new Double_t[n];
62 fPositionX = new Double_t[n];
63 fPositionY = new Double_t[n];
64 fResolution = 1;
65 fHistogram = 0;
66 fNPeaks = 0;
67}
68
69////////////////////////////////////////////////////////////////////////////////
70/// - maxpositions: maximum number of peaks
71/// - resolution: *NOT USED* determines resolution of the neighbouring peaks
72/// default value is 1 correspond to 3 sigma distance
73/// between peaks. Higher values allow higher resolution
74/// (smaller distance between peaks.
75/// May be set later through SetResolution.
76
77TSpectrum2::TSpectrum2(Int_t maxpositions, Double_t resolution) :TNamed("Spectrum", "Miroslav Morhac peak finder")
78{
79 Int_t n = maxpositions;
80 fMaxPeaks = n;
81 fPosition = new Double_t[n];
82 fPositionX = new Double_t[n];
83 fPositionY = new Double_t[n];
84 fHistogram = 0;
85 fNPeaks = 0;
86 SetResolution(resolution);
87}
88
89////////////////////////////////////////////////////////////////////////////////
90/// Destructor.
91
92TSpectrum2::~TSpectrum2()
93{
94 delete [] fPosition;
95 delete [] fPositionX;
96 delete [] fPositionY;
97 delete fHistogram;
98}
99
100////////////////////////////////////////////////////////////////////////////////
101/// static function: Set average window of searched peaks
102/// see TSpectrum2::SearchHighRes
103
104void TSpectrum2::SetAverageWindow(Int_t w)
105{
106 fgAverageWindow = w;
107}
108
109////////////////////////////////////////////////////////////////////////////////
110/// static function: Set max number of decon iterations in deconvolution operation
111/// see TSpectrum2::SearchHighRes
112
113void TSpectrum2::SetDeconIterations(Int_t n)
114{
115 fgIterations = n;
116}
117
118////////////////////////////////////////////////////////////////////////////////
119/// This function calculates the background spectrum in the input histogram h.
120/// The background is returned as a histogram.
121///
122/// Function parameters:
123/// - h: input 2-d histogram
124/// - numberIterations, (default value = 20)
125/// Increasing numberIterations make the result smoother and lower.
126/// - option: may contain one of the following options
127/// - to set the direction parameter
128/// "BackIncreasingWindow". By default the direction is BackDecreasingWindow
129/// - filterOrder-order of clipping filter. Possible values:
130/// - "BackOrder2" (default)
131/// - "BackOrder4"
132/// - "BackOrder6"
133/// - "BackOrder8"
134/// - "nosmoothing"- if selected, the background is not smoothed
135/// By default the background is smoothed.
136/// - smoothWindow-width of smoothing window. Possible values:
137/// - "BackSmoothing3" (default)
138/// - "BackSmoothing5"
139/// - "BackSmoothing7"
140/// - "BackSmoothing9"
141/// - "BackSmoothing11"
142/// - "BackSmoothing13"
143/// - "BackSmoothing15"
144/// - "Compton" if selected the estimation of Compton edge
145/// will be included.
146/// - "same" : if this option is specified, the resulting background
147/// histogram is superimposed on the picture in the current pad.
148///
149/// NOTE that the background is only evaluated in the current range of h.
150/// ie, if h has a bin range (set via h->GetXaxis()->SetRange(binmin,binmax),
151/// the returned histogram will be created with the same number of bins
152/// as the input histogram h, but only bins from binmin to binmax will be filled
153/// with the estimated background.
154
155TH1 *TSpectrum2::Background(const TH1 * h, Int_t number_of_iterations,
156 Option_t * option)
157{
158 Error("Background","function not yet implemented: h=%s, iter=%d, option=%sn"
159 , h->GetName(), number_of_iterations, option);
160 return 0;
161}
162
163////////////////////////////////////////////////////////////////////////////////
164/// Print the array of positions.
165
166void TSpectrum2::Print(Option_t *) const
167{
168 printf("\nNumber of positions = %d\n",fNPeaks);
169 for (Int_t i=0;i<fNPeaks;i++) {
170 printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
171 }
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// This function searches for peaks in source spectrum in hin
176/// The number of found peaks and their positions are written into
177/// the members fNpeaks and fPositionX.
178/// The search is performed in the current histogram range.
179///
180/// Function parameters:
181/// - hin: pointer to the histogram of source spectrum
182/// - sigma: sigma of searched peaks, for details we refer to manual
183/// - threshold: (default=0.05) peaks with amplitude less than
184/// threshold*highest_peak are discarded. 0<threshold<1
185///
186/// By default, the background is removed before deconvolution.
187/// Specify the option "nobackground" to not remove the background.
188///
189/// By default the "Markov" chain algorithm is used.
190/// Specify the option "noMarkov" to disable this algorithm
191/// Note that by default the source spectrum is replaced by a new spectrum
192///
193/// By default a polymarker object is created and added to the list of
194/// functions of the histogram. The histogram is drawn with the specified
195/// option and the polymarker object drawn on top of the histogram.
196/// The polymarker coordinates correspond to the npeaks peaks found in
197/// the histogram.
198/// A pointer to the polymarker object can be retrieved later via:
199/// ~~~ {.cpp}
200/// TList *functions = hin->GetListOfFunctions();
201/// TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker")
202/// ~~~
203/// Specify the option "goff" to disable the storage and drawing of the
204/// polymarker.
205
206Int_t TSpectrum2::Search(const TH1 * hin, Double_t sigma,
207 Option_t * option, Double_t threshold)
208{
209 if (hin == 0)
210 return 0;
211 Int_t dimension = hin->GetDimension();
212 if (dimension != 2) {
213 Error("Search", "Must be a 2-d histogram");
214 return 0;
215 }
216
217 TString opt = option;
218 opt.ToLower();
219 Bool_t background = kTRUE;
220 if (opt.Contains("nobackground")) {
221 background = kFALSE;
222 opt.ReplaceAll("nobackground","");
223 }
224 Bool_t markov = kTRUE;
225 if (opt.Contains("nomarkov")) {
226 markov = kFALSE;
227 opt.ReplaceAll("nomarkov","");
228 }
229
230 Int_t sizex = hin->GetXaxis()->GetNbins();
231 Int_t sizey = hin->GetYaxis()->GetNbins();
232 Int_t i, j, binx,biny, npeaks;
233 Double_t ** source = new Double_t*[sizex];
234 Double_t ** dest = new Double_t*[sizex];
235 for (i = 0; i < sizex; i++) {
236 source[i] = new Double_t[sizey];
237 dest[i] = new Double_t[sizey];
238 for (j = 0; j < sizey; j++) {
239 source[i][j] = hin->GetBinContent(i + 1, j + 1);
240 }
241 }
242 //npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, kTRUE, 3, kTRUE, 10);
243 //the smoothing option is used for 1-d but not for 2-d histograms
244 npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, background, fgIterations, markov, fgAverageWindow);
245
246 //The logic in the loop should be improved to use the fact
247 //that fPositionX,Y give a precise position inside a bin.
248 //The current algorithm takes the center of the bin only.
249 for (i = 0; i < npeaks; i++) {
250 binx = 1 + Int_t(fPositionX[i] + 0.5);
251 biny = 1 + Int_t(fPositionY[i] + 0.5);
252 fPositionX[i] = hin->GetXaxis()->GetBinCenter(binx);
253 fPositionY[i] = hin->GetYaxis()->GetBinCenter(biny);
254 }
255 for (i = 0; i < sizex; i++) {
256 delete [] source[i];
257 delete [] dest[i];
258 }
259 delete [] source;
260 delete [] dest;
261
262 if (opt.Contains("goff"))
263 return npeaks;
264 if (!npeaks) return 0;
265 TPolyMarker * pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
266 if (pm) {
267 hin->GetListOfFunctions()->Remove(pm);
268 delete pm;
269 }
270 pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
271 hin->GetListOfFunctions()->Add(pm);
272 pm->SetMarkerStyle(23);
273 pm->SetMarkerColor(kRed);
274 pm->SetMarkerSize(1.3);
275 ((TH1*)hin)->Draw(option);
276 return npeaks;
277}
278
279////////////////////////////////////////////////////////////////////////////////
280/// *NOT USED*
281/// resolution: determines resolution of the neighboring peaks
282/// default value is 1 correspond to 3 sigma distance
283/// between peaks. Higher values allow higher resolution
284/// (smaller distance between peaks.
285/// May be set later through SetResolution.
286
287void TSpectrum2::SetResolution(Double_t resolution)
288{
289 if (resolution > 1)
290 fResolution = resolution;
291 else
292 fResolution = 1;
293}
294
295////////////////////////////////////////////////////////////////////////////////
296/// This function calculates background spectrum from source spectrum.
297/// The result is placed to the array pointed by spectrum pointer.
298///
299/// Function parameters:
300/// - spectrum-pointer to the array of source spectrum
301/// - ssizex-x length of spectrum
302/// - ssizey-y length of spectrum
303/// - numberIterationsX-maximal x width of clipping window
304/// - numberIterationsY-maximal y width of clipping window
305/// for details we refer to manual
306/// - direction- direction of change of clipping window
307/// - possible values=kBackIncreasingWindow
308/// kBackDecreasingWindow
309/// - filterType-determines the algorithm of the filtering
310/// - possible values=kBackSuccessiveFiltering
311/// kBackOneStepFiltering
312///
313/// ### Background estimation
314///
315/// Goal: Separation of useful information (peaks) from useless information (background)
316///
317/// - method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping algorithm [1]
318///
319/// - there exist two algorithms for the estimation of new value in the channel \f$ i_1,i_2\f$
320///
321/// Algorithm based on Successive Comparisons:
322///
323/// It is an extension of one-dimensional SNIP algorithm to another dimension. For
324/// details we refer to [2].
325///
326/// Algorithm based on One Step Filtering:
327///
328/// New value in the estimated channel is calculated as:
329/// \f[ a = \nu_{p-1}(i_1,i_2)\f]
330/// \f[
331/// b = \frac{\nu_{p-1}(i_1-p,i_2-p) - 2\nu_{p-1}(i_1-p,i_2) + \nu_{p-1}(i_1-p,i_2+p) - 2\nu_{p-1}(i_1,i_2-p)}{4} +
332/// \f]
333/// \f[
334/// \frac{-2\nu_{p-1}(i_1,i_2+p) + \nu_{p-1}(i_1+p,i_2-p) - 2\nu_{p-1}(i_1+p,i_2) + \nu_{p-1}(i_1+p,i_2+p)}{4}
335/// \f]
336/// \f[ \nu_{p-1}(i_1,i_2) = min(a,b)\f]
337/// where p = 1, 2, ..., number_of_iterations.
338///
339/// #### References:
340///
341/// [1] C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the
342/// quantitative analysis of PIXE spectra in geoscience applications. NIM, B34 (1988), 396-402.
343///
344/// [2] M. Morhac;, J. Kliman, V.
345/// Matouoek, M. Veselsky, I. Turzo.:
346/// Background elimination methods for multidimensional gamma-ray spectra. NIM,
347/// A401 (1997) 113-132.
348///
349/// ### Example 1 - script Back_gamma64.c
350///
351/// \image html TSpectrum2_Background1.jpg Fig. 1 Original two-dimensional gamma-gamma-ray spectrum
352/// \image html TSpectrum2_Background2.jpg Fig. 2 Background estimated from data from Fig. 1 using decreasing clipping window with widths 4, 4 and algorithm based on successive comparisons. The estimate includes not only continuously changing background but also one-dimensional ridges.
353/// \image html TSpectrum2_Background3.jpg Fig. 3 Resulting peaks after subtraction of the estimated background (Fig. 2) from original two-dimensional gamma-gamma-ray spectrum (Fig. 1).
354///
355/// #### Script:
356///
357/// Example to illustrate the background estimator (class TSpectrum). To execute this example, do:
358///
359/// `root > .x Back_gamma64.C`
360///
361/// ~~~ {.cpp}
362/// #include <TSpectrum>
363/// void Back_gamma64() {
364/// Int_t i, j;
365/// Double_t nbinsx = 64;
366/// Double_t nbinsy = 64;
367/// Double_t xmin = 0;
368/// Double_t xmax = (Double_t)nbinsx;
369/// Double_t ymin = 0;
370/// Double_t ymax = (Double_t)nbinsy;
371/// Double_t ** source = new Double_t*[nbinsx];
372/// for (i=0;i<nbinsx;i++)
373/// source[i]=new Double_t[nbinsy];
374/// TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
375/// TFile *f = new TFile("spectra2/TSpectrum2.root");
376/// back=(TH2F*) f->Get("back1;1");
377/// TCanvas *Background = new TCanvas("Background","Estimation of background with increasing window",10,10,1000,700);
378/// TSpectrum *s = new TSpectrum();
379/// for (i = 0; i < nbinsx; i++){
380/// for (j = 0; j < nbinsy; j++){
381/// source[i][j] = back->GetBinContent(i + 1,j + 1);
382/// }
383/// }
384/// s->Background(source,nbinsx,nbinsy,4,4,kBackDecreasingWindow,kBackSuccessiveFiltering);
385/// for (i = 0; i < nbinsx; i++){
386/// for (j = 0; j < nbinsy; j++)
387/// back->SetBinContent(i + 1,j + 1, source[i][j]);
388/// }
389/// back->Draw("SURF");
390/// }
391/// ~~~
392///
393/// ### Example 2- script Back_gamma256.c
394///
395/// \image html TSpectrum2_Background4.jpg Fig. 4 Original two-dimensional gamma-gamma-ray spectrum 256x256 channels
396/// \image html TSpectrum2_Background5.jpg Fig. 5 Peaks after subtraction of the estimated background (increasing clipping window with widths 8, 8 and algorithm based on successive comparisons) from original two-dimensional gamma-gamma-ray spectrum (Fig. 4).
397///
398/// #### Script:
399///
400// Example to illustrate the background estimator (class TSpectrum).
401/// To execute this example, do
402///
403/// `root > .x Back_gamma256.C`
404///
405/// ~~~ {.cpp}
406/// #include <TSpectrum>
407/// void Back_gamma256() {
408/// Int_t i, j;
409/// Double_t nbinsx = 64;
410/// Double_t nbinsy = 64;
411/// Double_t xmin = 0;
412/// Double_t xmax = (Double_t)nbinsx;
413/// Double_t ymin = 0;
414/// Double_t ymax = (Double_t)nbinsy;
415/// Double_t** source = new Double_t*[nbinsx];
416/// for (i=0;i<nbinsx;i++)
417/// source[i]=new Double_t[nbinsy];
418/// TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
419/// TFile *f = new TFile("spectra2/TSpectrum2.root");
420/// back=(TH2F*) f->Get("back2;1");
421/// TCanvas *Background = new TCanvas("Background","Estimation of background with increasing window",10,10,1000,700);
422/// TSpectrum *s = new TSpectrum();
423/// for (i = 0; i < nbinsx; i++){
424/// for (j = 0; j < nbinsy; j++){
425/// source[i][j] = back->GetBinContent(i + 1,j + 1);
426/// }
427/// }
428/// s->Background(source,nbinsx,nbinsy,8,8,kBackIncreasingWindow,kBackSuccessiveFiltering);
429/// for (i = 0; i < nbinsx; i++){
430/// for (j = 0; j < nbinsy; j++)
431/// back->SetBinContent(i + 1,j + 1, source[i][j]);
432/// }
433/// back->Draw("SURF");
434/// }
435/// ~~~
436///
437/// Example 3- script Back_synt256.c
438///
439/// \image html TSpectrum2_Background6.jpg Fig. 6 Original two-dimensional synthetic spectrum 256x256 channels
440/// \image html TSpectrum2_Background7.jpg Fig. 7 Peaks after subtraction of the estimated background (increasing clipping window with widths 8, 8 and algorithm based on successive comparisons) from original two-dimensional gamma-gamma-ray spectrum (Fig. 6). One can observe artefacts (ridges) around the peaks due to the employed algorithm.
441/// \image html TSpectrum2_Background8.jpg Fig. 8 Peaks after subtraction of the estimated background (increasing clipping window with widths 8, 8 and algorithm based on one step filtering) from original two-dimensional gamma-gamma-ray spectrum (Fig. 6). The artefacts from the above given Fig. 7 disappeared.
442///
443/// #### Script:
444///
445/// Example to illustrate the background estimator (class TSpectrum).
446/// To execute this example, do
447///
448/// `root > .x Back_synt256.C`
449///
450/// ~~~ {.cpp}
451/// #include <TSpectrum>
452/// void Back_synt256() {
453/// Int_t i, j;
454/// Double_t nbinsx = 64;
455/// Double_t nbinsy = 64;
456/// Double_t xmin = 0;
457/// Double_t xmax = (Double_t)nbinsx;
458/// Double_t ymin = 0;
459/// Double_t ymax = (Double_t)nbinsy;
460/// Double_t** source = new Double_t*[nbinsx];
461/// for (i=0;i<nbinsx;i++)
462/// source[i]=new Double_t[nbinsy];
463/// TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
464/// TFile *f = new TFile("spectra2/TSpectrum2.root");
465/// back=(TH2F*) f->Get("back3;1");
466/// TCanvas *Background = new TCanvas("Background","Estimation of background with increasing window",10,10,1000,700);
467/// TSpectrum *s = new TSpectrum();
468/// for (i = 0; i < nbinsx; i++){
469/// for (j = 0; j < nbinsy; j++){
470/// source[i][j] = back->GetBinContent(i + 1,j + 1);
471/// }
472/// }
473/// s->Background(source,nbinsx,nbinsy,8,8,
474/// kBackIncreasingWindow,kBackSuccessiveFiltering);//kBackOneStepFiltering
475/// for (i = 0; i < nbinsx; i++){
476/// for (j = 0; j < nbinsy; j++)
477/// back->SetBinContent(i + 1,j + 1, source[i][j]);
478/// }
479/// back->Draw("SURF");
480/// }
481/// ~~~
482
483const char *TSpectrum2::Background(Double_t **spectrum,
484 Int_t ssizex, Int_t ssizey,
485 Int_t numberIterationsX,
486 Int_t numberIterationsY,
487 Int_t direction,
488 Int_t filterType)
489{
490 Int_t i, x, y, sampling, r1, r2;
491 Double_t a, b, p1, p2, p3, p4, s1, s2, s3, s4;
492 if (ssizex <= 0 || ssizey <= 0)
493 return "Wrong parameters";
494 if (numberIterationsX < 1 || numberIterationsY < 1)
495 return "Width of Clipping Window Must Be Positive";
496 if (ssizex < 2 * numberIterationsX + 1
497 || ssizey < 2 * numberIterationsY + 1)
498 return ("Too Large Clipping Window");
499 Double_t **working_space = new Double_t*[ssizex];
500 for (i = 0; i < ssizex; i++)
501 working_space[i] = new Double_t[ssizey];
502 sampling =
503 (Int_t) TMath::Max(numberIterationsX, numberIterationsY);
504 if (direction == kBackIncreasingWindow) {
505 if (filterType == kBackSuccessiveFiltering) {
506 for (i = 1; i <= sampling; i++) {
507 r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
508 (Int_t) TMath::Min(i, numberIterationsY);
509 for (y = r2; y < ssizey - r2; y++) {
510 for (x = r1; x < ssizex - r1; x++) {
511 a = spectrum[x][y];
512 p1 = spectrum[x - r1][y - r2];
513 p2 = spectrum[x - r1][y + r2];
514 p3 = spectrum[x + r1][y - r2];
515 p4 = spectrum[x + r1][y + r2];
516 s1 = spectrum[x][y - r2];
517 s2 = spectrum[x - r1][y];
518 s3 = spectrum[x + r1][y];
519 s4 = spectrum[x][y + r2];
520 b = (p1 + p2) / 2.0;
521 if (b > s2)
522 s2 = b;
523 b = (p1 + p3) / 2.0;
524 if (b > s1)
525 s1 = b;
526 b = (p2 + p4) / 2.0;
527 if (b > s4)
528 s4 = b;
529 b = (p3 + p4) / 2.0;
530 if (b > s3)
531 s3 = b;
532 s1 = s1 - (p1 + p3) / 2.0;
533 s2 = s2 - (p1 + p2) / 2.0;
534 s3 = s3 - (p3 + p4) / 2.0;
535 s4 = s4 - (p2 + p4) / 2.0;
536 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
537 p3 +
538 p4) / 4.0;
539 if (b < a && b > 0)
540 a = b;
541 working_space[x][y] = a;
542 }
543 }
544 for (y = r2; y < ssizey - r2; y++) {
545 for (x = r1; x < ssizex - r1; x++) {
546 spectrum[x][y] = working_space[x][y];
547 }
548 }
549 }
550 }
551
552 else if (filterType == kBackOneStepFiltering) {
553 for (i = 1; i <= sampling; i++) {
554 r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
555 (Int_t) TMath::Min(i, numberIterationsY);
556 for (y = r2; y < ssizey - r2; y++) {
557 for (x = r1; x < ssizex - r1; x++) {
558 a = spectrum[x][y];
559 b = -(spectrum[x - r1][y - r2] +
560 spectrum[x - r1][y + r2] + spectrum[x + r1][y -
561 r2]
562 + spectrum[x + r1][y + r2]) / 4 +
563 (spectrum[x][y - r2] + spectrum[x - r1][y] +
564 spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
565 if (b < a && b > 0)
566 a = b;
567 working_space[x][y] = a;
568 }
569 }
570 for (y = i; y < ssizey - i; y++) {
571 for (x = i; x < ssizex - i; x++) {
572 spectrum[x][y] = working_space[x][y];
573 }
574 }
575 }
576 }
577 }
578
579 else if (direction == kBackDecreasingWindow) {
580 if (filterType == kBackSuccessiveFiltering) {
581 for (i = sampling; i >= 1; i--) {
582 r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
583 (Int_t) TMath::Min(i, numberIterationsY);
584 for (y = r2; y < ssizey - r2; y++) {
585 for (x = r1; x < ssizex - r1; x++) {
586 a = spectrum[x][y];
587 p1 = spectrum[x - r1][y - r2];
588 p2 = spectrum[x - r1][y + r2];
589 p3 = spectrum[x + r1][y - r2];
590 p4 = spectrum[x + r1][y + r2];
591 s1 = spectrum[x][y - r2];
592 s2 = spectrum[x - r1][y];
593 s3 = spectrum[x + r1][y];
594 s4 = spectrum[x][y + r2];
595 b = (p1 + p2) / 2.0;
596 if (b > s2)
597 s2 = b;
598 b = (p1 + p3) / 2.0;
599 if (b > s1)
600 s1 = b;
601 b = (p2 + p4) / 2.0;
602 if (b > s4)
603 s4 = b;
604 b = (p3 + p4) / 2.0;
605 if (b > s3)
606 s3 = b;
607 s1 = s1 - (p1 + p3) / 2.0;
608 s2 = s2 - (p1 + p2) / 2.0;
609 s3 = s3 - (p3 + p4) / 2.0;
610 s4 = s4 - (p2 + p4) / 2.0;
611 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
612 p3 +
613 p4) / 4.0;
614 if (b < a && b > 0)
615 a = b;
616 working_space[x][y] = a;
617 }
618 }
619 for (y = r2; y < ssizey - r2; y++) {
620 for (x = r1; x < ssizex - r1; x++) {
621 spectrum[x][y] = working_space[x][y];
622 }
623 }
624 }
625 }
626
627 else if (filterType == kBackOneStepFiltering) {
628 for (i = sampling; i >= 1; i--) {
629 r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
630 (Int_t) TMath::Min(i, numberIterationsY);
631 for (y = r2; y < ssizey - r2; y++) {
632 for (x = r1; x < ssizex - r1; x++) {
633 a = spectrum[x][y];
634 b = -(spectrum[x - r1][y - r2] +
635 spectrum[x - r1][y + r2] + spectrum[x + r1][y -
636 r2]
637 + spectrum[x + r1][y + r2]) / 4 +
638 (spectrum[x][y - r2] + spectrum[x - r1][y] +
639 spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
640 if (b < a && b > 0)
641 a = b;
642 working_space[x][y] = a;
643 }
644 }
645 for (y = i; y < ssizey - i; y++) {
646 for (x = i; x < ssizex - i; x++) {
647 spectrum[x][y] = working_space[x][y];
648 }
649 }
650 }
651 }
652 }
653 for (i = 0; i < ssizex; i++)
654 delete[]working_space[i];
655 delete[]working_space;
656 return 0;
657}
658
659////////////////////////////////////////////////////////////////////////////////
660/// This function calculates smoothed spectrum from source spectrum
661/// based on Markov chain method.
662/// The result is placed in the array pointed by source pointer.
663///
664/// Function parameters:
665/// - source-pointer to the array of source spectrum
666/// - ssizex-x length of source
667/// - ssizey-y length of source
668/// - averWindow-width of averaging smoothing window
669///
670/// ### Smoothing
671///
672/// Goal: Suppression of statistical fluctuations the algorithm is based on discrete
673/// Markov chain, which has very simple invariant distribution
674/// \f[
675/// 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
676/// \f]
677/// \f$U_1\f$ being defined from the normalization condition \f$ \sum_{i=1}^{n} U_i = 1\f$
678/// n is the length of the smoothed spectrum and
679/// \f[
680/// 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]
681/// \f]
682/// is the probability of the change of the peak position from channel i to the channel i+1.
683/// \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
684/// of smoothing window. We have extended this algorithm to two dimensions.
685///
686/// #### Reference:
687///
688/// [1] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 (1996), 451.
689///
690/// ### Example 4 - script Smooth.c
691///
692/// \image html TSpectrum2_Smoothing1.jpg Fig. 9 Original noisy spectrum.
693/// \image html TSpectrum2_Smoothing2.jpg Fig. 10 Smoothed spectrum m=3 Peaks can hardly be observed. Peaks become apparent.
694/// \image html TSpectrum2_Smoothing3.jpg Fig. 11 Smoothed spectrum m=5
695/// \image html TSpectrum2_Smoothing4.jpg Fig.12 Smoothed spectrum m=7
696///
697/// #### Script:
698///
699/// Example to illustrate the Markov smoothing (class TSpectrum).
700/// To execute this example, do
701///
702/// `root > .x Smooth.C`
703///
704/// ~~~ {.cpp}
705/// #include <TSpectrum>
706/// void Smooth() {
707/// Int_t i, j;
708/// Double_t nbinsx = 256;
709/// Double_t nbinsy = 256;
710/// Double_t xmin = 0;
711/// Double_t xmax = (Double_t)nbinsx;
712/// Double_t ymin = 0;
713/// Double_t ymax = (Double_t)nbinsy;
714/// Double_t** source = new Double_t*[nbinsx];
715/// for (i=0;i<nbinsx;i++)
716/// source[i] = new Double_t[nbinsy];
717/// TH2F *smooth = new TH2F("smooth","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
718/// TFile *f = new TFile("spectra2/TSpectrum2.root");
719/// smooth=(TH2F*) f->Get("smooth1;1");
720/// TCanvas *Smoothing = new TCanvas("Smoothing","Markov smoothing",10,10,1000,700);
721/// TSpectrum *s = new TSpectrum();
722/// for (i = 0; i < nbinsx; i++){
723/// for (j = 0; j < nbinsy; j++){
724/// source[i][j] = smooth->GetBinContent(i + 1,j + 1);
725/// }
726/// }
727/// s->SmoothMarkov(source,nbinsx,nbinsx,3); //5,7
728/// for (i = 0; i < nbinsx; i++){
729/// for (j = 0; j < nbinsy; j++)
730/// smooth->SetBinContent(i + 1,j + 1, source[i][j]);
731/// }
732/// smooth->Draw("SURF");
733/// }
734/// ~~~
735
736const char* TSpectrum2::SmoothMarkov(Double_t **source, Int_t ssizex, Int_t ssizey, Int_t averWindow)
737{
738
739 Int_t xmin, xmax, ymin, ymax, i, j, l;
740 Double_t a, b, maxch;
741 Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy, plocha = 0;
742 if(averWindow <= 0)
743 return "Averaging Window must be positive";
744 Double_t **working_space = new Double_t*[ssizex];
745 for(i = 0; i < ssizex; i++)
746 working_space[i] = new Double_t[ssizey];
747 xmin = 0;
748 xmax = ssizex - 1;
749 ymin = 0;
750 ymax = ssizey - 1;
751 for(i = 0, maxch = 0; i < ssizex; i++){
752 for(j = 0; j < ssizey; j++){
753 working_space[i][j] = 0;
754 if(maxch < source[i][j])
755 maxch = source[i][j];
756
757 plocha += source[i][j];
758 }
759 }
760 if(maxch == 0) {
761 delete [] working_space;
762 return 0;
763 }
764
765 nom = 0;
766 working_space[xmin][ymin] = 1;
767 for(i = xmin; i < xmax; i++){
768 nip = source[i][ymin] / maxch;
769 nim = source[i + 1][ymin] / maxch;
770 sp = 0,sm = 0;
771 for(l = 1; l <= averWindow; l++){
772 if((i + l) > xmax)
773 a = source[xmax][ymin] / maxch;
774
775 else
776 a = source[i + l][ymin] / maxch;
777 b = a - nip;
778 if(a + nip <= 0)
779 a = 1;
780
781 else
782 a = TMath::Sqrt(a + nip);
783 b = b / a;
784 b = TMath::Exp(b);
785 sp = sp + b;
786 if(i - l + 1 < xmin)
787 a = source[xmin][ymin] / maxch;
788
789 else
790 a = source[i - l + 1][ymin] / maxch;
791 b = a - nim;
792 if(a + nim <= 0)
793 a = 1;
794
795 else
796 a = TMath::Sqrt(a + nim);
797 b = b / a;
798 b = TMath::Exp(b);
799 sm = sm + b;
800 }
801 a = sp / sm;
802 a = working_space[i + 1][ymin] = a * working_space[i][ymin];
803 nom = nom + a;
804 }
805 for(i = ymin; i < ymax; i++){
806 nip = source[xmin][i] / maxch;
807 nim = source[xmin][i + 1] / maxch;
808 sp = 0,sm = 0;
809 for(l = 1; l <= averWindow; l++){
810 if((i + l) > ymax)
811 a = source[xmin][ymax] / maxch;
812
813 else
814 a = source[xmin][i + l] / maxch;
815 b = a - nip;
816 if(a + nip <= 0)
817 a = 1;
818
819 else
820 a = TMath::Sqrt(a + nip);
821 b = b / a;
822 b = TMath::Exp(b);
823 sp = sp + b;
824 if(i - l + 1 < ymin)
825 a = source[xmin][ymin] / maxch;
826
827 else
828 a = source[xmin][i - l + 1] / maxch;
829 b = a - nim;
830 if(a + nim <= 0)
831 a = 1;
832
833 else
834 a = TMath::Sqrt(a + nim);
835 b = b / a;
836 b = TMath::Exp(b);
837 sm = sm + b;
838 }
839 a = sp / sm;
840 a = working_space[xmin][i + 1] = a * working_space[xmin][i];
841 nom = nom + a;
842 }
843 for(i = xmin; i < xmax; i++){
844 for(j = ymin; j < ymax; j++){
845 nip = source[i][j + 1] / maxch;
846 nim = source[i + 1][j + 1] / maxch;
847 spx = 0,smx = 0;
848 for(l = 1; l <= averWindow; l++){
849 if(i + l > xmax)
850 a = source[xmax][j] / maxch;
851
852 else
853 a = source[i + l][j] / maxch;
854 b = a - nip;
855 if(a + nip <= 0)
856 a = 1;
857
858 else
859 a = TMath::Sqrt(a + nip);
860 b = b / a;
861 b = TMath::Exp(b);
862 spx = spx + b;
863 if(i - l + 1 < xmin)
864 a = source[xmin][j] / maxch;
865
866 else
867 a = source[i - l + 1][j] / maxch;
868 b = a - nim;
869 if(a + nim <= 0)
870 a = 1;
871
872 else
873 a = TMath::Sqrt(a + nim);
874 b = b / a;
875 b = TMath::Exp(b);
876 smx = smx + b;
877 }
878 spy = 0,smy = 0;
879 nip = source[i + 1][j] / maxch;
880 nim = source[i + 1][j + 1] / maxch;
881 for (l = 1; l <= averWindow; l++) {
882 if (j + l > ymax) a = source[i][ymax]/maxch;
883 else a = source[i][j + l] / maxch;
884 b = a - nip;
885 if (a + nip <= 0) a = 1;
886 else a = TMath::Sqrt(a + nip);
887 b = b / a;
888 b = TMath::Exp(b);
889 spy = spy + b;
890 if (j - l + 1 < ymin) a = source[i][ymin] / maxch;
891 else a = source[i][j - l + 1] / maxch;
892 b = a - nim;
893 if (a + nim <= 0) a = 1;
894 else a = TMath::Sqrt(a + nim);
895 b = b / a;
896 b = TMath::Exp(b);
897 smy = smy + b;
898 }
899 a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx +smy);
900 working_space[i + 1][j + 1] = a;
901 nom = nom + a;
902 }
903 }
904 for(i = xmin; i <= xmax; i++){
905 for(j = ymin; j <= ymax; j++){
906 working_space[i][j] = working_space[i][j] / nom;
907 }
908 }
909 for(i = 0;i < ssizex; i++){
910 for(j = 0; j < ssizey; j++){
911 source[i][j] = plocha * working_space[i][j];
912 }
913 }
914 for (i = 0; i < ssizex; i++)
915 delete[]working_space[i];
916 delete[]working_space;
917 return 0;
918}
919
920////////////////////////////////////////////////////////////////////////////////
921/// This function calculates deconvolution from source spectrum
922/// according to response spectrum
923/// The result is placed in the matrix pointed by source pointer.
924///
925/// Function parameters:
926/// - source-pointer to the matrix of source spectrum
927/// - resp-pointer to the matrix of response spectrum
928/// - ssizex-x length of source and response spectra
929/// - ssizey-y length of source and response spectra
930/// - numberIterations, for details we refer to manual
931/// - numberRepetitions, for details we refer to manual
932/// - boost, boosting factor, for details we refer to manual
933///
934/// ### Deconvolution
935///
936/// Goal: Improvement of the resolution in spectra, decomposition of multiplets
937///
938/// Mathematical formulation of the 2-dimensional convolution system is
939///
940///\f[
941/// y(i_1,i_2) = \sum_{k_1=0}^{N_1-1} \sum_{k_2=0}^{N_2-1} h(i_1-k_1,i_2-k_2)x(k_1,k_2)
942///\f]
943///\f[
944/// i_1 = 0,1,2, ... ,N_1-1, i_2 = 0,1,2, ... ,N_2-1
945///\f]
946///
947/// where h(i,j) is the impulse response function, x, y are input and output
948/// matrices, respectively, \f$ N_1, N_2\f$ are the lengths of x and h matrices
949///
950/// - let us assume that we know the response and the output matrices (spectra) of the above given system.
951/// - the deconvolution represents solution of the overdetermined system of linear equations, i.e., the
952/// calculation of the matrix x.
953/// - from numerical stability point of view the operation of deconvolution is
954/// extremely critical (ill-posed problem) as well as time consuming operation.
955/// - the Gold deconvolution algorithm proves to work very well even for 2-dimensional
956/// systems. Generalisation of the algorithm for 2-dimensional systems was presented in [1], [2].
957/// - for Gold deconvolution algorithm as well as for boosted deconvolution algorithm we refer also to TSpectrum
958///
959/// #### References:
960///
961/// [1] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:
962/// Efficient one- and two-dimensional Gold deconvolution and its application to
963/// gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
964///
965/// [2] Morhac; M., Matouoek V., Kliman J., Efficient algorithm of multidimensional
966/// deconvolution and its application to nuclear data processing, Digital Signal
967/// Processing 13 (2003) 144.
968///
969/// ### Example 5 - script Decon.c
970///
971/// response function (usually peak) should be shifted to the beginning of the coordinate system (see Fig. 13)
972///
973/// \image html TSpectrum2_Deconvolution1.jpg Fig. 13 2-dimensional response spectrum
974/// \image html TSpectrum2_Deconvolution2.jpg Fig. 14 2-dimensional gamma-gamma-ray input spectrum (before deconvolution)
975/// \image html TSpectrum2_Deconvolution3.jpg Fig. 15 Spectrum from Fig. 14 after deconvolution (1000 iterations)
976///
977/// #### Script:
978///
979/// Example to illustrate the Gold deconvolution (class TSpectrum2).
980/// To execute this example, do
981///
982/// `root > .x Decon.C`
983///
984/// ~~~ {.cpp}
985/// #include <TSpectrum2>
986/// void Decon() {
987/// Int_t i, j;
988/// Double_t nbinsx = 256;
989/// Double_t nbinsy = 256;
990/// Double_t xmin = 0;
991/// Double_t xmax = (Double_t)nbinsx;
992/// Double_t ymin = 0;
993/// Double_t ymax = (Double_t)nbinsy;
994/// Double_t** source = new Double_t*[nbinsx];
995/// for (i=0;i<nbinsx;i++)
996/// source[i]=new Double_t[nbinsy];
997/// TH2F *decon = new TH2F("decon","Gold deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
998/// TFile *f = new TFile("spectra2/TSpectrum2.root");
999/// decon=(TH2F*) f->Get("decon1;1");
1000/// Double_t** response = new Double_t*[nbinsx];
1001/// for (i=0;i<nbinsx;i++)
1002/// response[i]=new Double_t[nbinsy];
1003/// TH2F *resp = new TH2F("resp","Response matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1004/// resp=(TH2F*) f->Get("resp1;1");
1005/// TCanvas *Deconvol = new TCanvas("Deconvolution","Gold deconvolution",10,10,1000,700);
1006/// TSpectrum *s = new TSpectrum();
1007/// for (i = 0; i < nbinsx; i++){
1008/// for (j = 0; j < nbinsy; j++){
1009/// source[i][j] = decon->GetBinContent(i + 1,j + 1);
1010/// }
1011/// }
1012/// for (i = 0; i < nbinsx; i++){
1013/// for (j = 0; j < nbinsy; j++){
1014/// response[i][j] = resp->GetBinContent(i + 1,j + 1);
1015/// }
1016/// }
1017/// s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
1018/// for (i = 0; i < nbinsx; i++){
1019/// for (j = 0; j < nbinsy; j++)
1020/// decon->SetBinContent(i + 1,j + 1, source[i][j]);
1021/// }
1022/// decon->Draw("SURF");
1023/// }
1024/// ~~~
1025///
1026/// ### Example 6 - script Decon2.c
1027///
1028/// \image html TSpectrum2_Deconvolution4.jpg Fig. 16 Response spectrum
1029/// \image html TSpectrum2_Deconvolution5.jpg Fig. 17 Original synthetic input spectrum (before deconvolution). It is composed of 17 peaks. 5 peaks are overlapping in the outlined multiplet and two peaks in doublet.
1030/// \image html TSpectrum2_Deconvolution6.jpg Fig. 18 Spectrum from Fig. 17 after deconvolution (1000 iterations). Resolution is improved but the peaks in multiplet remained unresolved.
1031///
1032/// #### Script:
1033///
1034/// Example to illustrate the Gold deconvolution (class TSpectrum2).
1035/// To execute this example, do
1036///
1037/// `root > .x Decon2.C`
1038///
1039/// ~~~ {.cpp}
1040/// #include <TSpectrum2>
1041/// void Decon2() {
1042/// Int_t i, j;
1043/// Double_t nbinsx = 64;
1044/// Double_t nbinsy = 64;
1045/// Double_t xmin = 0;
1046/// Double_t xmax = (Double_t)nbinsx;
1047/// Double_t ymin = 0;
1048/// Double_t ymax = (Double_t)nbinsy;
1049/// Double_t** source = new Double_t*[nbinsx];
1050/// for (i=0;i<nbinsx;i++)
1051/// source[i]=new Double_t[nbinsy];
1052/// TH2F *decon = new TH2F("decon","Gold deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1053/// TFile *f = new TFile("spectra2/TSpectrum2.root");
1054/// decon=(TH2F*) f->Get("decon2;1");
1055/// Double_t** response = new Double_t*[nbinsx];
1056/// for (i=0;i<nbinsx;i++)
1057/// response[i]=new Double_t[nbinsy];
1058/// TH2F *resp = new TH2F("resp","Response matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1059/// resp=(TH2F*) f->Get("resp2;1");
1060/// TCanvas *Deconvol = new TCanvas("Deconvolution","Gold deconvolution",10,10,1000,700);
1061/// TSpectrum *s = new TSpectrum();
1062/// for (i = 0; i < nbinsx; i++){
1063/// for (j = 0; j < nbinsy; j++){
1064/// source[i][j] = decon->GetBinContent(i + 1,j + 1);
1065/// }
1066/// }
1067/// for (i = 0; i < nbinsx; i++){
1068/// for (j = 0; j < nbinsy; j++){
1069/// response[i][j] = resp->GetBinContent(i + 1,j + 1);
1070/// }
1071/// }
1072/// s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
1073/// for (i = 0; i < nbinsx; i++){
1074/// for (j = 0; j < nbinsy; j++)
1075/// decon->SetBinContent(i + 1,j + 1, source[i][j]);
1076/// }
1077/// decon->Draw("SURF");
1078/// }
1079/// ~~~
1080///
1081/// ### Example 7 - script Decon2HR.c
1082///
1083/// \image html TSpectrum2_Deconvolution7.jpg Fig. 19 Spectrum from Fig. 17 after boosted deconvolution (50 iterations repeated 20 times, boosting coefficient was 1.2). All the peaks in multiplet as well as in doublet are completely decomposed.
1084///
1085/// #### Script:
1086///
1087/// Example to illustrate boosted Gold deconvolution (class TSpectrum2).
1088/// To execute this example, do
1089///
1090/// `root > .x Decon2HR.C`
1091///
1092/// ~~~ {.cpp}
1093/// #include <TSpectrum2>
1094/// void Decon2HR() {
1095/// Int_t i, j;
1096/// Double_t nbinsx = 64;
1097/// Double_t nbinsy = 64;
1098/// Double_t xmin = 0;
1099/// Double_t xmax = (Double_t)nbinsx;
1100/// Double_t ymin = 0;
1101/// Double_t ymax = (Double_t)nbinsy;
1102/// Double_t** source = new Double_t*[nbinsx];
1103/// for (i=0;i<nbinsx;i++)
1104/// source[i]=new Double_t[nbinsy];
1105/// TH2F *decon = new TH2F("decon","Boosted Gold deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1106/// TFile *f = new TFile("spectra2/TSpectrum2.root");
1107/// decon=(TH2F*) f->Get("decon2;1");
1108/// Double_t** response = new Double_t*[nbinsx];
1109/// for (i=0;i<nbinsx;i++)
1110/// response[i]=new Double_t[nbinsy];
1111/// TH2F *resp = new TH2F("resp","Response matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1112/// resp=(TH2F*) f->Get("resp2;1");
1113/// TCanvas *Deconvol = new TCanvas("Deconvolution","Gold deconvolution",10,10,1000,700);
1114/// TSpectrum *s = new TSpectrum();
1115/// for (i = 0; i < nbinsx; i++){
1116/// for (j = 0; j < nbinsy; j++){
1117/// source[i][j] = decon->GetBinContent(i + 1,j + 1);
1118/// }
1119/// }
1120/// for (i = 0; i < nbinsx; i++){
1121/// for (j = 0; j < nbinsy; j++){
1122/// response[i][j] = resp->GetBinContent(i + 1,j + 1);
1123/// }
1124/// }
1125/// s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
1126/// for (i = 0; i < nbinsx; i++){
1127/// for (j = 0; j < nbinsy; j++)
1128/// dec on->SetBinContent(i + 1,j + 1, source[i][j]);
1129/// }
1130/// decon->Draw("SURF");
1131/// }
1132/// ~~~
1133
1134const char *TSpectrum2::Deconvolution(Double_t **source, Double_t **resp,
1135 Int_t ssizex, Int_t ssizey,
1136 Int_t numberIterations,
1137 Int_t numberRepetitions,
1138 Double_t boost)
1139{
1140 Int_t i, j, lhx, lhy, i1, i2, j1, j2, k1, k2, lindex, i1min, i1max,
1141 i2min, i2max, j1min, j1max, j2min, j2max, positx = 0, posity = 0, repet;
1142 Double_t lda, ldb, ldc, area, maximum = 0;
1143 if (ssizex <= 0 || ssizey <= 0)
1144 return "Wrong parameters";
1145 if (numberIterations <= 0)
1146 return "Number of iterations must be positive";
1147 if (numberRepetitions <= 0)
1148 return "Number of repetitions must be positive";
1149 Double_t **working_space = new Double_t *[ssizex];
1150 for (i = 0; i < ssizex; i++)
1151 working_space[i] = new Double_t[5 * ssizey];
1152 area = 0;
1153 lhx = -1, lhy = -1;
1154 for (i = 0; i < ssizex; i++) {
1155 for (j = 0; j < ssizey; j++) {
1156 lda = resp[i][j];
1157 if (lda != 0) {
1158 if ((i + 1) > lhx)
1159 lhx = i + 1;
1160 if ((j + 1) > lhy)
1161 lhy = j + 1;
1162 }
1163 working_space[i][j] = lda;
1164 area = area + lda;
1165 if (lda > maximum) {
1166 maximum = lda;
1167 positx = i, posity = j;
1168 }
1169 }
1170 }
1171 if (lhx == -1 || lhy == -1) {
1172 delete [] working_space;
1173 return ("Zero response data");
1174 }
1175
1176//calculate ht*y and write into p
1177 for (i2 = 0; i2 < ssizey; i2++) {
1178 for (i1 = 0; i1 < ssizex; i1++) {
1179 ldc = 0;
1180 for (j2 = 0; j2 <= (lhy - 1); j2++) {
1181 for (j1 = 0; j1 <= (lhx - 1); j1++) {
1182 k2 = i2 + j2, k1 = i1 + j1;
1183 if (k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
1184 lda = working_space[j1][j2];
1185 ldb = source[k1][k2];
1186 ldc = ldc + lda * ldb;
1187 }
1188 }
1189 }
1190 working_space[i1][i2 + ssizey] = ldc;
1191 }
1192 }
1193
1194//calculate matrix b=ht*h
1195 i1min = -(lhx - 1), i1max = lhx - 1;
1196 i2min = -(lhy - 1), i2max = lhy - 1;
1197 for (i2 = i2min; i2 <= i2max; i2++) {
1198 for (i1 = i1min; i1 <= i1max; i1++) {
1199 ldc = 0;
1200 j2min = -i2;
1201 if (j2min < 0)
1202 j2min = 0;
1203 j2max = lhy - 1 - i2;
1204 if (j2max > lhy - 1)
1205 j2max = lhy - 1;
1206 for (j2 = j2min; j2 <= j2max; j2++) {
1207 j1min = -i1;
1208 if (j1min < 0)
1209 j1min = 0;
1210 j1max = lhx - 1 - i1;
1211 if (j1max > lhx - 1)
1212 j1max = lhx - 1;
1213 for (j1 = j1min; j1 <= j1max; j1++) {
1214 lda = working_space[j1][j2];
1215 if (i1 + j1 < ssizex && i2 + j2 < ssizey)
1216 ldb = working_space[i1 + j1][i2 + j2];
1217 else
1218 ldb = 0;
1219 ldc = ldc + lda * ldb;
1220 }
1221 }
1222 working_space[i1 - i1min][i2 - i2min + 2 * ssizey ] = ldc;
1223 }
1224 }
1225
1226//initialization in x1 matrix
1227 for (i2 = 0; i2 < ssizey; i2++) {
1228 for (i1 = 0; i1 < ssizex; i1++) {
1229 working_space[i1][i2 + 3 * ssizey] = 1;
1230 working_space[i1][i2 + 4 * ssizey] = 0;
1231 }
1232 }
1233
1234 //START OF ITERATIONS
1235 for (repet = 0; repet < numberRepetitions; repet++) {
1236 if (repet != 0) {
1237 for (i = 0; i < ssizex; i++) {
1238 for (j = 0; j < ssizey; j++) {
1239 working_space[i][j + 3 * ssizey] =
1240 TMath::Power(working_space[i][j + 3 * ssizey], boost);
1241 }
1242 }
1243 }
1244 for (lindex = 0; lindex < numberIterations; lindex++) {
1245 for (i2 = 0; i2 < ssizey; i2++) {
1246 for (i1 = 0; i1 < ssizex; i1++) {
1247 ldb = 0;
1248 j2min = i2;
1249 if (j2min > lhy - 1)
1250 j2min = lhy - 1;
1251 j2min = -j2min;
1252 j2max = ssizey - i2 - 1;
1253 if (j2max > lhy - 1)
1254 j2max = lhy - 1;
1255 j1min = i1;
1256 if (j1min > lhx - 1)
1257 j1min = lhx - 1;
1258 j1min = -j1min;
1259 j1max = ssizex - i1 - 1;
1260 if (j1max > lhx - 1)
1261 j1max = lhx - 1;
1262 for (j2 = j2min; j2 <= j2max; j2++) {
1263 for (j1 = j1min; j1 <= j1max; j1++) {
1264 ldc = working_space[j1 - i1min][j2 - i2min + 2 * ssizey];
1265 lda = working_space[i1 + j1][i2 + j2 + 3 * ssizey];
1266 ldb = ldb + lda * ldc;
1267 }
1268 }
1269 lda = working_space[i1][i2 + 3 * ssizey];
1270 ldc = working_space[i1][i2 + 1 * ssizey];
1271 if (ldc * lda != 0 && ldb != 0) {
1272 lda = lda * ldc / ldb;
1273 }
1274
1275 else
1276 lda = 0;
1277 working_space[i1][i2 + 4 * ssizey] = lda;
1278 }
1279 }
1280 for (i2 = 0; i2 < ssizey; i2++) {
1281 for (i1 = 0; i1 < ssizex; i1++)
1282 working_space[i1][i2 + 3 * ssizey] =
1283 working_space[i1][i2 + 4 * ssizey];
1284 }
1285 }
1286 }
1287 for (i = 0; i < ssizex; i++) {
1288 for (j = 0; j < ssizey; j++)
1289 source[(i + positx) % ssizex][(j + posity) % ssizey] =
1290 area * working_space[i][j + 3 * ssizey];
1291 }
1292 for (i = 0; i < ssizex; i++)
1293 delete[]working_space[i];
1294 delete[]working_space;
1295 return 0;
1296}
1297
1298////////////////////////////////////////////////////////////////////////////////
1299/// This function searches for peaks in source spectrum
1300/// It is based on deconvolution method. First the background is
1301/// removed (if desired), then Markov spectrum is calculated
1302/// (if desired), then the response function is generated
1303/// according to given sigma and deconvolution is carried out.
1304///
1305/// Function parameters:
1306/// - source-pointer to the matrix of source spectrum
1307/// - dest-pointer to the matrix of resulting deconvolved spectrum
1308/// - ssizex-x length of source spectrum
1309/// - ssizey-y length of source spectrum
1310/// - sigma-sigma of searched peaks, for details we refer to manual
1311/// - threshold-threshold value in % for selected peaks, peaks with
1312/// amplitude less than threshold*highest_peak/100
1313/// are ignored, see manual
1314/// - backgroundRemove-logical variable, set if the removal of
1315/// background before deconvolution is desired
1316/// - deconIterations-number of iterations in deconvolution operation
1317/// - markov-logical variable, if it is true, first the source spectrum
1318/// is replaced by new spectrum calculated using Markov
1319/// chains method.
1320/// - averWindow-averaging window of searched peaks, for details
1321/// we refer to manual (applies only for Markov method)
1322///
1323/// ### Peaks searching
1324///
1325/// Goal: to identify automatically the peaks in spectrum with the presence of the
1326/// continuous background, one-fold coincidences (ridges) and statistical
1327/// fluctuations - noise.
1328///
1329/// The common problems connected with correct peak identification in two-dimensional coincidence spectra are
1330///
1331/// - non-sensitivity to noise, i.e., only statistically relevant peaks should be identified
1332/// - non-sensitivity of the algorithm to continuous background
1333/// - non-sensitivity to one-fold coincidences (coincidences peak - background in both dimensions) and their crossings
1334/// - ability to identify peaks close to the edges of the spectrum region. Usually peak finders fail to detect them
1335/// - resolution, decomposition of doublets and multiplets. The algorithm should be able to recognise close positioned peaks.
1336/// - ability to identify peaks with different sigma
1337///
1338/// #### References:
1339///
1340/// [1] M.A. Mariscotti: A method for identification of peaks in the presence of
1341/// background and its application to spectrum analysis. NIM 50 (1967), 309-320.
1342///
1343/// [2] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:Identification
1344/// of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
1345/// 108-125.
1346///
1347/// [3] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
1348/// (1996), 451.
1349///
1350/// ### Examples of peak searching method
1351///
1352/// SearchHighRes function provides users with the possibility
1353/// to vary the input parameters and with the access to the output deconvolved data
1354/// in the destination spectrum. Based on the output data one can tune the
1355/// parameters.
1356///
1357/// ### Example 8 - script Src.c
1358///
1359/// \image html TSpectrum2_Searching1.jpg Fig. 20 Two-dimensional spectrum with found peaks denoted by markers (sigma=2, threshold=5%, 3 iterations steps in the deconvolution)
1360/// \image html TSpectrum2_Searching3.jpg Fig. 21 Spectrum from Fig. 20 after background elimination and deconvolution
1361///
1362/// #### Script:
1363///
1364/// Example to illustrate high resolution peak searching function (class TSpectrum).
1365/// To execute this example, do
1366///
1367/// `root > .x Src.C
1368///
1369/// ~~~ {.cpp}
1370/// #include <TSpectrum2>
1371/// void Src() {
1372/// Int_t i, j, nfound;
1373/// Double_t nbinsx = 64;
1374/// Double_t nbinsy = 64;
1375/// Double_t xmin = 0;
1376/// Double_t xmax = (Double_t)nbinsx;
1377/// Double_t ymin = 0;
1378/// Double_t ymax = (Double_t)nbinsy;
1379/// Double_t** source = new Double_t*[nbinsx];
1380/// for (i=0;i<nbinsx;i++)
1381/// source[i]=new Double_t[nbinsy];
1382/// Double_t** dest = new Double_t*[nbinsx];
1383/// for (i=0;i<nbinsx;i++)
1384/// dest[i]=new Double_t[nbinsy];
1385/// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);/
1386/// TFile *f = new TFile("spectra2/TSpectrum2.root");
1387/// search=(TH2F*) f->Get("search4;1");
1388/// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1389/// TSpectrum2 *s = new TSpectrum2();
1390/// for (i = 0; i < nbinsx; i++){
1391/// for (j = 0; j < nbinsy; j++){
1392/// source[i][j] = search->GetBinContent(i + 1,j + 1);
1393/// }
1394/// }
1395/// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);
1396/// printf("Found %d candidate peaks\n",nfound);
1397/// for(i=0;i<nfound;i++)
1398/// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1399/// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1400/// }
1401/// ~~~
1402///
1403/// ### Example 9 - script Src2.c
1404///
1405/// \image html TSpectrum2_Searching4.jpg Fig. 22 Two-dimensional noisy spectrum with found peaks denoted by markers (sigma=2, threshold=10%, 10 iterations steps in the deconvolution). One can observe that the algorithm is insensitive to the crossings of one-dimensional ridges. It identifies only two-coincidence peaks.
1406/// \image html TSpectrum2_Searching5.jpg Fig. 23 Spectrum from Fig. 22 after background elimination and deconvolution
1407///
1408/// #### Script:
1409///
1410/// Example to illustrate high resolution peak searching function (class TSpectrum).
1411/// To execute this example, do
1412///
1413/// `root > .x Src2.C`
1414///
1415/// ~~~ {.cpp}
1416/// #include <TSpectrum2>
1417/// void Src2() {
1418/// Int_t i, j, nfound;
1419/// Double_t nbinsx = 256;
1420/// Double_t nbinsy = 256;
1421/// Double_t xmin = 0;
1422/// Double_t xmax = (Double_t)nbinsx;
1423/// Double_t ymin = 0;
1424/// Double_t ymax = (Double_t)nbinsy;
1425/// Double_t** source = new Double_t*[nbinsx];
1426/// for (i=0;i<nbinsx;i++)
1427/// source[i]=new Double_t[nbinsy];
1428/// Double_t** dest = new Double_t*[nbinsx];
1429/// for (i=0;i<nbinsx;i++)
1430/// dest[i]=new Double_t[nbinsy];
1431/// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1432/// TFile *f = new TFile("spectra2/TSpectrum2.root");
1433/// search=(TH2F*) f->Get("back3;1");
1434/// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1435/// TSpectrum2 *s = new TSpectrum2();
1436/// for (i = 0; i < nbinsx; i++){
1437/// for (j = 0; j < nbinsy; j++){
1438/// source[i][j] = search->GetBinContent(i + 1,j + 1);
1439/// }
1440/// }
1441/// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 10, kTRUE, 10, kFALSE, 3);
1442/// printf("Found %d candidate peaks\n",nfound);
1443/// for(i=0;i<nfound;i++)
1444/// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1445/// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1446/// }
1447/// ~~~
1448///
1449/// ### Example 10 - script Src3.c
1450///
1451/// \image html TSpectrum2_Searching6.jpg Fig. 24 Two-dimensional spectrum with 15 found peaks denoted by markers. Some peaks are positioned close to each other. It is necessary to increase number of iterations in the deconvolution. In next 3 Figs. we shall study the influence of this parameter.
1452/// \image html TSpectrum2_Searching7.jpg Fig. 25 Spectrum from Fig. 24 after deconvolution (# of iterations = 3). Number of identified peaks = 13.
1453/// \image html TSpectrum2_Searching8.jpg Fig. 26 Spectrum from Fig. 24 after deconvolution (# of iterations = 10). Number of identified peaks = 13.
1454/// \image html TSpectrum2_Searching9.jpg Fig. 27 Spectrum from Fig. 24 after deconvolution (# of iterations = 100). Number of identified peaks = 15. Now the algorithm is able to decompose two doublets in the spectrum.
1455///
1456/// #### Script:
1457///
1458/// Example to illustrate high resolution peak searching function (class TSpectrum).
1459/// To execute this example, do
1460///
1461/// `root > .x Src3.C`
1462///
1463/// ~~~ {.cpp}
1464/// #include <TSpectrum2>
1465/// void Src3() {
1466/// Int_t i, j, nfound;
1467/// Double_t nbinsx = 64;
1468/// Double_t nbinsy = 64;
1469/// Double_t xmin = 0;
1470/// Double_t xmax = (Double_t)nbinsx;
1471/// Double_t ymin = 0;
1472/// Double_t ymax = (Double_t)nbinsy;
1473/// Double_t** source = new Double_t*[nbinsx];
1474/// for (i=0;i<nbinsx;i++)
1475/// source[i]=new Double_t[nbinsy];
1476/// Double_t** dest = new Double_t*[nbinsx];
1477/// for (i=0;i<nbinsx;i++)
1478/// dest[i]=new Double_t[nbinsy];
1479/// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1480/// TFile *f = new TFile("spectra2/TSpectrum2.root");
1481/// search=(TH2F*) f->Get("search1;1");
1482/// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1483/// TSpectrum2 *s = new TSpectrum2();
1484/// for (i = 0; i < nbinsx; i++){
1485/// for (j = 0; j < nbinsy; j++){
1486/// source[i][j] = search->GetBinContent(i + 1,j + 1);
1487/// }
1488/// }
1489/// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 2, kFALSE, 3, kFALSE, 1);//3, 10, 100
1490/// printf("Found %d candidate peaks\n",nfound);
1491/// for(i=0;i<nfound;i++)
1492/// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1493/// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1494/// }
1495/// ~~~
1496///
1497/// ### Example 11 - script Src4.c
1498///
1499/// \image html TSpectrum2_Searching10.jpg Fig. 28 Two-dimensional spectrum with peaks with different sigma denoted by markers (sigma=3, threshold=5%, 10 iterations steps in the deconvolution, Markov smoothing with window=3)
1500/// \image html TSpectrum2_Searching12.jpg Fig. 29 Spectrum from Fig. 28 after smoothing and deconvolution.
1501///
1502/// #### Script:
1503///
1504/// Example to illustrate high resolution peak searching function (class TSpectrum).
1505/// To execute this example, do
1506///
1507/// `root > .x Src4.C`
1508///
1509/// ~~~ {.cpp}
1510/// #include <TSpectrum2>
1511/// void Src4() {
1512/// Int_t i, j, nfound;
1513/// Double_t nbinsx = 64;
1514/// Double_t nbinsy = 64;
1515/// Double_t xmin = 0;
1516/// Double_t xmax = (Double_t)nbinsx;
1517/// Double_t ymin = 0;
1518/// Double_t ymax = (Double_t)nbinsy;
1519/// Double_t** source = new Double_t*[nbinsx];
1520/// for (i=0;i<nbinsx;i++)
1521/// source[i]=new Double_t[nbinsy];
1522/// Double_t** dest = new Double_t*[nbinsx];
1523/// for (i=0;i<nbinsx;i++)
1524/// dest[i]=new Double_t[nbinsy];
1525/// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1526/// TFile *f = new TFile("spectra2/TSpectrum2.root");
1527/// search=(TH2F*) f->Get("search2;1");
1528/// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1529/// TSpectrum2 *s = new TSpectrum2();
1530/// for (i = 0; i < nbinsx; i++){
1531/// for (j = 0; j < nbinsy; j++){
1532/// source[i][j] = search->GetBinContent(i + 1,j + 1);
1533/// }
1534/// }
1535/// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 3, 5, kFALSE, 10, kTRUE, 3);
1536/// printf("Found %d candidate peaks\n",nfound);
1537/// for(i=0;i<nfound;i++)
1538/// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1539/// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1540/// }
1541/// ~~~
1542///
1543/// ### Example 12 - script Src5.c`
1544///
1545/// \image html TSpectrum2_Searching13.jpg Fig. 30 Two-dimensional spectrum with peaks positioned close to the edges denoted by markers (sigma=2, threshold=5%, 10 iterations steps in the deconvolution)
1546/// \image html TSpectrum2_Searching14.jpg Fig. 31 Spectrum from Fig. 30 after deconvolution.
1547///
1548/// #### Script:
1549///
1550/// Example to illustrate high resolution peak searching function (class TSpectrum).
1551/// To execute this example, do
1552///
1553/// `root > .x Src5.C`
1554///
1555/// ~~~ {.cpp}
1556/// #include <TSpectrum2>
1557/// void Src5() {
1558/// Int_t i, j, nfound;
1559/// Double_t nbinsx = 64;
1560/// Double_t nbinsy = 64;
1561/// Double_t xmin = 0;
1562/// Double_t xmax = (Double_t)nbinsx;
1563/// Double_t ymin = 0;
1564/// Double_t ymax = (Double_t)nbinsy;
1565/// Double_t** source = new Double_t*[nbinsx];
1566/// for (i=0;i<nbinsx;i++)
1567/// source[i]=new Double_t[nbinsy];
1568/// Double_t** dest = new Double_t*[nbinsx];
1569/// for (i=0;i<nbinsx;i++)
1570/// dest[i]=new Double_t[nbinsy];
1571/// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1572/// TFile *f = new TFile("spectra2/TSpectrum2.root");
1573/// search=(TH2F*) f->Get("search3;1");
1574/// TCanvas *Searching = new TCanvas("Searching","High resolution peak searching",10,10,1000,700);
1575/// TSpectrum2 *s = new TSpectrum2();
1576/// for (i = 0; i < nbinsx; i++){
1577/// for (j = 0; j < nbinsy; j++){
1578/// source[i][j] = search->GetBinContent(i + 1,j + 1);
1579/// }
1580/// }
1581/// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kFALSE, 10, kFALSE, 1);
1582/// printf("Found %d candidate peaks\n",nfound);
1583/// for(i=0;i<nfound;i++)
1584/// printf("posx= %d, posy= %d, value=%d\n",(Int_t)(fPositionX[i]+0.5), (Int_t)(fPositionY[i]+0.5),
1585/// (Int_t)source[(Int_t)(fPositionX[i]+0.5)][(Int_t)(fPositionY[i]+0.5)]);
1586/// }
1587/// ~~~
1588
1589Int_t TSpectrum2::SearchHighRes(Double_t **source, Double_t **dest, Int_t ssizex, Int_t ssizey,
1590 Double_t sigma, Double_t threshold,
1591 Bool_t backgroundRemove,Int_t deconIterations,
1592 Bool_t markov, Int_t averWindow)
1593
1594{
1595 Int_t number_of_iterations = (Int_t)(4 * sigma + 0.5);
1596 Int_t k, lindex, priz;
1597 Double_t lda, ldb, ldc, area, maximum;
1598 Int_t xmin, xmax, l, peak_index = 0, ssizex_ext = ssizex + 4 * number_of_iterations, ssizey_ext = ssizey + 4 * number_of_iterations, shift = 2 * number_of_iterations;
1599 Int_t ymin, ymax, i, j;
1600 Double_t a, b, ax, ay, maxch, plocha = 0;
1601 Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy;
1602 Double_t p1, p2, p3, p4, s1, s2, s3, s4;
1603 Int_t x, y;
1604 Int_t lhx, lhy, i1, i2, j1, j2, k1, k2, i1min, i1max, i2min, i2max, j1min, j1max, j2min, j2max, positx, posity;
1605 if (sigma < 1) {
1606 Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
1607 return 0;
1608 }
1609
1610 if(threshold<=0||threshold>=100){
1611 Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
1612 return 0;
1613 }
1614
1615 j = (Int_t) (5.0 * sigma + 0.5);
1616 if (j >= PEAK_WINDOW / 2) {
1617 Error("SearchHighRes", "Too large sigma");
1618 return 0;
1619 }
1620
1621 if (markov == true) {
1622 if (averWindow <= 0) {
1623 Error("SearchHighRes", "Averaging window must be positive");
1624 return 0;
1625 }
1626 }
1627 if(backgroundRemove == true){
1628 if(ssizex_ext < 2 * number_of_iterations + 1 || ssizey_ext < 2 * number_of_iterations + 1){
1629 Error("SearchHighRes", "Too large clipping window");
1630 return 0;
1631 }
1632 }
1633 i = (Int_t)(4 * sigma + 0.5);
1634 i = 4 * i;
1635 Double_t **working_space = new Double_t *[ssizex + i];
1636 for (j = 0; j < ssizex + i; j++) {
1637 Double_t *wsk = working_space[j] = new Double_t[16 * (ssizey + i)];
1638 for (k=0;k<16 * (ssizey + i);k++) wsk[k] = 0;
1639 }
1640 for(j = 0; j < ssizey_ext; j++){
1641 for(i = 0; i < ssizex_ext; i++){
1642 if(i < shift){
1643 if(j < shift)
1644 working_space[i][j + ssizey_ext] = source[0][0];
1645
1646 else if(j >= ssizey + shift)
1647 working_space[i][j + ssizey_ext] = source[0][ssizey - 1];
1648
1649 else
1650 working_space[i][j + ssizey_ext] = source[0][j - shift];
1651 }
1652
1653 else if(i >= ssizex + shift){
1654 if(j < shift)
1655 working_space[i][j + ssizey_ext] = source[ssizex - 1][0];
1656
1657 else if(j >= ssizey + shift)
1658 working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1];
1659
1660 else
1661 working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift];
1662 }
1663
1664 else{
1665 if(j < shift)
1666 working_space[i][j + ssizey_ext] = source[i - shift][0];
1667
1668 else if(j >= ssizey + shift)
1669 working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1];
1670
1671 else
1672 working_space[i][j + ssizey_ext] = source[i - shift][j - shift];
1673 }
1674 }
1675 }
1676 if(backgroundRemove == true){
1677 for(i = 1; i <= number_of_iterations; i++){
1678 for(y = i; y < ssizey_ext - i; y++){
1679 for(x = i; x < ssizex_ext - i; x++){
1680 a = working_space[x][y + ssizey_ext];
1681 p1 = working_space[x - i][y + ssizey_ext - i];
1682 p2 = working_space[x - i][y + ssizey_ext + i];
1683 p3 = working_space[x + i][y + ssizey_ext - i];
1684 p4 = working_space[x + i][y + ssizey_ext + i];
1685 s1 = working_space[x][y + ssizey_ext - i];
1686 s2 = working_space[x - i][y + ssizey_ext];
1687 s3 = working_space[x + i][y + ssizey_ext];
1688 s4 = working_space[x][y + ssizey_ext + i];
1689 b = (p1 + p2) / 2.0;
1690 if(b > s2)
1691 s2 = b;
1692 b = (p1 + p3) / 2.0;
1693 if(b > s1)
1694 s1 = b;
1695 b = (p2 + p4) / 2.0;
1696 if(b > s4)
1697 s4 = b;
1698 b = (p3 + p4) / 2.0;
1699 if(b > s3)
1700 s3 = b;
1701 s1 = s1 - (p1 + p3) / 2.0;
1702 s2 = s2 - (p1 + p2) / 2.0;
1703 s3 = s3 - (p3 + p4) / 2.0;
1704 s4 = s4 - (p2 + p4) / 2.0;
1705 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
1706 if(b < a)
1707 a = b;
1708 working_space[x][y] = a;
1709 }
1710 }
1711 for(y = i;y < ssizey_ext - i; y++){
1712 for(x = i; x < ssizex_ext - i; x++){
1713 working_space[x][y + ssizey_ext] = working_space[x][y];
1714 }
1715 }
1716 }
1717 for(j = 0;j < ssizey_ext; j++){
1718 for(i = 0; i < ssizex_ext; i++){
1719 if(i < shift){
1720 if(j < shift)
1721 working_space[i][j + ssizey_ext] = source[0][0] - working_space[i][j + ssizey_ext];
1722
1723 else if(j >= ssizey + shift)
1724 working_space[i][j + ssizey_ext] = source[0][ssizey - 1] - working_space[i][j + ssizey_ext];
1725
1726 else
1727 working_space[i][j + ssizey_ext] = source[0][j - shift] - working_space[i][j + ssizey_ext];
1728 }
1729
1730 else if(i >= ssizex + shift){
1731 if(j < shift)
1732 working_space[i][j + ssizey_ext] = source[ssizex - 1][0] - working_space[i][j + ssizey_ext];
1733
1734 else if(j >= ssizey + shift)
1735 working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1] - working_space[i][j + ssizey_ext];
1736
1737 else
1738 working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift] - working_space[i][j + ssizey_ext];
1739 }
1740
1741 else{
1742 if(j < shift)
1743 working_space[i][j + ssizey_ext] = source[i - shift][0] - working_space[i][j + ssizey_ext];
1744
1745 else if(j >= ssizey + shift)
1746 working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1] - working_space[i][j + ssizey_ext];
1747
1748 else
1749 working_space[i][j + ssizey_ext] = source[i - shift][j - shift] - working_space[i][j + ssizey_ext];
1750 }
1751 }
1752 }
1753 }
1754 for(j = 0; j < ssizey_ext; j++){
1755 for(i = 0; i < ssizex_ext; i++){
1756 working_space[i][j + 15*ssizey_ext] = working_space[i][j + ssizey_ext];
1757 }
1758 }
1759 if(markov == true){
1760 for(i = 0;i < ssizex_ext; i++){
1761 for(j = 0; j < ssizey_ext; j++)
1762 working_space[i][j + 2 * ssizex_ext] = working_space[i][ssizey_ext + j];
1763 }
1764 xmin = 0;
1765 xmax = ssizex_ext - 1;
1766 ymin = 0;
1767 ymax = ssizey_ext - 1;
1768 for(i = 0, maxch = 0; i < ssizex_ext; i++){
1769 for(j = 0; j < ssizey_ext; j++){
1770 working_space[i][j] = 0;
1771 if(maxch < working_space[i][j + 2 * ssizey_ext])
1772 maxch = working_space[i][j + 2 * ssizey_ext];
1773 plocha += working_space[i][j + 2 * ssizey_ext];
1774 }
1775 }
1776 if(maxch == 0) {
1777 delete [] working_space;
1778 return 0;
1779 }
1780
1781 nom=0;
1782 working_space[xmin][ymin] = 1;
1783 for(i = xmin; i < xmax; i++){
1784 nip = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1785 nim = working_space[i + 1][ymin + 2 * ssizey_ext] / maxch;
1786 sp = 0,sm = 0;
1787 for(l = 1;l <= averWindow; l++){
1788 if((i + l) > xmax)
1789 a = working_space[xmax][ymin + 2 * ssizey_ext] / maxch;
1790
1791 else
1792 a = working_space[i + l][ymin + 2 * ssizey_ext] / maxch;
1793
1794 b = a - nip;
1795 if(a + nip <= 0)
1796 a = 1;
1797
1798 else
1799 a=TMath::Sqrt(a + nip);
1800 b = b / a;
1801 b = TMath::Exp(b);
1802 sp = sp + b;
1803 if(i - l + 1 < xmin)
1804 a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1805
1806 else
1807 a = working_space[i - l + 1][ymin + 2 * ssizey_ext] / maxch;
1808 b = a - nim;
1809 if(a + nim <= 0)
1810 a = 1;
1811
1812 else
1813 a=TMath::Sqrt(a + nim);
1814 b = b / a;
1815 b = TMath::Exp(b);
1816 sm = sm + b;
1817 }
1818 a = sp / sm;
1819 a = working_space[i + 1][ymin] = a * working_space[i][ymin];
1820 nom = nom + a;
1821 }
1822 for(i = ymin; i < ymax; i++){
1823 nip = working_space[xmin][i + 2 * ssizey_ext] / maxch;
1824 nim = working_space[xmin][i + 1 + 2 * ssizey_ext] / maxch;
1825 sp = 0,sm = 0;
1826 for(l = 1; l <= averWindow; l++){
1827 if((i + l) > ymax)
1828 a = working_space[xmin][ymax + 2 * ssizey_ext] / maxch;
1829
1830 else
1831 a = working_space[xmin][i + l + 2 * ssizey_ext] / maxch;
1832 b = a - nip;
1833 if(a + nip <= 0)
1834 a=1;
1835
1836 else
1837 a=TMath::Sqrt(a + nip);
1838 b = b / a;
1839 b = TMath::Exp(b);
1840 sp = sp + b;
1841 if(i - l + 1 < ymin)
1842 a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1843
1844 else
1845 a = working_space[xmin][i - l + 1 + 2 * ssizey_ext] / maxch;
1846 b = a - nim;
1847 if(a + nim <= 0)
1848 a = 1;
1849
1850 else
1851 a=TMath::Sqrt(a + nim);
1852 b = b / a;
1853 b = TMath::Exp(b);
1854 sm = sm + b;
1855 }
1856 a = sp / sm;
1857 a = working_space[xmin][i + 1] = a * working_space[xmin][i];
1858 nom = nom + a;
1859 }
1860 for(i = xmin; i < xmax; i++){
1861 for(j = ymin; j < ymax; j++){
1862 nip = working_space[i][j + 1 + 2 * ssizey_ext] / maxch;
1863 nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1864 spx = 0,smx = 0;
1865 for(l = 1; l <= averWindow; l++){
1866 if(i + l > xmax)
1867 a = working_space[xmax][j + 2 * ssizey_ext] / maxch;
1868
1869 else
1870 a = working_space[i + l][j + 2 * ssizey_ext] / maxch;
1871 b = a - nip;
1872 if(a + nip <= 0)
1873 a = 1;
1874
1875 else
1876 a=TMath::Sqrt(a + nip);
1877 b = b / a;
1878 b = TMath::Exp(b);
1879 spx = spx + b;
1880 if(i - l + 1 < xmin)
1881 a = working_space[xmin][j + 2 * ssizey_ext] / maxch;
1882
1883 else
1884 a = working_space[i - l + 1][j + 2 * ssizey_ext] / maxch;
1885 b = a - nim;
1886 if(a + nim <= 0)
1887 a=1;
1888
1889 else
1890 a=TMath::Sqrt(a + nim);
1891 b = b / a;
1892 b = TMath::Exp(b);
1893 smx = smx + b;
1894 }
1895 spy = 0,smy = 0;
1896 nip = working_space[i + 1][j + 2 * ssizey_ext] / maxch;
1897 nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1898 for(l = 1; l <= averWindow; l++){
1899 if(j + l > ymax)
1900 a = working_space[i][ymax + 2 * ssizey_ext] / maxch;
1901
1902 else
1903 a = working_space[i][j + l + 2 * ssizey_ext] / maxch;
1904 b = a - nip;
1905 if(a + nip <= 0)
1906 a = 1;
1907
1908 else
1909 a=TMath::Sqrt(a + nip);
1910 b = b / a;
1911 b = TMath::Exp(b);
1912 spy = spy + b;
1913 if(j - l + 1 < ymin)
1914 a = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1915
1916 else
1917 a = working_space[i][j - l + 1 + 2 * ssizey_ext] / maxch;
1918 b=a-nim;
1919 if(a + nim <= 0)
1920 a = 1;
1921 else
1922 a=TMath::Sqrt(a + nim);
1923 b = b / a;
1924 b = TMath::Exp(b);
1925 smy = smy + b;
1926 }
1927 a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx + smy);
1928 working_space[i + 1][j + 1] = a;
1929 nom = nom + a;
1930 }
1931 }
1932 for(i = xmin; i <= xmax; i++){
1933 for(j = ymin; j <= ymax; j++){
1934 working_space[i][j] = working_space[i][j] / nom;
1935 }
1936 }
1937 for(i = 0; i < ssizex_ext; i++){
1938 for(j = 0; j < ssizey_ext; j++){
1939 working_space[i][j + ssizey_ext] = working_space[i][j] * plocha;
1940 working_space[i][2 * ssizey_ext + j] = working_space[i][ssizey_ext + j];
1941 }
1942 }
1943 }
1944 //deconvolution starts
1945 area = 0;
1946 lhx = -1,lhy = -1;
1947 positx = 0,posity = 0;
1948 maximum = 0;
1949 //generate response matrix
1950 for(i = 0; i < ssizex_ext; i++){
1951 for(j = 0; j < ssizey_ext; j++){
1952 lda = (Double_t)i - 3 * sigma;
1953 ldb = (Double_t)j - 3 * sigma;
1954 lda = (lda * lda + ldb * ldb) / (2 * sigma * sigma);
1955 k=(Int_t)(1000 * TMath::Exp(-lda));
1956 lda = k;
1957 if(lda != 0){
1958 if((i + 1) > lhx)
1959 lhx = i + 1;
1960
1961 if((j + 1) > lhy)
1962 lhy = j + 1;
1963 }
1964 working_space[i][j] = lda;
1965 area = area + lda;
1966 if(lda > maximum){
1967 maximum = lda;
1968 positx = i,posity = j;
1969 }
1970 }
1971 }
1972 //read source matrix
1973 for(i = 0;i < ssizex_ext; i++){
1974 for(j = 0;j < ssizey_ext; j++){
1975 working_space[i][j + 14 * ssizey_ext] = TMath::Abs(working_space[i][j + ssizey_ext]);
1976 }
1977 }
1978 //calculate matrix b=ht*h
1979 i = lhx - 1;
1980 if(i > ssizex_ext)
1981 i = ssizex_ext;
1982
1983 j = lhy - 1;
1984 if(j>ssizey_ext)
1985 j = ssizey_ext;
1986
1987 i1min = -i,i1max = i;
1988 i2min = -j,i2max = j;
1989 for(i2 = i2min; i2 <= i2max; i2++){
1990 for(i1 = i1min; i1 <= i1max; i1++){
1991 ldc = 0;
1992 j2min = -i2;
1993 if(j2min < 0)
1994 j2min = 0;
1995
1996 j2max = lhy - 1 - i2;
1997 if(j2max > lhy - 1)
1998 j2max = lhy - 1;
1999
2000 for(j2 = j2min; j2 <= j2max; j2++){
2001 j1min = -i1;
2002 if(j1min < 0)
2003 j1min = 0;
2004
2005 j1max = lhx - 1 - i1;
2006 if(j1max > lhx - 1)
2007 j1max = lhx - 1;
2008
2009 for(j1 = j1min; j1 <= j1max; j1++){
2010 lda = working_space[j1][j2];
2011 ldb = working_space[i1 + j1][i2 + j2];
2012 ldc = ldc + lda * ldb;
2013 }
2014 }
2015 k = (i1 + ssizex_ext) / ssizex_ext;
2016 working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext] = ldc;
2017 }
2018 }
2019 //calculate at*y and write into p
2020 i = lhx - 1;
2021 if(i > ssizex_ext)
2022 i = ssizex_ext;
2023
2024 j = lhy - 1;
2025 if(j > ssizey_ext)
2026 j = ssizey_ext;
2027
2028 i2min = -j,i2max = ssizey_ext + j - 1;
2029 i1min = -i,i1max = ssizex_ext + i - 1;
2030 for(i2 = i2min; i2 <= i2max; i2++){
2031 for(i1=i1min;i1<=i1max;i1++){
2032 ldc=0;
2033 for(j2 = 0; j2 <= (lhy - 1); j2++){
2034 for(j1 = 0; j1 <= (lhx - 1); j1++){
2035 k2 = i2 + j2,k1 = i1 + j1;
2036 if(k2 >= 0 && k2 < ssizey_ext && k1 >= 0 && k1 < ssizex_ext){
2037 lda = working_space[j1][j2];
2038 ldb = working_space[k1][k2 + 14 * ssizey_ext];
2039 ldc = ldc + lda * ldb;
2040 }
2041 }
2042 }
2043 k = (i1 + ssizex_ext) / ssizex_ext;
2044 working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext] = ldc;
2045 }
2046 }
2047 //move matrix p
2048 for(i2 = 0; i2 < ssizey_ext; i2++){
2049 for(i1 = 0; i1 < ssizex_ext; i1++){
2050 k = (i1 + ssizex_ext) / ssizex_ext;
2051 ldb = working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext];
2052 working_space[i1][i2 + 14 * ssizey_ext] = ldb;
2053 }
2054 }
2055 //initialization in x1 matrix
2056 for(i2 = 0; i2 < ssizey_ext; i2++){
2057 for(i1 = 0; i1 < ssizex_ext; i1++){
2058 working_space[i1][i2 + ssizey_ext] = 1;
2059 working_space[i1][i2 + 2 * ssizey_ext] = 0;
2060 }
2061 }
2062 //START OF ITERATIONS
2063 for(lindex = 0; lindex < deconIterations; lindex++){
2064 for(i2 = 0; i2 < ssizey_ext; i2++){
2065 for(i1 = 0; i1 < ssizex_ext; i1++){
2066 lda = working_space[i1][i2 + ssizey_ext];
2067 ldc = working_space[i1][i2 + 14 * ssizey_ext];
2068 if(lda > 0.000001 && ldc > 0.000001){
2069 ldb=0;
2070 j2min=i2;
2071 if(j2min > lhy - 1)
2072 j2min = lhy - 1;
2073
2074 j2min = -j2min;
2075 j2max = ssizey_ext - i2 - 1;
2076 if(j2max > lhy - 1)
2077 j2max = lhy - 1;
2078
2079 j1min = i1;
2080 if(j1min > lhx - 1)
2081 j1min = lhx - 1;
2082
2083 j1min = -j1min;
2084 j1max = ssizex_ext - i1 - 1;
2085 if(j1max > lhx - 1)
2086 j1max = lhx - 1;
2087
2088 for(j2 = j2min; j2 <= j2max; j2++){
2089 for(j1 = j1min; j1 <= j1max; j1++){
2090 k = (j1 + ssizex_ext) / ssizex_ext;
2091 ldc = working_space[(j1 + ssizex_ext) % ssizex_ext][j2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext];
2092 lda = working_space[i1 + j1][i2 + j2 + ssizey_ext];
2093 ldb = ldb + lda * ldc;
2094 }
2095 }
2096 lda = working_space[i1][i2 + ssizey_ext];
2097 ldc = working_space[i1][i2 + 14 * ssizey_ext];
2098 if(ldc * lda != 0 && ldb != 0){
2099 lda =lda * ldc / ldb;
2100 }
2101
2102 else
2103 lda=0;
2104 working_space[i1][i2 + 2 * ssizey_ext] = lda;
2105 }
2106 }
2107 }
2108 for(i2 = 0; i2 < ssizey_ext; i2++){
2109 for(i1 = 0; i1 < ssizex_ext; i1++)
2110 working_space[i1][i2 + ssizey_ext] = working_space[i1][i2 + 2 * ssizey_ext];
2111 }
2112 }
2113 //looking for maximum
2114 maximum=0;
2115 for(i = 0; i < ssizex_ext; i++){
2116 for(j = 0; j < ssizey_ext; j++){
2117 working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext] = area * working_space[i][j + ssizey_ext];
2118 if(maximum < working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext])
2119 maximum = working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext];
2120
2121 }
2122 }
2123 //searching for peaks in deconvolved spectrum
2124 for(i = 1; i < ssizex_ext - 1; i++){
2125 for(j = 1; j < ssizey_ext - 1; j++){
2126 if(working_space[i][j] > working_space[i - 1][j] && working_space[i][j] > working_space[i - 1][j - 1] && working_space[i][j] > working_space[i][j - 1] && working_space[i][j] > working_space[i + 1][j - 1] && working_space[i][j] > working_space[i + 1][j] && working_space[i][j] > working_space[i + 1][j + 1] && working_space[i][j] > working_space[i][j + 1] && working_space[i][j] > working_space[i - 1][j + 1]){
2127 if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift){
2128 if(working_space[i][j] > threshold * maximum / 100.0){
2129 if(peak_index < fMaxPeaks){
2130 for(k = i - 1,a = 0,b = 0; k <= i + 1; k++){
2131 a += (Double_t)(k - shift) * working_space[k][j];
2132 b += working_space[k][j];
2133 }
2134 ax=a/b;
2135 if(ax < 0)
2136 ax = 0;
2137
2138 if(ax >= ssizex)
2139 ax = ssizex - 1;
2140
2141 for(k = j - 1,a = 0,b = 0; k <= j + 1; k++){
2142 a += (Double_t)(k - shift) * working_space[i][k];
2143 b += working_space[i][k];
2144 }
2145 ay=a/b;
2146 if(ay < 0)
2147 ay = 0;
2148
2149 if(ay >= ssizey)
2150 ay = ssizey - 1;
2151
2152 if(peak_index == 0){
2153 fPositionX[0] = ax;
2154 fPositionY[0] = ay;
2155 peak_index = 1;
2156 }
2157
2158 else{
2159 for(k = 0, priz = 0; k < peak_index && priz == 0; k++){
2160 if(working_space[shift+(Int_t)(ax+0.5)][15 * ssizey_ext + shift + (Int_t)(ay+0.5)] > working_space[shift+(Int_t)(fPositionX[k]+0.5)][15 * ssizey_ext + shift + (Int_t)(fPositionY[k]+0.5)])
2161 priz = 1;
2162 }
2163 if(priz == 0){
2164 if(k < fMaxPeaks){
2165 fPositionX[k] = ax;
2166 fPositionY[k] = ay;
2167 }
2168 }
2169
2170 else{
2171 for(l = peak_index; l >= k; l--){
2172 if(l < fMaxPeaks){
2173 fPositionX[l] = fPositionX[l - 1];
2174 fPositionY[l] = fPositionY[l - 1];
2175 }
2176 }
2177 fPositionX[k - 1] = ax;
2178 fPositionY[k - 1] = ay;
2179 }
2180 if(peak_index < fMaxPeaks)
2181 peak_index += 1;
2182 }
2183 }
2184 }
2185 }
2186 }
2187 }
2188 }
2189 //writing back deconvolved spectrum
2190 for(i = 0; i < ssizex; i++){
2191 for(j = 0; j < ssizey; j++){
2192 dest[i][j] = working_space[i + shift][j + shift];
2193 }
2194 }
2195 i = (Int_t)(4 * sigma + 0.5);
2196 i = 4 * i;
2197 for (j = 0; j < ssizex + i; j++)
2198 delete[]working_space[j];
2199 delete[]working_space;
2200 fNPeaks = peak_index;
2201 return fNPeaks;
2202}
2203
2204////////////////////////////////////////////////////////////////////////////////
2205/// static function (called by TH1), interface to TSpectrum2::Search
2206
2207Int_t TSpectrum2::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
2208{
2209 TSpectrum2 s;
2210 return s.Search(hist,sigma,option,threshold);
2211}
2212
2213////////////////////////////////////////////////////////////////////////////////
2214/// static function (called by TH1), interface to TSpectrum2::Background
2215
2216TH1 *TSpectrum2::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
2217{
2218 TSpectrum2 s;
2219 return s.Background(hist,niter,option);
2220}
b
#define b(i)
Definition: RSha256.hxx:100
s1
#define s1(x)
Definition: RSha256.hxx:91
h
#define h(i)
Definition: RSha256.hxx:106
Int_t
int Int_t
Definition: RtypesCore.h:41
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:88
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
Double_t
double Double_t
Definition: RtypesCore.h:55
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:87
Option_t
const char Option_t
Definition: RtypesCore.h:62
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:365
kRed
@ kRed
Definition: Rtypes.h:64
xmin
float xmin
Definition: THbookFile.cxx:93
ymin
float ymin
Definition: THbookFile.cxx:93
xmax
float xmax
Definition: THbookFile.cxx:93
ymax
float ymax
Definition: THbookFile.cxx:93
PEAK_WINDOW
#define PEAK_WINDOW
Definition: TSpectrum2.cxx:47
TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
TAxis::GetBinCenter
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
TAxis::GetNbins
Int_t GetNbins() const
Definition: TAxis.h:121
TH1
The TH1 histogram class.
Definition: TH1.h:56
TH1::GetDimension
virtual Int_t GetDimension() const
Definition: TH1.h:278
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
TH1::GetYaxis
TAxis * GetYaxis()
Definition: TH1.h:317
TH1::GetListOfFunctions
TList * GetListOfFunctions() const
Definition: TH1.h:239
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4882
TList::Add
virtual void Add(TObject *obj)
Definition: TList.h:87
TList::Remove
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:819
TList::FindObject
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition: TList.cxx:575
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
TPolyMarker
A PolyMarker is defined by an array on N points in a 2-D space.
Definition: TPolyMarker.h:31
TSpectrum2
Advanced 2-dimensional spectra processing.
Definition: TSpectrum2.h:18
TSpectrum2::fResolution
Double_t fResolution
NOT USED resolution of the neighboring peaks
Definition: TSpectrum2.h:25
TSpectrum2::fHistogram
TH1 * fHistogram
resulting histogram
Definition: TSpectrum2.h:26
TSpectrum2::fgIterations
static Int_t fgIterations
Maximum number of decon iterations (default=3)
Definition: TSpectrum2.h:28
TSpectrum2::SetDeconIterations
static void SetDeconIterations(Int_t n=3)
static function: Set max number of decon iterations in deconvolution operation see TSpectrum2::Search...
Definition: TSpectrum2.cxx:113
TSpectrum2::StaticBackground
static TH1 * StaticBackground(const TH1 *hist, Int_t niter=20, Option_t *option="")
static function (called by TH1), interface to TSpectrum2::Background
Definition: TSpectrum2.cxx:2216
TSpectrum2::SmoothMarkov
const char * SmoothMarkov(Double_t **source, Int_t ssizex, Int_t ssizey, Int_t averWindow)
This function calculates smoothed spectrum from source spectrum based on Markov chain method.
Definition: TSpectrum2.cxx:736
TSpectrum2::SetResolution
void SetResolution(Double_t resolution=1)
NOT USED resolution: determines resolution of the neighboring peaks default value is 1 correspond to ...
Definition: TSpectrum2.cxx:287
TSpectrum2::~TSpectrum2
virtual ~TSpectrum2()
Destructor.
Definition: TSpectrum2.cxx:92
TSpectrum2::fgAverageWindow
static Int_t fgAverageWindow
Average window of searched peaks.
Definition: TSpectrum2.h:27
TSpectrum2::fMaxPeaks
Int_t fMaxPeaks
Maximum number of peaks to be found.
Definition: TSpectrum2.h:20
TSpectrum2::TSpectrum2
TSpectrum2()
Constructor.
Definition: TSpectrum2.cxx:57
TSpectrum2::SearchHighRes
Int_t SearchHighRes(Double_t **source, Double_t **dest, Int_t ssizex, Int_t ssizey, 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: TSpectrum2.cxx:1589
TSpectrum2::fNPeaks
Int_t fNPeaks
number of peaks found
Definition: TSpectrum2.h:21
TSpectrum2::Background
virtual TH1 * Background(const TH1 *hist, Int_t niter=20, Option_t *option="")
This function calculates the background spectrum in the input histogram h.
Definition: TSpectrum2.cxx:155
TSpectrum2::Print
virtual void Print(Option_t *option="") const
Print the array of positions.
Definition: TSpectrum2.cxx:166
TSpectrum2::SetAverageWindow
static void SetAverageWindow(Int_t w=3)
static function: Set average window of searched peaks see TSpectrum2::SearchHighRes
Definition: TSpectrum2.cxx:104
TSpectrum2::Search
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05)
This function searches for peaks in source spectrum in hin The number of found peaks and their positi...
Definition: TSpectrum2.cxx:206
TSpectrum2::fPositionY
Double_t * fPositionY
[fNPeaks] Y position of peaks
Definition: TSpectrum2.h:24
TSpectrum2::StaticSearch
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
static function (called by TH1), interface to TSpectrum2::Search
Definition: TSpectrum2.cxx:2207
TSpectrum2::fPosition
Double_t * fPosition
[fNPeaks] array of current peak positions
Definition: TSpectrum2.h:22
TSpectrum2::Deconvolution
const char * Deconvolution(Double_t **source, Double_t **resp, Int_t ssizex, Int_t ssizey, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
This function calculates deconvolution from source spectrum according to response spectrum The result...
Definition: TSpectrum2.cxx:1134
TSpectrum2::fPositionX
Double_t * fPositionX
[fNPeaks] X position of peaks
Definition: TSpectrum2.h:23
TSpectrum2::kBackSuccessiveFiltering
@ kBackSuccessiveFiltering
Definition: TSpectrum2.h:34
TSpectrum2::kBackOneStepFiltering
@ kBackOneStepFiltering
Definition: TSpectrum2.h:35
TSpectrum2::kBackDecreasingWindow
@ kBackDecreasingWindow
Definition: TSpectrum2.h:33
TSpectrum2::kBackIncreasingWindow
@ kBackIncreasingWindow
Definition: TSpectrum2.h:32
TString
Basic string class.
Definition: TString.h:131
TString::ToLower
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
TString::ReplaceAll
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
y
Double_t y[n]
Definition: legend1.C:17
x
Double_t x[n]
Definition: legend1.C:17
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::VectorUtil::boost
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only re...
Definition: VectorUtil.h:327
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:715
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:723
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
l
auto * l
Definition: textangle.C:4
a
auto * a
Definition: textangle.C:12
dest
#define dest(otri, vertexptr)
Definition: triangle.c:1040