Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
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 \legacy{TSpectrum2, 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 - One-dimensional background estimation functions
14 - Two-dimensional background estimation functions
15 - One-dimensional smoothing functions
16 - Two-dimensional smoothing functions
17 - One-dimensional deconvolution functions
18 - Two-dimensional deconvolution functions
19 - One-dimensional peak search functions
20 - Two-dimensional peak search functions
21
22 The algorithms in this class have been published in the following references:
23
24 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.
25 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.
26 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.
27
28 These NIM papers are also available as doc or ps files from:
29
30 - [SpectrumDec.ps.gz](ftp://root.cern/root/SpectrumDec.ps.gz)
31 - [SpectrumSrc.ps.gz](ftp://root.cern/root/SpectrumSrc.ps.gz)
32 - [SpectrumBck.ps.gz](ftp://root.cern/root/SpectrumBck.ps.gz)
33
34 See also the
35 [online documentation](https://root.cern/guides/tspectrum-manual) and
36 [tutorials](https://root.cern/doc/master/group__tutorial__spectrum.html).
37
38 All the figures in this page were prepared using the DaqProVis
39 system, Data Acquisition, Processing and Visualization system,
40 developed at the Institute of Physics, Slovak Academy of Sciences, Bratislava,
41 Slovakia.
42*/
43
44#include "TSpectrum2.h"
45#include "TPolyMarker.h"
46#include "TList.h"
47#include "TH2.h"
48#include "TMath.h"
49#include "TVirtualPad.h"
50#define PEAK_WINDOW 1024
51
52Int_t TSpectrum2::fgIterations = 3;
53Int_t TSpectrum2::fgAverageWindow = 3;
54
55
56////////////////////////////////////////////////////////////////////////////////
57/// Constructor.
58
59TSpectrum2::TSpectrum2() :TNamed("Spectrum", "Miroslav Morhac peak finder")
60{
61 Int_t n = 100;
62 fMaxPeaks = n;
63 fPosition = new Double_t[n];
64 fPositionX = new Double_t[n];
65 fPositionY = new Double_t[n];
66 fResolution = 1;
67 fHistogram = nullptr;
68 fNPeaks = 0;
69}
70
71////////////////////////////////////////////////////////////////////////////////
72/// - maxpositions: maximum number of peaks
73/// - resolution: *NOT USED* determines resolution of the neighbouring peaks
74/// default value is 1 correspond to 3 sigma distance
75/// between peaks. Higher values allow higher resolution
76/// (smaller distance between peaks.
77/// May be set later through SetResolution.
78
79TSpectrum2::TSpectrum2(Int_t maxpositions, Double_t resolution) :TNamed("Spectrum", "Miroslav Morhac peak finder")
80{
81 Int_t n = maxpositions;
82 fMaxPeaks = n;
83 fPosition = new Double_t[n];
84 fPositionX = new Double_t[n];
85 fPositionY = new Double_t[n];
86 fHistogram = nullptr;
87 fNPeaks = 0;
88 SetResolution(resolution);
89}
90
91////////////////////////////////////////////////////////////////////////////////
92/// Destructor.
93
94TSpectrum2::~TSpectrum2()
95{
96 delete [] fPosition;
97 delete [] fPositionX;
98 delete [] fPositionY;
99 delete fHistogram;
100}
101
102////////////////////////////////////////////////////////////////////////////////
103/// static function: Set average window of searched peaks
104/// see TSpectrum2::SearchHighRes
105
110
111////////////////////////////////////////////////////////////////////////////////
112/// static function: Set max number of decon iterations in deconvolution operation
113/// see TSpectrum2::SearchHighRes
114
119
120////////////////////////////////////////////////////////////////////////////////
121/// This function calculates the background spectrum in the input histogram h.
122/// The background is returned as a histogram.
123///
124/// Function parameters:
125/// - h: input 2-d histogram
126/// - nIterX, nIterY, (default value = 20), iterations for X and Y
127/// Increasing number of iterations make the result smoother and lower.
128/// - option: may contain one of the following options
129/// - direction of change of clipping window
130/// - possible values=kBackIncreasingWindow
131/// kBackDecreasingWindow
132/// - filterType-determines the algorithm of the filtering
133/// - possible values=kBackSuccessiveFiltering
134/// kBackOneStepFiltering
135/// - "same" : if this option is specified, the resulting background
136/// histogram is superimposed on the picture in the current pad.
137///
138/// NOTE that the background is only evaluated in the current range of h.
139/// ie, if h has a bin range (set via h->GetXaxis()->SetRange(binmin,binmax),
140/// the returned histogram will be created with the same number of bins
141/// as the input histogram h, but only bins from binmin to binmax will be filled
142/// with the estimated background.
143
144TH1 *TSpectrum2::Background(const TH1 *h, Int_t nIterX, Int_t nIterY, Option_t *option)
145{
146 if (h == nullptr)
147 return nullptr;
148 Int_t dimension = h->GetDimension();
149 if (dimension != 2) {
150 Error("Background", "Only implemented for 2-d histograms");
151 return nullptr;
152 }
153 TString opt = option;
154 opt.ToLower();
155
156 // set options
157 Int_t direction = kBackDecreasingWindow;
158 if (opt.Contains("backincreasingwindow"))
159 direction = kBackIncreasingWindow;
160 Int_t filterType = kBackSuccessiveFiltering;
161 if (opt.Contains("backonestepfiltering"))
162 filterType = kBackOneStepFiltering;
163 Int_t firstX = h->GetXaxis()->GetFirst();
164 Int_t lastX = h->GetXaxis()->GetLast();
165 Int_t sizeX = lastX - firstX + 1;
166 Int_t firstY = h->GetYaxis()->GetFirst();
167 Int_t lastY = h->GetYaxis()->GetLast();
168 Int_t sizeY = lastY - firstY + 1;
169 Int_t i, j;
170 Double_t **source = new Double_t *[sizeX];
171 for (i = 0; i < sizeX; i++) {
172 source[i] = new Double_t[sizeY];
173 for (j = 0; j < sizeY; j++)
174 source[i][j] = h->GetBinContent(i + firstX, j + firstY);
175 }
176
177 // find background (source is input and in output contains the background
178 Background(source, sizeX, sizeY, nIterX, nIterY, direction, filterType);
179
180 // create output histogram containing background
181 // only bins in the range of the input histogram are filled
182 Int_t nch = strlen(h->GetName());
183 char *hbname = new char[nch + 20];
184 snprintf(hbname, nch + 20, "%s_background", h->GetName());
185 TH2 *hb = (TH2 *)h->Clone(hbname);
186 hb->Reset();
187 hb->GetListOfFunctions()->Delete();
188 for (i = 0; i < sizeX; i++)
189 for (j = 0; j < sizeY; j++)
190 hb->SetBinContent(i + firstX, j + firstY, source[i][j]);
191 hb->SetEntries(sizeX * sizeY);
192
193 // if option "same is specified, draw the result in the pad
194 if (opt.Contains("same")) {
195 if (gPad)
196 delete gPad->GetPrimitive(hbname);
197 hb->Draw("same");
198 }
199 for (i = 0; i < sizeX; i++) {
200 delete[] source[i];
201 }
202 delete[] source;
203 delete[] hbname;
204 return hb;
205}
206
207////////////////////////////////////////////////////////////////////////////////
208/// Print the array of positions.
209
210void TSpectrum2::Print(Option_t *) const
211{
212 printf("\nNumber of positions = %d\n",fNPeaks);
213 for (Int_t i=0;i<fNPeaks;i++) {
214 printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
215 }
216}
217
218////////////////////////////////////////////////////////////////////////////////
219/// This function searches for peaks in source spectrum in hin
220/// The number of found peaks and their positions are written into
221/// the members fNpeaks and fPositionX.
222/// The search is performed in the current histogram range.
223///
224/// Function parameters:
225/// - hin: pointer to the histogram of source spectrum
226/// - sigma: sigma of searched peaks, for details we refer to manual
227/// - threshold: (default=0.05) peaks with amplitude less than
228/// threshold*highest_peak are discarded. 0<threshold<1
229///
230/// By default, the background is removed before deconvolution.
231/// Specify the option "nobackground" to not remove the background.
232///
233/// By default the "Markov" chain algorithm is used.
234/// Specify the option "noMarkov" to disable this algorithm
235/// Note that by default the source spectrum is replaced by a new spectrum
236///
237/// By default a polymarker object is created and added to the list of
238/// functions of the histogram. The histogram is drawn with the specified
239/// option and the polymarker object drawn on top of the histogram.
240/// The polymarker coordinates correspond to the npeaks peaks found in
241/// the histogram.
242/// A pointer to the polymarker object can be retrieved later via:
243/// ~~~ {.cpp}
244/// TList *functions = hin->GetListOfFunctions();
245/// TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker")
246/// ~~~
247/// Specify the option "goff" to disable the storage and drawing of the
248/// polymarker.
249
250Int_t TSpectrum2::Search(const TH1 * hin, Double_t sigma,
251 Option_t * option, Double_t threshold)
252{
253 if (hin == nullptr)
254 return 0;
255 Int_t dimension = hin->GetDimension();
256 if (dimension != 2) {
257 Error("Search", "Must be a 2-d histogram");
258 return 0;
259 }
260
261 TString opt = option;
262 opt.ToLower();
263 Bool_t background = kTRUE;
264 if (opt.Contains("nobackground")) {
265 background = kFALSE;
266 opt.ReplaceAll("nobackground","");
267 }
268 Bool_t markov = kTRUE;
269 if (opt.Contains("nomarkov")) {
270 markov = kFALSE;
271 opt.ReplaceAll("nomarkov","");
272 }
273
274 Int_t sizex = hin->GetXaxis()->GetNbins();
275 Int_t sizey = hin->GetYaxis()->GetNbins();
276 Int_t i, j, binx,biny, npeaks;
277 Double_t ** source = new Double_t*[sizex];
278 Double_t ** dest = new Double_t*[sizex];
279 for (i = 0; i < sizex; i++) {
280 source[i] = new Double_t[sizey];
281 dest[i] = new Double_t[sizey];
282 for (j = 0; j < sizey; j++) {
283 source[i][j] = hin->GetBinContent(i + 1, j + 1);
284 }
285 }
286 //npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, kTRUE, 3, kTRUE, 10);
287 //the smoothing option is used for 1-d but not for 2-d histograms
288 npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, background, fgIterations, markov, fgAverageWindow);
289
290 //The logic in the loop should be improved to use the fact
291 //that fPositionX,Y give a precise position inside a bin.
292 //The current algorithm takes the center of the bin only.
293 for (i = 0; i < npeaks; i++) {
294 binx = 1 + Int_t(fPositionX[i] + 0.5);
295 biny = 1 + Int_t(fPositionY[i] + 0.5);
296 fPositionX[i] = hin->GetXaxis()->GetBinCenter(binx);
297 fPositionY[i] = hin->GetYaxis()->GetBinCenter(biny);
298 }
299 for (i = 0; i < sizex; i++) {
300 delete [] source[i];
301 delete [] dest[i];
302 }
303 delete [] source;
304 delete [] dest;
305
306 if (opt.Contains("goff"))
307 return npeaks;
308 if (!npeaks) return 0;
309 TPolyMarker * pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
310 if (pm) {
311 hin->GetListOfFunctions()->Remove(pm);
312 delete pm;
313 }
314 pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
315 hin->GetListOfFunctions()->Add(pm);
316 pm->SetMarkerStyle(23);
317 pm->SetMarkerColor(kRed);
318 pm->SetMarkerSize(1.3);
319 ((TH1*)hin)->Draw(option);
320 return npeaks;
321}
322
323////////////////////////////////////////////////////////////////////////////////
324/// *NOT USED*
325/// resolution: determines resolution of the neighboring peaks
326/// default value is 1 correspond to 3 sigma distance
327/// between peaks. Higher values allow higher resolution
328/// (smaller distance between peaks.
329/// May be set later through SetResolution.
330
331void TSpectrum2::SetResolution(Double_t resolution)
332{
333 if (resolution > 1)
334 fResolution = resolution;
335 else
336 fResolution = 1;
337}
338
339////////////////////////////////////////////////////////////////////////////////
340/// This function calculates background spectrum from source spectrum.
341/// The result is placed to the array pointed by spectrum pointer.
342///
343/// Function parameters:
344/// - spectrum-pointer to the array of source spectrum
345/// - ssizex-x length of spectrum
346/// - ssizey-y length of spectrum
347/// - numberIterationsX-maximal x width of clipping window
348/// - numberIterationsY-maximal y width of clipping window
349/// for details we refer to manual
350/// - direction- direction of change of clipping window
351/// - possible values=kBackIncreasingWindow
352/// kBackDecreasingWindow
353/// - filterType-determines the algorithm of the filtering
354/// - possible values=kBackSuccessiveFiltering
355/// kBackOneStepFiltering
356///
357/// ### Background estimation
358///
359/// Goal: Separation of useful information (peaks) from useless information (background)
360///
361/// - method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping algorithm [1]
362///
363/// - there exist two algorithms for the estimation of new value in the channel \f$ i_1,i_2\f$
364///
365/// Algorithm based on Successive Comparisons:
366///
367/// It is an extension of one-dimensional SNIP algorithm to another dimension. For
368/// details we refer to [2].
369///
370/// Algorithm based on One Step Filtering:
371///
372/// New value in the estimated channel is calculated as:
373/// \f[ a = \nu_{p-1}(i_1,i_2)\f]
374/// \f[
375/// 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} +
376/// \f]
377/// \f[
378/// \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}
379/// \f]
380/// \f[ \nu_{p-1}(i_1,i_2) = min(a,b)\f]
381/// where p = 1, 2, ..., number_of_iterations.
382///
383/// #### References:
384///
385/// [1] C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the
386/// quantitative analysis of PIXE spectra in geoscience applications. NIM, B34 (1988), 396-402.
387///
388/// [2] M. Morhac;, J. Kliman, V.
389/// Matouoek, M. Veselsky, I. Turzo.:
390/// Background elimination methods for multidimensional gamma-ray spectra. NIM,
391/// A401 (1997) 113-132.
392///
393/// ### Example 1 - Background_gamma64.C
394///
395/// Begin_Macro(source)
396/// ../../../tutorials/legacy/spectrum/Background_gamma64.C
397/// End_Macro
398///
399/// ### Example 2- Background_gamma256.C
400///
401/// Begin_Macro(source)
402/// ../../../tutorials/legacy/spectrum/Background_gamma256.C
403/// End_Macro
404///
405/// ### Example 3- Background_synt256.C
406///
407/// Begin_Macro(source)
408/// ../../../tutorials/legacy/spectrum/Background_synt256.C
409/// End_Macro
410
411const char *TSpectrum2::Background(Double_t **spectrum,
412 Int_t ssizex, Int_t ssizey,
413 Int_t numberIterationsX,
414 Int_t numberIterationsY,
415 Int_t direction,
416 Int_t filterType)
417{
418 Int_t i, x, y, sampling, r1, r2;
419 Double_t a, b, p1, p2, p3, p4, s1, s2, s3, s4;
420 if (ssizex <= 0 || ssizey <= 0)
421 return "Wrong parameters";
422 if (numberIterationsX < 1 || numberIterationsY < 1)
423 return "Width of Clipping Window Must Be Positive";
424 if (ssizex < 2 * numberIterationsX + 1
425 || ssizey < 2 * numberIterationsY + 1)
426 return ("Too Large Clipping Window");
427 Double_t **working_space = new Double_t*[ssizex];
428 for (i = 0; i < ssizex; i++)
429 working_space[i] = new Double_t[ssizey];
430 sampling =
431 (Int_t) TMath::Max(numberIterationsX, numberIterationsY);
432 if (direction == kBackIncreasingWindow) {
433 if (filterType == kBackSuccessiveFiltering) {
434 for (i = 1; i <= sampling; i++) {
435 r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
436 (Int_t) TMath::Min(i, numberIterationsY);
437 for (y = r2; y < ssizey - r2; y++) {
438 for (x = r1; x < ssizex - r1; x++) {
439 a = spectrum[x][y];
440 p1 = spectrum[x - r1][y - r2];
441 p2 = spectrum[x - r1][y + r2];
442 p3 = spectrum[x + r1][y - r2];
443 p4 = spectrum[x + r1][y + r2];
444 s1 = spectrum[x][y - r2];
445 s2 = spectrum[x - r1][y];
446 s3 = spectrum[x + r1][y];
447 s4 = spectrum[x][y + r2];
448 b = (p1 + p2) / 2.0;
449 if (b > s2)
450 s2 = b;
451 b = (p1 + p3) / 2.0;
452 if (b > s1)
453 s1 = b;
454 b = (p2 + p4) / 2.0;
455 if (b > s4)
456 s4 = b;
457 b = (p3 + p4) / 2.0;
458 if (b > s3)
459 s3 = b;
460 s1 = s1 - (p1 + p3) / 2.0;
461 s2 = s2 - (p1 + p2) / 2.0;
462 s3 = s3 - (p3 + p4) / 2.0;
463 s4 = s4 - (p2 + p4) / 2.0;
464 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
465 p3 +
466 p4) / 4.0;
467 if (b < a && b > 0)
468 a = b;
469 working_space[x][y] = a;
470 }
471 }
472 for (y = r2; y < ssizey - r2; y++) {
473 for (x = r1; x < ssizex - r1; x++) {
474 spectrum[x][y] = working_space[x][y];
475 }
476 }
477 }
478 }
479
480 else if (filterType == kBackOneStepFiltering) {
481 for (i = 1; i <= sampling; i++) {
482 r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
483 (Int_t) TMath::Min(i, numberIterationsY);
484 for (y = r2; y < ssizey - r2; y++) {
485 for (x = r1; x < ssizex - r1; x++) {
486 a = spectrum[x][y];
487 b = -(spectrum[x - r1][y - r2] +
488 spectrum[x - r1][y + r2] + spectrum[x + r1][y -
489 r2]
490 + spectrum[x + r1][y + r2]) / 4 +
491 (spectrum[x][y - r2] + spectrum[x - r1][y] +
492 spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
493 if (b < a && b > 0)
494 a = b;
495 working_space[x][y] = a;
496 }
497 }
498 for (y = i; y < ssizey - i; y++) {
499 for (x = i; x < ssizex - i; x++) {
500 spectrum[x][y] = working_space[x][y];
501 }
502 }
503 }
504 }
505 }
506
507 else if (direction == kBackDecreasingWindow) {
508 if (filterType == kBackSuccessiveFiltering) {
509 for (i = sampling; i >= 1; i--) {
510 r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
511 (Int_t) TMath::Min(i, numberIterationsY);
512 for (y = r2; y < ssizey - r2; y++) {
513 for (x = r1; x < ssizex - r1; x++) {
514 a = spectrum[x][y];
515 p1 = spectrum[x - r1][y - r2];
516 p2 = spectrum[x - r1][y + r2];
517 p3 = spectrum[x + r1][y - r2];
518 p4 = spectrum[x + r1][y + r2];
519 s1 = spectrum[x][y - r2];
520 s2 = spectrum[x - r1][y];
521 s3 = spectrum[x + r1][y];
522 s4 = spectrum[x][y + r2];
523 b = (p1 + p2) / 2.0;
524 if (b > s2)
525 s2 = b;
526 b = (p1 + p3) / 2.0;
527 if (b > s1)
528 s1 = b;
529 b = (p2 + p4) / 2.0;
530 if (b > s4)
531 s4 = b;
532 b = (p3 + p4) / 2.0;
533 if (b > s3)
534 s3 = b;
535 s1 = s1 - (p1 + p3) / 2.0;
536 s2 = s2 - (p1 + p2) / 2.0;
537 s3 = s3 - (p3 + p4) / 2.0;
538 s4 = s4 - (p2 + p4) / 2.0;
539 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
540 p3 +
541 p4) / 4.0;
542 if (b < a && b > 0)
543 a = b;
544 working_space[x][y] = a;
545 }
546 }
547 for (y = r2; y < ssizey - r2; y++) {
548 for (x = r1; x < ssizex - r1; x++) {
549 spectrum[x][y] = working_space[x][y];
550 }
551 }
552 }
553 }
554
555 else if (filterType == kBackOneStepFiltering) {
556 for (i = sampling; i >= 1; i--) {
557 r1 = (Int_t) TMath::Min(i, numberIterationsX), r2 =
558 (Int_t) TMath::Min(i, numberIterationsY);
559 for (y = r2; y < ssizey - r2; y++) {
560 for (x = r1; x < ssizex - r1; x++) {
561 a = spectrum[x][y];
562 b = -(spectrum[x - r1][y - r2] +
563 spectrum[x - r1][y + r2] + spectrum[x + r1][y -
564 r2]
565 + spectrum[x + r1][y + r2]) / 4 +
566 (spectrum[x][y - r2] + spectrum[x - r1][y] +
567 spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
568 if (b < a && b > 0)
569 a = b;
570 working_space[x][y] = a;
571 }
572 }
573 for (y = i; y < ssizey - i; y++) {
574 for (x = i; x < ssizex - i; x++) {
575 spectrum[x][y] = working_space[x][y];
576 }
577 }
578 }
579 }
580 }
581 for (i = 0; i < ssizex; i++)
582 delete[]working_space[i];
583 delete[]working_space;
584 return nullptr;
585}
586
587////////////////////////////////////////////////////////////////////////////////
588/// This function calculates smoothed spectrum from source spectrum
589/// based on Markov chain method.
590/// The result is placed in the array pointed by source pointer.
591///
592/// Function parameters:
593/// - source-pointer to the array of source spectrum
594/// - ssizex-x length of source
595/// - ssizey-y length of source
596/// - averWindow-width of averaging smoothing window
597///
598/// ### Smoothing
599///
600/// Goal: Suppression of statistical fluctuations the algorithm is based on discrete
601/// Markov chain, which has very simple invariant distribution
602/// \f[
603/// 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
604/// \f]
605/// \f$U_1\f$ being defined from the normalization condition \f$ \sum_{i=1}^{n} U_i = 1\f$
606/// n is the length of the smoothed spectrum and
607/// \f[
608/// 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]
609/// \f]
610/// is the probability of the change of the peak position from channel i to the channel i+1.
611/// \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
612/// of smoothing window. We have extended this algorithm to two dimensions.
613///
614/// #### Reference:
615///
616/// [1] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 (1996), 451.
617///
618/// ### Example 4 - Smooth.C
619///
620/// Begin_Macro(source)
621/// ../../../tutorials/legacy/spectrum/Smooth.C
622/// End_Macro
623
624const char* TSpectrum2::SmoothMarkov(Double_t **source, Int_t ssizex, Int_t ssizey, Int_t averWindow)
625{
626
627 Int_t xmin, xmax, ymin, ymax, i, j, l;
628 Double_t a, b, maxch;
629 Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy, plocha = 0;
630 if(averWindow <= 0)
631 return "Averaging Window must be positive";
632 Double_t **working_space = new Double_t*[ssizex];
633 for(i = 0; i < ssizex; i++)
634 working_space[i] = new Double_t[ssizey];
635 xmin = 0;
636 xmax = ssizex - 1;
637 ymin = 0;
638 ymax = ssizey - 1;
639 for(i = 0, maxch = 0; i < ssizex; i++){
640 for(j = 0; j < ssizey; j++){
641 working_space[i][j] = 0;
642 if(maxch < source[i][j])
643 maxch = source[i][j];
644
645 plocha += source[i][j];
646 }
647 }
648 if(maxch == 0) {
649 for (i = 0; i < ssizex; i++)
650 delete[]working_space[i];
651 delete [] working_space;
652 return nullptr;
653 }
654
655 nom = 0;
656 working_space[xmin][ymin] = 1;
657 for(i = xmin; i < xmax; i++){
658 nip = source[i][ymin] / maxch;
659 nim = source[i + 1][ymin] / maxch;
660 sp = 0,sm = 0;
661 for(l = 1; l <= averWindow; l++){
662 if((i + l) > xmax)
663 a = source[xmax][ymin] / maxch;
664
665 else
666 a = source[i + l][ymin] / maxch;
667 b = a - nip;
668 if(a + nip <= 0)
669 a = 1;
670
671 else
672 a = TMath::Sqrt(a + nip);
673 b = b / a;
674 b = TMath::Exp(b);
675 sp = sp + b;
676 if(i - l + 1 < xmin)
677 a = source[xmin][ymin] / maxch;
678
679 else
680 a = source[i - l + 1][ymin] / maxch;
681 b = a - nim;
682 if(a + nim <= 0)
683 a = 1;
684
685 else
686 a = TMath::Sqrt(a + nim);
687 b = b / a;
688 b = TMath::Exp(b);
689 sm = sm + b;
690 }
691 a = sp / sm;
692 a = working_space[i + 1][ymin] = a * working_space[i][ymin];
693 nom = nom + a;
694 }
695 for(i = ymin; i < ymax; i++){
696 nip = source[xmin][i] / maxch;
697 nim = source[xmin][i + 1] / maxch;
698 sp = 0,sm = 0;
699 for(l = 1; l <= averWindow; l++){
700 if((i + l) > ymax)
701 a = source[xmin][ymax] / maxch;
702
703 else
704 a = source[xmin][i + l] / maxch;
705 b = a - nip;
706 if(a + nip <= 0)
707 a = 1;
708
709 else
710 a = TMath::Sqrt(a + nip);
711 b = b / a;
712 b = TMath::Exp(b);
713 sp = sp + b;
714 if(i - l + 1 < ymin)
715 a = source[xmin][ymin] / maxch;
716
717 else
718 a = source[xmin][i - l + 1] / maxch;
719 b = a - nim;
720 if(a + nim <= 0)
721 a = 1;
722
723 else
724 a = TMath::Sqrt(a + nim);
725 b = b / a;
726 b = TMath::Exp(b);
727 sm = sm + b;
728 }
729 a = sp / sm;
730 a = working_space[xmin][i + 1] = a * working_space[xmin][i];
731 nom = nom + a;
732 }
733 for(i = xmin; i < xmax; i++){
734 for(j = ymin; j < ymax; j++){
735 nip = source[i][j + 1] / maxch;
736 nim = source[i + 1][j + 1] / maxch;
737 spx = 0,smx = 0;
738 for(l = 1; l <= averWindow; l++){
739 if(i + l > xmax)
740 a = source[xmax][j] / maxch;
741
742 else
743 a = source[i + l][j] / maxch;
744 b = a - nip;
745 if(a + nip <= 0)
746 a = 1;
747
748 else
749 a = TMath::Sqrt(a + nip);
750 b = b / a;
751 b = TMath::Exp(b);
752 spx = spx + b;
753 if(i - l + 1 < xmin)
754 a = source[xmin][j] / maxch;
755
756 else
757 a = source[i - l + 1][j] / maxch;
758 b = a - nim;
759 if(a + nim <= 0)
760 a = 1;
761
762 else
763 a = TMath::Sqrt(a + nim);
764 b = b / a;
765 b = TMath::Exp(b);
766 smx = smx + b;
767 }
768 spy = 0,smy = 0;
769 nip = source[i + 1][j] / maxch;
770 nim = source[i + 1][j + 1] / maxch;
771 for (l = 1; l <= averWindow; l++) {
772 if (j + l > ymax) a = source[i][ymax]/maxch;
773 else a = source[i][j + l] / maxch;
774 b = a - nip;
775 if (a + nip <= 0) a = 1;
776 else a = TMath::Sqrt(a + nip);
777 b = b / a;
778 b = TMath::Exp(b);
779 spy = spy + b;
780 if (j - l + 1 < ymin) a = source[i][ymin] / maxch;
781 else a = source[i][j - l + 1] / maxch;
782 b = a - nim;
783 if (a + nim <= 0) a = 1;
784 else a = TMath::Sqrt(a + nim);
785 b = b / a;
786 b = TMath::Exp(b);
787 smy = smy + b;
788 }
789 a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx +smy);
790 working_space[i + 1][j + 1] = a;
791 nom = nom + a;
792 }
793 }
794 for(i = xmin; i <= xmax; i++){
795 for(j = ymin; j <= ymax; j++){
796 working_space[i][j] = working_space[i][j] / nom;
797 }
798 }
799 for(i = 0;i < ssizex; i++){
800 for(j = 0; j < ssizey; j++){
801 source[i][j] = plocha * working_space[i][j];
802 }
803 }
804 for (i = 0; i < ssizex; i++)
805 delete[]working_space[i];
806 delete[]working_space;
807 return nullptr;
808}
809
810////////////////////////////////////////////////////////////////////////////////
811/// This function calculates deconvolution from source spectrum
812/// according to response spectrum
813/// The result is placed in the matrix pointed by source pointer.
814///
815/// Function parameters:
816/// - source-pointer to the matrix of source spectrum
817/// - resp-pointer to the matrix of response spectrum
818/// - ssizex-x length of source and response spectra
819/// - ssizey-y length of source and response spectra
820/// - numberIterations, for details we refer to manual
821/// - numberRepetitions, for details we refer to manual
822/// - boost, boosting factor, for details we refer to manual
823///
824/// ### Deconvolution
825///
826/// Goal: Improvement of the resolution in spectra, decomposition of multiplets
827///
828/// Mathematical formulation of the 2-dimensional convolution system is
829///
830///\f[
831/// 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)
832///\f]
833///\f[
834/// i_1 = 0,1,2, ... ,N_1-1, i_2 = 0,1,2, ... ,N_2-1
835///\f]
836///
837/// where h(i,j) is the impulse response function, x, y are input and output
838/// matrices, respectively, \f$ N_1, N_2\f$ are the lengths of x and h matrices
839///
840/// - let us assume that we know the response and the output matrices (spectra) of the above given system.
841/// - the deconvolution represents solution of the overdetermined system of linear equations, i.e., the
842/// calculation of the matrix x.
843/// - from numerical stability point of view the operation of deconvolution is
844/// extremely critical (ill-posed problem) as well as time consuming operation.
845/// - the Gold deconvolution algorithm proves to work very well even for 2-dimensional
846/// systems. Generalisation of the algorithm for 2-dimensional systems was presented in [1], [2].
847/// - for Gold deconvolution algorithm as well as for boosted deconvolution algorithm we refer also to TSpectrum
848///
849/// #### References:
850///
851/// [1] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:
852/// Efficient one- and two-dimensional Gold deconvolution and its application to
853/// gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
854///
855/// [2] Morhac; M., Matouoek V., Kliman J., Efficient algorithm of multidimensional
856/// deconvolution and its application to nuclear data processing, Digital Signal
857/// Processing 13 (2003) 144.
858///
859/// ### Example 5 - Deconvolution2_1.c
860///
861/// Begin_Macro(source)
862/// ../../../tutorials/legacy/spectrum/Deconvolution2_1.C
863/// End_Macro
864///
865/// ### Example 6 - Deconvolution2_2.C
866///
867/// Begin_Macro(source)
868/// ../../../tutorials/legacy/spectrum/Deconvolution2_2.C
869/// End_Macro
870///
871/// ### Example 7 - Deconvolution2_HR.C
872///
873/// Begin_Macro(source)
874/// ../../../tutorials/legacy/spectrum/Deconvolution2_HR.C
875/// End_Macro
876
877const char *TSpectrum2::Deconvolution(Double_t **source, Double_t **resp,
878 Int_t ssizex, Int_t ssizey,
879 Int_t numberIterations,
880 Int_t numberRepetitions,
881 Double_t boost)
882{
883 Int_t i, j, lhx, lhy, i1, i2, j1, j2, k1, k2, lindex, i1min, i1max,
884 i2min, i2max, j1min, j1max, j2min, j2max, positx = 0, posity = 0, repet;
885 Double_t lda, ldb, ldc, area, maximum = 0;
886 if (ssizex <= 0 || ssizey <= 0)
887 return "Wrong parameters";
888 if (numberIterations <= 0)
889 return "Number of iterations must be positive";
890 if (numberRepetitions <= 0)
891 return "Number of repetitions must be positive";
892 Double_t **working_space = new Double_t *[ssizex];
893 for (i = 0; i < ssizex; i++)
894 working_space[i] = new Double_t[5 * ssizey];
895 area = 0;
896 lhx = -1, lhy = -1;
897 for (i = 0; i < ssizex; i++) {
898 for (j = 0; j < ssizey; j++) {
899 lda = resp[i][j];
900 if (lda != 0) {
901 if ((i + 1) > lhx)
902 lhx = i + 1;
903 if ((j + 1) > lhy)
904 lhy = j + 1;
905 }
906 working_space[i][j] = lda;
907 area = area + lda;
908 if (lda > maximum) {
909 maximum = lda;
910 positx = i, posity = j;
911 }
912 }
913 }
914 if (lhx == -1 || lhy == -1) {
915 for (i = 0; i < ssizex; i++)
916 delete[]working_space[i];
917 delete [] working_space;
918 return ("Zero response data");
919 }
920
921//calculate ht*y and write into p
922 for (i2 = 0; i2 < ssizey; i2++) {
923 for (i1 = 0; i1 < ssizex; i1++) {
924 ldc = 0;
925 for (j2 = 0; j2 <= (lhy - 1); j2++) {
926 for (j1 = 0; j1 <= (lhx - 1); j1++) {
927 k2 = i2 + j2, k1 = i1 + j1;
928 if (k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
929 lda = working_space[j1][j2];
930 ldb = source[k1][k2];
931 ldc = ldc + lda * ldb;
932 }
933 }
934 }
935 working_space[i1][i2 + ssizey] = ldc;
936 }
937 }
938
939//calculate matrix b=ht*h
940 i1min = -(lhx - 1), i1max = lhx - 1;
941 i2min = -(lhy - 1), i2max = lhy - 1;
942 for (i2 = i2min; i2 <= i2max; i2++) {
943 for (i1 = i1min; i1 <= i1max; i1++) {
944 ldc = 0;
945 j2min = -i2;
946 if (j2min < 0)
947 j2min = 0;
948 j2max = lhy - 1 - i2;
949 if (j2max > lhy - 1)
950 j2max = lhy - 1;
951 for (j2 = j2min; j2 <= j2max; j2++) {
952 j1min = -i1;
953 if (j1min < 0)
954 j1min = 0;
955 j1max = lhx - 1 - i1;
956 if (j1max > lhx - 1)
957 j1max = lhx - 1;
958 for (j1 = j1min; j1 <= j1max; j1++) {
959 lda = working_space[j1][j2];
960 if (i1 + j1 < ssizex && i2 + j2 < ssizey)
961 ldb = working_space[i1 + j1][i2 + j2];
962 else
963 ldb = 0;
964 ldc = ldc + lda * ldb;
965 }
966 }
967 working_space[i1 - i1min][i2 - i2min + 2 * ssizey ] = ldc;
968 }
969 }
970
971//initialization in x1 matrix
972 for (i2 = 0; i2 < ssizey; i2++) {
973 for (i1 = 0; i1 < ssizex; i1++) {
974 working_space[i1][i2 + 3 * ssizey] = 1;
975 working_space[i1][i2 + 4 * ssizey] = 0;
976 }
977 }
978
979 //START OF ITERATIONS
980 for (repet = 0; repet < numberRepetitions; repet++) {
981 if (repet != 0) {
982 for (i = 0; i < ssizex; i++) {
983 for (j = 0; j < ssizey; j++) {
984 working_space[i][j + 3 * ssizey] =
985 TMath::Power(working_space[i][j + 3 * ssizey], boost);
986 }
987 }
988 }
989 for (lindex = 0; lindex < numberIterations; lindex++) {
990 for (i2 = 0; i2 < ssizey; i2++) {
991 for (i1 = 0; i1 < ssizex; i1++) {
992 ldb = 0;
993 j2min = i2;
994 if (j2min > lhy - 1)
995 j2min = lhy - 1;
996 j2min = -j2min;
997 j2max = ssizey - i2 - 1;
998 if (j2max > lhy - 1)
999 j2max = lhy - 1;
1000 j1min = i1;
1001 if (j1min > lhx - 1)
1002 j1min = lhx - 1;
1003 j1min = -j1min;
1004 j1max = ssizex - i1 - 1;
1005 if (j1max > lhx - 1)
1006 j1max = lhx - 1;
1007 for (j2 = j2min; j2 <= j2max; j2++) {
1008 for (j1 = j1min; j1 <= j1max; j1++) {
1009 ldc = working_space[j1 - i1min][j2 - i2min + 2 * ssizey];
1010 lda = working_space[i1 + j1][i2 + j2 + 3 * ssizey];
1011 ldb = ldb + lda * ldc;
1012 }
1013 }
1014 lda = working_space[i1][i2 + 3 * ssizey];
1015 ldc = working_space[i1][i2 + 1 * ssizey];
1016 if (ldc * lda != 0 && ldb != 0) {
1017 lda = lda * ldc / ldb;
1018 }
1019
1020 else
1021 lda = 0;
1022 working_space[i1][i2 + 4 * ssizey] = lda;
1023 }
1024 }
1025 for (i2 = 0; i2 < ssizey; i2++) {
1026 for (i1 = 0; i1 < ssizex; i1++)
1027 working_space[i1][i2 + 3 * ssizey] =
1028 working_space[i1][i2 + 4 * ssizey];
1029 }
1030 }
1031 }
1032 for (i = 0; i < ssizex; i++) {
1033 for (j = 0; j < ssizey; j++)
1034 source[(i + positx) % ssizex][(j + posity) % ssizey] =
1035 area * working_space[i][j + 3 * ssizey];
1036 }
1037 for (i = 0; i < ssizex; i++)
1038 delete[]working_space[i];
1039 delete[]working_space;
1040 return nullptr;
1041}
1042
1043////////////////////////////////////////////////////////////////////////////////
1044/// This function searches for peaks in source spectrum
1045/// It is based on deconvolution method. First the background is
1046/// removed (if desired), then Markov spectrum is calculated
1047/// (if desired), then the response function is generated
1048/// according to given sigma and deconvolution is carried out.
1049///
1050/// Function parameters:
1051/// - source-pointer to the matrix of source spectrum
1052/// - dest-pointer to the matrix of resulting deconvolved spectrum
1053/// - ssizex-x length of source spectrum
1054/// - ssizey-y length of source spectrum
1055/// - sigma-sigma of searched peaks, for details we refer to manual
1056/// - threshold-threshold value in % for selected peaks, peaks with
1057/// amplitude less than threshold*highest_peak/100
1058/// are ignored, see manual
1059/// - backgroundRemove-logical variable, set if the removal of
1060/// background before deconvolution is desired
1061/// - deconIterations-number of iterations in deconvolution operation
1062/// - markov-logical variable, if it is true, first the source spectrum
1063/// is replaced by new spectrum calculated using Markov
1064/// chains method.
1065/// - averWindow-averaging window of searched peaks, for details
1066/// we refer to manual (applies only for Markov method)
1067///
1068/// ### Peaks searching
1069///
1070/// Goal: to identify automatically the peaks in spectrum with the presence of the
1071/// continuous background, one-fold coincidences (ridges) and statistical
1072/// fluctuations - noise.
1073///
1074/// The common problems connected with correct peak identification in two-dimensional coincidence spectra are
1075///
1076/// - non-sensitivity to noise, i.e., only statistically relevant peaks should be identified
1077/// - non-sensitivity of the algorithm to continuous background
1078/// - non-sensitivity to one-fold coincidences (coincidences peak - background in both dimensions) and their crossings
1079/// - ability to identify peaks close to the edges of the spectrum region. Usually peak finders fail to detect them
1080/// - resolution, decomposition of doublets and multiplets. The algorithm should be able to recognise close positioned peaks.
1081/// - ability to identify peaks with different sigma
1082///
1083/// #### References:
1084///
1085/// [1] M.A. Mariscotti: A method for identification of peaks in the presence of
1086/// background and its application to spectrum analysis. NIM 50 (1967), 309-320.
1087///
1088/// [2] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:Identification
1089/// of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
1090/// 108-125.
1091///
1092/// [3] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
1093/// (1996), 451.
1094///
1095/// ### Examples of peak searching method
1096///
1097/// SearchHighRes function provides users with the possibility
1098/// to vary the input parameters and with the access to the output deconvolved data
1099/// in the destination spectrum. Based on the output data one can tune the
1100/// parameters.
1101///
1102/// ### Example 8 - Src.C
1103///
1104/// Begin_Macro(source)
1105/// ../../../tutorials/legacy/spectrum/Src.C
1106/// End_Macro
1107///
1108/// ### Example 9 - Src2.C
1109///
1110/// Begin_Macro(source)
1111/// ../../../tutorials/legacy/spectrum/Src2.C
1112/// End_Macro
1113///
1114/// ### Example 10 - Src3.C
1115///
1116/// Begin_Macro(source)
1117/// ../../../tutorials/legacy/spectrum/Src3.C
1118/// End_Macro
1119///
1120/// ### Example 11 - Src4.C
1121///
1122/// Begin_Macro(source)
1123/// ../../../tutorials/legacy/spectrum/Src4.C
1124/// End_Macro
1125///
1126/// ### Example 12 - Src5.C
1127///
1128/// Begin_Macro(source)
1129/// ../../../tutorials/legacy/spectrum/Src5.C
1130/// End_Macro
1131
1132Int_t TSpectrum2::SearchHighRes(Double_t **source, Double_t **dest, Int_t ssizex, Int_t ssizey,
1133 Double_t sigma, Double_t threshold,
1134 Bool_t backgroundRemove,Int_t deconIterations,
1135 Bool_t markov, Int_t averWindow)
1136
1137{
1138 Int_t number_of_iterations = (Int_t)(4 * sigma + 0.5);
1139 Int_t k, lindex, priz;
1140 Double_t lda, ldb, ldc, area, maximum;
1141 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;
1142 Int_t ymin, ymax, i, j;
1143 Double_t a, b, ax, ay, maxch, plocha = 0;
1144 Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy;
1145 Double_t p1, p2, p3, p4, s1, s2, s3, s4;
1146 Int_t x, y;
1147 Int_t lhx, lhy, i1, i2, j1, j2, k1, k2, i1min, i1max, i2min, i2max, j1min, j1max, j2min, j2max, positx, posity;
1148 if (sigma < 1) {
1149 Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
1150 return 0;
1151 }
1152
1153 if(threshold<=0||threshold>=100){
1154 Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
1155 return 0;
1156 }
1157
1158 j = (Int_t) (5.0 * sigma + 0.5);
1159 if (j >= PEAK_WINDOW / 2) {
1160 Error("SearchHighRes", "Too large sigma");
1161 return 0;
1162 }
1163
1164 if (markov == true) {
1165 if (averWindow <= 0) {
1166 Error("SearchHighRes", "Averaging window must be positive");
1167 return 0;
1168 }
1169 }
1170 if(backgroundRemove == true){
1171 if(ssizex_ext < 2 * number_of_iterations + 1 || ssizey_ext < 2 * number_of_iterations + 1){
1172 Error("SearchHighRes", "Too large clipping window");
1173 return 0;
1174 }
1175 }
1176 i = (Int_t)(4 * sigma + 0.5);
1177 i = 4 * i;
1178 Double_t **working_space = new Double_t *[ssizex + i];
1179 for (j = 0; j < ssizex + i; j++) {
1180 Double_t *wsk = working_space[j] = new Double_t[16 * (ssizey + i)];
1181 for (k=0;k<16 * (ssizey + i);k++) wsk[k] = 0;
1182 }
1183 for(j = 0; j < ssizey_ext; j++){
1184 for(i = 0; i < ssizex_ext; i++){
1185 if(i < shift){
1186 if(j < shift)
1187 working_space[i][j + ssizey_ext] = source[0][0];
1188
1189 else if(j >= ssizey + shift)
1190 working_space[i][j + ssizey_ext] = source[0][ssizey - 1];
1191
1192 else
1193 working_space[i][j + ssizey_ext] = source[0][j - shift];
1194 }
1195
1196 else if(i >= ssizex + shift){
1197 if(j < shift)
1198 working_space[i][j + ssizey_ext] = source[ssizex - 1][0];
1199
1200 else if(j >= ssizey + shift)
1201 working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1];
1202
1203 else
1204 working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift];
1205 }
1206
1207 else{
1208 if(j < shift)
1209 working_space[i][j + ssizey_ext] = source[i - shift][0];
1210
1211 else if(j >= ssizey + shift)
1212 working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1];
1213
1214 else
1215 working_space[i][j + ssizey_ext] = source[i - shift][j - shift];
1216 }
1217 }
1218 }
1219 if(backgroundRemove == true){
1220 for(i = 1; i <= number_of_iterations; i++){
1221 for(y = i; y < ssizey_ext - i; y++){
1222 for(x = i; x < ssizex_ext - i; x++){
1223 a = working_space[x][y + ssizey_ext];
1224 p1 = working_space[x - i][y + ssizey_ext - i];
1225 p2 = working_space[x - i][y + ssizey_ext + i];
1226 p3 = working_space[x + i][y + ssizey_ext - i];
1227 p4 = working_space[x + i][y + ssizey_ext + i];
1228 s1 = working_space[x][y + ssizey_ext - i];
1229 s2 = working_space[x - i][y + ssizey_ext];
1230 s3 = working_space[x + i][y + ssizey_ext];
1231 s4 = working_space[x][y + ssizey_ext + i];
1232 b = (p1 + p2) / 2.0;
1233 if(b > s2)
1234 s2 = b;
1235 b = (p1 + p3) / 2.0;
1236 if(b > s1)
1237 s1 = b;
1238 b = (p2 + p4) / 2.0;
1239 if(b > s4)
1240 s4 = b;
1241 b = (p3 + p4) / 2.0;
1242 if(b > s3)
1243 s3 = b;
1244 s1 = s1 - (p1 + p3) / 2.0;
1245 s2 = s2 - (p1 + p2) / 2.0;
1246 s3 = s3 - (p3 + p4) / 2.0;
1247 s4 = s4 - (p2 + p4) / 2.0;
1248 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
1249 if(b < a)
1250 a = b;
1251 working_space[x][y] = a;
1252 }
1253 }
1254 for(y = i;y < ssizey_ext - i; y++){
1255 for(x = i; x < ssizex_ext - i; x++){
1256 working_space[x][y + ssizey_ext] = working_space[x][y];
1257 }
1258 }
1259 }
1260 for(j = 0;j < ssizey_ext; j++){
1261 for(i = 0; i < ssizex_ext; i++){
1262 if(i < shift){
1263 if(j < shift)
1264 working_space[i][j + ssizey_ext] = source[0][0] - working_space[i][j + ssizey_ext];
1265
1266 else if(j >= ssizey + shift)
1267 working_space[i][j + ssizey_ext] = source[0][ssizey - 1] - working_space[i][j + ssizey_ext];
1268
1269 else
1270 working_space[i][j + ssizey_ext] = source[0][j - shift] - working_space[i][j + ssizey_ext];
1271 }
1272
1273 else if(i >= ssizex + shift){
1274 if(j < shift)
1275 working_space[i][j + ssizey_ext] = source[ssizex - 1][0] - working_space[i][j + ssizey_ext];
1276
1277 else if(j >= ssizey + shift)
1278 working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1] - working_space[i][j + ssizey_ext];
1279
1280 else
1281 working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift] - working_space[i][j + ssizey_ext];
1282 }
1283
1284 else{
1285 if(j < shift)
1286 working_space[i][j + ssizey_ext] = source[i - shift][0] - working_space[i][j + ssizey_ext];
1287
1288 else if(j >= ssizey + shift)
1289 working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1] - working_space[i][j + ssizey_ext];
1290
1291 else
1292 working_space[i][j + ssizey_ext] = source[i - shift][j - shift] - working_space[i][j + ssizey_ext];
1293 }
1294 }
1295 }
1296 }
1297 for(j = 0; j < ssizey_ext; j++){
1298 for(i = 0; i < ssizex_ext; i++){
1299 working_space[i][j + 15*ssizey_ext] = working_space[i][j + ssizey_ext];
1300 }
1301 }
1302 if(markov == true){
1303 for(i = 0;i < ssizex_ext; i++){
1304 for(j = 0; j < ssizey_ext; j++)
1305 working_space[i][j + 2 * ssizex_ext] = working_space[i][ssizey_ext + j];
1306 }
1307 xmin = 0;
1308 xmax = ssizex_ext - 1;
1309 ymin = 0;
1310 ymax = ssizey_ext - 1;
1311 for(i = 0, maxch = 0; i < ssizex_ext; i++){
1312 for(j = 0; j < ssizey_ext; j++){
1313 working_space[i][j] = 0;
1314 if(maxch < working_space[i][j + 2 * ssizey_ext])
1315 maxch = working_space[i][j + 2 * ssizey_ext];
1316 plocha += working_space[i][j + 2 * ssizey_ext];
1317 }
1318 }
1319 if(maxch == 0) {
1320 i = (Int_t)(4 * sigma + 0.5);
1321 i = 4 * i;
1322 for (j = 0; j < ssizex + i; j++)
1323 delete[]working_space[j];
1324 delete [] working_space;
1325 return 0;
1326 }
1327
1328 nom=0;
1329 working_space[xmin][ymin] = 1;
1330 for(i = xmin; i < xmax; i++){
1331 nip = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1332 nim = working_space[i + 1][ymin + 2 * ssizey_ext] / maxch;
1333 sp = 0,sm = 0;
1334 for(l = 1;l <= averWindow; l++){
1335 if((i + l) > xmax)
1336 a = working_space[xmax][ymin + 2 * ssizey_ext] / maxch;
1337
1338 else
1339 a = working_space[i + l][ymin + 2 * ssizey_ext] / maxch;
1340
1341 b = a - nip;
1342 if(a + nip <= 0)
1343 a = 1;
1344
1345 else
1346 a=TMath::Sqrt(a + nip);
1347 b = b / a;
1348 b = TMath::Exp(b);
1349 sp = sp + b;
1350 if(i - l + 1 < xmin)
1351 a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1352
1353 else
1354 a = working_space[i - l + 1][ymin + 2 * ssizey_ext] / maxch;
1355 b = a - nim;
1356 if(a + nim <= 0)
1357 a = 1;
1358
1359 else
1360 a=TMath::Sqrt(a + nim);
1361 b = b / a;
1362 b = TMath::Exp(b);
1363 sm = sm + b;
1364 }
1365 a = sp / sm;
1366 a = working_space[i + 1][ymin] = a * working_space[i][ymin];
1367 nom = nom + a;
1368 }
1369 for(i = ymin; i < ymax; i++){
1370 nip = working_space[xmin][i + 2 * ssizey_ext] / maxch;
1371 nim = working_space[xmin][i + 1 + 2 * ssizey_ext] / maxch;
1372 sp = 0,sm = 0;
1373 for(l = 1; l <= averWindow; l++){
1374 if((i + l) > ymax)
1375 a = working_space[xmin][ymax + 2 * ssizey_ext] / maxch;
1376
1377 else
1378 a = working_space[xmin][i + l + 2 * ssizey_ext] / maxch;
1379 b = a - nip;
1380 if(a + nip <= 0)
1381 a=1;
1382
1383 else
1384 a=TMath::Sqrt(a + nip);
1385 b = b / a;
1386 b = TMath::Exp(b);
1387 sp = sp + b;
1388 if(i - l + 1 < ymin)
1389 a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1390
1391 else
1392 a = working_space[xmin][i - l + 1 + 2 * ssizey_ext] / maxch;
1393 b = a - nim;
1394 if(a + nim <= 0)
1395 a = 1;
1396
1397 else
1398 a=TMath::Sqrt(a + nim);
1399 b = b / a;
1400 b = TMath::Exp(b);
1401 sm = sm + b;
1402 }
1403 a = sp / sm;
1404 a = working_space[xmin][i + 1] = a * working_space[xmin][i];
1405 nom = nom + a;
1406 }
1407 for(i = xmin; i < xmax; i++){
1408 for(j = ymin; j < ymax; j++){
1409 nip = working_space[i][j + 1 + 2 * ssizey_ext] / maxch;
1410 nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1411 spx = 0,smx = 0;
1412 for(l = 1; l <= averWindow; l++){
1413 if(i + l > xmax)
1414 a = working_space[xmax][j + 2 * ssizey_ext] / maxch;
1415
1416 else
1417 a = working_space[i + l][j + 2 * ssizey_ext] / maxch;
1418 b = a - nip;
1419 if(a + nip <= 0)
1420 a = 1;
1421
1422 else
1423 a=TMath::Sqrt(a + nip);
1424 b = b / a;
1425 b = TMath::Exp(b);
1426 spx = spx + b;
1427 if(i - l + 1 < xmin)
1428 a = working_space[xmin][j + 2 * ssizey_ext] / maxch;
1429
1430 else
1431 a = working_space[i - l + 1][j + 2 * ssizey_ext] / maxch;
1432 b = a - nim;
1433 if(a + nim <= 0)
1434 a=1;
1435
1436 else
1437 a=TMath::Sqrt(a + nim);
1438 b = b / a;
1439 b = TMath::Exp(b);
1440 smx = smx + b;
1441 }
1442 spy = 0,smy = 0;
1443 nip = working_space[i + 1][j + 2 * ssizey_ext] / maxch;
1444 nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1445 for(l = 1; l <= averWindow; l++){
1446 if(j + l > ymax)
1447 a = working_space[i][ymax + 2 * ssizey_ext] / maxch;
1448
1449 else
1450 a = working_space[i][j + l + 2 * ssizey_ext] / maxch;
1451 b = a - nip;
1452 if(a + nip <= 0)
1453 a = 1;
1454
1455 else
1456 a=TMath::Sqrt(a + nip);
1457 b = b / a;
1458 b = TMath::Exp(b);
1459 spy = spy + b;
1460 if(j - l + 1 < ymin)
1461 a = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1462
1463 else
1464 a = working_space[i][j - l + 1 + 2 * ssizey_ext] / maxch;
1465 b=a-nim;
1466 if(a + nim <= 0)
1467 a = 1;
1468 else
1469 a=TMath::Sqrt(a + nim);
1470 b = b / a;
1471 b = TMath::Exp(b);
1472 smy = smy + b;
1473 }
1474 a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx + smy);
1475 working_space[i + 1][j + 1] = a;
1476 nom = nom + a;
1477 }
1478 }
1479 for(i = xmin; i <= xmax; i++){
1480 for(j = ymin; j <= ymax; j++){
1481 working_space[i][j] = working_space[i][j] / nom;
1482 }
1483 }
1484 for(i = 0; i < ssizex_ext; i++){
1485 for(j = 0; j < ssizey_ext; j++){
1486 working_space[i][j + ssizey_ext] = working_space[i][j] * plocha;
1487 working_space[i][2 * ssizey_ext + j] = working_space[i][ssizey_ext + j];
1488 }
1489 }
1490 }
1491 //deconvolution starts
1492 area = 0;
1493 lhx = -1,lhy = -1;
1494 positx = 0,posity = 0;
1495 maximum = 0;
1496 //generate response matrix
1497 for(i = 0; i < ssizex_ext; i++){
1498 for(j = 0; j < ssizey_ext; j++){
1499 lda = (Double_t)i - 3 * sigma;
1500 ldb = (Double_t)j - 3 * sigma;
1501 lda = (lda * lda + ldb * ldb) / (2 * sigma * sigma);
1502 k=(Int_t)(1000 * TMath::Exp(-lda));
1503 lda = k;
1504 if(lda != 0){
1505 if((i + 1) > lhx)
1506 lhx = i + 1;
1507
1508 if((j + 1) > lhy)
1509 lhy = j + 1;
1510 }
1511 working_space[i][j] = lda;
1512 area = area + lda;
1513 if(lda > maximum){
1514 maximum = lda;
1515 positx = i,posity = j;
1516 }
1517 }
1518 }
1519 //read source matrix
1520 for(i = 0;i < ssizex_ext; i++){
1521 for(j = 0;j < ssizey_ext; j++){
1522 working_space[i][j + 14 * ssizey_ext] = TMath::Abs(working_space[i][j + ssizey_ext]);
1523 }
1524 }
1525 //calculate matrix b=ht*h
1526 i = lhx - 1;
1527 if(i > ssizex_ext)
1528 i = ssizex_ext;
1529
1530 j = lhy - 1;
1531 if(j>ssizey_ext)
1532 j = ssizey_ext;
1533
1534 i1min = -i,i1max = i;
1535 i2min = -j,i2max = j;
1536 for(i2 = i2min; i2 <= i2max; i2++){
1537 for(i1 = i1min; i1 <= i1max; i1++){
1538 ldc = 0;
1539 j2min = -i2;
1540 if(j2min < 0)
1541 j2min = 0;
1542
1543 j2max = lhy - 1 - i2;
1544 if(j2max > lhy - 1)
1545 j2max = lhy - 1;
1546
1547 for(j2 = j2min; j2 <= j2max; j2++){
1548 j1min = -i1;
1549 if(j1min < 0)
1550 j1min = 0;
1551
1552 j1max = lhx - 1 - i1;
1553 if(j1max > lhx - 1)
1554 j1max = lhx - 1;
1555
1556 for(j1 = j1min; j1 <= j1max; j1++){
1557 lda = working_space[j1][j2];
1558 ldb = working_space[i1 + j1][i2 + j2];
1559 ldc = ldc + lda * ldb;
1560 }
1561 }
1562 k = (i1 + ssizex_ext) / ssizex_ext;
1563 working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext] = ldc;
1564 }
1565 }
1566 //calculate at*y and write into p
1567 i = lhx - 1;
1568 if(i > ssizex_ext)
1569 i = ssizex_ext;
1570
1571 j = lhy - 1;
1572 if(j > ssizey_ext)
1573 j = ssizey_ext;
1574
1575 i2min = -j,i2max = ssizey_ext + j - 1;
1576 i1min = -i,i1max = ssizex_ext + i - 1;
1577 for(i2 = i2min; i2 <= i2max; i2++){
1578 for(i1=i1min;i1<=i1max;i1++){
1579 ldc=0;
1580 for(j2 = 0; j2 <= (lhy - 1); j2++){
1581 for(j1 = 0; j1 <= (lhx - 1); j1++){
1582 k2 = i2 + j2,k1 = i1 + j1;
1583 if(k2 >= 0 && k2 < ssizey_ext && k1 >= 0 && k1 < ssizex_ext){
1584 lda = working_space[j1][j2];
1585 ldb = working_space[k1][k2 + 14 * ssizey_ext];
1586 ldc = ldc + lda * ldb;
1587 }
1588 }
1589 }
1590 k = (i1 + ssizex_ext) / ssizex_ext;
1591 working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext] = ldc;
1592 }
1593 }
1594 //move matrix p
1595 for(i2 = 0; i2 < ssizey_ext; i2++){
1596 for(i1 = 0; i1 < ssizex_ext; i1++){
1597 k = (i1 + ssizex_ext) / ssizex_ext;
1598 ldb = working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext];
1599 working_space[i1][i2 + 14 * ssizey_ext] = ldb;
1600 }
1601 }
1602 //initialization in x1 matrix
1603 for(i2 = 0; i2 < ssizey_ext; i2++){
1604 for(i1 = 0; i1 < ssizex_ext; i1++){
1605 working_space[i1][i2 + ssizey_ext] = 1;
1606 working_space[i1][i2 + 2 * ssizey_ext] = 0;
1607 }
1608 }
1609 //START OF ITERATIONS
1610 for(lindex = 0; lindex < deconIterations; lindex++){
1611 for(i2 = 0; i2 < ssizey_ext; i2++){
1612 for(i1 = 0; i1 < ssizex_ext; i1++){
1613 lda = working_space[i1][i2 + ssizey_ext];
1614 ldc = working_space[i1][i2 + 14 * ssizey_ext];
1615 if(lda > 0.000001 && ldc > 0.000001){
1616 ldb=0;
1617 j2min=i2;
1618 if(j2min > lhy - 1)
1619 j2min = lhy - 1;
1620
1621 j2min = -j2min;
1622 j2max = ssizey_ext - i2 - 1;
1623 if(j2max > lhy - 1)
1624 j2max = lhy - 1;
1625
1626 j1min = i1;
1627 if(j1min > lhx - 1)
1628 j1min = lhx - 1;
1629
1630 j1min = -j1min;
1631 j1max = ssizex_ext - i1 - 1;
1632 if(j1max > lhx - 1)
1633 j1max = lhx - 1;
1634
1635 for(j2 = j2min; j2 <= j2max; j2++){
1636 for(j1 = j1min; j1 <= j1max; j1++){
1637 k = (j1 + ssizex_ext) / ssizex_ext;
1638 ldc = working_space[(j1 + ssizex_ext) % ssizex_ext][j2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext];
1639 lda = working_space[i1 + j1][i2 + j2 + ssizey_ext];
1640 ldb = ldb + lda * ldc;
1641 }
1642 }
1643 lda = working_space[i1][i2 + ssizey_ext];
1644 ldc = working_space[i1][i2 + 14 * ssizey_ext];
1645 if(ldc * lda != 0 && ldb != 0){
1646 lda =lda * ldc / ldb;
1647 }
1648
1649 else
1650 lda=0;
1651 working_space[i1][i2 + 2 * ssizey_ext] = lda;
1652 }
1653 }
1654 }
1655 for(i2 = 0; i2 < ssizey_ext; i2++){
1656 for(i1 = 0; i1 < ssizex_ext; i1++)
1657 working_space[i1][i2 + ssizey_ext] = working_space[i1][i2 + 2 * ssizey_ext];
1658 }
1659 }
1660 //looking for maximum
1661 maximum=0;
1662 for(i = 0; i < ssizex_ext; i++){
1663 for(j = 0; j < ssizey_ext; j++){
1664 working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext] = area * working_space[i][j + ssizey_ext];
1665 if(maximum < working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext])
1666 maximum = working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext];
1667
1668 }
1669 }
1670 //searching for peaks in deconvolved spectrum
1671 for(i = 1; i < ssizex_ext - 1; i++){
1672 for(j = 1; j < ssizey_ext - 1; j++){
1673 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]){
1674 if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift){
1675 if(working_space[i][j] > threshold * maximum / 100.0){
1676 if(peak_index < fMaxPeaks){
1677 for(k = i - 1,a = 0,b = 0; k <= i + 1; k++){
1678 a += (Double_t)(k - shift) * working_space[k][j];
1679 b += working_space[k][j];
1680 }
1681 ax=a/b;
1682 if(ax < 0)
1683 ax = 0;
1684
1685 if(ax >= ssizex)
1686 ax = ssizex - 1;
1687
1688 for(k = j - 1,a = 0,b = 0; k <= j + 1; k++){
1689 a += (Double_t)(k - shift) * working_space[i][k];
1690 b += working_space[i][k];
1691 }
1692 ay=a/b;
1693 if(ay < 0)
1694 ay = 0;
1695
1696 if(ay >= ssizey)
1697 ay = ssizey - 1;
1698
1699 if(peak_index == 0){
1700 fPositionX[0] = ax;
1701 fPositionY[0] = ay;
1702 peak_index = 1;
1703 }
1704
1705 else{
1706 for(k = 0, priz = 0; k < peak_index && priz == 0; k++){
1707 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)])
1708 priz = 1;
1709 }
1710 if(priz == 0){
1711 if(k < fMaxPeaks){
1712 fPositionX[k] = ax;
1713 fPositionY[k] = ay;
1714 }
1715 }
1716
1717 else{
1718 for(l = peak_index; l >= k; l--){
1719 if(l < fMaxPeaks){
1720 fPositionX[l] = fPositionX[l - 1];
1721 fPositionY[l] = fPositionY[l - 1];
1722 }
1723 }
1724 fPositionX[k - 1] = ax;
1725 fPositionY[k - 1] = ay;
1726 }
1727 if(peak_index < fMaxPeaks)
1728 peak_index += 1;
1729 }
1730 }
1731 }
1732 }
1733 }
1734 }
1735 }
1736 //writing back deconvolved spectrum
1737 for(i = 0; i < ssizex; i++){
1738 for(j = 0; j < ssizey; j++){
1739 dest[i][j] = working_space[i + shift][j + shift];
1740 }
1741 }
1742 i = (Int_t)(4 * sigma + 0.5);
1743 i = 4 * i;
1744 for (j = 0; j < ssizex + i; j++)
1745 delete[]working_space[j];
1746 delete[]working_space;
1747 fNPeaks = peak_index;
1748 return fNPeaks;
1749}
1750
1751////////////////////////////////////////////////////////////////////////////////
1752/// static function (called by TH2), interface to TSpectrum2::Search
1753
1759
1760////////////////////////////////////////////////////////////////////////////////
1761/// static function (called by TH2), interface to TSpectrum2::Background
1762
1763TH1 *TSpectrum2::StaticBackground(const TH1 *hist, Int_t nIterX, Int_t nIterY, Option_t *option)
1764{
1765 TSpectrum2 s;
1766 return s.Background(hist, nIterX, nIterY, option);
1767}
b
#define b(i)
Definition RSha256.hxx:100
a
#define a(i)
Definition RSha256.hxx:99
s1
#define s1(x)
Definition RSha256.hxx:91
h
#define h(i)
Definition RSha256.hxx:106
Int_t
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
Double_t
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
kTRUE
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
Option_t
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
kRed
@ kRed
Definition Rtypes.h:67
TRangeDynCast
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Definition TCollection.h:359
option
Option_t Option_t option
Definition TGWin32VirtualXProxy.cxx:44
dest
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t dest
Definition TGWin32VirtualXProxy.cxx:164
xmin
float xmin
Definition THbookFile.cxx:95
ymin
float ymin
Definition THbookFile.cxx:95
xmax
float xmax
Definition THbookFile.cxx:95
ymax
float ymax
Definition THbookFile.cxx:95
hbname
#define hbname
Definition THbookFile.cxx:125
PEAK_WINDOW
#define PEAK_WINDOW
Definition TSpectrum.cxx:37
gPad
#define gPad
Definition TVirtualPad.h:315
snprintf
#define snprintf
Definition civetweb.c:1579
TH1
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:109
TH2
Service class for 2-D histogram classes.
Definition TH2.h:30
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:1088
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:115
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:624
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:331
TSpectrum2::Background
virtual TH1 * Background(const TH1 *hist, Int_t nIterX=20, Int_t nIterY=20, Option_t *option="")
This function calculates the background spectrum in the input histogram h.
Definition TSpectrum2.cxx:144
TSpectrum2::StaticBackground
static TH1 * StaticBackground(const TH1 *hist, Int_t nIterX=20, Int_t nIterY=20, Option_t *option="")
static function (called by TH2), interface to TSpectrum2::Background
Definition TSpectrum2.cxx:1763
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:59
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
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:1132
TSpectrum2::Print
void Print(Option_t *option="") const override
Print the array of positions.
Definition TSpectrum2.cxx:210
TSpectrum2::fNPeaks
Int_t fNPeaks
number of peaks found
Definition TSpectrum2.h:21
TSpectrum2::SetAverageWindow
static void SetAverageWindow(Int_t w=3)
static function: Set average window of searched peaks see TSpectrum2::SearchHighRes
Definition TSpectrum2.cxx:106
TSpectrum2::~TSpectrum2
~TSpectrum2() override
Destructor.
Definition TSpectrum2.cxx:94
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:250
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 TH2), interface to TSpectrum2::Search
Definition TSpectrum2.cxx:1754
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:877
TSpectrum2::fPositionX
Double_t * fPositionX
[fNPeaks] X position of peaks
Definition TSpectrum2.h:23
TString
Basic string class.
Definition TString.h:138
TString::ToLower
void ToLower()
Change string to lower-case.
Definition TString.cxx:1189
TString::ReplaceAll
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition TString.h:713
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