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