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