Logo ROOT   6.16/01
Reference Guide
TSpectrum.cxx
Go to the documentation of this file.
1// @(#)root/spectrum:$Id$
2// Author: Miroslav Morhac 27/05/99
3
4#include "TSpectrum.h"
5#include "TPolyMarker.h"
6#include "TVirtualPad.h"
7#include "TList.h"
8#include "TH1.h"
9#include "TMath.h"
10
11/** \class TSpectrum
12 \ingroup Spectrum
13 \brief Advanced Spectra Processing
14 \author Miroslav Morhac
15
16 This class contains advanced spectra processing functions for:
17
18 - One-dimensional background estimation
19 - One-dimensional smoothing
20 - One-dimensional deconvolution
21 - One-dimensional peak search
22
23 The algorithms in this class have been published in the following references:
24
25 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.
26 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.
27 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.
28
29 These NIM papers are also available as doc or ps files from:
30
31 - [Spectrum.doc](https://root.cern.ch/download/Spectrum.doc)
32 - [SpectrumDec.ps.gz](https://root.cern.ch/download/SpectrumDec.ps.gz)
33 - [SpectrumSrc.ps.gz](https://root.cern.ch/download/SpectrumSrc.ps.gz)
34 - [SpectrumBck.ps.gz](https://root.cern.ch/download/SpectrumBck.ps.gz)
35
36 See also the
37 [online documentation](https://root.cern.ch/guides/tspectrum-manual) and
38 [tutorials](https://root.cern.ch/doc/master/group__tutorial__spectrum.html).
39*/
40
41Int_t TSpectrum::fgIterations = 3;
42Int_t TSpectrum::fgAverageWindow = 3;
43
44#define PEAK_WINDOW 1024
45ClassImp(TSpectrum);
46
47////////////////////////////////////////////////////////////////////////////////
48/// Constructor.
49
50TSpectrum::TSpectrum() :TNamed("Spectrum", "Miroslav Morhac peak finder")
51{
52 Int_t n = 100;
53 fMaxPeaks = n;
54 fPosition = new Double_t[n];
55 fPositionX = new Double_t[n];
56 fPositionY = new Double_t[n];
57 fResolution = 1;
58 fHistogram = 0;
59 fNPeaks = 0;
60}
61
62////////////////////////////////////////////////////////////////////////////////
63/// - maxpositions: maximum number of peaks
64/// - resolution: *NOT USED* determines resolution of the neighbouring peaks
65/// default value is 1 correspond to 3 sigma distance
66/// between peaks. Higher values allow higher resolution
67/// (smaller distance between peaks.
68/// May be set later through SetResolution.
69
70TSpectrum::TSpectrum(Int_t maxpositions, Double_t resolution)
71 :TNamed("Spectrum", "Miroslav Morhac peak finder")
72{
73 Int_t n = maxpositions;
74 if (n <= 0) n = 1;
75 fMaxPeaks = n;
76 fPosition = new Double_t[n];
77 fPositionX = new Double_t[n];
78 fPositionY = new Double_t[n];
79 fHistogram = 0;
80 fNPeaks = 0;
81 SetResolution(resolution);
82}
83
84////////////////////////////////////////////////////////////////////////////////
85/// Destructor.
86
87TSpectrum::~TSpectrum()
88{
89 delete [] fPosition;
90 delete [] fPositionX;
91 delete [] fPositionY;
92 delete fHistogram;
93}
94
95////////////////////////////////////////////////////////////////////////////////
96/// Static function: Set average window of searched peaks
97/// (see TSpectrum::SearchHighRes).
98
99void TSpectrum::SetAverageWindow(Int_t w)
100{
101 fgAverageWindow = w;
102}
103
104////////////////////////////////////////////////////////////////////////////////
105/// Static function: Set max number of decon iterations in deconvolution
106/// operation (see TSpectrum::SearchHighRes).
107
108void TSpectrum::SetDeconIterations(Int_t n)
109{
110 fgIterations = n;
111}
112
113////////////////////////////////////////////////////////////////////////////////
114/// One-dimensional background estimation function.
115///
116/// This function calculates the background spectrum in the input histogram h.
117/// The background is returned as a histogram.
118///
119/// #### Parameters:
120///
121/// - h: input 1-d histogram
122/// - numberIterations, (default value = 20)
123/// Increasing numberIterations make the result smoother and lower.
124/// - option: may contain one of the following options:
125///
126/// - to set the direction parameter
127/// "BackIncreasingWindow". By default the direction is BackDecreasingWindow
128/// - filterOrder-order of clipping filter, (default "BackOrder2")
129/// -possible values= "BackOrder4"
130/// "BackOrder6"
131/// "BackOrder8"
132/// - "nosmoothing"- if selected, the background is not smoothed
133/// By default the background is smoothed.
134/// - smoothWindow-width of smoothing window, (default is "BackSmoothing3")
135/// -possible values= "BackSmoothing5"
136/// "BackSmoothing7"
137/// "BackSmoothing9"
138/// "BackSmoothing11"
139/// "BackSmoothing13"
140/// "BackSmoothing15"
141/// - "Compton" if selected the estimation of Compton edge
142/// will be included.
143/// - "same" : if this option is specified, the resulting background
144/// histogram is superimposed on the picture in the current pad.
145///
146/// NOTE that the background is only evaluated in the current range of h.
147/// ie, if h has a bin range (set via `h->GetXaxis()->SetRange(binmin,binmax)`,
148/// the returned histogram will be created with the same number of bins
149/// as the input histogram h, but only bins from `binmin` to `binmax` will be filled
150/// with the estimated background.
151
152TH1 *TSpectrum::Background(const TH1 * h, int numberIterations,
153 Option_t * option)
154{
155 if (h == 0) return 0;
156 Int_t dimension = h->GetDimension();
157 if (dimension > 1) {
158 Error("Search", "Only implemented for 1-d histograms");
159 return 0;
160 }
161 TString opt = option;
162 opt.ToLower();
163
164 //set options
165 Int_t direction = kBackDecreasingWindow;
166 if (opt.Contains("backincreasingwindow")) direction = kBackIncreasingWindow;
167 Int_t filterOrder = kBackOrder2;
168 if (opt.Contains("backorder4")) filterOrder = kBackOrder4;
169 if (opt.Contains("backorder6")) filterOrder = kBackOrder6;
170 if (opt.Contains("backorder8")) filterOrder = kBackOrder8;
171 Bool_t smoothing = kTRUE;
172 if (opt.Contains("nosmoothing")) smoothing = kFALSE;
173 Int_t smoothWindow = kBackSmoothing3;
174 if (opt.Contains("backsmoothing5")) smoothWindow = kBackSmoothing5;
175 if (opt.Contains("backsmoothing7")) smoothWindow = kBackSmoothing7;
176 if (opt.Contains("backsmoothing9")) smoothWindow = kBackSmoothing9;
177 if (opt.Contains("backsmoothing11")) smoothWindow = kBackSmoothing11;
178 if (opt.Contains("backsmoothing13")) smoothWindow = kBackSmoothing13;
179 if (opt.Contains("backsmoothing15")) smoothWindow = kBackSmoothing15;
180 Bool_t compton = kFALSE;
181 if (opt.Contains("compton")) compton = kTRUE;
182
183 Int_t first = h->GetXaxis()->GetFirst();
184 Int_t last = h->GetXaxis()->GetLast();
185 Int_t size = last-first+1;
186 Int_t i;
187 Double_t * source = new Double_t[size];
188 for (i = 0; i < size; i++) source[i] = h->GetBinContent(i + first);
189
190 //find background (source is input and in output contains the background
191 Background(source,size,numberIterations, direction, filterOrder,smoothing,
192 smoothWindow,compton);
193
194 //create output histogram containing background
195 //only bins in the range of the input histogram are filled
196 Int_t nch = strlen(h->GetName());
197 char *hbname = new char[nch+20];
198 snprintf(hbname,nch+20,"%s_background",h->GetName());
199 TH1 *hb = (TH1*)h->Clone(hbname);
200 hb->Reset();
201 hb->GetListOfFunctions()->Delete();
202 hb->SetLineColor(2);
203 for (i=0; i< size; i++) hb->SetBinContent(i+first,source[i]);
204 hb->SetEntries(size);
205
206 //if option "same is specified, draw the result in the pad
207 if (opt.Contains("same")) {
208 if (gPad) delete gPad->GetPrimitive(hbname);
209 hb->Draw("same");
210 }
211 delete [] source;
212 delete [] hbname;
213 return hb;
214}
215
216////////////////////////////////////////////////////////////////////////////////
217/// Print the array of positions.
218
219void TSpectrum::Print(Option_t *) const
220{
221 printf("\nNumber of positions = %d\n",fNPeaks);
222 for (Int_t i=0;i<fNPeaks;i++) {
223 printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
224 }
225}
226
227////////////////////////////////////////////////////////////////////////////////
228/// One-dimensional peak search function
229///
230/// This function searches for peaks in source spectrum in hin
231/// The number of found peaks and their positions are written into
232/// the members fNpeaks and fPositionX.
233/// The search is performed in the current histogram range.
234///
235/// #### Parameters:
236///
237/// - hin: pointer to the histogram of source spectrum
238/// - sigma: sigma of searched peaks, for details we refer to manual
239/// - threshold: (default=0.05) peaks with amplitude less than
240/// threshold*highest_peak are discarded. 0<threshold<1
241///
242/// By default, the background is removed before deconvolution.
243/// Specify the option "nobackground" to not remove the background.
244///
245/// By default the "Markov" chain algorithm is used.
246/// Specify the option "noMarkov" to disable this algorithm
247/// Note that by default the source spectrum is replaced by a new spectrum
248///
249/// By default a polymarker object is created and added to the list of
250/// functions of the histogram. The histogram is drawn with the specified
251/// option and the polymarker object drawn on top of the histogram.
252/// The polymarker coordinates correspond to the npeaks peaks found in
253/// the histogram.
254///
255/// A pointer to the polymarker object can be retrieved later via:
256/// ~~~ {.cpp}
257/// TList *functions = hin->GetListOfFunctions();
258/// TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
259/// ~~~
260/// Specify the option "goff" to disable the storage and drawing of the
261/// polymarker.
262///
263/// To disable the final drawing of the histogram with the search results (in case
264/// you want to draw it yourself) specify "nodraw" in the options parameter.
265
266Int_t TSpectrum::Search(const TH1 * hin, Double_t sigma, Option_t * option,
267 Double_t threshold)
268{
269 if (hin == 0) return 0;
270 Int_t dimension = hin->GetDimension();
271 if (dimension > 2) {
272 Error("Search", "Only implemented for 1-d and 2-d histograms");
273 return 0;
274 }
275 if (threshold <=0 || threshold >= 1) {
276 Warning("Search","threshold must 0<threshold<1, threshold=0.05 assumed");
277 threshold = 0.05;
278 }
279 TString opt = option;
280 opt.ToLower();
281 Bool_t background = kTRUE;
282 if (opt.Contains("nobackground")) {
283 background = kFALSE;
284 opt.ReplaceAll("nobackground","");
285 }
286 Bool_t markov = kTRUE;
287 if (opt.Contains("nomarkov")) {
288 markov = kFALSE;
289 opt.ReplaceAll("nomarkov","");
290 }
291 Bool_t draw = kTRUE;
292 if (opt.Contains("nodraw")) {
293 draw = kFALSE;
294 opt.ReplaceAll("nodraw","");
295 }
296 if (dimension == 1) {
297 Int_t first = hin->GetXaxis()->GetFirst();
298 Int_t last = hin->GetXaxis()->GetLast();
299 Int_t size = last-first+1;
300 Int_t i, bin, npeaks;
301 Double_t * source = new Double_t[size];
302 Double_t * dest = new Double_t[size];
303 for (i = 0; i < size; i++) source[i] = hin->GetBinContent(i + first);
304 if (sigma < 1) {
305 sigma = size/fMaxPeaks;
306 if (sigma < 1) sigma = 1;
307 if (sigma > 8) sigma = 8;
308 }
309 npeaks = SearchHighRes(source, dest, size, sigma, 100*threshold,
310 background, fgIterations, markov, fgAverageWindow);
311
312 for (i = 0; i < npeaks; i++) {
313 bin = first + Int_t(fPositionX[i] + 0.5);
314 fPositionX[i] = hin->GetBinCenter(bin);
315 fPositionY[i] = hin->GetBinContent(bin);
316 }
317 delete [] source;
318 delete [] dest;
319
320 if (opt.Contains("goff"))
321 return npeaks;
322 if (!npeaks) return 0;
323 TPolyMarker * pm =
324 (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
325 if (pm) {
326 hin->GetListOfFunctions()->Remove(pm);
327 delete pm;
328 }
329 pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
330 hin->GetListOfFunctions()->Add(pm);
331 pm->SetMarkerStyle(23);
332 pm->SetMarkerColor(kRed);
333 pm->SetMarkerSize(1.3);
334 opt.ReplaceAll(" ","");
335 opt.ReplaceAll(",","");
336 if (draw)
337 ((TH1*)hin)->Draw(opt.Data());
338 return npeaks;
339 }
340 return 0;
341}
342
343////////////////////////////////////////////////////////////////////////////////
344/// *NOT USED*
345/// resolution: determines resolution of the neighbouring peaks
346/// default value is 1 correspond to 3 sigma distance
347/// between peaks. Higher values allow higher resolution
348/// (smaller distance between peaks.
349/// May be set later through SetResolution.
350
351void TSpectrum::SetResolution(Double_t resolution)
352{
353 if (resolution > 1)
354 fResolution = resolution;
355 else
356 fResolution = 1;
357}
358
359////////////////////////////////////////////////////////////////////////////////
360/// This function calculates background spectrum from source spectrum.
361/// The result is placed in the vector pointed by spe1945ctrum pointer.
362/// The goal is to separate the useful information (peaks) from useless
363/// information (background).
364///
365/// - method is based on Sensitive Nonlinear Iterative Peak (SNIP) clipping
366/// algorithm.
367/// - new value in the channel "i" is calculated
368///
369/// \f[
370/// v_p(i) = min \left\{ v_{p-1}(i)^{\frac{\left[v_{p-1}(i+p)+v_{p-1}(i-p)\right]}{2}} \right\}
371/// \f]
372///
373/// where p = 1, 2, ..., numberIterations. In fact it represents second order
374/// difference filter (-1,2,-1).
375///
376/// One can also change the
377/// direction of the change of the clipping window, the order of the clipping
378/// filter, to include smoothing, to set width of smoothing window and to include
379/// the estimation of Compton edges. On successful completion it returns 0. On
380/// error it returns pointer to the string describing error.
381///
382/// #### Parameters:
383///
384/// - spectrum: pointer to the vector of source spectrum
385/// - ssize: length of the spectrum vector
386/// - numberIterations: maximal width of clipping window,
387/// - direction: direction of change of clipping window.
388/// Possible values: kBackIncreasingWindow, kBackDecreasingWindow
389/// - filterOrder: order of clipping filter.
390/// Possible values: kBackOrder2, kBackOrder4, kBackOrder6, kBackOrder8
391/// - smoothing: logical variable whether the smoothing operation in the
392/// estimation of background will be included.
393/// Possible values: kFALSE, kTRUE
394/// - smoothWindow: width of smoothing window.
395/// Possible values: kBackSmoothing3, kBackSmoothing5, kBackSmoothing7,
396/// kBackSmoothing9, kBackSmoothing11, kBackSmoothing13, kBackSmoothing15.
397/// - compton: logical variable whether the estimation of Compton edge will be
398/// included. Possible values: kFALSE, kTRUE.
399///
400/// #### References:
401///
402/// 1. C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the
403/// quantitative analysis of PIXE spectra in geoscience applications. NIM, B34
404/// (1988), 396-402.
405///
406/// 2. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo:
407/// Background elimination methods for multidimensional gamma-ray spectra. NIM,
408/// A401 (1997) 113-132.
409///
410/// 3. D. D. Burgess, R. J. Tervo: Background estimation for gamma-ray
411/// spectroscopy. NIM 214 (1983), 431-434.
412///
413/// ### Example 1 script Background_incr.C:
414///
415/// Example of the estimation of background for number of iterations=6.
416/// Original spectrum is shown in black color, estimated background in red color.
417///
418/// Begin_Macro(source)
419/// ../../../tutorials/spectrum/Background_incr.C
420/// End_Macro
421///
422/// ### Example 2 script Background_decr.C:
423///
424/// In Example 1. one can notice that at the edges of the peaks the estimated
425/// background goes under the peaks. An alternative approach is to decrease the
426/// clipping window from a given value numberIterations to the value of one, which
427/// is presented in this example.
428///
429/// Example of the estimation of background for numberIterations=6 using
430/// decreasing clipping window algorithm. Original spectrum is shown in black
431/// color, estimated background in red color.
432///
433/// Begin_Macro(source)
434/// ../../../tutorials/spectrum/Background_decr.C
435/// End_Macro
436///
437/// ### Example 3 script Background_width.C:
438///
439/// The question is how to choose the width of the clipping window, i.e.,
440/// numberIterations parameter. The influence of this parameter on the estimated
441/// background is illustrated in Example 3.
442///
443/// Example of the influence of clipping window width on the estimated background
444/// for numberIterations=4 (red line), 6 (orange line) 8 (green line) using decreasing
445/// clipping window algorithm.
446///
447/// in general one should set this parameter so that the value
448/// 2*numberIterations+1 was greater than the widths of preserved objects (peaks).
449///
450/// Begin_Macro(source)
451/// ../../../tutorials/spectrum/Background_width.C
452/// End_Macro
453///
454/// ### Example 4 script Background_width2.C:
455///
456/// another example for very complex spectrum is given here.
457///
458/// Example of the influence of clipping window width on the estimated background
459/// for numberIterations=10 (red line), 20 (blue line), 30 (green line) and
460/// 40 (magenta line) using decreasing clipping window algorithm.
461///
462/// Begin_Macro(source)
463/// ../../../tutorials/spectrum/Background_width2.C
464/// End_Macro
465///
466/// ### Example 5 script Background_order.C:
467///
468/// Second order difference filter removes linear (quasi-linear) background and
469/// preserves symmetrical peaks. However if the shape of the background is more
470/// complex one can employ higher-order clipping filters.
471///
472/// Example of the influence of clipping filter difference order on the estimated
473/// background for fNnumberIterations=40, 2-nd order red line, 4-th order blue line,
474/// 6-th order green line and 8-th order magenta line, and using decreasing
475/// clipping window algorithm.
476///
477/// Begin_Macro(source)
478/// ../../../tutorials/spectrum/Background_order.C
479/// End_Macro
480///
481/// ### Example 6 script Background_smooth.C:
482///
483/// The estimate of the background can be influenced by noise present in the
484/// spectrum. We proposed the algorithm of the background estimate with
485/// simultaneous smoothing. In the original algorithm without smoothing, the
486/// estimated background snatches the lower spikes in the noise. Consequently,
487/// the areas of peaks are biased by this error.
488///
489/// \image html TSpectrum_Background_smooth1.jpg Principle of background estimation algorithm with simultaneous smoothing.
490///
491/// Begin_Macro(source)
492/// ../../../tutorials/spectrum/Background_smooth.C
493/// End_Macro
494///
495/// ### Example 8 script Background_compton.C:
496///
497/// Sometimes it is necessary to include also the Compton edges into the estimate of
498/// the background. This example presents the synthetic spectrum
499/// with Compton edges. The background was estimated using the 8-th order filter
500/// with the estimation of the Compton edges using decreasing
501/// clipping window algorithm (numberIterations=10) with smoothing
502/// (smoothingWindow=5).
503///
504/// Example of the estimate of the background with Compton edges (red line) for
505/// numberIterations=10, 8-th order difference filter, using decreasing clipping
506/// window algorithm and smoothing (smoothingWindow=5).
507///
508/// Begin_Macro(source)
509/// ../../../tutorials/spectrum/Background_compton.C
510/// End_Macro
511
512const char *TSpectrum::Background(Double_t *spectrum, int ssize,
513 int numberIterations,
514 int direction, int filterOrder,
515 bool smoothing,int smoothWindow,
516 bool compton)
517{
518 int i, j, w, bw, b1, b2, priz;
519 Double_t a, b, c, d, e, yb1, yb2, ai, av, men, b4, c4, d4, e4, b6, c6, d6, e6, f6, g6, b8, c8, d8, e8, f8, g8, h8, i8;
520 if (ssize <= 0)
521 return "Wrong Parameters";
522 if (numberIterations < 1)
523 return "Width of Clipping Window Must Be Positive";
524 if (ssize < 2 * numberIterations + 1)
525 return "Too Large Clipping Window";
526 if (smoothing == kTRUE && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15)
527 return "Incorrect width of smoothing window";
528 Double_t *working_space = new Double_t[2 * ssize];
529 for (i = 0; i < ssize; i++){
530 working_space[i] = spectrum[i];
531 working_space[i + ssize] = spectrum[i];
532 }
533 bw=(smoothWindow-1)/2;
534 if (direction == kBackIncreasingWindow)
535 i = 1;
536 else if(direction == kBackDecreasingWindow)
537 i = numberIterations;
538 if (filterOrder == kBackOrder2) {
539 do{
540 for (j = i; j < ssize - i; j++) {
541 if (smoothing == kFALSE){
542 a = working_space[ssize + j];
543 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
544 if (b < a)
545 a = b;
546 working_space[j] = a;
547 }
548
549 else if (smoothing == kTRUE){
550 a = working_space[ssize + j];
551 av = 0;
552 men = 0;
553 for (w = j - bw; w <= j + bw; w++){
554 if ( w >= 0 && w < ssize){
555 av += working_space[ssize + w];
556 men +=1;
557 }
558 }
559 av = av / men;
560 b = 0;
561 men = 0;
562 for (w = j - i - bw; w <= j - i + bw; w++){
563 if ( w >= 0 && w < ssize){
564 b += working_space[ssize + w];
565 men +=1;
566 }
567 }
568 b = b / men;
569 c = 0;
570 men = 0;
571 for (w = j + i - bw; w <= j + i + bw; w++){
572 if ( w >= 0 && w < ssize){
573 c += working_space[ssize + w];
574 men +=1;
575 }
576 }
577 c = c / men;
578 b = (b + c) / 2;
579 if (b < a)
580 av = b;
581 working_space[j]=av;
582 }
583 }
584 for (j = i; j < ssize - i; j++)
585 working_space[ssize + j] = working_space[j];
586 if (direction == kBackIncreasingWindow)
587 i+=1;
588 else if(direction == kBackDecreasingWindow)
589 i-=1;
590 }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
591 }
592
593 else if (filterOrder == kBackOrder4) {
594 do{
595 for (j = i; j < ssize - i; j++) {
596 if (smoothing == kFALSE){
597 a = working_space[ssize + j];
598 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
599 c = 0;
600 ai = i / 2;
601 c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
602 c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
603 c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
604 c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
605 if (b < c)
606 b = c;
607 if (b < a)
608 a = b;
609 working_space[j] = a;
610 }
611
612 else if (smoothing == kTRUE){
613 a = working_space[ssize + j];
614 av = 0;
615 men = 0;
616 for (w = j - bw; w <= j + bw; w++){
617 if ( w >= 0 && w < ssize){
618 av += working_space[ssize + w];
619 men +=1;
620 }
621 }
622 av = av / men;
623 b = 0;
624 men = 0;
625 for (w = j - i - bw; w <= j - i + bw; w++){
626 if ( w >= 0 && w < ssize){
627 b += working_space[ssize + w];
628 men +=1;
629 }
630 }
631 b = b / men;
632 c = 0;
633 men = 0;
634 for (w = j + i - bw; w <= j + i + bw; w++){
635 if ( w >= 0 && w < ssize){
636 c += working_space[ssize + w];
637 men +=1;
638 }
639 }
640 c = c / men;
641 b = (b + c) / 2;
642 ai = i / 2;
643 b4 = 0, men = 0;
644 for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
645 if (w >= 0 && w < ssize){
646 b4 += working_space[ssize + w];
647 men +=1;
648 }
649 }
650 b4 = b4 / men;
651 c4 = 0, men = 0;
652 for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
653 if (w >= 0 && w < ssize){
654 c4 += working_space[ssize + w];
655 men +=1;
656 }
657 }
658 c4 = c4 / men;
659 d4 = 0, men = 0;
660 for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
661 if (w >= 0 && w < ssize){
662 d4 += working_space[ssize + w];
663 men +=1;
664 }
665 }
666 d4 = d4 / men;
667 e4 = 0, men = 0;
668 for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
669 if (w >= 0 && w < ssize){
670 e4 += working_space[ssize + w];
671 men +=1;
672 }
673 }
674 e4 = e4 / men;
675 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
676 if (b < b4)
677 b = b4;
678 if (b < a)
679 av = b;
680 working_space[j]=av;
681 }
682 }
683 for (j = i; j < ssize - i; j++)
684 working_space[ssize + j] = working_space[j];
685 if (direction == kBackIncreasingWindow)
686 i+=1;
687 else if(direction == kBackDecreasingWindow)
688 i-=1;
689 }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
690 }
691
692 else if (filterOrder == kBackOrder6) {
693 do{
694 for (j = i; j < ssize - i; j++) {
695 if (smoothing == kFALSE){
696 a = working_space[ssize + j];
697 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
698 c = 0;
699 ai = i / 2;
700 c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
701 c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
702 c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
703 c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
704 d = 0;
705 ai = i / 3;
706 d += working_space[ssize + j - (Int_t) (3 * ai)] / 20;
707 d -= 6 * working_space[ssize + j - (Int_t) (2 * ai)] / 20;
708 d += 15 * working_space[ssize + j - (Int_t) ai] / 20;
709 d += 15 * working_space[ssize + j + (Int_t) ai] / 20;
710 d -= 6 * working_space[ssize + j + (Int_t) (2 * ai)] / 20;
711 d += working_space[ssize + j + (Int_t) (3 * ai)] / 20;
712 if (b < d)
713 b = d;
714 if (b < c)
715 b = c;
716 if (b < a)
717 a = b;
718 working_space[j] = a;
719 }
720
721 else if (smoothing == kTRUE){
722 a = working_space[ssize + j];
723 av = 0;
724 men = 0;
725 for (w = j - bw; w <= j + bw; w++){
726 if ( w >= 0 && w < ssize){
727 av += working_space[ssize + w];
728 men +=1;
729 }
730 }
731 av = av / men;
732 b = 0;
733 men = 0;
734 for (w = j - i - bw; w <= j - i + bw; w++){
735 if ( w >= 0 && w < ssize){
736 b += working_space[ssize + w];
737 men +=1;
738 }
739 }
740 b = b / men;
741 c = 0;
742 men = 0;
743 for (w = j + i - bw; w <= j + i + bw; w++){
744 if ( w >= 0 && w < ssize){
745 c += working_space[ssize + w];
746 men +=1;
747 }
748 }
749 c = c / men;
750 b = (b + c) / 2;
751 ai = i / 2;
752 b4 = 0, men = 0;
753 for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
754 if (w >= 0 && w < ssize){
755 b4 += working_space[ssize + w];
756 men +=1;
757 }
758 }
759 b4 = b4 / men;
760 c4 = 0, men = 0;
761 for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
762 if (w >= 0 && w < ssize){
763 c4 += working_space[ssize + w];
764 men +=1;
765 }
766 }
767 c4 = c4 / men;
768 d4 = 0, men = 0;
769 for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
770 if (w >= 0 && w < ssize){
771 d4 += working_space[ssize + w];
772 men +=1;
773 }
774 }
775 d4 = d4 / men;
776 e4 = 0, men = 0;
777 for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
778 if (w >= 0 && w < ssize){
779 e4 += working_space[ssize + w];
780 men +=1;
781 }
782 }
783 e4 = e4 / men;
784 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
785 ai = i / 3;
786 b6 = 0, men = 0;
787 for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
788 if (w >= 0 && w < ssize){
789 b6 += working_space[ssize + w];
790 men +=1;
791 }
792 }
793 b6 = b6 / men;
794 c6 = 0, men = 0;
795 for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
796 if (w >= 0 && w < ssize){
797 c6 += working_space[ssize + w];
798 men +=1;
799 }
800 }
801 c6 = c6 / men;
802 d6 = 0, men = 0;
803 for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
804 if (w >= 0 && w < ssize){
805 d6 += working_space[ssize + w];
806 men +=1;
807 }
808 }
809 d6 = d6 / men;
810 e6 = 0, men = 0;
811 for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
812 if (w >= 0 && w < ssize){
813 e6 += working_space[ssize + w];
814 men +=1;
815 }
816 }
817 e6 = e6 / men;
818 f6 = 0, men = 0;
819 for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
820 if (w >= 0 && w < ssize){
821 f6 += working_space[ssize + w];
822 men +=1;
823 }
824 }
825 f6 = f6 / men;
826 g6 = 0, men = 0;
827 for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
828 if (w >= 0 && w < ssize){
829 g6 += working_space[ssize + w];
830 men +=1;
831 }
832 }
833 g6 = g6 / men;
834 b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
835 if (b < b6)
836 b = b6;
837 if (b < b4)
838 b = b4;
839 if (b < a)
840 av = b;
841 working_space[j]=av;
842 }
843 }
844 for (j = i; j < ssize - i; j++)
845 working_space[ssize + j] = working_space[j];
846 if (direction == kBackIncreasingWindow)
847 i+=1;
848 else if(direction == kBackDecreasingWindow)
849 i-=1;
850 }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
851 }
852
853 else if (filterOrder == kBackOrder8) {
854 do{
855 for (j = i; j < ssize - i; j++) {
856 if (smoothing == kFALSE){
857 a = working_space[ssize + j];
858 b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
859 c = 0;
860 ai = i / 2;
861 c -= working_space[ssize + j - (Int_t) (2 * ai)] / 6;
862 c += 4 * working_space[ssize + j - (Int_t) ai] / 6;
863 c += 4 * working_space[ssize + j + (Int_t) ai] / 6;
864 c -= working_space[ssize + j + (Int_t) (2 * ai)] / 6;
865 d = 0;
866 ai = i / 3;
867 d += working_space[ssize + j - (Int_t) (3 * ai)] / 20;
868 d -= 6 * working_space[ssize + j - (Int_t) (2 * ai)] / 20;
869 d += 15 * working_space[ssize + j - (Int_t) ai] / 20;
870 d += 15 * working_space[ssize + j + (Int_t) ai] / 20;
871 d -= 6 * working_space[ssize + j + (Int_t) (2 * ai)] / 20;
872 d += working_space[ssize + j + (Int_t) (3 * ai)] / 20;
873 e = 0;
874 ai = i / 4;
875 e -= working_space[ssize + j - (Int_t) (4 * ai)] / 70;
876 e += 8 * working_space[ssize + j - (Int_t) (3 * ai)] / 70;
877 e -= 28 * working_space[ssize + j - (Int_t) (2 * ai)] / 70;
878 e += 56 * working_space[ssize + j - (Int_t) ai] / 70;
879 e += 56 * working_space[ssize + j + (Int_t) ai] / 70;
880 e -= 28 * working_space[ssize + j + (Int_t) (2 * ai)] / 70;
881 e += 8 * working_space[ssize + j + (Int_t) (3 * ai)] / 70;
882 e -= working_space[ssize + j + (Int_t) (4 * ai)] / 70;
883 if (b < e)
884 b = e;
885 if (b < d)
886 b = d;
887 if (b < c)
888 b = c;
889 if (b < a)
890 a = b;
891 working_space[j] = a;
892 }
893
894 else if (smoothing == kTRUE){
895 a = working_space[ssize + j];
896 av = 0;
897 men = 0;
898 for (w = j - bw; w <= j + bw; w++){
899 if ( w >= 0 && w < ssize){
900 av += working_space[ssize + w];
901 men +=1;
902 }
903 }
904 av = av / men;
905 b = 0;
906 men = 0;
907 for (w = j - i - bw; w <= j - i + bw; w++){
908 if ( w >= 0 && w < ssize){
909 b += working_space[ssize + w];
910 men +=1;
911 }
912 }
913 b = b / men;
914 c = 0;
915 men = 0;
916 for (w = j + i - bw; w <= j + i + bw; w++){
917 if ( w >= 0 && w < ssize){
918 c += working_space[ssize + w];
919 men +=1;
920 }
921 }
922 c = c / men;
923 b = (b + c) / 2;
924 ai = i / 2;
925 b4 = 0, men = 0;
926 for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
927 if (w >= 0 && w < ssize){
928 b4 += working_space[ssize + w];
929 men +=1;
930 }
931 }
932 b4 = b4 / men;
933 c4 = 0, men = 0;
934 for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
935 if (w >= 0 && w < ssize){
936 c4 += working_space[ssize + w];
937 men +=1;
938 }
939 }
940 c4 = c4 / men;
941 d4 = 0, men = 0;
942 for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
943 if (w >= 0 && w < ssize){
944 d4 += working_space[ssize + w];
945 men +=1;
946 }
947 }
948 d4 = d4 / men;
949 e4 = 0, men = 0;
950 for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
951 if (w >= 0 && w < ssize){
952 e4 += working_space[ssize + w];
953 men +=1;
954 }
955 }
956 e4 = e4 / men;
957 b4 = (-b4 + 4 * c4 + 4 * d4 - e4) / 6;
958 ai = i / 3;
959 b6 = 0, men = 0;
960 for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
961 if (w >= 0 && w < ssize){
962 b6 += working_space[ssize + w];
963 men +=1;
964 }
965 }
966 b6 = b6 / men;
967 c6 = 0, men = 0;
968 for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
969 if (w >= 0 && w < ssize){
970 c6 += working_space[ssize + w];
971 men +=1;
972 }
973 }
974 c6 = c6 / men;
975 d6 = 0, men = 0;
976 for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
977 if (w >= 0 && w < ssize){
978 d6 += working_space[ssize + w];
979 men +=1;
980 }
981 }
982 d6 = d6 / men;
983 e6 = 0, men = 0;
984 for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
985 if (w >= 0 && w < ssize){
986 e6 += working_space[ssize + w];
987 men +=1;
988 }
989 }
990 e6 = e6 / men;
991 f6 = 0, men = 0;
992 for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
993 if (w >= 0 && w < ssize){
994 f6 += working_space[ssize + w];
995 men +=1;
996 }
997 }
998 f6 = f6 / men;
999 g6 = 0, men = 0;
1000 for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1001 if (w >= 0 && w < ssize){
1002 g6 += working_space[ssize + w];
1003 men +=1;
1004 }
1005 }
1006 g6 = g6 / men;
1007 b6 = (b6 - 6 * c6 + 15 * d6 + 15 * e6 - 6 * f6 + g6) / 20;
1008 ai = i / 4;
1009 b8 = 0, men = 0;
1010 for (w = j - (Int_t)(4 * ai) - bw; w <= j - (Int_t)(4 * ai) + bw; w++){
1011 if (w >= 0 && w < ssize){
1012 b8 += working_space[ssize + w];
1013 men +=1;
1014 }
1015 }
1016 b8 = b8 / men;
1017 c8 = 0, men = 0;
1018 for (w = j - (Int_t)(3 * ai) - bw; w <= j - (Int_t)(3 * ai) + bw; w++){
1019 if (w >= 0 && w < ssize){
1020 c8 += working_space[ssize + w];
1021 men +=1;
1022 }
1023 }
1024 c8 = c8 / men;
1025 d8 = 0, men = 0;
1026 for (w = j - (Int_t)(2 * ai) - bw; w <= j - (Int_t)(2 * ai) + bw; w++){
1027 if (w >= 0 && w < ssize){
1028 d8 += working_space[ssize + w];
1029 men +=1;
1030 }
1031 }
1032 d8 = d8 / men;
1033 e8 = 0, men = 0;
1034 for (w = j - (Int_t)ai - bw; w <= j - (Int_t)ai + bw; w++){
1035 if (w >= 0 && w < ssize){
1036 e8 += working_space[ssize + w];
1037 men +=1;
1038 }
1039 }
1040 e8 = e8 / men;
1041 f8 = 0, men = 0;
1042 for (w = j + (Int_t)ai - bw; w <= j + (Int_t)ai + bw; w++){
1043 if (w >= 0 && w < ssize){
1044 f8 += working_space[ssize + w];
1045 men +=1;
1046 }
1047 }
1048 f8 = f8 / men;
1049 g8 = 0, men = 0;
1050 for (w = j + (Int_t)(2 * ai) - bw; w <= j + (Int_t)(2 * ai) + bw; w++){
1051 if (w >= 0 && w < ssize){
1052 g8 += working_space[ssize + w];
1053 men +=1;
1054 }
1055 }
1056 g8 = g8 / men;
1057 h8 = 0, men = 0;
1058 for (w = j + (Int_t)(3 * ai) - bw; w <= j + (Int_t)(3 * ai) + bw; w++){
1059 if (w >= 0 && w < ssize){
1060 h8 += working_space[ssize + w];
1061 men +=1;
1062 }
1063 }
1064 h8 = h8 / men;
1065 i8 = 0, men = 0;
1066 for (w = j + (Int_t)(4 * ai) - bw; w <= j + (Int_t)(4 * ai) + bw; w++){
1067 if (w >= 0 && w < ssize){
1068 i8 += working_space[ssize + w];
1069 men +=1;
1070 }
1071 }
1072 i8 = i8 / men;
1073 b8 = ( -b8 + 8 * c8 - 28 * d8 + 56 * e8 - 56 * f8 - 28 * g8 + 8 * h8 - i8)/70;
1074 if (b < b8)
1075 b = b8;
1076 if (b < b6)
1077 b = b6;
1078 if (b < b4)
1079 b = b4;
1080 if (b < a)
1081 av = b;
1082 working_space[j]=av;
1083 }
1084 }
1085 for (j = i; j < ssize - i; j++)
1086 working_space[ssize + j] = working_space[j];
1087 if (direction == kBackIncreasingWindow)
1088 i += 1;
1089 else if(direction == kBackDecreasingWindow)
1090 i -= 1;
1091 }while((direction == kBackIncreasingWindow && i <= numberIterations) || (direction == kBackDecreasingWindow && i >= 1));
1092 }
1093
1094 if (compton == kTRUE) {
1095 for (i = 0, b2 = 0; i < ssize; i++){
1096 b1 = b2;
1097 a = working_space[i], b = spectrum[i];
1098 j = i;
1099 if (TMath::Abs(a - b) >= 1) {
1100 b1 = i - 1;
1101 if (b1 < 0)
1102 b1 = 0;
1103 yb1 = working_space[b1];
1104 for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){
1105 a = working_space[b2], b = spectrum[b2];
1106 c = c + b - yb1;
1107 if (TMath::Abs(a - b) < 1) {
1108 priz = 1;
1109 yb2 = b;
1110 }
1111 }
1112 if (b2 == ssize)
1113 b2 -= 1;
1114 yb2 = working_space[b2];
1115 if (yb1 <= yb2){
1116 for (j = b1, c = 0; j <= b2; j++){
1117 b = spectrum[j];
1118 c = c + b - yb1;
1119 }
1120 if (c > 1){
1121 c = (yb2 - yb1) / c;
1122 for (j = b1, d = 0; j <= b2 && j < ssize; j++){
1123 b = spectrum[j];
1124 d = d + b - yb1;
1125 a = c * d + yb1;
1126 working_space[ssize + j] = a;
1127 }
1128 }
1129 }
1130
1131 else{
1132 for (j = b2, c = 0; j >= b1; j--){
1133 b = spectrum[j];
1134 c = c + b - yb2;
1135 }
1136 if (c > 1){
1137 c = (yb1 - yb2) / c;
1138 for (j = b2, d = 0;j >= b1 && j >= 0; j--){
1139 b = spectrum[j];
1140 d = d + b - yb2;
1141 a = c * d + yb2;
1142 working_space[ssize + j] = a;
1143 }
1144 }
1145 }
1146 i=b2;
1147 }
1148 }
1149 }
1150
1151 for (j = 0; j < ssize; j++)
1152 spectrum[j] = working_space[ssize + j];
1153 delete[]working_space;
1154 return 0;
1155}
1156
1157////////////////////////////////////////////////////////////////////////////////
1158/// One-dimensional markov spectrum smoothing function
1159///
1160/// This function calculates smoothed spectrum from source spectrum based on
1161/// Markov chain method. The result is placed in the array pointed by source
1162/// pointer. On successful completion it returns 0. On error it returns pointer
1163/// to the string describing error.
1164///
1165/// #### Parameters:
1166///
1167/// - source: pointer to the array of source spectrum
1168/// - ssize: length of source array
1169/// - averWindow: width of averaging smoothing window
1170///
1171/// The goal of this function is the suppression of the statistical fluctuations.
1172/// The algorithm is based on discrete Markov chain, which has very simple
1173/// invariant distribution:
1174///
1175/// \f[
1176/// U_2 = \frac{p_{1,2}}{p_{2,1}}U_1, U_3 = \frac{p_{2,3}}{p_{3,2}}U_2U_1, ... , U_n = \frac{p_{n-1,n}}{p_{n,n-1}}U_{n-1}...U_2U_1
1177/// \f]
1178/// \f$ U_1\f$ being defined from the normalization condition
1179/// \f$ \sum_{i=1}^{n} U_i=1\f$. \f$ n \f$ is the length of the smoothed spectrum and
1180/// \f[
1181/// p_{i,i\pm 1} = A_i\sum_{k=1}^{m} exp\left[ \frac{y(i\pm k)-y(i)}{y(i\pm k)+y(i)}\right]
1182/// \f]
1183///
1184/// #### Reference:
1185///
1186/// 1. Z.K. Silagadze, A new algorithm for automatic photopeak searches.
1187/// NIM A 376 (1996), 451.
1188///
1189/// ### Example 14 - script Smoothing.C
1190///
1191/// Begin_Macro(source)
1192/// ../../../tutorials/spectrum/Smoothing.C
1193/// End_Macro
1194
1195const char* TSpectrum::SmoothMarkov(Double_t *source, int ssize, int averWindow)
1196{
1197 int xmin, xmax, i, l;
1198 Double_t a, b, maxch;
1199 Double_t nom, nip, nim, sp, sm, area = 0;
1200 if(averWindow <= 0)
1201 return "Averaging Window must be positive";
1202 Double_t *working_space = new Double_t[ssize];
1203 xmin = 0,xmax = ssize - 1;
1204 for(i = 0, maxch = 0; i < ssize; i++){
1205 working_space[i]=0;
1206 if(maxch < source[i])
1207 maxch = source[i];
1208
1209 area += source[i];
1210 }
1211 if(maxch == 0) {
1212 delete [] working_space;
1213 return 0 ;
1214 }
1215
1216 nom = 1;
1217 working_space[xmin] = 1;
1218 for(i = xmin; i < xmax; i++){
1219 nip = source[i] / maxch;
1220 nim = source[i + 1] / maxch;
1221 sp = 0,sm = 0;
1222 for(l = 1; l <= averWindow; l++){
1223 if((i + l) > xmax)
1224 a = source[xmax] / maxch;
1225
1226 else
1227 a = source[i + l] / maxch;
1228 b = a - nip;
1229 if(a + nip <= 0)
1230 a = 1;
1231
1232 else
1233 a = TMath::Sqrt(a + nip);
1234 b = b / a;
1235 b = TMath::Exp(b);
1236 sp = sp + b;
1237 if((i - l + 1) < xmin)
1238 a = source[xmin] / maxch;
1239
1240 else
1241 a = source[i - l + 1] / maxch;
1242 b = a - nim;
1243 if(a + nim <= 0)
1244 a = 1;
1245 else
1246 a = TMath::Sqrt(a + nim);
1247 b = b / a;
1248 b = TMath::Exp(b);
1249 sm = sm + b;
1250 }
1251 a = sp / sm;
1252 a = working_space[i + 1] = working_space[i] * a;
1253 nom = nom + a;
1254 }
1255 for(i = xmin; i <= xmax; i++){
1256 working_space[i] = working_space[i] / nom;
1257 }
1258 for(i = 0; i < ssize; i++)
1259 source[i] = working_space[i] * area;
1260 delete[]working_space;
1261 return 0;
1262}
1263
1264////////////////////////////////////////////////////////////////////////////////
1265/// One-dimensional deconvolution function
1266///
1267/// This function calculates deconvolution from source spectrum according to
1268/// response spectrum using Gold deconvolution algorithm. The result is placed
1269/// in the vector pointed by source pointer. On successful completion it
1270/// returns 0. On error it returns pointer to the string describing error. If
1271/// desired after every numberIterations one can apply boosting operation
1272/// (exponential function with exponent given by boost coefficient) and repeat
1273/// it numberRepetitions times.
1274///
1275/// #### Parameters:
1276///
1277/// - source: pointer to the vector of source spectrum
1278/// - response: pointer to the vector of response spectrum
1279/// - ssize: length of source and response spectra
1280/// - numberIterations, for details we refer to the reference given below
1281/// - numberRepetitions, for repeated boosted deconvolution
1282/// - boost, boosting coefficient
1283///
1284/// The goal of this function is the improvement of the resolution in spectra,
1285/// decomposition of multiplets. The mathematical formulation of
1286/// the convolution system is:
1287///
1288/// \f[
1289/// y(i) = \sum_{k=0}^{N-1} h(i-k)x(k), i=0,1,2,...,N-1
1290/// \f]
1291///
1292/// where h(i) is the impulse response function, x, y are input and output
1293/// vectors, respectively, N is the length of x and h vectors. In matrix form
1294/// we have:
1295/**
1296 \f[
1297 \begin{bmatrix}
1298 y(0) \\
1299 y(1) \\
1300 \dots \\
1301 y(2N-2) \\
1302 y(2N-1)
1303 \end{bmatrix}
1304 =
1305 \begin{bmatrix}
1306 h(0) & 0 & 0 & \dots & 0 \\
1307 h(1) & h(0) & 0 & \dots & \dots \\
1308 \dots & h(1) & h(0) & \dots & \dots \\
1309 \dots & \dots & h(1) & \dots & \dots \\
1310 \dots & \dots & \dots & \dots & \dots \\
1311 h(N-1) & \dots & \dots &\dots & 0 \\
1312 0 & h(N-1) & \dots & \dots & h(0) \\
1313 0 & 0 & h(N-1) & \dots & h(1) \\
1314 \dots & \dots & \dots & \dots & \dots \\
1315 0 & 0 & 0 & \dots & h(N-1)
1316 \end{bmatrix}
1317 \begin{bmatrix}
1318 x(0) \\
1319 x(1) \\
1320 \dots \\
1321 x(N-1)
1322 \end{bmatrix}
1323 \f]
1324*/
1325/// Let us assume that we know the response and the output vector (spectrum) of
1326/// the above given system. The deconvolution represents solution of the
1327/// overdetermined system of linear equations, i.e., the calculation of the
1328/// vector x. From numerical stability point of view the operation of
1329/// deconvolution is extremely critical (ill-posed problem) as well as time
1330/// consuming operation. The Gold deconvolution algorithm proves to work very
1331/// well, other methods (Fourier, VanCittert etc) oscillate. It is suitable to
1332/// process positive definite data (e.g. histograms).
1333///
1334/// #### Gold deconvolution algorithm:
1335/**
1336 \f[
1337 y = Hx \\
1338 H^{T}y = H^{T}Hx \\
1339 y^{'} = H^{'}x \\
1340 x_{i}^{(k+1)} = \frac{y_{i}^{'}}{\sum_{m=0}^{N-1}H_{im}^{'}x_{m}{(k)}}x_{i}{(k)}, i=0,1,2,...,N-1 \\
1341 k = 1,2,3,...,L
1342 x^{0} = [1,1, ..., 1]^T
1343 \f]
1344*/
1345/// Where L is given number of iterations (numberIterations parameter).
1346///
1347/// #### Boosted deconvolution:
1348///
1349/// 1. Set the initial solution:
1350/// \f$ x^{(0)} = [1,1,...,1]^{T} \f$
1351/// 2. Set required number of repetitions R and iterations L.
1352/// 3. Set r = 1.
1353/// 4. Using Gold deconvolution algorithm for k=1,2,...,L find
1354/// \f$ x^{(L)} \f$
1355/// 5. If r = R stop calculation, else
1356///
1357/// 1. Apply boosting operation, i.e., set
1358/// \f$ x^{(0)}(i) = [x^{(L)}(i)]^{p} \f$
1359/// i=0,1,...N-1 and p is boosting coefficient >0.
1360/// 2. r = r + 1
1361/// 3. continue in 4.
1362///
1363/// #### References:
1364///
1365/// 1. Gold R., ANL-6984, Argonne National Laboratories, Argonne Ill, 1964.
1366/// 2. Coote G.E., Iterative smoothing and deconvolution of one- and two-dimensional
1367/// elemental distribution data, NIM B 130 (1997) 118.
1368/// 3. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky,
1369/// I. Turzo: Efficient one- and two-dimensional Gold deconvolution and
1370/// its application to gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
1371/// 4. Morhac; M., Matouoek V., Kliman J., Efficient algorithm of multidimensional
1372/// deconvolution and its application to nuclear data processing, Digital Signal
1373/// Processing 13 (2003) 144.
1374///
1375/// ### Example 8 - script Deconvolution.C :
1376///
1377/// response function (usually peak) should be shifted left to the first
1378/// non-zero channel (bin).
1379///
1380/// \image html TSpectrum_Deconvolution2.jpg Principle how the response matrix is composed inside of the Deconvolution function.
1381///
1382/// Begin_Macro(source)
1383/// ../../../tutorials/spectrum/Deconvolution.C
1384/// End_Macro
1385///
1386/// ### Examples of Gold deconvolution method:
1387///
1388/// First let us study the influence of the number of iterations on the
1389/// deconvolved spectrum (Fig. 12).
1390///
1391/// \image html TSpectrum_Deconvolution_wide1.jpg Fig. 12 Study of Gold deconvolution algorithm.The original source spectrum is drawn with black color, spectrum after 100 iterations with red color, spectrum after 1000 iterations with blue color, spectrum after 10000 iterations with green color and spectrum after 100000 iterations with magenta color.
1392///
1393/// For relatively narrow peaks in the above given example the Gold
1394/// deconvolution method is able to decompose overlapping peaks practically to
1395/// delta - functions. In the next example we have chosen a synthetic data
1396/// (spectrum, 256 channels) consisting of 5 very closely positioned, relatively
1397/// wide peaks (sigma =5), with added noise (Fig. 13). Thin lines represent
1398/// pure Gaussians (see Table 1); thick line is a resulting spectrum with
1399/// additive noise (10% of the amplitude of small peaks).
1400///
1401/// \image html TSpectrum_Deconvolution_wide2.jpg Fig. 13 Testing example of synthetic spectrum composed of 5 Gaussians with added noise.
1402///
1403/// | Peak # | Position | Height | Area |
1404/// |----------|----------|--------|--------|
1405/// | 1 | 50 | 500 | 10159 |
1406/// | 2 | 70 | 3000 | 60957 |
1407/// | 3 | 80 | 1000 | 20319 |
1408/// | 4 | 100 | 5000 | 101596 |
1409/// | 5 | 110 | 500 | 10159 |
1410///
1411/// Table 1 Positions, heights and areas of peaks in the spectrum shown in Fig. 13.
1412///
1413/// In ideal case, we should obtain the result given in Fig. 14. The areas of
1414/// the Gaussian components of the spectrum are concentrated completely to
1415/// delta-functions. When solving the overdetermined system of linear equations
1416/// with data from Fig. 13 in the sense of minimum least squares criterion
1417/// without any regularisation we obtain the result with large oscillations
1418/// (Fig. 15). From mathematical point of view, it is the optimal solution in
1419/// the unconstrained space of independent variables. From physical point of
1420/// view we are interested only in a meaningful solution. Therefore, we have to
1421/// employ regularisation techniques (e.g. Gold deconvolution) and/or to
1422/// confine the space of allowed solutions to subspace of positive solutions.
1423///
1424/// \image html TSpectrum_Deconvolution_wide3.jpg Fig. 14 The same spectrum like in Fig. 13, outlined bars show the contents of present components (peaks).
1425/// \image html TSpectrum_Deconvolution_wide4.jpg Fig. 15 Least squares solution of the system of linear equations without regularisation.
1426///
1427/// ### Example 9 - script Deconvolution_wide.C
1428///
1429/// When we employ Gold deconvolution algorithm we obtain the result given in
1430/// Fig. 16. One can observe that the resulting spectrum is smooth. On the
1431/// other hand the method is not able to decompose completely the peaks in the
1432/// spectrum.
1433///
1434/// Example of Gold deconvolution for closely positioned wide peaks. The original
1435/// source spectrum is drawn with black color, the spectrum after the deconvolution
1436/// (10000 iterations) with red color.
1437///
1438/// Begin_Macro(source)
1439/// ../../../tutorials/spectrum/Deconvolution_wide.C
1440/// End_Macro
1441///
1442/// ### Example 10 - script Deconvolution_wide_boost.C :
1443///
1444/// Further let us employ boosting operation into deconvolution (Fig. 17).
1445///
1446/// The original source spectrum is drawn with black color, the spectrum after
1447/// the deconvolution with red color. Number of iterations = 200, number of
1448/// repetitions = 50 and boosting coefficient = 1.2.
1449///
1450/// One can observe that peaks are decomposed practically to delta functions.
1451/// Number of peaks is correct, positions of big peaks as well as their areas
1452/// are relatively well estimated. However there is a considerable error in
1453/// the estimation of the position of small right hand peak.
1454///
1455/// Begin_Macro(source)
1456/// ../../../tutorials/spectrum/Deconvolution_wide_boost.C
1457/// End_Macro
1458
1459const char *TSpectrum::Deconvolution(Double_t *source, const Double_t *response,
1460 int ssize, int numberIterations,
1461 int numberRepetitions, Double_t boost )
1462{
1463 if (ssize <= 0)
1464 return "Wrong Parameters";
1465
1466 if (numberRepetitions <= 0)
1467 return "Wrong Parameters";
1468
1469 // working_space-pointer to the working vector
1470 // (its size must be 4*ssize of source spectrum)
1471 Double_t *working_space = new Double_t[4 * ssize];
1472 int i, j, k, lindex, posit, lh_gold, l, repet;
1473 Double_t lda, ldb, ldc, area, maximum;
1474 area = 0;
1475 lh_gold = -1;
1476 posit = 0;
1477 maximum = 0;
1478
1479//read response vector
1480 for (i = 0; i < ssize; i++) {
1481 lda = response[i];
1482 if (lda != 0)
1483 lh_gold = i + 1;
1484 working_space[i] = lda;
1485 area += lda;
1486 if (lda > maximum) {
1487 maximum = lda;
1488 posit = i;
1489 }
1490 }
1491 if (lh_gold == -1) {
1492 delete [] working_space;
1493 return "ZERO RESPONSE VECTOR";
1494 }
1495
1496//read source vector
1497 for (i = 0; i < ssize; i++)
1498 working_space[2 * ssize + i] = source[i];
1499
1500// create matrix at*a and vector at*y
1501 for (i = 0; i < ssize; i++){
1502 lda = 0;
1503 for (j = 0; j < ssize; j++){
1504 ldb = working_space[j];
1505 k = i + j;
1506 if (k < ssize){
1507 ldc = working_space[k];
1508 lda = lda + ldb * ldc;
1509 }
1510 }
1511 working_space[ssize + i] = lda;
1512 lda = 0;
1513 for (k = 0; k < ssize; k++){
1514 l = k - i;
1515 if (l >= 0){
1516 ldb = working_space[l];
1517 ldc = working_space[2 * ssize + k];
1518 lda = lda + ldb * ldc;
1519 }
1520 }
1521 working_space[3 * ssize + i]=lda;
1522 }
1523
1524// move vector at*y
1525 for (i = 0; i < ssize; i++){
1526 working_space[2 * ssize + i] = working_space[3 * ssize + i];
1527 }
1528
1529//initialization of resulting vector
1530 for (i = 0; i < ssize; i++)
1531 working_space[i] = 1;
1532
1533 //**START OF ITERATIONS**
1534 for (repet = 0; repet < numberRepetitions; repet++) {
1535 if (repet != 0) {
1536 for (i = 0; i < ssize; i++)
1537 working_space[i] = TMath::Power(working_space[i], boost);
1538 }
1539 for (lindex = 0; lindex < numberIterations; lindex++) {
1540 for (i = 0; i < ssize; i++) {
1541 if (working_space[2 * ssize + i] > 0.000001
1542 && working_space[i] > 0.000001) {
1543 lda = 0;
1544 for (j = 0; j < lh_gold; j++) {
1545 ldb = working_space[j + ssize];
1546 if (j != 0){
1547 k = i + j;
1548 ldc = 0;
1549 if (k < ssize)
1550 ldc = working_space[k];
1551 k = i - j;
1552 if (k >= 0)
1553 ldc += working_space[k];
1554 }
1555
1556 else
1557 ldc = working_space[i];
1558 lda = lda + ldb * ldc;
1559 }
1560 ldb = working_space[2 * ssize + i];
1561 if (lda != 0)
1562 lda = ldb / lda;
1563
1564 else
1565 lda = 0;
1566 ldb = working_space[i];
1567 lda = lda * ldb;
1568 working_space[3 * ssize + i] = lda;
1569 }
1570 }
1571 for (i = 0; i < ssize; i++)
1572 working_space[i] = working_space[3 * ssize + i];
1573 }
1574 }
1575
1576//shift resulting spectrum
1577 for (i = 0; i < ssize; i++) {
1578 lda = working_space[i];
1579 j = i + posit;
1580 j = j % ssize;
1581 working_space[ssize + j] = lda;
1582 }
1583
1584//write back resulting spectrum
1585 for (i = 0; i < ssize; i++)
1586 source[i] = area * working_space[ssize + i];
1587 delete[]working_space;
1588 return 0;
1589}
1590
1591
1592////////////////////////////////////////////////////////////////////////////////
1593/// One-dimensional deconvolution function.
1594///
1595/// This function calculates deconvolution from source spectrum according to
1596/// response spectrum using Richardson-Lucy deconvolution algorithm. The result
1597/// is placed in the vector pointed by source pointer. On successful completion
1598/// it returns 0. On error it returns pointer to the string describing error.
1599/// If desired after every numberIterations one can apply boosting operation
1600/// (exponential function with exponent given by boost coefficient) and repeat
1601/// it numberRepetitions times (see Gold deconvolution).
1602///
1603/// #### Parameters:
1604///
1605/// - source: pointer to the vector of source spectrum
1606/// - response: pointer to the vector of response spectrum
1607/// - ssize: length of source and response spectra
1608/// - numberIterations, for details we refer to the reference given above
1609/// - numberRepetitions, for repeated boosted deconvolution
1610/// - boost, boosting coefficient
1611///
1612/// ### Richardson-Lucy deconvolution algorithm:
1613///
1614/// For discrete systems it has the form:
1615/**
1616 \f[
1617 x^{(n)}(i) = x^{(n-1)}(i) \sum_{j=0}^{N-1}h(i,j)\frac{y(j)}{\sum_{k=0}^{M-1}h(j,k)x^{(n-1)}(k)} \\
1618 i \in \left<0,M-1\right>
1619 \f]
1620*/
1621/// for positive input data and response matrix this iterative method forces
1622/// the deconvoluted spectra to be non-negative. The Richardson-Lucy
1623/// iteration converges to the maximum likelihood solution for Poisson statistics
1624/// in the data.
1625///
1626/// #### References:
1627///
1628/// 1. Abreu M.C. et al., A four-dimensional deconvolution method to correct NA38
1629/// experimental data, NIM A 405 (1998) 139.
1630/// 2. Lucy L.B., A.J. 79 (1974) 745.
1631/// 3. Richardson W.H., J. Opt. Soc. Am. 62 (1972) 55.
1632///
1633/// ### Examples of Richardson-Lucy deconvolution method:
1634///
1635/// ### Example 11 - script DeconvolutionRL_wide.C :
1636///
1637/// When we employ Richardson-Lucy deconvolution algorithm to our data from
1638/// Fig. 13 we obtain the following result. One can observe improvements
1639/// as compared to the result achieved by Gold deconvolution. Nevertheless it is
1640/// unable to decompose the multiplet.
1641///
1642/// Example of Richardson-Lucy deconvolution for closely positioned wide peaks.
1643/// The original source spectrum is drawn with black color, the spectrum after
1644/// the deconvolution (10000 iterations) with red color.
1645///
1646/// Begin_Macro(source)
1647/// ../../../tutorials/spectrum/DeconvolutionRL_wide.C
1648/// End_Macro
1649///
1650/// ### Example 12 - script DeconvolutionRL_wide_boost.C :
1651///
1652/// Further let us employ boosting operation into deconvolution.
1653///
1654/// The original source spectrum is drawn with black color, the spectrum after
1655/// the deconvolution with red color. Number of iterations = 200, number of
1656/// repetitions = 50 and boosting coefficient = 1.2.
1657///
1658/// One can observe improvements in the estimation of peak positions as compared
1659/// to the results achieved by Gold deconvolution.
1660///
1661/// Begin_Macro(source)
1662/// ../../../tutorials/spectrum/DeconvolutionRL_wide_boost.C
1663/// End_Macro
1664
1665const char *TSpectrum::DeconvolutionRL(Double_t *source, const Double_t *response,
1666 int ssize, int numberIterations,
1667 int numberRepetitions, Double_t boost )
1668{
1669 if (ssize <= 0)
1670 return "Wrong Parameters";
1671
1672 if (numberRepetitions <= 0)
1673 return "Wrong Parameters";
1674
1675 // working_space-pointer to the working vector
1676 // (its size must be 4*ssize of source spectrum)
1677 Double_t *working_space = new Double_t[4 * ssize];
1678 int i, j, k, lindex, posit, lh_gold, repet, kmin, kmax;
1679 Double_t lda, ldb, ldc, maximum;
1680 lh_gold = -1;
1681 posit = 0;
1682 maximum = 0;
1683
1684//read response vector
1685 for (i = 0; i < ssize; i++) {
1686 lda = response[i];
1687 if (lda != 0)
1688 lh_gold = i + 1;
1689 working_space[ssize + i] = lda;
1690 if (lda > maximum) {
1691 maximum = lda;
1692 posit = i;
1693 }
1694 }
1695 if (lh_gold == -1) {
1696 delete [] working_space;
1697 return "ZERO RESPONSE VECTOR";
1698 }
1699
1700//read source vector
1701 for (i = 0; i < ssize; i++)
1702 working_space[2 * ssize + i] = source[i];
1703
1704//initialization of resulting vector
1705 for (i = 0; i < ssize; i++){
1706 if (i <= ssize - lh_gold)
1707 working_space[i] = 1;
1708
1709 else
1710 working_space[i] = 0;
1711
1712 }
1713 //**START OF ITERATIONS**
1714 for (repet = 0; repet < numberRepetitions; repet++) {
1715 if (repet != 0) {
1716 for (i = 0; i < ssize; i++)
1717 working_space[i] = TMath::Power(working_space[i], boost);
1718 }
1719 for (lindex = 0; lindex < numberIterations; lindex++) {
1720 for (i = 0; i <= ssize - lh_gold; i++){
1721 lda = 0;
1722 if (working_space[i] > 0){//x[i]
1723 for (j = i; j < i + lh_gold; j++){
1724 ldb = working_space[2 * ssize + j];//y[j]
1725 if (j < ssize){
1726 if (ldb > 0){//y[j]
1727 kmax = j;
1728 if (kmax > lh_gold - 1)
1729 kmax = lh_gold - 1;
1730 kmin = j + lh_gold - ssize;
1731 if (kmin < 0)
1732 kmin = 0;
1733 ldc = 0;
1734 for (k = kmax; k >= kmin; k--){
1735 ldc += working_space[ssize + k] * working_space[j - k];//h[k]*x[j-k]
1736 }
1737 if (ldc > 0)
1738 ldb = ldb / ldc;
1739
1740 else
1741 ldb = 0;
1742 }
1743 ldb = ldb * working_space[ssize + j - i];//y[j]*h[j-i]/suma(h[j][k]x[k])
1744 }
1745 lda += ldb;
1746 }
1747 lda = lda * working_space[i];
1748 }
1749 working_space[3 * ssize + i] = lda;
1750 }
1751 for (i = 0; i < ssize; i++)
1752 working_space[i] = working_space[3 * ssize + i];
1753 }
1754 }
1755
1756//shift resulting spectrum
1757 for (i = 0; i < ssize; i++) {
1758 lda = working_space[i];
1759 j = i + posit;
1760 j = j % ssize;
1761 working_space[ssize + j] = lda;
1762 }
1763
1764//write back resulting spectrum
1765 for (i = 0; i < ssize; i++)
1766 source[i] = working_space[ssize + i];
1767 delete[]working_space;
1768 return 0;
1769}
1770
1771
1772////////////////////////////////////////////////////////////////////////////////
1773/// One-dimensional unfolding function
1774///
1775/// This function unfolds source spectrum according to response matrix columns.
1776/// The result is placed in the vector pointed by source pointer.
1777/// The coefficients of the resulting vector represent contents of the columns
1778/// (weights) in the input vector. On successful completion it returns 0. On
1779/// error it returns pointer to the string describing error. If desired after
1780/// every numberIterations one can apply boosting operation (exponential
1781/// function with exponent given by boost coefficient) and repeat it
1782/// numberRepetitions times. For details we refer to [1].
1783///
1784/// #### Parameters:
1785///
1786/// - source: pointer to the vector of source spectrum
1787/// - respMatrix: pointer to the matrix of response spectra
1788/// - ssizex: length of source spectrum and # of rows of the response
1789/// matrix. ssizex must be >= ssizey.
1790/// - ssizey: length of destination coefficients and # of columns of the response
1791/// matrix.
1792/// - numberIterations: number of iterations
1793/// - numberRepetitions: number of repetitions for boosted deconvolution.
1794/// It must be greater or equal to one.
1795/// - boost: boosting coefficient, applies only if numberRepetitions is
1796/// greater than one.
1797///
1798/// ### Unfolding:
1799///
1800/// The goal is the decomposition of spectrum to a given set of component spectra.
1801///
1802/// The mathematical formulation of the discrete linear system is:
1803///
1804/// \f[
1805/// y(i) = \sum_{k=0}^{N_y-1} h(i,k)x(k), i = 0,1,2,...,N_x-1
1806/// \f]
1807/**
1808 \f[
1809 \begin{bmatrix}
1810 y(0) \\
1811 y(1) \\
1812 \dots \\
1813 y(N_x-1)
1814 \end{bmatrix}
1815 =
1816 \begin{bmatrix}
1817 h(0,0) & h(0,1) & \dots & h(0,N_y-1) \\
1818 h(1,0) & h(1,1) & \dots & h(1,N_y-1) \\
1819 \dots \\
1820 h(N_x-1,0) & h(N_x-1,1) & \dots & h(N_x-1,N_y-1)
1821 \end{bmatrix}
1822 \begin{bmatrix}
1823 x(0) \\
1824 x(1) \\
1825 \dots \\
1826 x(N_y-1)
1827 \end{bmatrix}
1828 \f]
1829*/
1830///
1831/// #### References:
1832///
1833/// 1. Jandel M., Morhac; M., Kliman J., Krupa L., Matouoek
1834/// V., Hamilton J. H., Ramaya A. V.:
1835/// Decomposition of continuum gamma-ray spectra using synthesised response matrix.
1836/// NIM A 516 (2004), 172-183.
1837///
1838/// ### Example of unfolding:
1839///
1840/// ### Example 13 - script Unfolding.C:
1841///
1842/// \image html TSpectrum_Unfolding3.gif Fig. 20 Response matrix composed of neutron spectra of pure chemical elements.
1843/// \image html TSpectrum_Unfolding2.jpg Fig. 21 Source neutron spectrum to be decomposed
1844/// \image html TSpectrum_Unfolding3.jpg Fig. 22 Spectrum after decomposition, contains 10 coefficients, which correspond to contents of chemical components (dominant 8-th and 10-th components, i.e. O, Si)
1845///
1846/// #### Script:
1847///
1848/// ~~~ {.cpp}
1849/// // Example to illustrate unfolding function (class TSpectrum).
1850/// // To execute this example, do
1851/// // root > .x Unfolding.C
1852///
1853/// void Unfolding() {
1854/// Int_t i, j;
1855/// Int_t nbinsx = 2048;
1856/// Int_t nbinsy = 10;
1857/// double xmin = 0;
1858/// double xmax = nbinsx;
1859/// double ymin = 0;
1860/// double ymax = nbinsy;
1861/// double *source = new double[nbinsx];
1862/// double **response = new double *[nbinsy];
1863/// for (i=0;i<nbinsy;i++) response[i] = new double[nbinsx];
1864/// TH1F *h = new TH1F("h","",nbinsx,xmin,xmax);
1865/// TH1F *d = new TH1F("d","Decomposition - unfolding",nbinsx,xmin,xmax);
1866/// TH2F *decon_unf_resp = new TH2F("decon_unf_resp","Root File",nbinsy,ymin,ymax,nbinsx,xmin,xmax);
1867/// TFile *f = new TFile("TSpectrum.root");
1868/// h = (TH1F*) f->Get("decon_unf_in;1");
1869/// TFile *fr = new TFile("TSpectrum.root");
1870/// decon_unf_resp = (TH2F*) fr->Get("decon_unf_resp;1");
1871/// for (i = 0; i < nbinsx; i++) source[i] = h->GetBinContent(i + 1);
1872/// for (i = 0; i < nbinsy; i++){
1873/// for (j = 0; j< nbinsx; j++){
1874/// response[i][j] = decon_unf_resp->GetBinContent(i + 1, j + 1);
1875/// }
1876/// }
1877/// TCanvas *Decon1 = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("Decon1");
1878/// if (!Decon1) Decon1 = new TCanvas("Decon1","Decon1",10,10,1000,700);
1879/// h->Draw("L");
1880/// TSpectrum *s = new TSpectrum();
1881/// s->Unfolding(source,(const double**)response,nbinsx,nbinsy,1000,1,1);
1882/// for (i = 0; i < nbinsy; i++) d->SetBinContent(i + 1,source[i]);
1883/// d->SetLineColor(kRed);
1884/// d->SetAxisRange(0,nbinsy);
1885/// d->Draw("");
1886/// }
1887/// ~~~
1888
1889const char *TSpectrum::Unfolding(Double_t *source,
1890 const Double_t **respMatrix,
1891 int ssizex, int ssizey,
1892 int numberIterations,
1893 int numberRepetitions, Double_t boost)
1894{
1895 int i, j, k, lindex, lhx = 0, repet;
1896 Double_t lda, ldb, ldc, area;
1897 if (ssizex <= 0 || ssizey <= 0)
1898 return "Wrong Parameters";
1899 if (ssizex < ssizey)
1900 return "Sizex must be greater than sizey)";
1901 if (numberIterations <= 0)
1902 return "Number of iterations must be positive";
1903 Double_t *working_space =
1904 new Double_t[ssizex * ssizey + 2 * ssizey * ssizey + 4 * ssizex];
1905
1906/*read response matrix*/
1907 for (j = 0; j < ssizey && lhx != -1; j++) {
1908 area = 0;
1909 lhx = -1;
1910 for (i = 0; i < ssizex; i++) {
1911 lda = respMatrix[j][i];
1912 if (lda != 0) {
1913 lhx = i + 1;
1914 }
1915 working_space[j * ssizex + i] = lda;
1916 area = area + lda;
1917 }
1918 if (lhx != -1) {
1919 for (i = 0; i < ssizex; i++)
1920 working_space[j * ssizex + i] /= area;
1921 }
1922 }
1923 if (lhx == -1) {
1924 delete [] working_space;
1925 return ("ZERO COLUMN IN RESPONSE MATRIX");
1926 }
1927
1928/*read source vector*/
1929 for (i = 0; i < ssizex; i++)
1930 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
1931 source[i];
1932
1933/*create matrix at*a + at*y */
1934 for (i = 0; i < ssizey; i++) {
1935 for (j = 0; j < ssizey; j++) {
1936 lda = 0;
1937 for (k = 0; k < ssizex; k++) {
1938 ldb = working_space[ssizex * i + k];
1939 ldc = working_space[ssizex * j + k];
1940 lda = lda + ldb * ldc;
1941 }
1942 working_space[ssizex * ssizey + ssizey * i + j] = lda;
1943 }
1944 lda = 0;
1945 for (k = 0; k < ssizex; k++) {
1946 ldb = working_space[ssizex * i + k];
1947 ldc =
1948 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
1949 k];
1950 lda = lda + ldb * ldc;
1951 }
1952 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
1953 lda;
1954 }
1955
1956/*move vector at*y*/
1957 for (i = 0; i < ssizey; i++)
1958 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
1959 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
1960
1961/*create matrix at*a*at*a + vector at*a*at*y */
1962 for (i = 0; i < ssizey; i++) {
1963 for (j = 0; j < ssizey; j++) {
1964 lda = 0;
1965 for (k = 0; k < ssizey; k++) {
1966 ldb = working_space[ssizex * ssizey + ssizey * i + k];
1967 ldc = working_space[ssizex * ssizey + ssizey * j + k];
1968 lda = lda + ldb * ldc;
1969 }
1970 working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j] =
1971 lda;
1972 }
1973 lda = 0;
1974 for (k = 0; k < ssizey; k++) {
1975 ldb = working_space[ssizex * ssizey + ssizey * i + k];
1976 ldc =
1977 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex +
1978 k];
1979 lda = lda + ldb * ldc;
1980 }
1981 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] =
1982 lda;
1983 }
1984
1985/*move at*a*at*y*/
1986 for (i = 0; i < ssizey; i++)
1987 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i] =
1988 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
1989
1990/*initialization in resulting vector */
1991 for (i = 0; i < ssizey; i++)
1992 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = 1;
1993
1994 /***START OF ITERATIONS***/
1995 for (repet = 0; repet < numberRepetitions; repet++) {
1996 if (repet != 0) {
1997 for (i = 0; i < ssizey; i++)
1998 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] = TMath::Power(working_space[ssizex * ssizey + 2 * ssizey * ssizey + i], boost);
1999 }
2000 for (lindex = 0; lindex < numberIterations; lindex++) {
2001 for (i = 0; i < ssizey; i++) {
2002 lda = 0;
2003 for (j = 0; j < ssizey; j++) {
2004 ldb =
2005 working_space[ssizex * ssizey + ssizey * ssizey + ssizey * i + j];
2006 ldc = working_space[ssizex * ssizey + 2 * ssizey * ssizey + j];
2007 lda = lda + ldb * ldc;
2008 }
2009 ldb =
2010 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 2 * ssizex + i];
2011 if (lda != 0) {
2012 lda = ldb / lda;
2013 }
2014
2015 else
2016 lda = 0;
2017 ldb = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2018 lda = lda * ldb;
2019 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i] = lda;
2020 }
2021 for (i = 0; i < ssizey; i++)
2022 working_space[ssizex * ssizey + 2 * ssizey * ssizey + i] =
2023 working_space[ssizex * ssizey + 2 * ssizey * ssizey + 3 * ssizex + i];
2024 }
2025 }
2026
2027/*write back resulting spectrum*/
2028 for (i = 0; i < ssizex; i++) {
2029 if (i < ssizey)
2030 source[i] = working_space[ssizex * ssizey + 2 * ssizey * ssizey + i];
2031
2032 else
2033 source[i] = 0;
2034 }
2035 delete[]working_space;
2036 return 0;
2037}
2038
2039
2040////////////////////////////////////////////////////////////////////////////////
2041/// One-dimensional high-resolution peak search function
2042///
2043/// This function searches for peaks in source spectrum. It is based on
2044/// deconvolution method. First the background is removed (if desired), then
2045/// Markov smoothed spectrum is calculated (if desired), then the response
2046/// function is generated according to given sigma and deconvolution is
2047/// carried out. The order of peaks is arranged according to their heights in
2048/// the spectrum after background elimination. The highest peak is the first in
2049/// the list. On success it returns number of found peaks.
2050///
2051/// #### Parameters:
2052///
2053/// - source: pointer to the vector of source spectrum.
2054/// - destVector: pointer to the vector of resulting deconvolved spectrum.
2055/// - ssize: length of source spectrum.
2056/// - sigma: sigma of searched peaks, for details we refer to manual.
2057/// - threshold: threshold value in % for selected peaks, peaks with
2058/// amplitude less than threshold*highest_peak/100
2059/// are ignored, see manual.
2060/// - backgroundRemove: logical variable, set if the removal of
2061/// background before deconvolution is desired.
2062/// - deconIterations-number of iterations in deconvolution operation.
2063/// - markov: logical variable, if it is true, first the source spectrum
2064/// is replaced by new spectrum calculated using Markov
2065/// chains method.
2066/// - averWindow: averaging window of searched peaks, for details
2067/// we refer to manual (applies only for Markov method).
2068///
2069/// ### Peaks searching:
2070///
2071/// The goal of this function is to identify automatically the peaks in spectrum
2072/// with the presence of the continuous background and statistical
2073/// fluctuations - noise.
2074///
2075/// The common problems connected with correct peak identification are:
2076///
2077/// - non-sensitivity to noise, i.e., only statistically
2078/// relevant peaks should be identified.
2079/// - non-sensitivity of the algorithm to continuous
2080/// background.
2081/// - ability to identify peaks close to the edges of the
2082/// spectrum region. Usually peak finders fail to detect them.
2083/// - resolution, decomposition of Double_tts and multiplets.
2084/// The algorithm should be able to recognise close positioned peaks.
2085/// - ability to identify peaks with different sigma.
2086///
2087/// \image html TSpectrum_Searching1.jpg Fig. 27 An example of one-dimensional synthetic spectrum with found peaks denoted by markers.
2088///
2089/// #### References:
2090///
2091/// 1. M.A. Mariscotti: A method for identification of peaks in the presence of
2092/// background and its application to spectrum analysis. NIM 50 (1967),
2093/// 309-320.
2094/// 2. M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky,
2095/// I. Turzo.:Identification of peaks in
2096/// multidimensional coincidence gamma-ray spectra. NIM, A443 (2000) 108-125.
2097/// 3. Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM
2098/// A 376 (1996), 451.
2099///
2100/// Examples of peak searching method:
2101///
2102/// The SearchHighRes function provides users with the possibility to vary the
2103/// input parameters and with the access to the output deconvolved data in the
2104/// destination spectrum. Based on the output data one can tune the parameters.
2105///
2106/// ### Example 15 - script SearchHR1.C:
2107///
2108/// One-dimensional spectrum with found peaks denoted by markers, 3 iterations
2109/// steps in the deconvolution.
2110///
2111/// #### Script:
2112///
2113/// Begin_Macro(source)
2114/// ../../../tutorials/spectrum/SearchHR1.C
2115/// End_Macro
2116///
2117/// ### Example 16 - script SearchHR3.C:
2118///
2119/// Influence of number of iterations (3-red, 10-blue, 100- green, 1000-magenta),
2120/// sigma=8, smoothing width=3.
2121///
2122/// Begin_Macro(source)
2123/// ../../../tutorials/spectrum/SearchHR3.C
2124/// End_Macro
2125
2126Int_t TSpectrum::SearchHighRes(Double_t *source,Double_t *destVector, int ssize,
2127 Double_t sigma, Double_t threshold,
2128 bool backgroundRemove,int deconIterations,
2129 bool markov, int averWindow)
2130{
2131 int i, j, numberIterations = (Int_t)(7 * sigma + 0.5);
2132 Double_t a, b, c;
2133 int k, lindex, posit, imin, imax, jmin, jmax, lh_gold, priz;
2134 Double_t lda, ldb, ldc, area, maximum, maximum_decon;
2135 int xmin, xmax, l, peak_index = 0, size_ext = ssize + 2 * numberIterations, shift = numberIterations, bw = 2, w;
2136 Double_t maxch;
2137 Double_t nom, nip, nim, sp, sm, plocha = 0;
2138 Double_t m0low=0,m1low=0,m2low=0,l0low=0,l1low=0,detlow,av,men;
2139 if (sigma < 1) {
2140 Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
2141 return 0;
2142 }
2143
2144 if(threshold<=0 || threshold>=100){
2145 Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
2146 return 0;
2147 }
2148
2149 j = (Int_t) (5.0 * sigma + 0.5);
2150 if (j >= PEAK_WINDOW / 2) {
2151 Error("SearchHighRes", "Too large sigma");
2152 return 0;
2153 }
2154
2155 if (markov == true) {
2156 if (averWindow <= 0) {
2157 Error("SearchHighRes", "Averaging window must be positive");
2158 return 0;
2159 }
2160 }
2161
2162 if(backgroundRemove == true){
2163 if(ssize < 2 * numberIterations + 1){
2164 Error("SearchHighRes", "Too large clipping window");
2165 return 0;
2166 }
2167 }
2168
2169 k = int(2 * sigma+0.5);
2170 if(k >= 2){
2171 for(i = 0;i < k;i++){
2172 a = i,b = source[i];
2173 m0low += 1,m1low += a,m2low += a * a,l0low += b,l1low += a * b;
2174 }
2175 detlow = m0low * m2low - m1low * m1low;
2176 if(detlow != 0)
2177 l1low = (-l0low * m1low + l1low * m0low) / detlow;
2178
2179 else
2180 l1low = 0;
2181 if(l1low > 0)
2182 l1low=0;
2183 }
2184
2185 else{
2186 l1low = 0;
2187 }
2188
2189 i = (Int_t)(7 * sigma + 0.5);
2190 i = 2 * i;
2191 Double_t *working_space = new Double_t [7 * (ssize + i)];
2192 for (j=0;j<7 * (ssize + i);j++) working_space[j] = 0;
2193 for(i = 0; i < size_ext; i++){
2194 if(i < shift){
2195 a = i - shift;
2196 working_space[i + size_ext] = source[0] + l1low * a;
2197 if(working_space[i + size_ext] < 0)
2198 working_space[i + size_ext]=0;
2199 }
2200
2201 else if(i >= ssize + shift){
2202 a = i - (ssize - 1 + shift);
2203 working_space[i + size_ext] = source[ssize - 1];
2204 if(working_space[i + size_ext] < 0)
2205 working_space[i + size_ext]=0;
2206 }
2207
2208 else
2209 working_space[i + size_ext] = source[i - shift];
2210 }
2211
2212 if(backgroundRemove == true){
2213 for(i = 1; i <= numberIterations; i++){
2214 for(j = i; j < size_ext - i; j++){
2215 if(markov == false){
2216 a = working_space[size_ext + j];
2217 b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2218 if(b < a)
2219 a = b;
2220
2221 working_space[j]=a;
2222 }
2223
2224 else{
2225 a = working_space[size_ext + j];
2226 av = 0;
2227 men = 0;
2228 for (w = j - bw; w <= j + bw; w++){
2229 if ( w >= 0 && w < size_ext){
2230 av += working_space[size_ext + w];
2231 men +=1;
2232 }
2233 }
2234 av = av / men;
2235 b = 0;
2236 men = 0;
2237 for (w = j - i - bw; w <= j - i + bw; w++){
2238 if ( w >= 0 && w < size_ext){
2239 b += working_space[size_ext + w];
2240 men +=1;
2241 }
2242 }
2243 b = b / men;
2244 c = 0;
2245 men = 0;
2246 for (w = j + i - bw; w <= j + i + bw; w++){
2247 if ( w >= 0 && w < size_ext){
2248 c += working_space[size_ext + w];
2249 men +=1;
2250 }
2251 }
2252 c = c / men;
2253 b = (b + c) / 2;
2254 if (b < a)
2255 av = b;
2256 working_space[j]=av;
2257 }
2258 }
2259 for(j = i; j < size_ext - i; j++)
2260 working_space[size_ext + j] = working_space[j];
2261 }
2262 for(j = 0;j < size_ext; j++){
2263 if(j < shift){
2264 a = j - shift;
2265 b = source[0] + l1low * a;
2266 if (b < 0) b = 0;
2267 working_space[size_ext + j] = b - working_space[size_ext + j];
2268 }
2269
2270 else if(j >= ssize + shift){
2271 a = j - (ssize - 1 + shift);
2272 b = source[ssize - 1];
2273 if (b < 0) b = 0;
2274 working_space[size_ext + j] = b - working_space[size_ext + j];
2275 }
2276
2277 else{
2278 working_space[size_ext + j] = source[j - shift] - working_space[size_ext + j];
2279 }
2280 }
2281 for(j = 0;j < size_ext; j++){
2282 if(working_space[size_ext + j] < 0) working_space[size_ext + j] = 0;
2283 }
2284 }
2285
2286 for(i = 0; i < size_ext; i++){
2287 working_space[i + 6*size_ext] = working_space[i + size_ext];
2288 }
2289
2290 if(markov == true){
2291 for(j = 0; j < size_ext; j++)
2292 working_space[2 * size_ext + j] = working_space[size_ext + j];
2293 xmin = 0,xmax = size_ext - 1;
2294 for(i = 0, maxch = 0; i < size_ext; i++){
2295 working_space[i] = 0;
2296 if(maxch < working_space[2 * size_ext + i])
2297 maxch = working_space[2 * size_ext + i];
2298 plocha += working_space[2 * size_ext + i];
2299 }
2300 if(maxch == 0) {
2301 delete [] working_space;
2302 return 0;
2303 }
2304
2305 nom = 1;
2306 working_space[xmin] = 1;
2307 for(i = xmin; i < xmax; i++){
2308 nip = working_space[2 * size_ext + i] / maxch;
2309 nim = working_space[2 * size_ext + i + 1] / maxch;
2310 sp = 0,sm = 0;
2311 for(l = 1; l <= averWindow; l++){
2312 if((i + l) > xmax)
2313 a = working_space[2 * size_ext + xmax] / maxch;
2314
2315 else
2316 a = working_space[2 * size_ext + i + l] / maxch;
2317
2318 b = a - nip;
2319 if(a + nip <= 0)
2320 a=1;
2321
2322 else
2323 a = TMath::Sqrt(a + nip);
2324
2325 b = b / a;
2326 b = TMath::Exp(b);
2327 sp = sp + b;
2328 if((i - l + 1) < xmin)
2329 a = working_space[2 * size_ext + xmin] / maxch;
2330
2331 else
2332 a = working_space[2 * size_ext + i - l + 1] / maxch;
2333
2334 b = a - nim;
2335 if(a + nim <= 0)
2336 a = 1;
2337
2338 else
2339 a = TMath::Sqrt(a + nim);
2340
2341 b = b / a;
2342 b = TMath::Exp(b);
2343 sm = sm + b;
2344 }
2345 a = sp / sm;
2346 a = working_space[i + 1] = working_space[i] * a;
2347 nom = nom + a;
2348 }
2349 for(i = xmin; i <= xmax; i++){
2350 working_space[i] = working_space[i] / nom;
2351 }
2352 for(j = 0; j < size_ext; j++)
2353 working_space[size_ext + j] = working_space[j] * plocha;
2354 for(j = 0; j < size_ext; j++){
2355 working_space[2 * size_ext + j] = working_space[size_ext + j];
2356 }
2357 if(backgroundRemove == true){
2358 for(i = 1; i <= numberIterations; i++){
2359 for(j = i; j < size_ext - i; j++){
2360 a = working_space[size_ext + j];
2361 b = (working_space[size_ext + j - i] + working_space[size_ext + j + i]) / 2.0;
2362 if(b < a)
2363 a = b;
2364 working_space[j] = a;
2365 }
2366 for(j = i; j < size_ext - i; j++)
2367 working_space[size_ext + j] = working_space[j];
2368 }
2369 for(j = 0; j < size_ext; j++){
2370 working_space[size_ext + j] = working_space[2 * size_ext + j] - working_space[size_ext + j];
2371 }
2372 }
2373 }
2374//deconvolution starts
2375 area = 0;
2376 lh_gold = -1;
2377 posit = 0;
2378 maximum = 0;
2379//generate response vector
2380 for(i = 0; i < size_ext; i++){
2381 lda = (Double_t)i - 3 * sigma;
2382 lda = lda * lda / (2 * sigma * sigma);
2383 j = (Int_t)(1000 * TMath::Exp(-lda));
2384 lda = j;
2385 if(lda != 0)
2386 lh_gold = i + 1;
2387
2388 working_space[i] = lda;
2389 area = area + lda;
2390 if(lda > maximum){
2391 maximum = lda;
2392 posit = i;
2393 }
2394 }
2395//read source vector
2396 for(i = 0; i < size_ext; i++)
2397 working_space[2 * size_ext + i] = TMath::Abs(working_space[size_ext + i]);
2398//create matrix at*a(vector b)
2399 i = lh_gold - 1;
2400 if(i > size_ext)
2401 i = size_ext;
2402
2403 imin = -i,imax = i;
2404 for(i = imin; i <= imax; i++){
2405 lda = 0;
2406 jmin = 0;
2407 if(i < 0)
2408 jmin = -i;
2409 jmax = lh_gold - 1 - i;
2410 if(jmax > (lh_gold - 1))
2411 jmax = lh_gold - 1;
2412
2413 for(j = jmin;j <= jmax; j++){
2414 ldb = working_space[j];
2415 ldc = working_space[i + j];
2416 lda = lda + ldb * ldc;
2417 }
2418 working_space[size_ext + i - imin] = lda;
2419 }
2420//create vector p
2421 i = lh_gold - 1;
2422 imin = -i,imax = size_ext + i - 1;
2423 for(i = imin; i <= imax; i++){
2424 lda = 0;
2425 for(j = 0; j <= (lh_gold - 1); j++){
2426 ldb = working_space[j];
2427 k = i + j;
2428 if(k >= 0 && k < size_ext){
2429 ldc = working_space[2 * size_ext + k];
2430 lda = lda + ldb * ldc;
2431 }
2432
2433 }
2434 working_space[4 * size_ext + i - imin] = lda;
2435 }
2436//move vector p
2437 for(i = imin; i <= imax; i++)
2438 working_space[2 * size_ext + i - imin] = working_space[4 * size_ext + i - imin];
2439//initialization of resulting vector
2440 for(i = 0; i < size_ext; i++)
2441 working_space[i] = 1;
2442//START OF ITERATIONS
2443 for(lindex = 0; lindex < deconIterations; lindex++){
2444 for(i = 0; i < size_ext; i++){
2445 if(TMath::Abs(working_space[2 * size_ext + i]) > 0.00001 && TMath::Abs(working_space[i]) > 0.00001){
2446 lda=0;
2447 jmin = lh_gold - 1;
2448 if(jmin > i)
2449 jmin = i;
2450
2451 jmin = -jmin;
2452 jmax = lh_gold - 1;
2453 if(jmax > (size_ext - 1 - i))
2454 jmax=size_ext-1-i;
2455
2456 for(j = jmin; j <= jmax; j++){
2457 ldb = working_space[j + lh_gold - 1 + size_ext];
2458 ldc = working_space[i + j];
2459 lda = lda + ldb * ldc;
2460 }
2461 ldb = working_space[2 * size_ext + i];
2462 if(lda != 0)
2463 lda = ldb / lda;
2464
2465 else
2466 lda = 0;
2467
2468 ldb = working_space[i];
2469 lda = lda * ldb;
2470 working_space[3 * size_ext + i] = lda;
2471 }
2472 }
2473 for(i = 0; i < size_ext; i++){
2474 working_space[i] = working_space[3 * size_ext + i];
2475 }
2476 }
2477//shift resulting spectrum
2478 for(i=0;i<size_ext;i++){
2479 lda = working_space[i];
2480 j = i + posit;
2481 j = j % size_ext;
2482 working_space[size_ext + j] = lda;
2483 }
2484//write back resulting spectrum
2485 maximum = 0, maximum_decon = 0;
2486 j = lh_gold - 1;
2487 for(i = 0; i < size_ext - j; i++){
2488 if(i >= shift && i < ssize + shift){
2489 working_space[i] = area * working_space[size_ext + i + j];
2490 if(maximum_decon < working_space[i])
2491 maximum_decon = working_space[i];
2492 if(maximum < working_space[6 * size_ext + i])
2493 maximum = working_space[6 * size_ext + i];
2494 }
2495
2496 else
2497 working_space[i] = 0;
2498 }
2499 lda=1;
2500 if(lda>threshold)
2501 lda=threshold;
2502 lda=lda/100;
2503
2504//searching for peaks in deconvolved spectrum
2505 for(i = 1; i < size_ext - 1; i++){
2506 if(working_space[i] > working_space[i - 1] && working_space[i] > working_space[i + 1]){
2507 if(i >= shift && i < ssize + shift){
2508 if(working_space[i] > lda*maximum_decon && working_space[6 * size_ext + i] > threshold * maximum / 100.0){
2509 for(j = i - 1, a = 0, b = 0; j <= i + 1; j++){
2510 a += (Double_t)(j - shift) * working_space[j];
2511 b += working_space[j];
2512 }
2513 a = a / b;
2514 if(a < 0)
2515 a = 0;
2516
2517 if(a >= ssize)
2518 a = ssize - 1;
2519 if(peak_index == 0){
2520 fPositionX[0] = a;
2521 peak_index = 1;
2522 }
2523
2524 else{
2525 for(j = 0, priz = 0; j < peak_index && priz == 0; j++){
2526 if(working_space[6 * size_ext + shift + (Int_t)a] > working_space[6 * size_ext + shift + (Int_t)fPositionX[j]])
2527 priz = 1;
2528 }
2529 if(priz == 0){
2530 if(j < fMaxPeaks){
2531 fPositionX[j] = a;
2532 }
2533 }
2534
2535 else{
2536 for(k = peak_index; k >= j; k--){
2537 if(k < fMaxPeaks){
2538 fPositionX[k] = fPositionX[k - 1];
2539 }
2540 }
2541 fPositionX[j - 1] = a;
2542 }
2543 if(peak_index < fMaxPeaks)
2544 peak_index += 1;
2545 }
2546 }
2547 }
2548 }
2549 }
2550
2551 for(i = 0; i < ssize; i++) destVector[i] = working_space[i + shift];
2552 delete [] working_space;
2553 fNPeaks = peak_index;
2554 if(peak_index == fMaxPeaks)
2555 Warning("SearchHighRes", "Peak buffer full");
2556 return fNPeaks;
2557}
2558
2559////////////////////////////////////////////////////////////////////////////////
2560/// Old name of SearcHighRes introduced for back compatibility.
2561/// This function will be removed after the June 2006 release
2562
2563Int_t TSpectrum::Search1HighRes(Double_t *source,Double_t *destVector, int ssize,
2564 Double_t sigma, Double_t threshold,
2565 bool backgroundRemove,int deconIterations,
2566 bool markov, int averWindow)
2567{
2568
2569
2570 return SearchHighRes(source,destVector,ssize,sigma,threshold,backgroundRemove,
2571 deconIterations,markov,averWindow);
2572}
2573
2574////////////////////////////////////////////////////////////////////////////////
2575/// Static function, interface to TSpectrum::Search.
2576
2577Int_t TSpectrum::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
2578{
2579
2580 TSpectrum s;
2581 return s.Search(hist,sigma,option,threshold);
2582}
2583
2584////////////////////////////////////////////////////////////////////////////////
2585/// Static function, interface to TSpectrum::Background.
2586
2587TH1 *TSpectrum::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
2588{
2589 TSpectrum s;
2590 return s.Background(hist,niter,option);
2591}
d
#define d(i)
Definition: RSha256.hxx:102
b
#define b(i)
Definition: RSha256.hxx:100
c
#define c(i)
Definition: RSha256.hxx:101
h
#define h(i)
Definition: RSha256.hxx:106
e
#define e(i)
Definition: RSha256.hxx:103
Int_t
int Int_t
Definition: RtypesCore.h:41
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:88
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
Double_t
double Double_t
Definition: RtypesCore.h:55
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:87
Option_t
const char Option_t
Definition: RtypesCore.h:62
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:363
kRed
@ kRed
Definition: Rtypes.h:63
xmin
float xmin
Definition: THbookFile.cxx:93
xmax
float xmax
Definition: THbookFile.cxx:93
hbname
#define hbname
Definition: THbookFile.cxx:123
PEAK_WINDOW
#define PEAK_WINDOW
Definition: TSpectrum.cxx:44
gPad
#define gPad
Definition: TVirtualPad.h:286
snprintf
#define snprintf
Definition: civetweb.c:1540
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
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::GetLast
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:455
TAxis::GetFirst
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:444
TH1
The TH1 histogram class.
Definition: TH1.h:56
TH1::GetBinCenter
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8462
TH1::GetDimension
virtual Int_t GetDimension() const
Definition: TH1.h:278
TH1::Reset
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH1.cxx:6608
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
TH1::SetBinContent
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8542
TH1::GetListOfFunctions
TList * GetListOfFunctions() const
Definition: TH1.h:239
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2974
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4790
TH1::SetEntries
virtual void SetEntries(Double_t n)
Definition: TH1.h:381
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:818
TList::FindObject
virtual TObject * FindObject(const char *name) const
Delete a TObjLink object.
Definition: TList.cxx:574
TList::Delete
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition: TList.cxx:467
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
TObject::Warning
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
TPolyMarker
A PolyMarker is defined by an array on N points in a 2-D space.
Definition: TPolyMarker.h:31
TSpectrum
Advanced Spectra Processing.
Definition: TSpectrum.h:18
TSpectrum::fResolution
Double_t fResolution
NOT USED resolution of the neighboring peaks
Definition: TSpectrum.h:30
TSpectrum::StaticBackground
static TH1 * StaticBackground(const TH1 *hist, Int_t niter=20, Option_t *option="")
Static function, interface to TSpectrum::Background.
Definition: TSpectrum.cxx:2587
TSpectrum::fMaxPeaks
Int_t fMaxPeaks
Maximum number of peaks to be found.
Definition: TSpectrum.h:25
TSpectrum::fPositionX
Double_t * fPositionX
[fNPeaks] X position of peaks
Definition: TSpectrum.h:28
TSpectrum::kBackOrder8
@ kBackOrder8
Definition: TSpectrum.h:40
TSpectrum::kBackSmoothing11
@ kBackSmoothing11
Definition: TSpectrum.h:47
TSpectrum::kBackSmoothing15
@ kBackSmoothing15
Definition: TSpectrum.h:49
TSpectrum::kBackOrder2
@ kBackOrder2
Definition: TSpectrum.h:37
TSpectrum::kBackOrder6
@ kBackOrder6
Definition: TSpectrum.h:39
TSpectrum::kBackSmoothing7
@ kBackSmoothing7
Definition: TSpectrum.h:45
TSpectrum::kBackIncreasingWindow
@ kBackIncreasingWindow
Definition: TSpectrum.h:41
TSpectrum::kBackDecreasingWindow
@ kBackDecreasingWindow
Definition: TSpectrum.h:42
TSpectrum::kBackSmoothing13
@ kBackSmoothing13
Definition: TSpectrum.h:48
TSpectrum::kBackSmoothing9
@ kBackSmoothing9
Definition: TSpectrum.h:46
TSpectrum::kBackOrder4
@ kBackOrder4
Definition: TSpectrum.h:38
TSpectrum::kBackSmoothing5
@ kBackSmoothing5
Definition: TSpectrum.h:44
TSpectrum::kBackSmoothing3
@ kBackSmoothing3
Definition: TSpectrum.h:43
TSpectrum::SmoothMarkov
const char * SmoothMarkov(Double_t *source, Int_t ssize, Int_t averWindow)
One-dimensional markov spectrum smoothing function.
Definition: TSpectrum.cxx:1195
TSpectrum::Unfolding
const char * Unfolding(Double_t *source, const Double_t **respMatrix, Int_t ssizex, Int_t ssizey, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional unfolding function.
Definition: TSpectrum.cxx:1889
TSpectrum::SetResolution
void SetResolution(Double_t resolution=1)
NOT USED resolution: determines resolution of the neighbouring peaks default value is 1 correspond to...
Definition: TSpectrum.cxx:351
TSpectrum::~TSpectrum
virtual ~TSpectrum()
Destructor.
Definition: TSpectrum.cxx:87
TSpectrum::TSpectrum
TSpectrum()
Constructor.
Definition: TSpectrum.cxx:50
TSpectrum::Background
virtual TH1 * Background(const TH1 *hist, Int_t niter=20, Option_t *option="")
One-dimensional background estimation function.
Definition: TSpectrum.cxx:152
TSpectrum::SearchHighRes
Int_t SearchHighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
One-dimensional high-resolution peak search function.
Definition: TSpectrum.cxx:2126
TSpectrum::Search
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05)
One-dimensional peak search function.
Definition: TSpectrum.cxx:266
TSpectrum::fNPeaks
Int_t fNPeaks
number of peaks found
Definition: TSpectrum.h:26
TSpectrum::Search1HighRes
Int_t Search1HighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
Old name of SearcHighRes introduced for back compatibility.
Definition: TSpectrum.cxx:2563
TSpectrum::Deconvolution
const char * Deconvolution(Double_t *source, const Double_t *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional deconvolution function.
Definition: TSpectrum.cxx:1459
TSpectrum::fgAverageWindow
static Int_t fgAverageWindow
Average window of searched peaks.
Definition: TSpectrum.h:32
TSpectrum::SetDeconIterations
static void SetDeconIterations(Int_t n=3)
Static function: Set max number of decon iterations in deconvolution operation (see TSpectrum::Search...
Definition: TSpectrum.cxx:108
TSpectrum::fPosition
Double_t * fPosition
[fNPeaks] array of current peak positions
Definition: TSpectrum.h:27
TSpectrum::SetAverageWindow
static void SetAverageWindow(Int_t w=3)
Static function: Set average window of searched peaks (see TSpectrum::SearchHighRes).
Definition: TSpectrum.cxx:99
TSpectrum::DeconvolutionRL
const char * DeconvolutionRL(Double_t *source, const Double_t *response, Int_t ssize, Int_t numberIterations, Int_t numberRepetitions, Double_t boost)
One-dimensional deconvolution function.
Definition: TSpectrum.cxx:1665
TSpectrum::fgIterations
static Int_t fgIterations
Maximum number of decon iterations (default=3)
Definition: TSpectrum.h:33
TSpectrum::fHistogram
TH1 * fHistogram
resulting histogram
Definition: TSpectrum.h:31
TSpectrum::Print
virtual void Print(Option_t *option="") const
Print the array of positions.
Definition: TSpectrum.cxx:219
TSpectrum::StaticSearch
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05)
Static function, interface to TSpectrum::Search.
Definition: TSpectrum.cxx:2577
TSpectrum::fPositionY
Double_t * fPositionY
[fNPeaks] Y position of peaks
Definition: TSpectrum.h:29
TString
Basic string class.
Definition: TString.h:131
TString::ToLower
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1100
TString::Data
const char * Data() const
Definition: TString.h:364
TString::ReplaceAll
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
TString::Contains
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::VectorUtil::boost
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only re...
Definition: VectorUtil.h:327
TGeoUnit::s
static constexpr double s
Definition: TGeoSystemOfUnits.h:162
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:715
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:723
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
first
Definition: first.py:1
l
auto * l
Definition: textangle.C:4
a
auto * a
Definition: textangle.C:12
dest
#define dest(otri, vertexptr)
Definition: triangle.c:1040