Logo ROOT  

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