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 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 for (i = 0; i < ssizex; i++)
606 delete[]working_space[i];
607 delete [] working_space;
608 return 0;
609 }
610
611 nom = 0;
612 working_space[xmin][ymin] = 1;
613 for(i = xmin; i < xmax; i++){
614 nip = source[i][ymin] / maxch;
615 nim = source[i + 1][ymin] / maxch;
616 sp = 0,sm = 0;
617 for(l = 1; l <= averWindow; l++){
618 if((i + l) > xmax)
619 a = source[xmax][ymin] / maxch;
620
621 else
622 a = source[i + l][ymin] / maxch;
623 b = a - nip;
624 if(a + nip <= 0)
625 a = 1;
626
627 else
628 a = TMath::Sqrt(a + nip);
629 b = b / a;
630 b = TMath::Exp(b);
631 sp = sp + b;
632 if(i - l + 1 < xmin)
633 a = source[xmin][ymin] / maxch;
634
635 else
636 a = source[i - l + 1][ymin] / maxch;
637 b = a - nim;
638 if(a + nim <= 0)
639 a = 1;
640
641 else
642 a = TMath::Sqrt(a + nim);
643 b = b / a;
644 b = TMath::Exp(b);
645 sm = sm + b;
646 }
647 a = sp / sm;
648 a = working_space[i + 1][ymin] = a * working_space[i][ymin];
649 nom = nom + a;
650 }
651 for(i = ymin; i < ymax; i++){
652 nip = source[xmin][i] / maxch;
653 nim = source[xmin][i + 1] / maxch;
654 sp = 0,sm = 0;
655 for(l = 1; l <= averWindow; l++){
656 if((i + l) > ymax)
657 a = source[xmin][ymax] / maxch;
658
659 else
660 a = source[xmin][i + l] / maxch;
661 b = a - nip;
662 if(a + nip <= 0)
663 a = 1;
664
665 else
666 a = TMath::Sqrt(a + nip);
667 b = b / a;
668 b = TMath::Exp(b);
669 sp = sp + b;
670 if(i - l + 1 < ymin)
671 a = source[xmin][ymin] / maxch;
672
673 else
674 a = source[xmin][i - l + 1] / maxch;
675 b = a - nim;
676 if(a + nim <= 0)
677 a = 1;
678
679 else
680 a = TMath::Sqrt(a + nim);
681 b = b / a;
682 b = TMath::Exp(b);
683 sm = sm + b;
684 }
685 a = sp / sm;
686 a = working_space[xmin][i + 1] = a * working_space[xmin][i];
687 nom = nom + a;
688 }
689 for(i = xmin; i < xmax; i++){
690 for(j = ymin; j < ymax; j++){
691 nip = source[i][j + 1] / maxch;
692 nim = source[i + 1][j + 1] / maxch;
693 spx = 0,smx = 0;
694 for(l = 1; l <= averWindow; l++){
695 if(i + l > xmax)
696 a = source[xmax][j] / maxch;
697
698 else
699 a = source[i + l][j] / maxch;
700 b = a - nip;
701 if(a + nip <= 0)
702 a = 1;
703
704 else
705 a = TMath::Sqrt(a + nip);
706 b = b / a;
707 b = TMath::Exp(b);
708 spx = spx + b;
709 if(i - l + 1 < xmin)
710 a = source[xmin][j] / maxch;
711
712 else
713 a = source[i - l + 1][j] / maxch;
714 b = a - nim;
715 if(a + nim <= 0)
716 a = 1;
717
718 else
719 a = TMath::Sqrt(a + nim);
720 b = b / a;
721 b = TMath::Exp(b);
722 smx = smx + b;
723 }
724 spy = 0,smy = 0;
725 nip = source[i + 1][j] / maxch;
726 nim = source[i + 1][j + 1] / maxch;
727 for (l = 1; l <= averWindow; l++) {
728 if (j + l > ymax) a = source[i][ymax]/maxch;
729 else a = source[i][j + l] / maxch;
730 b = a - nip;
731 if (a + nip <= 0) a = 1;
732 else a = TMath::Sqrt(a + nip);
733 b = b / a;
734 b = TMath::Exp(b);
735 spy = spy + b;
736 if (j - l + 1 < ymin) a = source[i][ymin] / maxch;
737 else a = source[i][j - l + 1] / maxch;
738 b = a - nim;
739 if (a + nim <= 0) a = 1;
740 else a = TMath::Sqrt(a + nim);
741 b = b / a;
742 b = TMath::Exp(b);
743 smy = smy + b;
744 }
745 a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx +smy);
746 working_space[i + 1][j + 1] = a;
747 nom = nom + a;
748 }
749 }
750 for(i = xmin; i <= xmax; i++){
751 for(j = ymin; j <= ymax; j++){
752 working_space[i][j] = working_space[i][j] / nom;
753 }
754 }
755 for(i = 0;i < ssizex; i++){
756 for(j = 0; j < ssizey; j++){
757 source[i][j] = plocha * working_space[i][j];
758 }
759 }
760 for (i = 0; i < ssizex; i++)
761 delete[]working_space[i];
762 delete[]working_space;
763 return 0;
764}
765
766////////////////////////////////////////////////////////////////////////////////
767/// This function calculates deconvolution from source spectrum
768/// according to response spectrum
769/// The result is placed in the matrix pointed by source pointer.
770///
771/// Function parameters:
772/// - source-pointer to the matrix of source spectrum
773/// - resp-pointer to the matrix of response spectrum
774/// - ssizex-x length of source and response spectra
775/// - ssizey-y length of source and response spectra
776/// - numberIterations, for details we refer to manual
777/// - numberRepetitions, for details we refer to manual
778/// - boost, boosting factor, for details we refer to manual
779///
780/// ### Deconvolution
781///
782/// Goal: Improvement of the resolution in spectra, decomposition of multiplets
783///
784/// Mathematical formulation of the 2-dimensional convolution system is
785///
786///\f[
787/// 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)
788///\f]
789///\f[
790/// i_1 = 0,1,2, ... ,N_1-1, i_2 = 0,1,2, ... ,N_2-1
791///\f]
792///
793/// where h(i,j) is the impulse response function, x, y are input and output
794/// matrices, respectively, \f$ N_1, N_2\f$ are the lengths of x and h matrices
795///
796/// - let us assume that we know the response and the output matrices (spectra) of the above given system.
797/// - the deconvolution represents solution of the overdetermined system of linear equations, i.e., the
798/// calculation of the matrix x.
799/// - from numerical stability point of view the operation of deconvolution is
800/// extremely critical (ill-posed problem) as well as time consuming operation.
801/// - the Gold deconvolution algorithm proves to work very well even for 2-dimensional
802/// systems. Generalisation of the algorithm for 2-dimensional systems was presented in [1], [2].
803/// - for Gold deconvolution algorithm as well as for boosted deconvolution algorithm we refer also to TSpectrum
804///
805/// #### References:
806///
807/// [1] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:
808/// Efficient one- and two-dimensional Gold deconvolution and its application to
809/// gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
810///
811/// [2] Morhac; M., Matouoek V., Kliman J., Efficient algorithm of multidimensional
812/// deconvolution and its application to nuclear data processing, Digital Signal
813/// Processing 13 (2003) 144.
814///
815/// ### Example 5 - Deconvolution2_1.c
816///
817/// Begin_Macro(source)
818/// ../../../tutorials/spectrum/Deconvolution2_1.C
819/// End_Macro
820///
821/// ### Example 6 - Deconvolution2_2.C
822///
823/// Begin_Macro(source)
824/// ../../../tutorials/spectrum/Deconvolution2_2.C
825/// End_Macro
826///
827/// ### Example 7 - Deconvolution2_HR.C
828///
829/// Begin_Macro(source)
830/// ../../../tutorials/spectrum/Deconvolution2_HR.C
831/// End_Macro
832
833const char *TSpectrum2::Deconvolution(Double_t **source, Double_t **resp,
834 Int_t ssizex, Int_t ssizey,
835 Int_t numberIterations,
836 Int_t numberRepetitions,
837 Double_t boost)
838{
839 Int_t i, j, lhx, lhy, i1, i2, j1, j2, k1, k2, lindex, i1min, i1max,
840 i2min, i2max, j1min, j1max, j2min, j2max, positx = 0, posity = 0, repet;
841 Double_t lda, ldb, ldc, area, maximum = 0;
842 if (ssizex <= 0 || ssizey <= 0)
843 return "Wrong parameters";
844 if (numberIterations <= 0)
845 return "Number of iterations must be positive";
846 if (numberRepetitions <= 0)
847 return "Number of repetitions must be positive";
848 Double_t **working_space = new Double_t *[ssizex];
849 for (i = 0; i < ssizex; i++)
850 working_space[i] = new Double_t[5 * ssizey];
851 area = 0;
852 lhx = -1, lhy = -1;
853 for (i = 0; i < ssizex; i++) {
854 for (j = 0; j < ssizey; j++) {
855 lda = resp[i][j];
856 if (lda != 0) {
857 if ((i + 1) > lhx)
858 lhx = i + 1;
859 if ((j + 1) > lhy)
860 lhy = j + 1;
861 }
862 working_space[i][j] = lda;
863 area = area + lda;
864 if (lda > maximum) {
865 maximum = lda;
866 positx = i, posity = j;
867 }
868 }
869 }
870 if (lhx == -1 || lhy == -1) {
871 for (i = 0; i < ssizex; i++)
872 delete[]working_space[i];
873 delete [] working_space;
874 return ("Zero response data");
875 }
876
877//calculate ht*y and write into p
878 for (i2 = 0; i2 < ssizey; i2++) {
879 for (i1 = 0; i1 < ssizex; i1++) {
880 ldc = 0;
881 for (j2 = 0; j2 <= (lhy - 1); j2++) {
882 for (j1 = 0; j1 <= (lhx - 1); j1++) {
883 k2 = i2 + j2, k1 = i1 + j1;
884 if (k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
885 lda = working_space[j1][j2];
886 ldb = source[k1][k2];
887 ldc = ldc + lda * ldb;
888 }
889 }
890 }
891 working_space[i1][i2 + ssizey] = ldc;
892 }
893 }
894
895//calculate matrix b=ht*h
896 i1min = -(lhx - 1), i1max = lhx - 1;
897 i2min = -(lhy - 1), i2max = lhy - 1;
898 for (i2 = i2min; i2 <= i2max; i2++) {
899 for (i1 = i1min; i1 <= i1max; i1++) {
900 ldc = 0;
901 j2min = -i2;
902 if (j2min < 0)
903 j2min = 0;
904 j2max = lhy - 1 - i2;
905 if (j2max > lhy - 1)
906 j2max = lhy - 1;
907 for (j2 = j2min; j2 <= j2max; j2++) {
908 j1min = -i1;
909 if (j1min < 0)
910 j1min = 0;
911 j1max = lhx - 1 - i1;
912 if (j1max > lhx - 1)
913 j1max = lhx - 1;
914 for (j1 = j1min; j1 <= j1max; j1++) {
915 lda = working_space[j1][j2];
916 if (i1 + j1 < ssizex && i2 + j2 < ssizey)
917 ldb = working_space[i1 + j1][i2 + j2];
918 else
919 ldb = 0;
920 ldc = ldc + lda * ldb;
921 }
922 }
923 working_space[i1 - i1min][i2 - i2min + 2 * ssizey ] = ldc;
924 }
925 }
926
927//initialization in x1 matrix
928 for (i2 = 0; i2 < ssizey; i2++) {
929 for (i1 = 0; i1 < ssizex; i1++) {
930 working_space[i1][i2 + 3 * ssizey] = 1;
931 working_space[i1][i2 + 4 * ssizey] = 0;
932 }
933 }
934
935 //START OF ITERATIONS
936 for (repet = 0; repet < numberRepetitions; repet++) {
937 if (repet != 0) {
938 for (i = 0; i < ssizex; i++) {
939 for (j = 0; j < ssizey; j++) {
940 working_space[i][j + 3 * ssizey] =
941 TMath::Power(working_space[i][j + 3 * ssizey], boost);
942 }
943 }
944 }
945 for (lindex = 0; lindex < numberIterations; lindex++) {
946 for (i2 = 0; i2 < ssizey; i2++) {
947 for (i1 = 0; i1 < ssizex; i1++) {
948 ldb = 0;
949 j2min = i2;
950 if (j2min > lhy - 1)
951 j2min = lhy - 1;
952 j2min = -j2min;
953 j2max = ssizey - i2 - 1;
954 if (j2max > lhy - 1)
955 j2max = lhy - 1;
956 j1min = i1;
957 if (j1min > lhx - 1)
958 j1min = lhx - 1;
959 j1min = -j1min;
960 j1max = ssizex - i1 - 1;
961 if (j1max > lhx - 1)
962 j1max = lhx - 1;
963 for (j2 = j2min; j2 <= j2max; j2++) {
964 for (j1 = j1min; j1 <= j1max; j1++) {
965 ldc = working_space[j1 - i1min][j2 - i2min + 2 * ssizey];
966 lda = working_space[i1 + j1][i2 + j2 + 3 * ssizey];
967 ldb = ldb + lda * ldc;
968 }
969 }
970 lda = working_space[i1][i2 + 3 * ssizey];
971 ldc = working_space[i1][i2 + 1 * ssizey];
972 if (ldc * lda != 0 && ldb != 0) {
973 lda = lda * ldc / ldb;
974 }
975
976 else
977 lda = 0;
978 working_space[i1][i2 + 4 * ssizey] = lda;
979 }
980 }
981 for (i2 = 0; i2 < ssizey; i2++) {
982 for (i1 = 0; i1 < ssizex; i1++)
983 working_space[i1][i2 + 3 * ssizey] =
984 working_space[i1][i2 + 4 * ssizey];
985 }
986 }
987 }
988 for (i = 0; i < ssizex; i++) {
989 for (j = 0; j < ssizey; j++)
990 source[(i + positx) % ssizex][(j + posity) % ssizey] =
991 area * working_space[i][j + 3 * ssizey];
992 }
993 for (i = 0; i < ssizex; i++)
994 delete[]working_space[i];
995 delete[]working_space;
996 return 0;
997}
998
999////////////////////////////////////////////////////////////////////////////////
1000/// This function searches for peaks in source spectrum
1001/// It is based on deconvolution method. First the background is
1002/// removed (if desired), then Markov spectrum is calculated
1003/// (if desired), then the response function is generated
1004/// according to given sigma and deconvolution is carried out.
1005///
1006/// Function parameters:
1007/// - source-pointer to the matrix of source spectrum
1008/// - dest-pointer to the matrix of resulting deconvolved spectrum
1009/// - ssizex-x length of source spectrum
1010/// - ssizey-y length of source spectrum
1011/// - sigma-sigma of searched peaks, for details we refer to manual
1012/// - threshold-threshold value in % for selected peaks, peaks with
1013/// amplitude less than threshold*highest_peak/100
1014/// are ignored, see manual
1015/// - backgroundRemove-logical variable, set if the removal of
1016/// background before deconvolution is desired
1017/// - deconIterations-number of iterations in deconvolution operation
1018/// - markov-logical variable, if it is true, first the source spectrum
1019/// is replaced by new spectrum calculated using Markov
1020/// chains method.
1021/// - averWindow-averaging window of searched peaks, for details
1022/// we refer to manual (applies only for Markov method)
1023///
1024/// ### Peaks searching
1025///
1026/// Goal: to identify automatically the peaks in spectrum with the presence of the
1027/// continuous background, one-fold coincidences (ridges) and statistical
1028/// fluctuations - noise.
1029///
1030/// The common problems connected with correct peak identification in two-dimensional coincidence spectra are
1031///
1032/// - non-sensitivity to noise, i.e., only statistically relevant peaks should be identified
1033/// - non-sensitivity of the algorithm to continuous background
1034/// - non-sensitivity to one-fold coincidences (coincidences peak - background in both dimensions) and their crossings
1035/// - ability to identify peaks close to the edges of the spectrum region. Usually peak finders fail to detect them
1036/// - resolution, decomposition of doublets and multiplets. The algorithm should be able to recognise close positioned peaks.
1037/// - ability to identify peaks with different sigma
1038///
1039/// #### References:
1040///
1041/// [1] M.A. Mariscotti: A method for identification of peaks in the presence of
1042/// background and its application to spectrum analysis. NIM 50 (1967), 309-320.
1043///
1044/// [2] M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:Identification
1045/// of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
1046/// 108-125.
1047///
1048/// [3] Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
1049/// (1996), 451.
1050///
1051/// ### Examples of peak searching method
1052///
1053/// SearchHighRes function provides users with the possibility
1054/// to vary the input parameters and with the access to the output deconvolved data
1055/// in the destination spectrum. Based on the output data one can tune the
1056/// parameters.
1057///
1058/// ### Example 8 - Src.C
1059///
1060/// Begin_Macro(source)
1061/// ../../../tutorials/spectrum/Src.C
1062/// End_Macro
1063///
1064/// ### Example 9 - Src2.C
1065///
1066/// Begin_Macro(source)
1067/// ../../../tutorials/spectrum/Src2.C
1068/// End_Macro
1069///
1070/// ### Example 10 - Src3.C
1071///
1072/// Begin_Macro(source)
1073/// ../../../tutorials/spectrum/Src3.C
1074/// End_Macro
1075///
1076/// ### Example 11 - Src4.C
1077///
1078/// Begin_Macro(source)
1079/// ../../../tutorials/spectrum/Src4.C
1080/// End_Macro
1081///
1082/// ### Example 12 - Src5.C
1083///
1084/// Begin_Macro(source)
1085/// ../../../tutorials/spectrum/Src5.C
1086/// End_Macro
1087
1088Int_t TSpectrum2::SearchHighRes(Double_t **source, Double_t **dest, Int_t ssizex, Int_t ssizey,
1089 Double_t sigma, Double_t threshold,
1090 Bool_t backgroundRemove,Int_t deconIterations,
1091 Bool_t markov, Int_t averWindow)
1092
1093{
1094 Int_t number_of_iterations = (Int_t)(4 * sigma + 0.5);
1095 Int_t k, lindex, priz;
1096 Double_t lda, ldb, ldc, area, maximum;
1097 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;
1098 Int_t ymin, ymax, i, j;
1099 Double_t a, b, ax, ay, maxch, plocha = 0;
1100 Double_t nom, nip, nim, sp, sm, spx, spy, smx, smy;
1101 Double_t p1, p2, p3, p4, s1, s2, s3, s4;
1102 Int_t x, y;
1103 Int_t lhx, lhy, i1, i2, j1, j2, k1, k2, i1min, i1max, i2min, i2max, j1min, j1max, j2min, j2max, positx, posity;
1104 if (sigma < 1) {
1105 Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
1106 return 0;
1107 }
1108
1109 if(threshold<=0||threshold>=100){
1110 Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
1111 return 0;
1112 }
1113
1114 j = (Int_t) (5.0 * sigma + 0.5);
1115 if (j >= PEAK_WINDOW / 2) {
1116 Error("SearchHighRes", "Too large sigma");
1117 return 0;
1118 }
1119
1120 if (markov == true) {
1121 if (averWindow <= 0) {
1122 Error("SearchHighRes", "Averaging window must be positive");
1123 return 0;
1124 }
1125 }
1126 if(backgroundRemove == true){
1127 if(ssizex_ext < 2 * number_of_iterations + 1 || ssizey_ext < 2 * number_of_iterations + 1){
1128 Error("SearchHighRes", "Too large clipping window");
1129 return 0;
1130 }
1131 }
1132 i = (Int_t)(4 * sigma + 0.5);
1133 i = 4 * i;
1134 Double_t **working_space = new Double_t *[ssizex + i];
1135 for (j = 0; j < ssizex + i; j++) {
1136 Double_t *wsk = working_space[j] = new Double_t[16 * (ssizey + i)];
1137 for (k=0;k<16 * (ssizey + i);k++) wsk[k] = 0;
1138 }
1139 for(j = 0; j < ssizey_ext; j++){
1140 for(i = 0; i < ssizex_ext; i++){
1141 if(i < shift){
1142 if(j < shift)
1143 working_space[i][j + ssizey_ext] = source[0][0];
1144
1145 else if(j >= ssizey + shift)
1146 working_space[i][j + ssizey_ext] = source[0][ssizey - 1];
1147
1148 else
1149 working_space[i][j + ssizey_ext] = source[0][j - shift];
1150 }
1151
1152 else if(i >= ssizex + shift){
1153 if(j < shift)
1154 working_space[i][j + ssizey_ext] = source[ssizex - 1][0];
1155
1156 else if(j >= ssizey + shift)
1157 working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1];
1158
1159 else
1160 working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift];
1161 }
1162
1163 else{
1164 if(j < shift)
1165 working_space[i][j + ssizey_ext] = source[i - shift][0];
1166
1167 else if(j >= ssizey + shift)
1168 working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1];
1169
1170 else
1171 working_space[i][j + ssizey_ext] = source[i - shift][j - shift];
1172 }
1173 }
1174 }
1175 if(backgroundRemove == true){
1176 for(i = 1; i <= number_of_iterations; i++){
1177 for(y = i; y < ssizey_ext - i; y++){
1178 for(x = i; x < ssizex_ext - i; x++){
1179 a = working_space[x][y + ssizey_ext];
1180 p1 = working_space[x - i][y + ssizey_ext - i];
1181 p2 = working_space[x - i][y + ssizey_ext + i];
1182 p3 = working_space[x + i][y + ssizey_ext - i];
1183 p4 = working_space[x + i][y + ssizey_ext + i];
1184 s1 = working_space[x][y + ssizey_ext - i];
1185 s2 = working_space[x - i][y + ssizey_ext];
1186 s3 = working_space[x + i][y + ssizey_ext];
1187 s4 = working_space[x][y + ssizey_ext + i];
1188 b = (p1 + p2) / 2.0;
1189 if(b > s2)
1190 s2 = b;
1191 b = (p1 + p3) / 2.0;
1192 if(b > s1)
1193 s1 = b;
1194 b = (p2 + p4) / 2.0;
1195 if(b > s4)
1196 s4 = b;
1197 b = (p3 + p4) / 2.0;
1198 if(b > s3)
1199 s3 = b;
1200 s1 = s1 - (p1 + p3) / 2.0;
1201 s2 = s2 - (p1 + p2) / 2.0;
1202 s3 = s3 - (p3 + p4) / 2.0;
1203 s4 = s4 - (p2 + p4) / 2.0;
1204 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
1205 if(b < a)
1206 a = b;
1207 working_space[x][y] = a;
1208 }
1209 }
1210 for(y = i;y < ssizey_ext - i; y++){
1211 for(x = i; x < ssizex_ext - i; x++){
1212 working_space[x][y + ssizey_ext] = working_space[x][y];
1213 }
1214 }
1215 }
1216 for(j = 0;j < ssizey_ext; j++){
1217 for(i = 0; i < ssizex_ext; i++){
1218 if(i < shift){
1219 if(j < shift)
1220 working_space[i][j + ssizey_ext] = source[0][0] - working_space[i][j + ssizey_ext];
1221
1222 else if(j >= ssizey + shift)
1223 working_space[i][j + ssizey_ext] = source[0][ssizey - 1] - working_space[i][j + ssizey_ext];
1224
1225 else
1226 working_space[i][j + ssizey_ext] = source[0][j - shift] - working_space[i][j + ssizey_ext];
1227 }
1228
1229 else if(i >= ssizex + shift){
1230 if(j < shift)
1231 working_space[i][j + ssizey_ext] = source[ssizex - 1][0] - working_space[i][j + ssizey_ext];
1232
1233 else if(j >= ssizey + shift)
1234 working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1] - working_space[i][j + ssizey_ext];
1235
1236 else
1237 working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift] - working_space[i][j + ssizey_ext];
1238 }
1239
1240 else{
1241 if(j < shift)
1242 working_space[i][j + ssizey_ext] = source[i - shift][0] - working_space[i][j + ssizey_ext];
1243
1244 else if(j >= ssizey + shift)
1245 working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1] - working_space[i][j + ssizey_ext];
1246
1247 else
1248 working_space[i][j + ssizey_ext] = source[i - shift][j - shift] - working_space[i][j + ssizey_ext];
1249 }
1250 }
1251 }
1252 }
1253 for(j = 0; j < ssizey_ext; j++){
1254 for(i = 0; i < ssizex_ext; i++){
1255 working_space[i][j + 15*ssizey_ext] = working_space[i][j + ssizey_ext];
1256 }
1257 }
1258 if(markov == true){
1259 for(i = 0;i < ssizex_ext; i++){
1260 for(j = 0; j < ssizey_ext; j++)
1261 working_space[i][j + 2 * ssizex_ext] = working_space[i][ssizey_ext + j];
1262 }
1263 xmin = 0;
1264 xmax = ssizex_ext - 1;
1265 ymin = 0;
1266 ymax = ssizey_ext - 1;
1267 for(i = 0, maxch = 0; i < ssizex_ext; i++){
1268 for(j = 0; j < ssizey_ext; j++){
1269 working_space[i][j] = 0;
1270 if(maxch < working_space[i][j + 2 * ssizey_ext])
1271 maxch = working_space[i][j + 2 * ssizey_ext];
1272 plocha += working_space[i][j + 2 * ssizey_ext];
1273 }
1274 }
1275 if(maxch == 0) {
1276 i = (Int_t)(4 * sigma + 0.5);
1277 i = 4 * i;
1278 for (j = 0; j < ssizex + i; j++)
1279 delete[]working_space[j];
1280 delete [] working_space;
1281 return 0;
1282 }
1283
1284 nom=0;
1285 working_space[xmin][ymin] = 1;
1286 for(i = xmin; i < xmax; i++){
1287 nip = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1288 nim = working_space[i + 1][ymin + 2 * ssizey_ext] / maxch;
1289 sp = 0,sm = 0;
1290 for(l = 1;l <= averWindow; l++){
1291 if((i + l) > xmax)
1292 a = working_space[xmax][ymin + 2 * ssizey_ext] / maxch;
1293
1294 else
1295 a = working_space[i + l][ymin + 2 * ssizey_ext] / maxch;
1296
1297 b = a - nip;
1298 if(a + nip <= 0)
1299 a = 1;
1300
1301 else
1302 a=TMath::Sqrt(a + nip);
1303 b = b / a;
1304 b = TMath::Exp(b);
1305 sp = sp + b;
1306 if(i - l + 1 < xmin)
1307 a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1308
1309 else
1310 a = working_space[i - l + 1][ymin + 2 * ssizey_ext] / maxch;
1311 b = a - nim;
1312 if(a + nim <= 0)
1313 a = 1;
1314
1315 else
1316 a=TMath::Sqrt(a + nim);
1317 b = b / a;
1318 b = TMath::Exp(b);
1319 sm = sm + b;
1320 }
1321 a = sp / sm;
1322 a = working_space[i + 1][ymin] = a * working_space[i][ymin];
1323 nom = nom + a;
1324 }
1325 for(i = ymin; i < ymax; i++){
1326 nip = working_space[xmin][i + 2 * ssizey_ext] / maxch;
1327 nim = working_space[xmin][i + 1 + 2 * ssizey_ext] / maxch;
1328 sp = 0,sm = 0;
1329 for(l = 1; l <= averWindow; l++){
1330 if((i + l) > ymax)
1331 a = working_space[xmin][ymax + 2 * ssizey_ext] / maxch;
1332
1333 else
1334 a = working_space[xmin][i + l + 2 * ssizey_ext] / maxch;
1335 b = a - nip;
1336 if(a + nip <= 0)
1337 a=1;
1338
1339 else
1340 a=TMath::Sqrt(a + nip);
1341 b = b / a;
1342 b = TMath::Exp(b);
1343 sp = sp + b;
1344 if(i - l + 1 < ymin)
1345 a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
1346
1347 else
1348 a = working_space[xmin][i - l + 1 + 2 * ssizey_ext] / maxch;
1349 b = a - nim;
1350 if(a + nim <= 0)
1351 a = 1;
1352
1353 else
1354 a=TMath::Sqrt(a + nim);
1355 b = b / a;
1356 b = TMath::Exp(b);
1357 sm = sm + b;
1358 }
1359 a = sp / sm;
1360 a = working_space[xmin][i + 1] = a * working_space[xmin][i];
1361 nom = nom + a;
1362 }
1363 for(i = xmin; i < xmax; i++){
1364 for(j = ymin; j < ymax; j++){
1365 nip = working_space[i][j + 1 + 2 * ssizey_ext] / maxch;
1366 nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1367 spx = 0,smx = 0;
1368 for(l = 1; l <= averWindow; l++){
1369 if(i + l > xmax)
1370 a = working_space[xmax][j + 2 * ssizey_ext] / maxch;
1371
1372 else
1373 a = working_space[i + l][j + 2 * ssizey_ext] / maxch;
1374 b = a - nip;
1375 if(a + nip <= 0)
1376 a = 1;
1377
1378 else
1379 a=TMath::Sqrt(a + nip);
1380 b = b / a;
1381 b = TMath::Exp(b);
1382 spx = spx + b;
1383 if(i - l + 1 < xmin)
1384 a = working_space[xmin][j + 2 * ssizey_ext] / maxch;
1385
1386 else
1387 a = working_space[i - l + 1][j + 2 * ssizey_ext] / maxch;
1388 b = a - nim;
1389 if(a + nim <= 0)
1390 a=1;
1391
1392 else
1393 a=TMath::Sqrt(a + nim);
1394 b = b / a;
1395 b = TMath::Exp(b);
1396 smx = smx + b;
1397 }
1398 spy = 0,smy = 0;
1399 nip = working_space[i + 1][j + 2 * ssizey_ext] / maxch;
1400 nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
1401 for(l = 1; l <= averWindow; l++){
1402 if(j + l > ymax)
1403 a = working_space[i][ymax + 2 * ssizey_ext] / maxch;
1404
1405 else
1406 a = working_space[i][j + l + 2 * ssizey_ext] / maxch;
1407 b = a - nip;
1408 if(a + nip <= 0)
1409 a = 1;
1410
1411 else
1412 a=TMath::Sqrt(a + nip);
1413 b = b / a;
1414 b = TMath::Exp(b);
1415 spy = spy + b;
1416 if(j - l + 1 < ymin)
1417 a = working_space[i][ymin + 2 * ssizey_ext] / maxch;
1418
1419 else
1420 a = working_space[i][j - l + 1 + 2 * ssizey_ext] / maxch;
1421 b=a-nim;
1422 if(a + nim <= 0)
1423 a = 1;
1424 else
1425 a=TMath::Sqrt(a + nim);
1426 b = b / a;
1427 b = TMath::Exp(b);
1428 smy = smy + b;
1429 }
1430 a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx + smy);
1431 working_space[i + 1][j + 1] = a;
1432 nom = nom + a;
1433 }
1434 }
1435 for(i = xmin; i <= xmax; i++){
1436 for(j = ymin; j <= ymax; j++){
1437 working_space[i][j] = working_space[i][j] / nom;
1438 }
1439 }
1440 for(i = 0; i < ssizex_ext; i++){
1441 for(j = 0; j < ssizey_ext; j++){
1442 working_space[i][j + ssizey_ext] = working_space[i][j] * plocha;
1443 working_space[i][2 * ssizey_ext + j] = working_space[i][ssizey_ext + j];
1444 }
1445 }
1446 }
1447 //deconvolution starts
1448 area = 0;
1449 lhx = -1,lhy = -1;
1450 positx = 0,posity = 0;
1451 maximum = 0;
1452 //generate response matrix
1453 for(i = 0; i < ssizex_ext; i++){
1454 for(j = 0; j < ssizey_ext; j++){
1455 lda = (Double_t)i - 3 * sigma;
1456 ldb = (Double_t)j - 3 * sigma;
1457 lda = (lda * lda + ldb * ldb) / (2 * sigma * sigma);
1458 k=(Int_t)(1000 * TMath::Exp(-lda));
1459 lda = k;
1460 if(lda != 0){
1461 if((i + 1) > lhx)
1462 lhx = i + 1;
1463
1464 if((j + 1) > lhy)
1465 lhy = j + 1;
1466 }
1467 working_space[i][j] = lda;
1468 area = area + lda;
1469 if(lda > maximum){
1470 maximum = lda;
1471 positx = i,posity = j;
1472 }
1473 }
1474 }
1475 //read source matrix
1476 for(i = 0;i < ssizex_ext; i++){
1477 for(j = 0;j < ssizey_ext; j++){
1478 working_space[i][j + 14 * ssizey_ext] = TMath::Abs(working_space[i][j + ssizey_ext]);
1479 }
1480 }
1481 //calculate matrix b=ht*h
1482 i = lhx - 1;
1483 if(i > ssizex_ext)
1484 i = ssizex_ext;
1485
1486 j = lhy - 1;
1487 if(j>ssizey_ext)
1488 j = ssizey_ext;
1489
1490 i1min = -i,i1max = i;
1491 i2min = -j,i2max = j;
1492 for(i2 = i2min; i2 <= i2max; i2++){
1493 for(i1 = i1min; i1 <= i1max; i1++){
1494 ldc = 0;
1495 j2min = -i2;
1496 if(j2min < 0)
1497 j2min = 0;
1498
1499 j2max = lhy - 1 - i2;
1500 if(j2max > lhy - 1)
1501 j2max = lhy - 1;
1502
1503 for(j2 = j2min; j2 <= j2max; j2++){
1504 j1min = -i1;
1505 if(j1min < 0)
1506 j1min = 0;
1507
1508 j1max = lhx - 1 - i1;
1509 if(j1max > lhx - 1)
1510 j1max = lhx - 1;
1511
1512 for(j1 = j1min; j1 <= j1max; j1++){
1513 lda = working_space[j1][j2];
1514 ldb = working_space[i1 + j1][i2 + j2];
1515 ldc = ldc + lda * ldb;
1516 }
1517 }
1518 k = (i1 + ssizex_ext) / ssizex_ext;
1519 working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext] = ldc;
1520 }
1521 }
1522 //calculate at*y and write into p
1523 i = lhx - 1;
1524 if(i > ssizex_ext)
1525 i = ssizex_ext;
1526
1527 j = lhy - 1;
1528 if(j > ssizey_ext)
1529 j = ssizey_ext;
1530
1531 i2min = -j,i2max = ssizey_ext + j - 1;
1532 i1min = -i,i1max = ssizex_ext + i - 1;
1533 for(i2 = i2min; i2 <= i2max; i2++){
1534 for(i1=i1min;i1<=i1max;i1++){
1535 ldc=0;
1536 for(j2 = 0; j2 <= (lhy - 1); j2++){
1537 for(j1 = 0; j1 <= (lhx - 1); j1++){
1538 k2 = i2 + j2,k1 = i1 + j1;
1539 if(k2 >= 0 && k2 < ssizey_ext && k1 >= 0 && k1 < ssizex_ext){
1540 lda = working_space[j1][j2];
1541 ldb = working_space[k1][k2 + 14 * ssizey_ext];
1542 ldc = ldc + lda * ldb;
1543 }
1544 }
1545 }
1546 k = (i1 + ssizex_ext) / ssizex_ext;
1547 working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext] = ldc;
1548 }
1549 }
1550 //move matrix p
1551 for(i2 = 0; i2 < ssizey_ext; i2++){
1552 for(i1 = 0; i1 < ssizex_ext; i1++){
1553 k = (i1 + ssizex_ext) / ssizex_ext;
1554 ldb = working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext];
1555 working_space[i1][i2 + 14 * ssizey_ext] = ldb;
1556 }
1557 }
1558 //initialization in x1 matrix
1559 for(i2 = 0; i2 < ssizey_ext; i2++){
1560 for(i1 = 0; i1 < ssizex_ext; i1++){
1561 working_space[i1][i2 + ssizey_ext] = 1;
1562 working_space[i1][i2 + 2 * ssizey_ext] = 0;
1563 }
1564 }
1565 //START OF ITERATIONS
1566 for(lindex = 0; lindex < deconIterations; lindex++){
1567 for(i2 = 0; i2 < ssizey_ext; i2++){
1568 for(i1 = 0; i1 < ssizex_ext; i1++){
1569 lda = working_space[i1][i2 + ssizey_ext];
1570 ldc = working_space[i1][i2 + 14 * ssizey_ext];
1571 if(lda > 0.000001 && ldc > 0.000001){
1572 ldb=0;
1573 j2min=i2;
1574 if(j2min > lhy - 1)
1575 j2min = lhy - 1;
1576
1577 j2min = -j2min;
1578 j2max = ssizey_ext - i2 - 1;
1579 if(j2max > lhy - 1)
1580 j2max = lhy - 1;
1581
1582 j1min = i1;
1583 if(j1min > lhx - 1)
1584 j1min = lhx - 1;
1585
1586 j1min = -j1min;
1587 j1max = ssizex_ext - i1 - 1;
1588 if(j1max > lhx - 1)
1589 j1max = lhx - 1;
1590
1591 for(j2 = j2min; j2 <= j2max; j2++){
1592 for(j1 = j1min; j1 <= j1max; j1++){
1593 k = (j1 + ssizex_ext) / ssizex_ext;
1594 ldc = working_space[(j1 + ssizex_ext) % ssizex_ext][j2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext];
1595 lda = working_space[i1 + j1][i2 + j2 + ssizey_ext];
1596 ldb = ldb + lda * ldc;
1597 }
1598 }
1599 lda = working_space[i1][i2 + ssizey_ext];
1600 ldc = working_space[i1][i2 + 14 * ssizey_ext];
1601 if(ldc * lda != 0 && ldb != 0){
1602 lda =lda * ldc / ldb;
1603 }
1604
1605 else
1606 lda=0;
1607 working_space[i1][i2 + 2 * ssizey_ext] = lda;
1608 }
1609 }
1610 }
1611 for(i2 = 0; i2 < ssizey_ext; i2++){
1612 for(i1 = 0; i1 < ssizex_ext; i1++)
1613 working_space[i1][i2 + ssizey_ext] = working_space[i1][i2 + 2 * ssizey_ext];
1614 }
1615 }
1616 //looking for maximum
1617 maximum=0;
1618 for(i = 0; i < ssizex_ext; i++){
1619 for(j = 0; j < ssizey_ext; j++){
1620 working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext] = area * working_space[i][j + ssizey_ext];
1621 if(maximum < working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext])
1622 maximum = working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext];
1623
1624 }
1625 }
1626 //searching for peaks in deconvolved spectrum
1627 for(i = 1; i < ssizex_ext - 1; i++){
1628 for(j = 1; j < ssizey_ext - 1; j++){
1629 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]){
1630 if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift){
1631 if(working_space[i][j] > threshold * maximum / 100.0){
1632 if(peak_index < fMaxPeaks){
1633 for(k = i - 1,a = 0,b = 0; k <= i + 1; k++){
1634 a += (Double_t)(k - shift) * working_space[k][j];
1635 b += working_space[k][j];
1636 }
1637 ax=a/b;
1638 if(ax < 0)
1639 ax = 0;
1640
1641 if(ax >= ssizex)
1642 ax = ssizex - 1;
1643
1644 for(k = j - 1,a = 0,b = 0; k <= j + 1; k++){
1645 a += (Double_t)(k - shift) * working_space[i][k];
1646 b += working_space[i][k];
1647 }
1648 ay=a/b;
1649 if(ay < 0)
1650 ay = 0;
1651
1652 if(ay >= ssizey)
1653 ay = ssizey - 1;
1654
1655 if(peak_index == 0){
1656 fPositionX[0] = ax;
1657 fPositionY[0] = ay;
1658 peak_index = 1;
1659 }
1660
1661 else{
1662 for(k = 0, priz = 0; k < peak_index && priz == 0; k++){
1663 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)])
1664 priz = 1;
1665 }
1666 if(priz == 0){
1667 if(k < fMaxPeaks){
1668 fPositionX[k] = ax;
1669 fPositionY[k] = ay;
1670 }
1671 }
1672
1673 else{
1674 for(l = peak_index; l >= k; l--){
1675 if(l < fMaxPeaks){
1676 fPositionX[l] = fPositionX[l - 1];
1677 fPositionY[l] = fPositionY[l - 1];
1678 }
1679 }
1680 fPositionX[k - 1] = ax;
1681 fPositionY[k - 1] = ay;
1682 }
1683 if(peak_index < fMaxPeaks)
1684 peak_index += 1;
1685 }
1686 }
1687 }
1688 }
1689 }
1690 }
1691 }
1692 //writing back deconvolved spectrum
1693 for(i = 0; i < ssizex; i++){
1694 for(j = 0; j < ssizey; j++){
1695 dest[i][j] = working_space[i + shift][j + shift];
1696 }
1697 }
1698 i = (Int_t)(4 * sigma + 0.5);
1699 i = 4 * i;
1700 for (j = 0; j < ssizex + i; j++)
1701 delete[]working_space[j];
1702 delete[]working_space;
1703 fNPeaks = peak_index;
1704 return fNPeaks;
1705}
1706
1707////////////////////////////////////////////////////////////////////////////////
1708/// static function (called by TH1), interface to TSpectrum2::Search
1709
1710Int_t TSpectrum2::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
1711{
1712 TSpectrum2 s;
1713 return s.Search(hist,sigma,option,threshold);
1714}
1715
1716////////////////////////////////////////////////////////////////////////////////
1717/// static function (called by TH1), interface to TSpectrum2::Background
1718
1719TH1 *TSpectrum2::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
1720{
1721 TSpectrum2 s;
1722 return s.Background(hist,niter,option);
1723}
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
const Bool_t kFALSE
Definition RtypesCore.h:92
Double_t
double Double_t
Definition RtypesCore.h:59
kTRUE
const Bool_t kTRUE
Definition RtypesCore.h:91
Option_t
const char Option_t
Definition RtypesCore.h:66
ClassImp
#define ClassImp(name)
Definition Rtypes.h:364
kRed
@ kRed
Definition Rtypes.h:66
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:45
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:478
TAxis::GetNbins
Int_t GetNbins() const
Definition TAxis.h:121
TH1
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
TH1::GetDimension
virtual Int_t GetDimension() const
Definition TH1.h:282
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition TH1.h:320
TH1::GetYaxis
TAxis * GetYaxis()
Definition TH1.h:321
TH1::GetListOfFunctions
TList * GetListOfFunctions() const
Definition TH1.h:243
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:4993
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:822
TList::FindObject
virtual TObject * FindObject(const char *name) const
Find an object in this list using its name.
Definition TList.cxx:578
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:893
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:1719
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:1088
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:1710
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:833
TSpectrum2::fPositionX
Double_t * fPositionX
[fNPeaks] X position of peaks
Definition TSpectrum2.h:23
TSpectrum2::kBackSuccessiveFiltering
@ kBackSuccessiveFiltering
Definition TSpectrum2.h:34
TSpectrum2::kBackOneStepFiltering
@ kBackOneStepFiltering
Definition TSpectrum2.h:35
TSpectrum2::kBackDecreasingWindow
@ kBackDecreasingWindow
Definition TSpectrum2.h:33
TSpectrum2::kBackIncreasingWindow
@ kBackIncreasingWindow
Definition TSpectrum2.h:32
TString
Basic string class.
Definition TString.h:136
TString::ToLower
void ToLower()
Change string to lower-case.
Definition TString.cxx:1145
TString::ReplaceAll
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition TString.h:692
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:624
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)
Definition TMathBase.h:212
TMath::Exp
Double_t Exp(Double_t x)
Definition TMath.h:727
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition TMath.h:691
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition TMath.h:735
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
dest
#define dest(otri, vertexptr)
Definition triangle.c:1040