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