Logo ROOT  

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