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