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