Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
TSpectrumFit.cxx
Go to the documentation of this file.
1// @(#)root/spectrum:$Id$
2// Author: Miroslav Morhac 25/09/06
3
4/** \class TSpectrumFit
5 \ingroup Spectrum
6 \brief Advanced 1-dimensional spectra fitting functions
7 \author Miroslav Morhac
8
9 \legacy{TSpectrumFit, For modeling a spectrum fitting and estimating the background one can use RooFit while for deconvolution and unfolding one can use TUnfold.}
10
11 Class for fitting 1D spectra using AWMI (algorithm without matrix
12 inversion) and conjugate gradient algorithms for symmetrical
13 matrices (Stiefel-Hestens method). AWMI method allows to fit
14 simultaneously 100s up to 1000s peaks. Stiefel method is very stable,
15 it converges faster, but is more time consuming
16
17 The algorithms in this class have been published in the following
18 references:
19 1. M. Morhac et al.: Efficient fitting algorithms applied to
20 analysis of coincidence gamma-ray spectra. Computer Physics
21 Communications, Vol 172/1 (2005) pp. 19-41.
22
23 2. M. Morhac et al.: Study of fitting algorithms applied to
24 simultaneous analysis of large number of peaks in gamma-ray spectra.
25 Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003.
26
27*/
28
29#include "TSpectrumFit.h"
30#include "TMath.h"
31
32
33////////////////////////////////////////////////////////////////////////////////
34/// Default constructor
35
36TSpectrumFit::TSpectrumFit() :TNamed("SpectrumFit", "Miroslav Morhac peak fitter")
37{
38 fNPeaks = 0;
39 fNumberIterations = 1;
40 fXmin = 0;
41 fXmax = 100;
42 fStatisticType = kFitOptimChiCounts;
43 fAlphaOptim = kFitAlphaHalving;
44 fPower = kFitPower2;
45 fFitTaylor = kFitTaylorOrderFirst;
46 fAlpha =1;
47 fChi = 0;
48 fPositionInit = nullptr;
49 fPositionCalc = nullptr;
50 fPositionErr = nullptr;
51 fFixPosition = nullptr;
52 fAmpInit = nullptr;
53 fAmpCalc = nullptr;
54 fAmpErr = nullptr;
55 fFixAmp = nullptr;
56 fArea = nullptr;
57 fAreaErr = nullptr;
58 fSigmaInit = 2;
59 fSigmaCalc = 1;
60 fSigmaErr = 0;
61 fTInit = 0;
62 fTCalc = 0;
63 fTErr = 0;
64 fBInit = 1;
65 fBCalc = 0;
66 fBErr = 0;
67 fSInit = 0;
68 fSCalc = 0;
69 fSErr = 0;
70 fA0Init = 0;
71 fA0Calc = 0;
72 fA0Err = 0;
73 fA1Init = 0;
74 fA1Calc = 0;
75 fA1Err = 0;
76 fA2Init = 0;
77 fA2Calc = 0;
78 fA2Err = 0;
79 fFixSigma = false;
80 fFixT = true;
81 fFixB = true;
82 fFixS = true;
83 fFixA0 = true;
84 fFixA1 = true;
85 fFixA2 = true;
86}
87
88////////////////////////////////////////////////////////////////////////////////
89/// numberPeaks: number of fitted peaks (must be greater than zero)
90///
91/// the constructor allocates arrays for all fitted parameters (peak positions, amplitudes etc) and sets the member
92/// variables to their default values. One can change these variables by member functions (setters) of TSpectrumFit class.
93///
94/// Shape function of the fitted peaks is
95///
96/// \image html spectrumfit_constructor_image001.gif
97///
98/// where a represents vector of
99/// fitted parameters (positions p(j), amplitudes A(j), sigma, relative amplitudes
100/// T, S and slope B).
101
102TSpectrumFit::TSpectrumFit(Int_t numberPeaks) :TNamed("SpectrumFit", "Miroslav Morhac peak fitter")
103{
104 if (numberPeaks <= 0){
105 Error ("TSpectrumFit","Invalid number of peaks, must be > than 0");
106 return;
107 }
108 fNPeaks = numberPeaks;
109 fNumberIterations = 1;
110 fXmin = 0;
111 fXmax = 100;
112 fStatisticType = kFitOptimChiCounts;
113 fAlphaOptim = kFitAlphaHalving;
114 fPower = kFitPower2;
115 fFitTaylor = kFitTaylorOrderFirst;
116 fAlpha =1;
117 fChi = 0;
118 fPositionInit = new Double_t[numberPeaks];
119 fPositionCalc = new Double_t[numberPeaks];
120 fPositionErr = new Double_t[numberPeaks];
121 fFixPosition = new Bool_t[numberPeaks];
122 fAmpInit = new Double_t[numberPeaks];
123 fAmpCalc = new Double_t[numberPeaks];
124 fAmpErr = new Double_t[numberPeaks];
125 fFixAmp = new Bool_t[numberPeaks];
126 fArea = new Double_t[numberPeaks];
127 fAreaErr = new Double_t[numberPeaks];
128 fSigmaInit = 2;
129 fSigmaCalc = 1;
130 fSigmaErr = 0;
131 fTInit = 0;
132 fTCalc = 0;
133 fTErr = 0;
134 fBInit = 1;
135 fBCalc = 0;
136 fBErr = 0;
137 fSInit = 0;
138 fSCalc = 0;
139 fSErr = 0;
140 fA0Init = 0;
141 fA0Calc = 0;
142 fA0Err = 0;
143 fA1Init = 0;
144 fA1Calc = 0;
145 fA1Err = 0;
146 fA2Init = 0;
147 fA2Calc = 0;
148 fA2Err = 0;
149 fFixSigma = false;
150 fFixT = true;
151 fFixB = true;
152 fFixS = true;
153 fFixA0 = true;
154 fFixA1 = true;
155 fFixA2 = true;
156}
157
158////////////////////////////////////////////////////////////////////////////////
159/// Destructor
160
161TSpectrumFit::~TSpectrumFit()
162{
163 delete [] fPositionInit;
164 delete [] fPositionCalc;
165 delete [] fPositionErr;
166 delete [] fFixPosition;
167 delete [] fAmpInit;
168 delete [] fAmpCalc;
169 delete [] fAmpErr;
170 delete [] fFixAmp;
171 delete [] fArea;
172 delete [] fAreaErr;
173}
174
175/////////////////////////////////////////////////////////////////////////////////
176// This function calculates error function of x.
177
178Double_t TSpectrumFit::Erfc(Double_t x)
179{
180 Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap = 0.47047;
181 Double_t a, t, c, w;
182 a = TMath::Abs(x);
183 w = 1. + dap * a;
184 t = 1. / w;
185 w = a * a;
186 if (w < 700)
187 c = exp(-w);
188
189 else {
190 c = 0;
191 }
192 c = c * t * (da1 + t * (da2 + t * da3));
193 if (x < 0)
194 c = 1. - c;
195 return (c);
196}
197
198////////////////////////////////////////////////////////////////////////////////
199/// This function calculates derivative of error function of x.
200
201Double_t TSpectrumFit::Derfc(Double_t x)
202{
203 Double_t a, t, c, w;
204 Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap = 0.47047;
205 a = TMath::Abs(x);
206 w = 1. + dap * a;
207 t = 1. / w;
208 w = a * a;
209 if (w < 700)
210 c = exp(-w);
211
212 else {
213 c = 0;
214 }
215 c = (-1.) * dap * c * t * t * (da1 + t * (2. * da2 + t * 3. * da3)) -
216 2. * a * Erfc(a);
217 return (c);
218}
219
220////////////////////////////////////////////////////////////////////////////////
221/// This function calculates derivative of peak shape function (see manual)
222/// according to amplitude of peak.
223/// Function parameters:
224/// - i-channel
225/// - i0-position of peak
226/// - sigma-sigma of peak
227/// - t, s-relative amplitudes
228/// - b-slope
229
230Double_t TSpectrumFit::Deramp(Double_t i, Double_t i0, Double_t sigma, Double_t t,
231 Double_t s, Double_t b)
232{
233 Double_t p, q, r, a;
234 p = (i - i0) / sigma;
235 if ((p * p) < 700)
236 q = exp(-p * p);
237
238 else {
239 q = 0;
240 }
241 r = 0;
242 if (t != 0) {
243 a = p / b;
244 if (a > 700)
245 a = 700;
246 r = t * exp(a) / 2.;
247 }
248 if (r != 0)
249 r = r * Erfc(p + 1. / (2. * b));
250 q = q + r;
251 if (s != 0)
252 q = q + s * Erfc(p) / 2.;
253 return (q);
254}
255
256////////////////////////////////////////////////////////////////////////////////
257/// This function calculates derivative of peak shape function (see manual)
258/// according to peak position.
259/// Function parameters:
260/// - i-channel
261/// - amp-amplitude of peak
262/// - i0-position of peak
263/// - sigma-sigma of peak
264/// - t, s-relative amplitudes
265/// - b-slope
266
267Double_t TSpectrumFit::Deri0(Double_t i, Double_t amp, Double_t i0, Double_t sigma,
268 Double_t t, Double_t s, Double_t b)
269{
270 Double_t p, r1, r2, r3, r4, c, d, e;
271 p = (i - i0) / sigma;
272 d = 2. * sigma;
273 if ((p * p) < 700)
274 r1 = 2. * p * exp(-p * p) / sigma;
275
276 else {
277 r1 = 0;
278 }
279 r2 = 0, r3 = 0;
280 if (t != 0) {
281 c = p + 1. / (2. * b);
282 e = p / b;
283 if (e > 700)
284 e = 700;
285 r2 = -t * exp(e) * Erfc(c) / (d * b);
286 r3 = -t * exp(e) * Derfc(c) / d;
287 }
288 r4 = 0;
289 if (s != 0)
290 r4 = -s * Derfc(p) / d;
291 r1 = amp * (r1 + r2 + r3 + r4);
292 return (r1);
293}
294
295////////////////////////////////////////////////////////////////////////////////
296/// This function calculates second derivative of peak shape function
297/// (see manual) according to peak position.
298/// Function parameters:
299/// - i-channel
300/// - amp-amplitude of peak
301/// - i0-position of peak
302/// - sigma-width of peak
303
304Double_t TSpectrumFit::Derderi0(Double_t i, Double_t amp, Double_t i0,
305 Double_t sigma)
306{
307 Double_t p, r1, r2, r3, r4;
308 p = (i - i0) / sigma;
309 if ((p * p) < 700)
310 r1 = exp(-p * p);
311
312 else {
313 r1 = 0;
314 }
315 r1 = r1 * (4 * p * p - 2) / (sigma * sigma);
316 r2 = 0, r3 = 0, r4 = 0;
317 r1 = amp * (r1 + r2 + r3 + r4);
318 return (r1);
319}
320
321////////////////////////////////////////////////////////////////////////////////
322/// This function calculates derivative of peaks shape function (see manual)
323/// according to sigma of peaks.
324/// Function parameters:
325/// - num_of_fitted_peaks-number of fitted peaks
326/// - i-channel
327/// - parameter-array of peaks parameters (amplitudes and positions)
328/// - sigma-sigma of peak
329/// - t, s-relative amplitudes
330/// - b-slope
331
332Double_t TSpectrumFit::Dersigma(Int_t num_of_fitted_peaks, Double_t i,
333 const Double_t *parameter, Double_t sigma,
334 Double_t t, Double_t s, Double_t b)
335{
336 Int_t j;
337 Double_t r, p, r1, r2, r3, r4, c, d, e;
338 r = 0;
339 d = 2. * sigma;
340 for (j = 0; j < num_of_fitted_peaks; j++) {
341 p = (i - parameter[2 * j + 1]) / sigma;
342 r1 = 0;
343 if (TMath::Abs(p) < 3) {
344 if ((p * p) < 700)
345 r1 = 2. * p * p * exp(-p * p) / sigma;
346
347 else {
348 r1 = 0;
349 }
350 }
351 r2 = 0, r3 = 0;
352 if (t != 0) {
353 c = p + 1. / (2. * b);
354 e = p / b;
355 if (e > 700)
356 e = 700;
357 r2 = -t * p * exp(e) * Erfc(c) / (d * b);
358 r3 = -t * p * exp(e) * Derfc(c) / d;
359 }
360 r4 = 0;
361 if (s != 0)
362 r4 = -s * p * Derfc(p) / d;
363 r = r + parameter[2 * j] * (r1 + r2 + r3 + r4);
364 }
365 return (r);
366}
367
368////////////////////////////////////////////////////////////////////////////////
369/// This function calculates second derivative of peaks shape function
370/// (see manual) according to sigma of peaks.
371/// Function parameters:
372/// - num_of_fitted_peaks-number of fitted peaks
373/// - i-channel
374/// - parameter-array of peaks parameters (amplitudes and positions)
375/// - sigma-sigma of peak
376
377Double_t TSpectrumFit::Derdersigma(Int_t num_of_fitted_peaks, Double_t i,
378 const Double_t *parameter, Double_t sigma)
379{
380 Int_t j;
381 Double_t r, p, r1, r2, r3, r4;
382 r = 0;
383 for (j = 0; j < num_of_fitted_peaks; j++) {
384 p = (i - parameter[2 * j + 1]) / sigma;
385 r1 = 0;
386 if (TMath::Abs(p) < 3) {
387 if ((p * p) < 700)
388 r1 = exp(-p * p) * p * p * (4. * p * p - 6) / (sigma * sigma);
389
390 else {
391 r1 = 0;
392 }
393 }
394 r2 = 0, r3 = 0, r4 = 0;
395 r = r + parameter[2 * j] * (r1 + r2 + r3 + r4);
396 }
397 return (r);
398}
399
400////////////////////////////////////////////////////////////////////////////////
401/// This function calculates derivative of peaks shape function (see manual)
402/// according to relative amplitude t.
403/// Function parameters:
404/// - num_of_fitted_peaks-number of fitted peaks
405/// - i-channel
406/// - parameter-array of peaks parameters (amplitudes and positions)
407/// - sigma-sigma of peak
408/// - b-slope
409
410Double_t TSpectrumFit::Dert(Int_t num_of_fitted_peaks, Double_t i,
411 const Double_t *parameter, Double_t sigma, Double_t b)
412{
413 Int_t j;
414 Double_t r, p, r1, c, e;
415 r = 0;
416 for (j = 0; j < num_of_fitted_peaks; j++) {
417 p = (i - parameter[2 * j + 1]) / sigma;
418 c = p + 1. / (2. * b);
419 e = p / b;
420 if (e > 700)
421 e = 700;
422 r1 = exp(e) * Erfc(c);
423 r = r + parameter[2 * j] * r1;
424 }
425 r = r / 2.;
426 return (r);
427}
428
429////////////////////////////////////////////////////////////////////////////////
430/// This function calculates derivative of peaks shape function (see manual)
431/// according to relative amplitude s.
432/// Function parameters:
433/// - num_of_fitted_peaks-number of fitted peaks
434/// - i-channel
435/// - parameter-array of peaks parameters (amplitudes and positions)
436/// - sigma-sigma of peak
437
438Double_t TSpectrumFit::Ders(Int_t num_of_fitted_peaks, Double_t i,
439 const Double_t *parameter, Double_t sigma)
440{
441 Int_t j;
442 Double_t r, p, r1;
443 r = 0;
444 for (j = 0; j < num_of_fitted_peaks; j++) {
445 p = (i - parameter[2 * j + 1]) / sigma;
446 r1 = Erfc(p);
447 r = r + parameter[2 * j] * r1;
448 }
449 r = r / 2.;
450 return (r);
451}
452
453////////////////////////////////////////////////////////////////////////////////
454/// This function calculates derivative of peaks shape function (see manual)
455/// according to slope b.
456/// Function parameters:
457/// - num_of_fitted_peaks-number of fitted peaks
458/// - i-channel
459/// - parameter-array of peaks parameters (amplitudes and positions)
460/// - sigma-sigma of peak
461/// - t-relative amplitude
462/// - b-slope
463
464Double_t TSpectrumFit::Derb(Int_t num_of_fitted_peaks, Double_t i,
465 const Double_t *parameter, Double_t sigma, Double_t t,
466 Double_t b)
467{
468 Int_t j;
469 Double_t r, p, r1, c, e;
470 r = 0;
471 for (j = 0; j < num_of_fitted_peaks && t != 0; j++) {
472 p = (i - parameter[2 * j + 1]) / sigma;
473 c = p + 1. / (2. * b);
474 e = p / b;
475 r1 = p * Erfc(c);
476 r1 = r1 + Derfc(c) / 2.;
477 if (e > 700)
478 e = 700;
479 if (e < -700)
480 r1 = 0;
481
482 else
483 r1 = r1 * exp(e);
484 r = r + parameter[2 * j] * r1;
485 }
486 r = -r * t / (2. * b * b);
487 return (r);
488}
489
490////////////////////////////////////////////////////////////////////////////////
491/// Derivative of background according to a1
492
493Double_t TSpectrumFit::Dera1(Double_t i)
494{
495 return (i);
496}
497
498////////////////////////////////////////////////////////////////////////////////
499/// Derivative of background according to a2
500
501Double_t TSpectrumFit::Dera2(Double_t i)
502{
503 return (i * i);
504}
505
506////////////////////////////////////////////////////////////////////////////////
507/// This function calculates peaks shape function (see manual)
508/// Function parameters:
509/// - num_of_fitted_peaks-number of fitted peaks
510/// - i-channel
511/// - parameter-array of peaks parameters (amplitudes and positions)
512/// - sigma-sigma of peak
513/// - t, s-relative amplitudes
514/// - b-slope
515/// - a0, a1, a2- background coefficients
516
517Double_t TSpectrumFit::Shape(Int_t num_of_fitted_peaks, Double_t i,
518 const Double_t *parameter, Double_t sigma, Double_t t,
519 Double_t s, Double_t b, Double_t a0, Double_t a1,
520 Double_t a2)
521{
522 Int_t j;
523 Double_t r, p, r1, r2, r3, c, e;
524 r = 0;
525 for (j = 0; j < num_of_fitted_peaks; j++) {
526 if (sigma > 0.0001)
527 p = (i - parameter[2 * j + 1]) / sigma;
528
529 else {
530 if (i == parameter[2 * j + 1])
531 p = 0;
532
533 else
534 p = 10;
535 }
536 r1 = 0;
537 if (TMath::Abs(p) < 3) {
538 if ((p * p) < 700)
539 r1 = exp(-p * p);
540
541 else {
542 r1 = 0;
543 }
544 }
545 r2 = 0;
546 if (t != 0) {
547 c = p + 1. / (2. * b);
548 e = p / b;
549 if (e > 700)
550 e = 700;
551 r2 = t * exp(e) * Erfc(c) / 2.;
552 }
553 r3 = 0;
554 if (s != 0)
555 r3 = s * Erfc(p) / 2.;
556 r = r + parameter[2 * j] * (r1 + r2 + r3);
557 }
558 r = r + a0 + a1 * i + a2 * i * i;
559 return (r);
560}
561
562////////////////////////////////////////////////////////////////////////////////
563/// This function calculates area of a peak
564/// Function parameters:
565/// - a-amplitude of the peak
566/// - sigma-sigma of peak
567/// - t-relative amplitude
568/// - b-slope
569
570Double_t TSpectrumFit::Area(Double_t a, Double_t sigma, Double_t t, Double_t b)
571{
572 Double_t odm_pi = 1.7724538, r = 0;
573 if (b != 0)
574 r = 0.5 / b;
575 r = (-1.) * r * r;
576 if (TMath::Abs(r) < 700)
577 r = a * sigma * (odm_pi + t * b * exp(r));
578
579 else {
580 r = a * sigma * odm_pi;
581 }
582 return (r);
583}
584
585////////////////////////////////////////////////////////////////////////////////
586/// This function calculates derivative of the area of peak
587/// according to its amplitude.
588/// Function parameters:
589/// - sigma-sigma of peak
590/// - t-relative amplitudes
591/// - b-slope
592
593Double_t TSpectrumFit::Derpa(Double_t sigma, Double_t t, Double_t b)
594{
595 Double_t odm_pi = 1.7724538, r;
596 r = 0.5 / b;
597 r = (-1.) * r * r;
598 if (TMath::Abs(r) < 700)
599 r = sigma * (odm_pi + t * b * exp(r));
600
601 else {
602 r = sigma * odm_pi;
603 }
604 return (r);
605}
606
607////////////////////////////////////////////////////////////////////////////////
608/// This function calculates derivative of the area of peak
609/// according to sigma of peaks.
610/// Function parameters:
611/// - a-amplitude of peak
612/// - t-relative amplitudes
613/// - b-slope
614
615Double_t TSpectrumFit::Derpsigma(Double_t a, Double_t t, Double_t b)
616{
617 Double_t odm_pi = 1.7724538, r;
618 r = 0.5 / b;
619 r = (-1.) * r * r;
620 if (TMath::Abs(r) < 700)
621 r = a * (odm_pi + t * b * exp(r));
622
623 else {
624 r = a * odm_pi;
625 }
626 return (r);
627}
628
629////////////////////////////////////////////////////////////////////////////////
630/// This function calculates derivative of the area of peak
631/// according to t parameter.
632/// Function parameters:
633/// - sigma-sigma of peak
634/// - t-relative amplitudes
635/// - b-slope
636
637Double_t TSpectrumFit::Derpt(Double_t a, Double_t sigma, Double_t b)
638{
639 Double_t r;
640 r = 0.5 / b;
641 r = (-1.) * r * r;
642 if (TMath::Abs(r) < 700)
643 r = a * sigma * b * exp(r);
644
645 else {
646 r = 0;
647 }
648 return (r);
649}
650
651////////////////////////////////////////////////////////////////////////////////
652/// This function calculates derivative of the area of peak
653/// according to b parameter.
654/// Function parameters:
655/// - sigma-sigma of peak
656/// - t-relative amplitudes
657/// - b-slope
658
659Double_t TSpectrumFit::Derpb(Double_t a, Double_t sigma, Double_t t, Double_t b)
660{
661 Double_t r;
662 r = (-1) * 0.25 / (b * b);
663 if (TMath::Abs(r) < 700)
664 r = a * sigma * t * exp(r) * (1 - 2 * r);
665
666 else {
667 r = 0;
668 }
669 return (r);
670}
671
672////////////////////////////////////////////////////////////////////////////////
673/// Power function
674
675Double_t TSpectrumFit::Ourpowl(Double_t a, Int_t pw)
676{
677 Double_t c;
678 Double_t a2 = a*a;
679 c = 1;
680 if (pw > 0) c *= a2;
681 if (pw > 2) c *= a2;
682 if (pw > 4) c *= a2;
683 if (pw > 6) c *= a2;
684 if (pw > 8) c *= a2;
685 if (pw > 10) c *= a2;
686 if (pw > 12) c *= a2;
687 return (c);
688}
689
690////////////////////////////////////////////////////////////////////////////////
691/// This function fits the source spectrum. The calling program should
692/// fill in input parameters of the TSpectrumFit class
693/// The fitted parameters are written into
694/// TSpectrumFit class output parameters and fitted data are written into
695/// source spectrum.
696///
697/// Function parameters:
698/// - source-pointer to the vector of source spectrum
699///
700/// ### Fitting
701///
702/// Goal: to estimate simultaneously peak shape parameters in spectra with large
703/// number of peaks
704///
705/// - peaks can be fitted separately, each peak (or multiplets) in a region or
706/// together all peaks in a spectrum. To fit separately each peak one needs to
707/// determine the fitted region. However it can happen that the regions of
708/// neighbouring peaks are overlapping. Then the results of fitting are very poor.
709/// On the other hand, when fitting together all peaks found in a spectrum, one
710/// needs to have a method that is stable (converges) and fast enough to carry out
711/// fitting in reasonable time
712///
713/// - we have implemented the non-symmetrical semi-empirical peak shape function [1]
714///
715/// - it contains the symmetrical Gaussian as well as non-symmetrical terms.
716///
717/// \image html spectrumfit_awmi_image001.gif
718///
719/// where T and S are relative amplitudes and B is slope.
720///
721/// - algorithm without matrix inversion (AWMI) allows fitting tens, hundreds
722/// of peaks simultaneously that represent sometimes thousands of parameters [2], [5].
723///
724/// #### References:
725///
726/// [1] Phillps G.W., Marlow K.W., NIM 137 (1976) 525.
727///
728/// [2] I. A. Slavic: Nonlinear least-squares fitting without matrix inversion
729/// applied to complex Gaussian spectra analysis. NIM 134 (1976) 285-289.
730///
731/// [3] T. Awaya: A new method for curve fitting to the data with low statistics
732/// not using chi-square method. NIM 165 (1979) 317-323.
733///
734/// [4] T. Hauschild, M. Jentschel: Comparison of maximum likelihood estimation
735/// and chi-square statistics applied to counting experiments. NIM A 457 (2001)
736/// 384-401.
737///
738/// [5] M. Morhac, J. Kliman, M. Jandel, L. Krupa, V. Matouoek: Study of fitting
739/// algorithms applied to simultaneous analysis of large number of peaks in -ray
740/// spectra. Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003
741///
742/// ### Example - script FitAwmi.c:
743///
744/// \image html spectrumfit_awmi_image002.jpg Fig. 1 Original spectrum (black line) and fitted spectrum using AWMI algorithm (red line) and number of iteration steps = 1000. Positions of fitted peaks are denoted by markers
745///
746/// #### Script:
747///
748/// Example to illustrate fitting function using AWMI algorithm.
749/// To execute this example, do:
750///
751/// `root > .x FitAwmi.C`
752///
753/// ~~~ {.cpp}
754/// void FitAwmi() {
755/// Double_t a;
756/// Int_t i,nfound=0,bin;
757/// Int_t nbins = 256;
758/// Int_t xmin = 0;
759/// Int_t xmax = nbins;
760/// Double_t * source = new Double_t[nbins];
761/// Double_t * dest = new Double_t[nbins];
762/// TH1F *h = new TH1F("h","Fitting using AWMI algorithm",nbins,xmin,xmax);
763/// TH1F *d = new TH1F("d","",nbins,xmin,xmax);
764/// TFile *f = new TFile("TSpectrum.root");
765/// h=(TH1F*) f->Get("fit;1");
766/// for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);
767/// TCanvas *Fit1 = gROOT->GetListOfCanvases()->FindObject("Fit1");
768/// if (!Fit1) Fit1 = new TCanvas("Fit1","Fit1",10,10,1000,700);
769/// h->Draw("L");
770/// TSpectrum *s = new TSpectrum();
771/// //searching for candidate peaks positions
772/// nfound = s->SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);
773/// Bool_t *FixPos =new Bool_t[nfound];
774/// Bool_t *FixAmp = new Bool_t[nfound];
775/// for(i = 0; i< nfound ; i++){
776/// FixPos[i] = kFALSE;
777/// FixAmp[i] = kFALSE;
778/// }
779/// //filling in the
780/// initial estimates of the input parameters
781/// Double_t *PosX = new Double_t[nfound];
782/// Double_t *PosY = new Double_t[nfound];
783/// PosX = s->GetPositionX();
784/// for (i = 0; i < nfound; i++) {
785/// a=PosX[i];
786/// bin = 1 + Int_t(a + 0.5);
787/// PosY[i] = h->GetBinContent(bin);
788/// }
789/// TSpectrumFit *pfit=new TSpectrumFit(nfound);
790/// pfit->SetFitParameters(xmin, xmax-1, 1000, 0.1, pfit->kFitOptimChiCounts,
791/// pfit->kFitAlphaHalving, pfit->kFitPower2,
792/// pfit->kFitTaylorOrderFirst);
793/// pfit->SetPeakParameters(2, kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);
794/// pfit->FitAwmi(source);
795/// Double_t *CalcPositions = new Double_t[nfound];
796/// Double_t *CalcAmplitudes = new Double_t[nfound];
797/// CalcPositions=pfit->GetPositions();
798/// CalcAmplitudes=pfit->GetAmplitudes();
799/// for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
800/// d->SetLineColor(kRed);
801/// d->Draw("SAME L");
802/// for (i = 0; i < nfound; i++) {
803/// a=CalcPositions[i];
804/// bin = 1 + Int_t(a + 0.5);
805/// PosX[i] = d->GetBinCenter(bin);
806/// PosY[i] = d->GetBinContent(bin);
807/// }
808/// TPolyMarker * pm = (TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");
809/// if (pm) {
810/// h->GetListOfFunctions()->Remove(pm);
811/// delete pm;
812/// }
813/// pm = new TPolyMarker(nfound, PosX, PosY);
814/// h->GetListOfFunctions()->Add(pm);
815/// pm->SetMarkerStyle(23);
816/// pm->SetMarkerColor(kRed);
817/// pm->SetMarkerSize(1);
818/// }
819/// ~~~
820
821void TSpectrumFit::FitAwmi(Double_t *source)
822{
823 Int_t i, j, k, shift =
824 2 * fNPeaks + 7, peak_vel, rozmer, iter, pw, regul_cycle,
825 flag;
826 Double_t a, b, c, d = 0, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
827 0, pi, pmin = 0, chi_cel = 0, chi_er;
828 Double_t *working_space = new Double_t[5 * (2 * fNPeaks + 7)];
829 for (i = 0, j = 0; i < fNPeaks; i++) {
830 working_space[2 * i] = fAmpInit[i]; //vector parameter
831 if (fFixAmp[i] == false) {
832 working_space[shift + j] = fAmpInit[i]; //vector xk
833 j += 1;
834 }
835 working_space[2 * i + 1] = fPositionInit[i]; //vector parameter
836 if (fFixPosition[i] == false) {
837 working_space[shift + j] = fPositionInit[i]; //vector xk
838 j += 1;
839 }
840 }
841 peak_vel = 2 * i;
842 working_space[2 * i] = fSigmaInit; //vector parameter
843 if (fFixSigma == false) {
844 working_space[shift + j] = fSigmaInit; //vector xk
845 j += 1;
846 }
847 working_space[2 * i + 1] = fTInit; //vector parameter
848 if (fFixT == false) {
849 working_space[shift + j] = fTInit; //vector xk
850 j += 1;
851 }
852 working_space[2 * i + 2] = fBInit; //vector parameter
853 if (fFixB == false) {
854 working_space[shift + j] = fBInit; //vector xk
855 j += 1;
856 }
857 working_space[2 * i + 3] = fSInit; //vector parameter
858 if (fFixS == false) {
859 working_space[shift + j] = fSInit; //vector xk
860 j += 1;
861 }
862 working_space[2 * i + 4] = fA0Init; //vector parameter
863 if (fFixA0 == false) {
864 working_space[shift + j] = fA0Init; //vector xk
865 j += 1;
866 }
867 working_space[2 * i + 5] = fA1Init; //vector parameter
868 if (fFixA1 == false) {
869 working_space[shift + j] = fA1Init; //vector xk
870 j += 1;
871 }
872 working_space[2 * i + 6] = fA2Init; //vector parameter
873 if (fFixA2 == false) {
874 working_space[shift + j] = fA2Init; //vector xk
875 j += 1;
876 }
877 rozmer = j;
878 if (rozmer == 0){
879 delete [] working_space;
880 Error ("FitAwmi","All parameters are fixed");
881 return;
882 }
883 if (rozmer >= fXmax - fXmin + 1){
884 delete [] working_space;
885 Error ("FitAwmi","Number of fitted parameters is larger than # of fitted points");
886 return;
887 }
888 for (iter = 0; iter < fNumberIterations; iter++) {
889 for (j = 0; j < rozmer; j++) {
890 working_space[2 * shift + j] = 0, working_space[3 * shift + j] = 0; //der,temp
891 }
892
893 //filling vectors
894 alpha = fAlpha;
895 chi_opt = 0, pw = fPower - 2;
896 for (i = fXmin; i <= fXmax; i++) {
897 yw = source[i];
898 ywm = yw;
899 f = Shape(fNPeaks, (Double_t) i, working_space,
900 working_space[peak_vel], working_space[peak_vel + 1],
901 working_space[peak_vel + 3],
902 working_space[peak_vel + 2],
903 working_space[peak_vel + 4],
904 working_space[peak_vel + 5],
905 working_space[peak_vel + 6]);
906 if (fStatisticType == kFitOptimMaxLikelihood) {
907 if (f > 0.00001)
908 chi_opt += yw * TMath::Log(f) - f;
909 }
910
911 else {
912 if (ywm != 0)
913 chi_opt += (yw - f) * (yw - f) / ywm;
914 }
915 if (fStatisticType == kFitOptimChiFuncValues) {
916 ywm = f;
917 if (f < 0.00001)
918 ywm = 0.00001;
919 }
920
921 else if (fStatisticType == kFitOptimMaxLikelihood) {
922 ywm = f;
923 if (f < 0.001)
924 ywm = 0.001;
925 }
926
927 else {
928 if (ywm == 0)
929 ywm = 1;
930 }
931
932 //calculation of gradient vector
933 for (j = 0, k = 0; j < fNPeaks; j++) {
934 if (fFixAmp[j] == false) {
935 a = Deramp((Double_t) i, working_space[2 * j + 1],
936 working_space[peak_vel],
937 working_space[peak_vel + 1],
938 working_space[peak_vel + 3],
939 working_space[peak_vel + 2]);
940 if (ywm != 0) {
941 c = Ourpowl(a, pw);
942 if (fStatisticType == kFitOptimChiFuncValues) {
943 b = a * (yw * yw - f * f) / (ywm * ywm);
944 working_space[2 * shift + k] += b * c; //der
945 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
946 working_space[3 * shift + k] += b * c; //temp
947 }
948
949 else {
950 b = a * (yw - f) / ywm;
951 working_space[2 * shift + k] += b * c; //der
952 b = a * a / ywm;
953 working_space[3 * shift + k] += b * c; //temp
954 }
955 }
956 k += 1;
957 }
958 if (fFixPosition[j] == false) {
959 a = Deri0((Double_t) i, working_space[2 * j],
960 working_space[2 * j + 1],
961 working_space[peak_vel],
962 working_space[peak_vel + 1],
963 working_space[peak_vel + 3],
964 working_space[peak_vel + 2]);
965 if (fFitTaylor == kFitTaylorOrderSecond)
966 d = Derderi0((Double_t) i, working_space[2 * j],
967 working_space[2 * j + 1],
968 working_space[peak_vel]);
969 if (ywm != 0) {
970 c = Ourpowl(a, pw);
971 if (TMath::Abs(a) > 0.00000001
972 && fFitTaylor == kFitTaylorOrderSecond) {
973 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
974 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
975 d = 0;
976 }
977
978 else
979 d = 0;
980 a = a + d;
981 if (fStatisticType == kFitOptimChiFuncValues) {
982 b = a * (yw * yw - f * f) / (ywm * ywm);
983 working_space[2 * shift + k] += b * c; //der
984 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
985 working_space[3 * shift + k] += b * c; //temp
986 }
987
988 else {
989 b = a * (yw - f) / ywm;
990 working_space[2 * shift + k] += b * c; //der
991 b = a * a / ywm;
992 working_space[3 * shift + k] += b * c; //temp
993 }
994 }
995 k += 1;
996 }
997 }
998 if (fFixSigma == false) {
999 a = Dersigma(fNPeaks, (Double_t) i, working_space,
1000 working_space[peak_vel],
1001 working_space[peak_vel + 1],
1002 working_space[peak_vel + 3],
1003 working_space[peak_vel + 2]);
1004 if (fFitTaylor == kFitTaylorOrderSecond)
1005 d = Derdersigma(fNPeaks, (Double_t) i,
1006 working_space, working_space[peak_vel]);
1007 if (ywm != 0) {
1008 c = Ourpowl(a, pw);
1009 if (TMath::Abs(a) > 0.00000001
1010 && fFitTaylor == kFitTaylorOrderSecond) {
1011 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
1012 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
1013 d = 0;
1014 }
1015
1016 else
1017 d = 0;
1018 a = a + d;
1019 if (fStatisticType == kFitOptimChiFuncValues) {
1020 b = a * (yw * yw - f * f) / (ywm * ywm);
1021 working_space[2 * shift + k] += b * c; //der
1022 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1023 working_space[3 * shift + k] += b * c; //temp
1024 }
1025
1026 else {
1027 b = a * (yw - f) / ywm;
1028 working_space[2 * shift + k] += b * c; //der
1029 b = a * a / ywm;
1030 working_space[3 * shift + k] += b * c; //temp
1031 }
1032 }
1033 k += 1;
1034 }
1035 if (fFixT == false) {
1036 a = Dert(fNPeaks, (Double_t) i, working_space,
1037 working_space[peak_vel],
1038 working_space[peak_vel + 2]);
1039 if (ywm != 0) {
1040 c = Ourpowl(a, pw);
1041 if (fStatisticType == kFitOptimChiFuncValues) {
1042 b = a * (yw * yw - f * f) / (ywm * ywm);
1043 working_space[2 * shift + k] += b * c; //der
1044 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1045 working_space[3 * shift + k] += b * c; //temp
1046 }
1047
1048 else {
1049 b = a * (yw - f) / ywm;
1050 working_space[2 * shift + k] += b * c; //der
1051 b = a * a / ywm;
1052 working_space[3 * shift + k] += b * c; //temp
1053 }
1054 }
1055 k += 1;
1056 }
1057 if (fFixB == false) {
1058 a = Derb(fNPeaks, (Double_t) i, working_space,
1059 working_space[peak_vel], working_space[peak_vel + 1],
1060 working_space[peak_vel + 2]);
1061 if (ywm != 0) {
1062 c = Ourpowl(a, pw);
1063 if (fStatisticType == kFitOptimChiFuncValues) {
1064 b = a * (yw * yw - f * f) / (ywm * ywm);
1065 working_space[2 * shift + k] += b * c; //der
1066 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1067 working_space[3 * shift + k] += b * c; //temp
1068 }
1069
1070 else {
1071 b = a * (yw - f) / ywm;
1072 working_space[2 * shift + k] += b * c; //der
1073 b = a * a / ywm;
1074 working_space[3 * shift + k] += b * c; //temp
1075 }
1076 }
1077 k += 1;
1078 }
1079 if (fFixS == false) {
1080 a = Ders(fNPeaks, (Double_t) i, working_space,
1081 working_space[peak_vel]);
1082 if (ywm != 0) {
1083 c = Ourpowl(a, pw);
1084 if (fStatisticType == kFitOptimChiFuncValues) {
1085 b = a * (yw * yw - f * f) / (ywm * ywm);
1086 working_space[2 * shift + k] += b * c; //der
1087 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1088 working_space[3 * shift + k] += b * c; //temp
1089 }
1090
1091 else {
1092 b = a * (yw - f) / ywm;
1093 working_space[2 * shift + k] += b * c; //der
1094 b = a * a / ywm;
1095 working_space[3 * shift + k] += b * c; //temp
1096 }
1097 }
1098 k += 1;
1099 }
1100 if (fFixA0 == false) {
1101 a = 1.;
1102 if (ywm != 0) {
1103 c = Ourpowl(a, pw);
1104 if (fStatisticType == kFitOptimChiFuncValues) {
1105 b = a * (yw * yw - f * f) / (ywm * ywm);
1106 working_space[2 * shift + k] += b * c; //der
1107 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1108 working_space[3 * shift + k] += b * c; //temp
1109 }
1110
1111 else {
1112 b = a * (yw - f) / ywm;
1113 working_space[2 * shift + k] += b * c; //der
1114 b = a * a / ywm;
1115 working_space[3 * shift + k] += b * c; //temp
1116 }
1117 }
1118 k += 1;
1119 }
1120 if (fFixA1 == false) {
1121 a = Dera1((Double_t) i);
1122 if (ywm != 0) {
1123 c = Ourpowl(a, pw);
1124 if (fStatisticType == kFitOptimChiFuncValues) {
1125 b = a * (yw * yw - f * f) / (ywm * ywm);
1126 working_space[2 * shift + k] += b * c; //der
1127 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1128 working_space[3 * shift + k] += b * c; //temp
1129 }
1130
1131 else {
1132 b = a * (yw - f) / ywm;
1133 working_space[2 * shift + k] += b * c; //der
1134 b = a * a / ywm;
1135 working_space[3 * shift + k] += b * c; //temp
1136 }
1137 }
1138 k += 1;
1139 }
1140 if (fFixA2 == false) {
1141 a = Dera2((Double_t) i);
1142 if (ywm != 0) {
1143 c = Ourpowl(a, pw);
1144 if (fStatisticType == kFitOptimChiFuncValues) {
1145 b = a * (yw * yw - f * f) / (ywm * ywm);
1146 working_space[2 * shift + k] += b * c; //der
1147 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1148 working_space[3 * shift + k] += b * c; //temp
1149 }
1150
1151 else {
1152 b = a * (yw - f) / ywm;
1153 working_space[2 * shift + k] += b * c; //der
1154 b = a * a / ywm;
1155 working_space[3 * shift + k] += b * c; //temp
1156 }
1157 }
1158 k += 1;
1159 }
1160 }
1161 for (j = 0; j < rozmer; j++) {
1162 if (TMath::Abs(working_space[3 * shift + j]) > 0.000001)
1163 working_space[2 * shift + j] = working_space[2 * shift + j] / TMath::Abs(working_space[3 * shift + j]); //der[j]=der[j]/temp[j]
1164 else
1165 working_space[2 * shift + j] = 0; //der[j]
1166 }
1167
1168 //calculate chi_opt
1169 chi2 = chi_opt;
1170 chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
1171
1172 //calculate new parameters
1173 regul_cycle = 0;
1174 for (j = 0; j < rozmer; j++) {
1175 working_space[4 * shift + j] = working_space[shift + j]; //temp_xk[j]=xk[j]
1176 }
1177
1178 do {
1179 if (fAlphaOptim == kFitAlphaOptimal) {
1180 if (fStatisticType != kFitOptimMaxLikelihood)
1181 chi_min = 10000 * chi2;
1182
1183 else
1184 chi_min = 0.1 * chi2;
1185 flag = 0;
1186 for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
1187 for (j = 0; j < rozmer; j++) {
1188 working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
1189 }
1190 for (i = 0, j = 0; i < fNPeaks; i++) {
1191 if (fFixAmp[i] == false) {
1192 if (working_space[shift + j] < 0) //xk[j]
1193 working_space[shift + j] = 0; //xk[j]
1194 working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
1195 j += 1;
1196 }
1197 if (fFixPosition[i] == false) {
1198 if (working_space[shift + j] < fXmin) //xk[j]
1199 working_space[shift + j] = fXmin; //xk[j]
1200 if (working_space[shift + j] > fXmax) //xk[j]
1201 working_space[shift + j] = fXmax; //xk[j]
1202 working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
1203 j += 1;
1204 }
1205 }
1206 if (fFixSigma == false) {
1207 if (working_space[shift + j] < 0.001) { //xk[j]
1208 working_space[shift + j] = 0.001; //xk[j]
1209 }
1210 working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
1211 j += 1;
1212 }
1213 if (fFixT == false) {
1214 working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
1215 j += 1;
1216 }
1217 if (fFixB == false) {
1218 if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
1219 if (working_space[shift + j] < 0) //xk[j]
1220 working_space[shift + j] = -0.001; //xk[j]
1221 else
1222 working_space[shift + j] = 0.001; //xk[j]
1223 }
1224 working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
1225 j += 1;
1226 }
1227 if (fFixS == false) {
1228 working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
1229 j += 1;
1230 }
1231 if (fFixA0 == false) {
1232 working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
1233 j += 1;
1234 }
1235 if (fFixA1 == false) {
1236 working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
1237 j += 1;
1238 }
1239 if (fFixA2 == false) {
1240 working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
1241 j += 1;
1242 }
1243 chi2 = 0;
1244 for (i = fXmin; i <= fXmax; i++) {
1245 yw = source[i];
1246 ywm = yw;
1247 f = Shape(fNPeaks, (Double_t) i, working_space,
1248 working_space[peak_vel],
1249 working_space[peak_vel + 1],
1250 working_space[peak_vel + 3],
1251 working_space[peak_vel + 2],
1252 working_space[peak_vel + 4],
1253 working_space[peak_vel + 5],
1254 working_space[peak_vel + 6]);
1255 if (fStatisticType == kFitOptimChiFuncValues) {
1256 ywm = f;
1257 if (f < 0.00001)
1258 ywm = 0.00001;
1259 }
1260 if (fStatisticType == kFitOptimMaxLikelihood) {
1261 if (f > 0.00001)
1262 chi2 += yw * TMath::Log(f) - f;
1263 }
1264
1265 else {
1266 if (ywm != 0)
1267 chi2 += (yw - f) * (yw - f) / ywm;
1268 }
1269 }
1270 if ((chi2 < chi_min
1271 && fStatisticType != kFitOptimMaxLikelihood)
1272 || (chi2 > chi_min
1273 && fStatisticType == kFitOptimMaxLikelihood)) {
1274 pmin = pi, chi_min = chi2;
1275 }
1276
1277 else
1278 flag = 1;
1279 if (pi == 0.1)
1280 chi_min = chi2;
1281 chi = chi_min;
1282 }
1283 if (pmin != 0.1) {
1284 for (j = 0; j < rozmer; j++) {
1285 working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pmin*alpha*der[j]
1286 }
1287 for (i = 0, j = 0; i < fNPeaks; i++) {
1288 if (fFixAmp[i] == false) {
1289 if (working_space[shift + j] < 0) //xk[j]
1290 working_space[shift + j] = 0; //xk[j]
1291 working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
1292 j += 1;
1293 }
1294 if (fFixPosition[i] == false) {
1295 if (working_space[shift + j] < fXmin) //xk[j]
1296 working_space[shift + j] = fXmin; //xk[j]
1297 if (working_space[shift + j] > fXmax) //xk[j]
1298 working_space[shift + j] = fXmax; //xk[j]
1299 working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
1300 j += 1;
1301 }
1302 }
1303 if (fFixSigma == false) {
1304 if (working_space[shift + j] < 0.001) { //xk[j]
1305 working_space[shift + j] = 0.001; //xk[j]
1306 }
1307 working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
1308 j += 1;
1309 }
1310 if (fFixT == false) {
1311 working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
1312 j += 1;
1313 }
1314 if (fFixB == false) {
1315 if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
1316 if (working_space[shift + j] < 0) //xk[j]
1317 working_space[shift + j] = -0.001; //xk[j]
1318 else
1319 working_space[shift + j] = 0.001; //xk[j]
1320 }
1321 working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
1322 j += 1;
1323 }
1324 if (fFixS == false) {
1325 working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
1326 j += 1;
1327 }
1328 if (fFixA0 == false) {
1329 working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
1330 j += 1;
1331 }
1332 if (fFixA1 == false) {
1333 working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
1334 j += 1;
1335 }
1336 if (fFixA2 == false) {
1337 working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
1338 j += 1;
1339 }
1340 chi = chi_min;
1341 }
1342 }
1343
1344 else {
1345 for (j = 0; j < rozmer; j++) {
1346 working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
1347 }
1348 for (i = 0, j = 0; i < fNPeaks; i++) {
1349 if (fFixAmp[i] == false) {
1350 if (working_space[shift + j] < 0) //xk[j]
1351 working_space[shift + j] = 0; //xk[j]
1352 working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
1353 j += 1;
1354 }
1355 if (fFixPosition[i] == false) {
1356 if (working_space[shift + j] < fXmin) //xk[j]
1357 working_space[shift + j] = fXmin; //xk[j]
1358 if (working_space[shift + j] > fXmax) //xk[j]
1359 working_space[shift + j] = fXmax; //xk[j]
1360 working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
1361 j += 1;
1362 }
1363 }
1364 if (fFixSigma == false) {
1365 if (working_space[shift + j] < 0.001) { //xk[j]
1366 working_space[shift + j] = 0.001; //xk[j]
1367 }
1368 working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
1369 j += 1;
1370 }
1371 if (fFixT == false) {
1372 working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
1373 j += 1;
1374 }
1375 if (fFixB == false) {
1376 if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
1377 if (working_space[shift + j] < 0) //xk[j]
1378 working_space[shift + j] = -0.001; //xk[j]
1379 else
1380 working_space[shift + j] = 0.001; //xk[j]
1381 }
1382 working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
1383 j += 1;
1384 }
1385 if (fFixS == false) {
1386 working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
1387 j += 1;
1388 }
1389 if (fFixA0 == false) {
1390 working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
1391 j += 1;
1392 }
1393 if (fFixA1 == false) {
1394 working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
1395 j += 1;
1396 }
1397 if (fFixA2 == false) {
1398 working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
1399 j += 1;
1400 }
1401 chi = 0;
1402 for (i = fXmin; i <= fXmax; i++) {
1403 yw = source[i];
1404 ywm = yw;
1405 f = Shape(fNPeaks, (Double_t) i, working_space,
1406 working_space[peak_vel],
1407 working_space[peak_vel + 1],
1408 working_space[peak_vel + 3],
1409 working_space[peak_vel + 2],
1410 working_space[peak_vel + 4],
1411 working_space[peak_vel + 5],
1412 working_space[peak_vel + 6]);
1413 if (fStatisticType == kFitOptimChiFuncValues) {
1414 ywm = f;
1415 if (f < 0.00001)
1416 ywm = 0.00001;
1417 }
1418 if (fStatisticType == kFitOptimMaxLikelihood) {
1419 if (f > 0.00001)
1420 chi += yw * TMath::Log(f) - f;
1421 }
1422
1423 else {
1424 if (ywm != 0)
1425 chi += (yw - f) * (yw - f) / ywm;
1426 }
1427 }
1428 }
1429 chi2 = chi;
1430 chi = TMath::Sqrt(TMath::Abs(chi));
1431 if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
1432 alpha = alpha * chi_opt / (2 * chi);
1433
1434 else if (fAlphaOptim == kFitAlphaOptimal)
1435 alpha = alpha / 10.0;
1436 iter += 1;
1437 regul_cycle += 1;
1438 } while (((chi > chi_opt
1439 && fStatisticType != kFitOptimMaxLikelihood)
1440 || (chi < chi_opt
1441 && fStatisticType == kFitOptimMaxLikelihood))
1442 && regul_cycle < kFitNumRegulCycles);
1443 for (j = 0; j < rozmer; j++) {
1444 working_space[4 * shift + j] = 0; //temp_xk[j]
1445 working_space[2 * shift + j] = 0; //der[j]
1446 }
1447 for (i = fXmin, chi_cel = 0; i <= fXmax; i++) {
1448 yw = source[i];
1449 if (yw == 0)
1450 yw = 1;
1451 f = Shape(fNPeaks, (Double_t) i, working_space,
1452 working_space[peak_vel], working_space[peak_vel + 1],
1453 working_space[peak_vel + 3],
1454 working_space[peak_vel + 2],
1455 working_space[peak_vel + 4],
1456 working_space[peak_vel + 5],
1457 working_space[peak_vel + 6]);
1458 chi_opt = (yw - f) * (yw - f) / yw;
1459 chi_cel += (yw - f) * (yw - f) / yw;
1460
1461 //calculate gradient vector
1462 for (j = 0, k = 0; j < fNPeaks; j++) {
1463 if (fFixAmp[j] == false) {
1464 a = Deramp((Double_t) i, working_space[2 * j + 1],
1465 working_space[peak_vel],
1466 working_space[peak_vel + 1],
1467 working_space[peak_vel + 3],
1468 working_space[peak_vel + 2]);
1469 if (yw != 0) {
1470 c = Ourpowl(a, pw);
1471 working_space[2 * shift + k] += chi_opt * c; //der[k]
1472 b = a * a / yw;
1473 working_space[4 * shift + k] += b * c; //temp_xk[k]
1474 }
1475 k += 1;
1476 }
1477 if (fFixPosition[j] == false) {
1478 a = Deri0((Double_t) i, working_space[2 * j],
1479 working_space[2 * j + 1],
1480 working_space[peak_vel],
1481 working_space[peak_vel + 1],
1482 working_space[peak_vel + 3],
1483 working_space[peak_vel + 2]);
1484 if (yw != 0) {
1485 c = Ourpowl(a, pw);
1486 working_space[2 * shift + k] += chi_opt * c; //der[k]
1487 b = a * a / yw;
1488 working_space[4 * shift + k] += b * c; //temp_xk[k]
1489 }
1490 k += 1;
1491 }
1492 }
1493 if (fFixSigma == false) {
1494 a = Dersigma(fNPeaks, (Double_t) i, working_space,
1495 working_space[peak_vel],
1496 working_space[peak_vel + 1],
1497 working_space[peak_vel + 3],
1498 working_space[peak_vel + 2]);
1499 if (yw != 0) {
1500 c = Ourpowl(a, pw);
1501 working_space[2 * shift + k] += chi_opt * c; //der[k]
1502 b = a * a / yw;
1503 working_space[4 * shift + k] += b * c; //temp_xk[k]
1504 }
1505 k += 1;
1506 }
1507 if (fFixT == false) {
1508 a = Dert(fNPeaks, (Double_t) i, working_space,
1509 working_space[peak_vel],
1510 working_space[peak_vel + 2]);
1511 if (yw != 0) {
1512 c = Ourpowl(a, pw);
1513 working_space[2 * shift + k] += chi_opt * c; //der[k]
1514 b = a * a / yw;
1515 working_space[4 * shift + k] += b * c; //temp_xk[k]
1516 }
1517 k += 1;
1518 }
1519 if (fFixB == false) {
1520 a = Derb(fNPeaks, (Double_t) i, working_space,
1521 working_space[peak_vel], working_space[peak_vel + 1],
1522 working_space[peak_vel + 2]);
1523 if (yw != 0) {
1524 c = Ourpowl(a, pw);
1525 working_space[2 * shift + k] += chi_opt * c; //der[k]
1526 b = a * a / yw;
1527 working_space[4 * shift + k] += b * c; //temp_xk[k]
1528 }
1529 k += 1;
1530 }
1531 if (fFixS == false) {
1532 a = Ders(fNPeaks, (Double_t) i, working_space,
1533 working_space[peak_vel]);
1534 if (yw != 0) {
1535 c = Ourpowl(a, pw);
1536 working_space[2 * shift + k] += chi_opt * c; //der[k]
1537 b = a * a / yw;
1538 working_space[4 * shift + k] += b * c; //tem_xk[k]
1539 }
1540 k += 1;
1541 }
1542 if (fFixA0 == false) {
1543 a = 1.0;
1544 if (yw != 0) {
1545 c = Ourpowl(a, pw);
1546 working_space[2 * shift + k] += chi_opt * c; //der[k]
1547 b = a * a / yw;
1548 working_space[4 * shift + k] += b * c; //temp_xk[k]
1549 }
1550 k += 1;
1551 }
1552 if (fFixA1 == false) {
1553 a = Dera1((Double_t) i);
1554 if (yw != 0) {
1555 c = Ourpowl(a, pw);
1556 working_space[2 * shift + k] += chi_opt * c; //der[k]
1557 b = a * a / yw;
1558 working_space[4 * shift + k] += b * c; //temp_xk[k]
1559 }
1560 k += 1;
1561 }
1562 if (fFixA2 == false) {
1563 a = Dera2((Double_t) i);
1564 if (yw != 0) {
1565 c = Ourpowl(a, pw);
1566 working_space[2 * shift + k] += chi_opt * c; //der[k]
1567 b = a * a / yw;
1568 working_space[4 * shift + k] += b * c; //temp_xk[k]
1569 }
1570 k += 1;
1571 }
1572 }
1573 }
1574 b = fXmax - fXmin + 1 - rozmer;
1575 chi_er = chi_cel / b;
1576 for (i = 0, j = 0; i < fNPeaks; i++) {
1577 fArea[i] =
1578 Area(working_space[2 * i], working_space[peak_vel],
1579 working_space[peak_vel + 1], working_space[peak_vel + 2]);
1580 if (fFixAmp[i] == false) {
1581 fAmpCalc[i] = working_space[shift + j]; //xk[j]
1582 if (working_space[3 * shift + j] != 0)
1583 fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1584 if (fArea[i] > 0) {
1585 a = Derpa(working_space[peak_vel],
1586 working_space[peak_vel + 1],
1587 working_space[peak_vel + 2]);
1588 b = working_space[4 * shift + j]; //temp_xk[j]
1589 if (b == 0)
1590 b = 1;
1591
1592 else
1593 b = 1 / b;
1594 fAreaErr[i] = TMath::Sqrt(TMath::Abs(a * a * b * chi_er));
1595 }
1596
1597 else
1598 fAreaErr[i] = 0;
1599 j += 1;
1600 }
1601
1602 else {
1603 fAmpCalc[i] = fAmpInit[i];
1604 fAmpErr[i] = 0;
1605 fAreaErr[i] = 0;
1606 }
1607 if (fFixPosition[i] == false) {
1608 fPositionCalc[i] = working_space[shift + j]; //xk[j]
1609 if (working_space[3 * shift + j] != 0) //temp[j]
1610 fPositionErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1611 j += 1;
1612 }
1613
1614 else {
1615 fPositionCalc[i] = fPositionInit[i];
1616 fPositionErr[i] = 0;
1617 }
1618 }
1619 if (fFixSigma == false) {
1620 fSigmaCalc = working_space[shift + j]; //xk[j]
1621 if (working_space[3 * shift + j] != 0) //temp[j]
1622 fSigmaErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1623 j += 1;
1624 }
1625
1626 else {
1627 fSigmaCalc = fSigmaInit;
1628 fSigmaErr = 0;
1629 }
1630 if (fFixT == false) {
1631 fTCalc = working_space[shift + j]; //xk[j]
1632 if (working_space[3 * shift + j] != 0) //temp[j]
1633 fTErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1634 j += 1;
1635 }
1636
1637 else {
1638 fTCalc = fTInit;
1639 fTErr = 0;
1640 }
1641 if (fFixB == false) {
1642 fBCalc = working_space[shift + j]; //xk[j]
1643 if (working_space[3 * shift + j] != 0) //temp[j]
1644 fBErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1645 j += 1;
1646 }
1647
1648 else {
1649 fBCalc = fBInit;
1650 fBErr = 0;
1651 }
1652 if (fFixS == false) {
1653 fSCalc = working_space[shift + j]; //xk[j]
1654 if (working_space[3 * shift + j] != 0) //temp[j]
1655 fSErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1656 j += 1;
1657 }
1658
1659 else {
1660 fSCalc = fSInit;
1661 fSErr = 0;
1662 }
1663 if (fFixA0 == false) {
1664 fA0Calc = working_space[shift + j]; //xk[j]
1665 if (working_space[3 * shift + j] != 0) //temp[j]
1666 fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1667 j += 1;
1668 }
1669
1670 else {
1671 fA0Calc = fA0Init;
1672 fA0Err = 0;
1673 }
1674 if (fFixA1 == false) {
1675 fA1Calc = working_space[shift + j]; //xk[j]
1676 if (working_space[3 * shift + j] != 0) //temp[j]
1677 fA1Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1678 j += 1;
1679 }
1680
1681 else {
1682 fA1Calc = fA1Init;
1683 fA1Err = 0;
1684 }
1685 if (fFixA2 == false) {
1686 fA2Calc = working_space[shift + j]; //xk[j]
1687 if (working_space[3 * shift + j] != 0) //temp[j]
1688 fA2Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
1689 j += 1;
1690 }
1691
1692 else {
1693 fA2Calc = fA2Init;
1694 fA2Err = 0;
1695 }
1696 b = fXmax - fXmin + 1 - rozmer;
1697 fChi = chi_cel / b;
1698 for (i = fXmin; i <= fXmax; i++) {
1699 f = Shape(fNPeaks, (Double_t) i, working_space,
1700 working_space[peak_vel], working_space[peak_vel + 1],
1701 working_space[peak_vel + 3], working_space[peak_vel + 2],
1702 working_space[peak_vel + 4], working_space[peak_vel + 5],
1703 working_space[peak_vel + 6]);
1704 source[i] = f;
1705 }
1706 delete[]working_space;
1707 return;
1708}
1709
1710////////////////////////////////////////////////////////////////////////////////
1711/// This function calculates solution of the system of linear equations.
1712/// The matrix a should have a dimension size*(size+4)
1713/// The calling function should fill in the matrix, the column size should
1714/// contain vector y (right side of the system of equations). The result is
1715/// placed into size+1 column of the matrix.
1716/// according to sigma of peaks.
1717///
1718/// Function parameters:
1719/// - a-matrix with dimension size*(size+4)
1720/// - size-number of rows of the matrix
1721
1722void TSpectrumFit::StiefelInversion(Double_t **a, Int_t size)
1723{
1724 Int_t i, j, k = 0;
1725 Double_t sk = 0, b, lambdak, normk, normk_old = 0;
1726
1727 do {
1728 normk = 0;
1729
1730 //calculation of rk and norm
1731 for (i = 0; i < size; i++) {
1732 a[i][size + 2] = -a[i][size]; //rk=-C
1733 for (j = 0; j < size; j++) {
1734 a[i][size + 2] += a[i][j] * a[j][size + 1]; //A*xk-C
1735 }
1736 normk += a[i][size + 2] * a[i][size + 2]; //calculation normk
1737 }
1738
1739 //calculation of sk
1740 if (k != 0) {
1741 sk = normk / normk_old;
1742 }
1743
1744 //calculation of uk
1745 for (i = 0; i < size; i++) {
1746 a[i][size + 3] = -a[i][size + 2] + sk * a[i][size + 3]; //uk=-rk+sk*uk-1
1747 }
1748
1749 //calculation of lambdak
1750 lambdak = 0;
1751 for (i = 0; i < size; i++) {
1752 for (j = 0, b = 0; j < size; j++) {
1753 b += a[i][j] * a[j][size + 3]; //A*uk
1754 }
1755 lambdak += b * a[i][size + 3];
1756 }
1757 if (TMath::Abs(lambdak) > 1e-50) //computer zero
1758 lambdak = normk / lambdak;
1759
1760 else
1761 lambdak = 0;
1762 for (i = 0; i < size; i++)
1763 a[i][size + 1] += lambdak * a[i][size + 3]; //xk+1=xk+lambdak*uk
1764 normk_old = normk;
1765 k += 1;
1766 } while (k < size && TMath::Abs(normk) > 1e-50); //computer zero
1767 return;
1768}
1769
1770////////////////////////////////////////////////////////////////////////////////
1771/// This function fits the source spectrum. The calling program should
1772/// fill in input parameters
1773/// The fitted parameters are written into
1774/// output parameters and fitted data are written into
1775/// source spectrum.
1776///
1777/// Function parameters:
1778/// - source-pointer to the vector of source spectrum
1779///
1780/// ### Example - script FitStiefel.c:
1781///
1782/// \image html spectrumfit_stiefel_image001.jpg Fig. 2 Original spectrum (black line) and fitted spectrum using Stiefel-Hestens method (red line) and number of iteration steps = 100. Positions of fitted peaks are denoted by markers
1783///
1784/// #### Script:
1785///
1786/// Example to illustrate fitting function using Stiefel-Hestens method.
1787/// To execute this example, do:
1788///
1789/// root > .x FitStiefel.C
1790///
1791/// ~~~ {.cpp}
1792/// void FitStiefel() {
1793/// Double_t a;
1794/// Int_t i,nfound=0,bin;
1795/// Int_t nbins = 256;
1796/// Int_t xmin = 0;
1797/// Int_t xmax = nbins;
1798/// Double_t * source = new Double_t[nbins];
1799/// Double_t * dest = new Double_t[nbins];
1800/// TH1F *h = new TH1F("h","Fitting using AWMI algorithm",nbins,xmin,xmax);
1801/// TH1F *d = new TH1F("d","",nbins,xmin,xmax);
1802/// TFile *f = new TFile("TSpectrum.root");
1803/// h=(TH1F*) f->Get("fit;1");
1804/// for (i = 0; i < nbins;i++) source[i]=h->GetBinContent(i + 1);
1805/// TCanvas *Fit1 = gROOT->GetListOfCanvases()->FindObject("Fit1");
1806/// if (!Fit1) Fit1 = new TCanvas("Fit1","Fit1",10,10,1000,700);
1807/// h->Draw("L");
1808/// TSpectrum *s = new TSpectrum();
1809/// //searching for candidate peaks positions
1810/// nfound = s->SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);
1811/// Bool_t *FixPos = new Bool_t[nfound];
1812/// Bool_t *FixAmp = new Bool_t[nfound];
1813/// for(i = 0; i< nfound ; i++){
1814/// FixPos[i] = kFALSE;
1815/// FixAmp[i] = kFALSE;
1816/// }
1817/// //filling in the initial estimates of the input parameters
1818/// Double_t *PosX = new Double_t[nfound];
1819/// Double_t *PosY = new Double_t[nfound];
1820/// PosX = s->GetPositionX();
1821/// for (i = 0; i < nfound; i++) {
1822/// a=PosX[i];
1823/// bin = 1 + Int_t(a + 0.5);
1824/// PosY[i] = h->GetBinContent(bin);
1825/// }
1826/// TSpectrumFit *pfit = new TSpectrumFit(nfound);
1827/// pfit->SetFitParameters(xmin, xmax-1, 1000, 0.1, pfit->kFitOptimChiCounts,
1828/// pfit->kFitAlphaHalving, pfit->kFitPower2,
1829/// pfit->kFitTaylorOrderFirst);
1830/// pfit->SetPeakParameters(2, kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);
1831/// pfit->FitStiefel(source);
1832/// Double_t *CalcPositions = new Double_t[nfound];
1833/// Double_t *CalcAmplitudes = new Double_t[nfound];
1834/// CalcPositions=pfit->GetPositions();
1835/// CalcAmplitudes=pfit->GetAmplitudes();
1836/// for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);
1837/// d->SetLineColor(kRed);
1838/// d->Draw("SAMEL");
1839/// for (i = 0; i < nfound; i++) {
1840/// a=CalcPositions[i];
1841/// bin = 1 + Int_t(a + 0.5);
1842/// PosX[i] = d->GetBinCenter(bin);
1843/// PosY[i] = d->GetBinContent(bin);
1844/// }
1845/// TPolyMarker * pm = (TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");
1846/// if (pm) {
1847/// h->GetListOfFunctions()->Remove(pm);
1848/// delete pm;
1849/// }
1850/// pm = new TPolyMarker(nfound, PosX, PosY);
1851/// h->GetListOfFunctions()->Add(pm);
1852/// pm->SetMarkerStyle(23);
1853/// pm->SetMarkerColor(kRed);
1854/// pm->SetMarkerSize(1);
1855/// }
1856/// ~~~
1857
1858void TSpectrumFit::FitStiefel(Double_t *source)
1859{
1860 Int_t i, j, k, shift =
1861 2 * fNPeaks + 7, peak_vel, rozmer, iter, regul_cycle,
1862 flag;
1863 Double_t a, b, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
1864 0, pi, pmin = 0, chi_cel = 0, chi_er;
1865 Double_t *working_space = new Double_t[5 * (2 * fNPeaks + 7)];
1866 for (i = 0, j = 0; i < fNPeaks; i++) {
1867 working_space[2 * i] = fAmpInit[i]; //vector parameter
1868 if (fFixAmp[i] == false) {
1869 working_space[shift + j] = fAmpInit[i]; //vector xk
1870 j += 1;
1871 }
1872 working_space[2 * i + 1] = fPositionInit[i]; //vector parameter
1873 if (fFixPosition[i] == false) {
1874 working_space[shift + j] = fPositionInit[i]; //vector xk
1875 j += 1;
1876 }
1877 }
1878 peak_vel = 2 * i;
1879 working_space[2 * i] = fSigmaInit; //vector parameter
1880 if (fFixSigma == false) {
1881 working_space[shift + j] = fSigmaInit; //vector xk
1882 j += 1;
1883 }
1884 working_space[2 * i + 1] = fTInit; //vector parameter
1885 if (fFixT == false) {
1886 working_space[shift + j] = fTInit; //vector xk
1887 j += 1;
1888 }
1889 working_space[2 * i + 2] = fBInit; //vector parameter
1890 if (fFixB == false) {
1891 working_space[shift + j] = fBInit; //vector xk
1892 j += 1;
1893 }
1894 working_space[2 * i + 3] = fSInit; //vector parameter
1895 if (fFixS == false) {
1896 working_space[shift + j] = fSInit; //vector xk
1897 j += 1;
1898 }
1899 working_space[2 * i + 4] = fA0Init; //vector parameter
1900 if (fFixA0 == false) {
1901 working_space[shift + j] = fA0Init; //vector xk
1902 j += 1;
1903 }
1904 working_space[2 * i + 5] = fA1Init; //vector parameter
1905 if (fFixA1 == false) {
1906 working_space[shift + j] = fA1Init; //vector xk
1907 j += 1;
1908 }
1909 working_space[2 * i + 6] = fA2Init; //vector parameter
1910 if (fFixA2 == false) {
1911 working_space[shift + j] = fA2Init; //vector xk
1912 j += 1;
1913 }
1914 rozmer = j;
1915 if (rozmer == 0){
1916 Error ("FitAwmi","All parameters are fixed");
1917 delete [] working_space;
1918 return;
1919 }
1920 if (rozmer >= fXmax - fXmin + 1){
1921 Error ("FitAwmi","Number of fitted parameters is larger than # of fitted points");
1922 delete [] working_space;
1923 return;
1924 }
1925 Double_t **working_matrix = new Double_t *[rozmer];
1926 for (i = 0; i < rozmer; i++)
1927 working_matrix[i] = new Double_t[rozmer + 4];
1928 for (iter = 0; iter < fNumberIterations; iter++) {
1929 for (j = 0; j < rozmer; j++) {
1930 working_space[3 * shift + j] = 0; //temp
1931 for (k = 0; k < (rozmer + 4); k++) {
1932 working_matrix[j][k] = 0;
1933 }
1934 }
1935
1936 //filling working matrix
1937 alpha = fAlpha;
1938 chi_opt = 0;
1939 for (i = fXmin; i <= fXmax; i++) {
1940
1941 //calculation of gradient vector
1942 for (j = 0, k = 0; j < fNPeaks; j++) {
1943 if (fFixAmp[j] == false) {
1944 working_space[2 * shift + k] =
1945 Deramp((Double_t) i, working_space[2 * j + 1],
1946 working_space[peak_vel],
1947 working_space[peak_vel + 1],
1948 working_space[peak_vel + 3],
1949 working_space[peak_vel + 2]);
1950 k += 1;
1951 }
1952 if (fFixPosition[j] == false) {
1953 working_space[2 * shift + k] =
1954 Deri0((Double_t) i, working_space[2 * j],
1955 working_space[2 * j + 1], working_space[peak_vel],
1956 working_space[peak_vel + 1],
1957 working_space[peak_vel + 3],
1958 working_space[peak_vel + 2]);
1959 k += 1;
1960 }
1961 } if (fFixSigma == false) {
1962 working_space[2 * shift + k] =
1963 Dersigma(fNPeaks, (Double_t) i, working_space,
1964 working_space[peak_vel],
1965 working_space[peak_vel + 1],
1966 working_space[peak_vel + 3],
1967 working_space[peak_vel + 2]);
1968 k += 1;
1969 }
1970 if (fFixT == false) {
1971 working_space[2 * shift + k] =
1972 Dert(fNPeaks, (Double_t) i, working_space,
1973 working_space[peak_vel], working_space[peak_vel + 2]);
1974 k += 1;
1975 }
1976 if (fFixB == false) {
1977 working_space[2 * shift + k] =
1978 Derb(fNPeaks, (Double_t) i, working_space,
1979 working_space[peak_vel], working_space[peak_vel + 1],
1980 working_space[peak_vel + 2]);
1981 k += 1;
1982 }
1983 if (fFixS == false) {
1984 working_space[2 * shift + k] =
1985 Ders(fNPeaks, (Double_t) i, working_space,
1986 working_space[peak_vel]);
1987 k += 1;
1988 }
1989 if (fFixA0 == false) {
1990 working_space[2 * shift + k] = 1.;
1991 k += 1;
1992 }
1993 if (fFixA1 == false) {
1994 working_space[2 * shift + k] = Dera1((Double_t) i);
1995 k += 1;
1996 }
1997 if (fFixA2 == false) {
1998 working_space[2 * shift + k] = Dera2((Double_t) i);
1999 k += 1;
2000 }
2001 yw = source[i];
2002 ywm = yw;
2003 f = Shape(fNPeaks, (Double_t) i, working_space,
2004 working_space[peak_vel], working_space[peak_vel + 1],
2005 working_space[peak_vel + 3],
2006 working_space[peak_vel + 2],
2007 working_space[peak_vel + 4],
2008 working_space[peak_vel + 5],
2009 working_space[peak_vel + 6]);
2010 if (fStatisticType == kFitOptimMaxLikelihood) {
2011 if (f > 0.00001)
2012 chi_opt += yw * TMath::Log(f) - f;
2013 }
2014
2015 else {
2016 if (ywm != 0)
2017 chi_opt += (yw - f) * (yw - f) / ywm;
2018 }
2019 if (fStatisticType == kFitOptimChiFuncValues) {
2020 ywm = f;
2021 if (f < 0.00001)
2022 ywm = 0.00001;
2023 }
2024
2025 else if (fStatisticType == kFitOptimMaxLikelihood) {
2026 ywm = f;
2027 if (f < 0.00001)
2028 ywm = 0.00001;
2029 }
2030
2031 else {
2032 if (ywm == 0)
2033 ywm = 1;
2034 }
2035 for (j = 0; j < rozmer; j++) {
2036 for (k = 0; k < rozmer; k++) {
2037 b = working_space[2 * shift +
2038 j] * working_space[2 * shift + k] / ywm;
2039 if (fStatisticType == kFitOptimChiFuncValues)
2040 b = b * (4 * yw - 2 * f) / ywm;
2041 working_matrix[j][k] += b;
2042 if (j == k)
2043 working_space[3 * shift + j] += b;
2044 }
2045 }
2046 if (fStatisticType == kFitOptimChiFuncValues)
2047 b = (f * f - yw * yw) / (ywm * ywm);
2048
2049 else
2050 b = (f - yw) / ywm;
2051 for (j = 0; j < rozmer; j++) {
2052 working_matrix[j][rozmer] -= b * working_space[2 * shift + j];
2053 }
2054 }
2055 for (i = 0; i < rozmer; i++) {
2056 working_matrix[i][rozmer + 1] = 0; //xk
2057 }
2058 StiefelInversion(working_matrix, rozmer);
2059 for (i = 0; i < rozmer; i++) {
2060 working_space[2 * shift + i] = working_matrix[i][rozmer + 1]; //der
2061 }
2062
2063 //calculate chi_opt
2064 chi2 = chi_opt;
2065 chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
2066
2067 //calculate new parameters
2068 regul_cycle = 0;
2069 for (j = 0; j < rozmer; j++) {
2070 working_space[4 * shift + j] = working_space[shift + j]; //temp_xk[j]=xk[j]
2071 }
2072
2073 do {
2074 if (fAlphaOptim == kFitAlphaOptimal) {
2075 if (fStatisticType != kFitOptimMaxLikelihood)
2076 chi_min = 10000 * chi2;
2077
2078 else
2079 chi_min = 0.1 * chi2;
2080 flag = 0;
2081 for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
2082 for (j = 0; j < rozmer; j++) {
2083 working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
2084 }
2085 for (i = 0, j = 0; i < fNPeaks; i++) {
2086 if (fFixAmp[i] == false) {
2087 if (working_space[shift + j] < 0) //xk[j]
2088 working_space[shift + j] = 0; //xk[j]
2089 working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
2090 j += 1;
2091 }
2092 if (fFixPosition[i] == false) {
2093 if (working_space[shift + j] < fXmin) //xk[j]
2094 working_space[shift + j] = fXmin; //xk[j]
2095 if (working_space[shift + j] > fXmax) //xk[j]
2096 working_space[shift + j] = fXmax; //xk[j]
2097 working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
2098 j += 1;
2099 }
2100 }
2101 if (fFixSigma == false) {
2102 if (working_space[shift + j] < 0.001) { //xk[j]
2103 working_space[shift + j] = 0.001; //xk[j]
2104 }
2105 working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
2106 j += 1;
2107 }
2108 if (fFixT == false) {
2109 working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
2110 j += 1;
2111 }
2112 if (fFixB == false) {
2113 if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2114 if (working_space[shift + j] < 0) //xk[j]
2115 working_space[shift + j] = -0.001; //xk[j]
2116 else
2117 working_space[shift + j] = 0.001; //xk[j]
2118 }
2119 working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
2120 j += 1;
2121 }
2122 if (fFixS == false) {
2123 working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
2124 j += 1;
2125 }
2126 if (fFixA0 == false) {
2127 working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
2128 j += 1;
2129 }
2130 if (fFixA1 == false) {
2131 working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
2132 j += 1;
2133 }
2134 if (fFixA2 == false) {
2135 working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
2136 j += 1;
2137 }
2138 chi2 = 0;
2139 for (i = fXmin; i <= fXmax; i++) {
2140 yw = source[i];
2141 ywm = yw;
2142 f = Shape(fNPeaks, (Double_t) i, working_space,
2143 working_space[peak_vel],
2144 working_space[peak_vel + 1],
2145 working_space[peak_vel + 3],
2146 working_space[peak_vel + 2],
2147 working_space[peak_vel + 4],
2148 working_space[peak_vel + 5],
2149 working_space[peak_vel + 6]);
2150 if (fStatisticType == kFitOptimChiFuncValues) {
2151 ywm = f;
2152 if (f < 0.00001)
2153 ywm = 0.00001;
2154 }
2155 if (fStatisticType == kFitOptimMaxLikelihood) {
2156 if (f > 0.00001)
2157 chi2 += yw * TMath::Log(f) - f;
2158 }
2159
2160 else {
2161 if (ywm != 0)
2162 chi2 += (yw - f) * (yw - f) / ywm;
2163 }
2164 }
2165 if ((chi2 < chi_min
2166 && fStatisticType != kFitOptimMaxLikelihood)
2167 || (chi2 > chi_min
2168 && fStatisticType == kFitOptimMaxLikelihood)) {
2169 pmin = pi, chi_min = chi2;
2170 }
2171
2172 else
2173 flag = 1;
2174 if (pi == 0.1)
2175 chi_min = chi2;
2176 chi = chi_min;
2177 }
2178 if (pmin != 0.1) {
2179 for (j = 0; j < rozmer; j++) {
2180 working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pmin*alpha*der[j]
2181 }
2182 for (i = 0, j = 0; i < fNPeaks; i++) {
2183 if (fFixAmp[i] == false) {
2184 if (working_space[shift + j] < 0) //xk[j]
2185 working_space[shift + j] = 0; //xk[j]
2186 working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
2187 j += 1;
2188 }
2189 if (fFixPosition[i] == false) {
2190 if (working_space[shift + j] < fXmin) //xk[j]
2191 working_space[shift + j] = fXmin; //xk[j]
2192 if (working_space[shift + j] > fXmax) //xk[j]
2193 working_space[shift + j] = fXmax; //xk[j]
2194 working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
2195 j += 1;
2196 }
2197 }
2198 if (fFixSigma == false) {
2199 if (working_space[shift + j] < 0.001) { //xk[j]
2200 working_space[shift + j] = 0.001; //xk[j]
2201 }
2202 working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
2203 j += 1;
2204 }
2205 if (fFixT == false) {
2206 working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
2207 j += 1;
2208 }
2209 if (fFixB == false) {
2210 if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2211 if (working_space[shift + j] < 0) //xk[j]
2212 working_space[shift + j] = -0.001; //xk[j]
2213 else
2214 working_space[shift + j] = 0.001; //xk[j]
2215 }
2216 working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
2217 j += 1;
2218 }
2219 if (fFixS == false) {
2220 working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
2221 j += 1;
2222 }
2223 if (fFixA0 == false) {
2224 working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
2225 j += 1;
2226 }
2227 if (fFixA1 == false) {
2228 working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
2229 j += 1;
2230 }
2231 if (fFixA2 == false) {
2232 working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
2233 j += 1;
2234 }
2235 chi = chi_min;
2236 }
2237 }
2238
2239 else {
2240 for (j = 0; j < rozmer; j++) {
2241 working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+alpha*der[j]
2242 }
2243 for (i = 0, j = 0; i < fNPeaks; i++) {
2244 if (fFixAmp[i] == false) {
2245 if (working_space[shift + j] < 0) //xk[j]
2246 working_space[shift + j] = 0; //xk[j]
2247 working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
2248 j += 1;
2249 }
2250 if (fFixPosition[i] == false) {
2251 if (working_space[shift + j] < fXmin) //xk[j]
2252 working_space[shift + j] = fXmin; //xk[j]
2253 if (working_space[shift + j] > fXmax) //xk[j]
2254 working_space[shift + j] = fXmax; //xk[j]
2255 working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
2256 j += 1;
2257 }
2258 }
2259 if (fFixSigma == false) {
2260 if (working_space[shift + j] < 0.001) { //xk[j]
2261 working_space[shift + j] = 0.001; //xk[j]
2262 }
2263 working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
2264 j += 1;
2265 }
2266 if (fFixT == false) {
2267 working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
2268 j += 1;
2269 }
2270 if (fFixB == false) {
2271 if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2272 if (working_space[shift + j] < 0) //xk[j]
2273 working_space[shift + j] = -0.001; //xk[j]
2274 else
2275 working_space[shift + j] = 0.001; //xk[j]
2276 }
2277 working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
2278 j += 1;
2279 }
2280 if (fFixS == false) {
2281 working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
2282 j += 1;
2283 }
2284 if (fFixA0 == false) {
2285 working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
2286 j += 1;
2287 }
2288 if (fFixA1 == false) {
2289 working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
2290 j += 1;
2291 }
2292 if (fFixA2 == false) {
2293 working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
2294 j += 1;
2295 }
2296 chi = 0;
2297 for (i = fXmin; i <= fXmax; i++) {
2298 yw = source[i];
2299 ywm = yw;
2300 f = Shape(fNPeaks, (Double_t) i, working_space,
2301 working_space[peak_vel],
2302 working_space[peak_vel + 1],
2303 working_space[peak_vel + 3],
2304 working_space[peak_vel + 2],
2305 working_space[peak_vel + 4],
2306 working_space[peak_vel + 5],
2307 working_space[peak_vel + 6]);
2308 if (fStatisticType == kFitOptimChiFuncValues) {
2309 ywm = f;
2310 if (f < 0.00001)
2311 ywm = 0.00001;
2312 }
2313 if (fStatisticType == kFitOptimMaxLikelihood) {
2314 if (f > 0.00001)
2315 chi += yw * TMath::Log(f) - f;
2316 }
2317
2318 else {
2319 if (ywm != 0)
2320 chi += (yw - f) * (yw - f) / ywm;
2321 }
2322 }
2323 }
2324 chi2 = chi;
2325 chi = TMath::Sqrt(TMath::Abs(chi));
2326 if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
2327 alpha = alpha * chi_opt / (2 * chi);
2328
2329 else if (fAlphaOptim == kFitAlphaOptimal)
2330 alpha = alpha / 10.0;
2331 iter += 1;
2332 regul_cycle += 1;
2333 } while (((chi > chi_opt
2334 && fStatisticType != kFitOptimMaxLikelihood)
2335 || (chi < chi_opt
2336 && fStatisticType == kFitOptimMaxLikelihood))
2337 && regul_cycle < kFitNumRegulCycles);
2338 for (j = 0; j < rozmer; j++) {
2339 working_space[4 * shift + j] = 0; //temp_xk[j]
2340 working_space[2 * shift + j] = 0; //der[j]
2341 }
2342 for (i = fXmin, chi_cel = 0; i <= fXmax; i++) {
2343 yw = source[i];
2344 if (yw == 0)
2345 yw = 1;
2346 f = Shape(fNPeaks, (Double_t) i, working_space,
2347 working_space[peak_vel], working_space[peak_vel + 1],
2348 working_space[peak_vel + 3],
2349 working_space[peak_vel + 2],
2350 working_space[peak_vel + 4],
2351 working_space[peak_vel + 5],
2352 working_space[peak_vel + 6]);
2353 chi_opt = (yw - f) * (yw - f) / yw;
2354 chi_cel += (yw - f) * (yw - f) / yw;
2355
2356 //calculate gradient vector
2357 for (j = 0, k = 0; j < fNPeaks; j++) {
2358 if (fFixAmp[j] == false) {
2359 a = Deramp((Double_t) i, working_space[2 * j + 1],
2360 working_space[peak_vel],
2361 working_space[peak_vel + 1],
2362 working_space[peak_vel + 3],
2363 working_space[peak_vel + 2]);
2364 if (yw != 0) {
2365 working_space[2 * shift + k] += chi_opt; //der[k]
2366 b = a * a / yw;
2367 working_space[4 * shift + k] += b; //temp_xk[k]
2368 }
2369 k += 1;
2370 }
2371 if (fFixPosition[j] == false) {
2372 a = Deri0((Double_t) i, working_space[2 * j],
2373 working_space[2 * j + 1],
2374 working_space[peak_vel],
2375 working_space[peak_vel + 1],
2376 working_space[peak_vel + 3],
2377 working_space[peak_vel + 2]);
2378 if (yw != 0) {
2379 working_space[2 * shift + k] += chi_opt; //der[k]
2380 b = a * a / yw;
2381 working_space[4 * shift + k] += b; //temp_xk[k]
2382 }
2383 k += 1;
2384 }
2385 }
2386 if (fFixSigma == false) {
2387 a = Dersigma(fNPeaks, (Double_t) i, working_space,
2388 working_space[peak_vel],
2389 working_space[peak_vel + 1],
2390 working_space[peak_vel + 3],
2391 working_space[peak_vel + 2]);
2392 if (yw != 0) {
2393 working_space[2 * shift + k] += chi_opt; //der[k]
2394 b = a * a / yw;
2395 working_space[4 * shift + k] += b; //temp_xk[k]
2396 }
2397 k += 1;
2398 }
2399 if (fFixT == false) {
2400 a = Dert(fNPeaks, (Double_t) i, working_space,
2401 working_space[peak_vel],
2402 working_space[peak_vel + 2]);
2403 if (yw != 0) {
2404 working_space[2 * shift + k] += chi_opt; //der[k]
2405 b = a * a / yw;
2406 working_space[4 * shift + k] += b; //temp_xk[k]
2407 }
2408 k += 1;
2409 }
2410 if (fFixB == false) {
2411 a = Derb(fNPeaks, (Double_t) i, working_space,
2412 working_space[peak_vel], working_space[peak_vel + 1],
2413 working_space[peak_vel + 2]);
2414 if (yw != 0) {
2415 working_space[2 * shift + k] += chi_opt; //der[k]
2416 b = a * a / yw;
2417 working_space[4 * shift + k] += b; //temp_xk[k]
2418 }
2419 k += 1;
2420 }
2421 if (fFixS == false) {
2422 a = Ders(fNPeaks, (Double_t) i, working_space,
2423 working_space[peak_vel]);
2424 if (yw != 0) {
2425 working_space[2 * shift + k] += chi_opt; //der[k]
2426 b = a * a / yw;
2427 working_space[4 * shift + k] += b; //tem_xk[k]
2428 }
2429 k += 1;
2430 }
2431 if (fFixA0 == false) {
2432 a = 1.0;
2433 if (yw != 0) {
2434 working_space[2 * shift + k] += chi_opt; //der[k]
2435 b = a * a / yw;
2436 working_space[4 * shift + k] += b; //temp_xk[k]
2437 }
2438 k += 1;
2439 }
2440 if (fFixA1 == false) {
2441 a = Dera1((Double_t) i);
2442 if (yw != 0) {
2443 working_space[2 * shift + k] += chi_opt; //der[k]
2444 b = a * a / yw;
2445 working_space[4 * shift + k] += b; //temp_xk[k]
2446 }
2447 k += 1;
2448 }
2449 if (fFixA2 == false) {
2450 a = Dera2((Double_t) i);
2451 if (yw != 0) {
2452 working_space[2 * shift + k] += chi_opt; //der[k]
2453 b = a * a / yw;
2454 working_space[4 * shift + k] += b; //temp_xk[k]
2455 }
2456 k += 1;
2457 }
2458 }
2459 }
2460 b = fXmax - fXmin + 1 - rozmer;
2461 chi_er = chi_cel / b;
2462 for (i = 0, j = 0; i < fNPeaks; i++) {
2463 fArea[i] =
2464 Area(working_space[2 * i], working_space[peak_vel],
2465 working_space[peak_vel + 1], working_space[peak_vel + 2]);
2466 if (fFixAmp[i] == false) {
2467 fAmpCalc[i] = working_space[shift + j]; //xk[j]
2468 if (working_space[3 * shift + j] != 0)
2469 fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2470 if (fArea[i] > 0) {
2471 a = Derpa(working_space[peak_vel],
2472 working_space[peak_vel + 1],
2473 working_space[peak_vel + 2]);
2474 b = working_space[4 * shift + j]; //temp_xk[j]
2475 if (b == 0)
2476 b = 1;
2477
2478 else
2479 b = 1 / b;
2480 fAreaErr[i] = TMath::Sqrt(TMath::Abs(a * a * b * chi_er));
2481 }
2482
2483 else
2484 fAreaErr[i] = 0;
2485 j += 1;
2486 }
2487
2488 else {
2489 fAmpCalc[i] = fAmpInit[i];
2490 fAmpErr[i] = 0;
2491 fAreaErr[i] = 0;
2492 }
2493 if (fFixPosition[i] == false) {
2494 fPositionCalc[i] = working_space[shift + j]; //xk[j]
2495 if (working_space[3 * shift + j] != 0) //temp[j]
2496 fPositionErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //Der[j]/temp[j]
2497 j += 1;
2498 }
2499
2500 else {
2501 fPositionCalc[i] = fPositionInit[i];
2502 fPositionErr[i] = 0;
2503 }
2504 }
2505 if (fFixSigma == false) {
2506 fSigmaCalc = working_space[shift + j]; //xk[j]
2507 if (working_space[3 * shift + j] != 0) //temp[j]
2508 fSigmaErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2509 j += 1;
2510 }
2511
2512 else {
2513 fSigmaCalc = fSigmaInit;
2514 fSigmaErr = 0;
2515 }
2516 if (fFixT == false) {
2517 fTCalc = working_space[shift + j]; //xk[j]
2518 if (working_space[3 * shift + j] != 0) //temp[j]
2519 fTErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2520 j += 1;
2521 }
2522
2523 else {
2524 fTCalc = fTInit;
2525 fTErr = 0;
2526 }
2527 if (fFixB == false) {
2528 fBCalc = working_space[shift + j]; //xk[j]
2529 if (working_space[3 * shift + j] != 0) //temp[j]
2530 fBErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2531 j += 1;
2532 }
2533
2534 else {
2535 fBCalc = fBInit;
2536 fBErr = 0;
2537 }
2538 if (fFixS == false) {
2539 fSCalc = working_space[shift + j]; //xk[j]
2540 if (working_space[3 * shift + j] != 0) //temp[j]
2541 fSErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2542 j += 1;
2543 }
2544
2545 else {
2546 fSCalc = fSInit;
2547 fSErr = 0;
2548 }
2549 if (fFixA0 == false) {
2550 fA0Calc = working_space[shift + j]; //xk[j]
2551 if (working_space[3 * shift + j] != 0) //temp[j]
2552 fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2553 j += 1;
2554 }
2555
2556 else {
2557 fA0Calc = fA0Init;
2558 fA0Err = 0;
2559 }
2560 if (fFixA1 == false) {
2561 fA1Calc = working_space[shift + j]; //xk[j]
2562 if (working_space[3 * shift + j] != 0) //temp[j]
2563 fA1Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2564 j += 1;
2565 }
2566
2567 else {
2568 fA1Calc = fA1Init;
2569 fA1Err = 0;
2570 }
2571 if (fFixA2 == false) {
2572 fA2Calc = working_space[shift + j]; //xk[j]
2573 if (working_space[3 * shift + j] != 0) //temp[j]
2574 fA2Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2575 j += 1;
2576 }
2577
2578 else {
2579 fA2Calc = fA2Init;
2580 fA2Err = 0;
2581 }
2582 b = fXmax - fXmin + 1 - rozmer;
2583 fChi = chi_cel / b;
2584 for (i = fXmin; i <= fXmax; i++) {
2585 f = Shape(fNPeaks, (Double_t) i, working_space,
2586 working_space[peak_vel], working_space[peak_vel + 1],
2587 working_space[peak_vel + 3], working_space[peak_vel + 2],
2588 working_space[peak_vel + 4], working_space[peak_vel + 5],
2589 working_space[peak_vel + 6]);
2590 source[i] = f;
2591 }
2592 for (i = 0; i < rozmer; i++)
2593 delete [] working_matrix[i];
2594 delete [] working_matrix;
2595 delete [] working_space;
2596 return;
2597}
2598
2599////////////////////////////////////////////////////////////////////////////////
2600/// This function sets the following fitting parameters:
2601/// - xmin, xmax - fitting region
2602/// - numberIterations - # of desired iterations in the fit
2603/// - alpha - convergence coefficient, it should be positive number and <=1, for details see references
2604/// - statisticType - type of statistics, possible values kFitOptimChiCounts (chi square statistics with counts as weighting coefficients), kFitOptimChiFuncValues (chi square statistics with function values as weighting coefficients),kFitOptimMaxLikelihood
2605/// - alphaOptim - optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal
2606/// - power - possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting function.
2607/// - fitTaylor - order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi fitting function.
2608
2609void TSpectrumFit::SetFitParameters(Int_t xmin,Int_t xmax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor)
2610{
2611 if(xmin<0 || xmax <= xmin){
2612 Error("SetFitParameters", "Wrong range");
2613 return;
2614 }
2615 if (numberIterations <= 0){
2616 Error("SetFitParameters","Invalid number of iterations, must be positive");
2617 return;
2618 }
2619 if (alpha <= 0 || alpha > 1){
2620 Error ("SetFitParameters","Invalid step coefficient alpha, must be > than 0 and <=1");
2621 return;
2622 }
2623 if (statisticType != kFitOptimChiCounts
2624 && statisticType != kFitOptimChiFuncValues
2625 && statisticType != kFitOptimMaxLikelihood){
2626 Error("SetFitParameters","Wrong type of statistic");
2627 return;
2628 }
2629 if (alphaOptim != kFitAlphaHalving
2630 && alphaOptim != kFitAlphaOptimal){
2631 Error("SetFitParameters","Wrong optimization algorithm");
2632 return;
2633 }
2634 if (power != kFitPower2 && power != kFitPower4
2635 && power != kFitPower6 && power != kFitPower8
2636 && power != kFitPower10 && power != kFitPower12){
2637 Error("SetFitParameters","Wrong power");
2638 return;
2639 }
2640 if (fitTaylor != kFitTaylorOrderFirst
2641 && fitTaylor != kFitTaylorOrderSecond){
2642 Error("SetFitParameters","Wrong order of Taylor development");
2643 return;
2644 }
2645 fXmin=xmin,fXmax=xmax,fNumberIterations=numberIterations,fAlpha=alpha,fStatisticType=statisticType,fAlphaOptim=alphaOptim,fPower=power,fFitTaylor=fitTaylor;
2646}
2647
2648////////////////////////////////////////////////////////////////////////////////
2649/// This function sets the following fitting parameters of peaks:
2650/// - sigma - initial value of sigma parameter
2651/// - fixSigma - logical value of sigma parameter, which allows to fix the parameter (not to fit)
2652/// - positionInit - array of initial values of peaks positions
2653/// - fixPosition - array of logical values which allow to fix appropriate positions (not fit). However they are present in the estimated functional.
2654/// - ampInit - array of initial values of peaks amplitudes
2655/// - fixAmp - array of logical values which allow to fix appropriate amplitudes (not fit). However they are present in the estimated functional
2656
2657void TSpectrumFit::SetPeakParameters(Double_t sigma, Bool_t fixSigma, const Double_t *positionInit, const Bool_t *fixPosition, const Double_t *ampInit, const Bool_t *fixAmp)
2658{
2659 Int_t i;
2660 if (sigma <= 0){
2661 Error ("SetPeakParameters","Invalid sigma, must be > than 0");
2662 return;
2663 }
2664 for(i=0; i < fNPeaks; i++){
2665 if((Int_t) positionInit[i] < fXmin || (Int_t) positionInit[i] > fXmax){
2666 Error ("SetPeakParameters","Invalid peak position, must be in the range fXmin, fXmax");
2667 return;
2668 }
2669 if(ampInit[i] < 0){
2670 Error ("SetPeakParameters","Invalid peak amplitude, must be > than 0");
2671 return;
2672 }
2673 }
2674 fSigmaInit = sigma,fFixSigma = fixSigma;
2675 for(i=0; i < fNPeaks; i++){
2676 fPositionInit[i] = (Double_t) positionInit[i];
2677 fFixPosition[i] = fixPosition[i];
2678 fAmpInit[i] = (Double_t) ampInit[i];
2679 fFixAmp[i] = fixAmp[i];
2680 }
2681}
2682
2683////////////////////////////////////////////////////////////////////////////////
2684/// This function sets the following fitting parameters of background:
2685/// - a0Init - initial value of a0 parameter (background is estimated as a0+a1*x+a2*x*x)
2686/// - fixA0 - logical value of a0 parameter, which allows to fix the parameter (not to fit)
2687/// - a1Init - initial value of a1 parameter
2688/// - fixA1 - logical value of a1 parameter, which allows to fix the parameter (not to fit)
2689/// - a2Init - initial value of a2 parameter
2690/// - fixA2 - logical value of a2 parameter, which allows to fix the parameter (not to fit)
2691
2701
2702////////////////////////////////////////////////////////////////////////////////
2703/// This function sets the following fitting parameters of tails of peaks
2704/// - tInit - initial value of t parameter
2705/// - fixT - logical value of t parameter, which allows to fix the parameter (not to fit)
2706/// - bInit - initial value of b parameter
2707/// - fixB - logical value of b parameter, which allows to fix the parameter (not to fit)
2708/// - sInit - initial value of s parameter
2709/// - fixS - logical value of s parameter, which allows to fix the parameter (not to fit)
2710
2711void TSpectrumFit::SetTailParameters(Double_t tInit, Bool_t fixT, Double_t bInit, Bool_t fixB, Double_t sInit, Bool_t fixS)
2712{
2713 fTInit = tInit;
2714 fFixT = fixT;
2715 fBInit = bInit;
2716 fFixB = fixB;
2717 fSInit = sInit;
2718 fFixS = fixS;
2719}
2720
2721////////////////////////////////////////////////////////////////////////////////
2722/// This function gets the sigma parameter and its error
2723/// - sigma - gets the fitted value of sigma parameter
2724/// - sigmaErr - gets error value of sigma parameter
2725
2731
2732////////////////////////////////////////////////////////////////////////////////
2733/// This function gets the background parameters and their errors
2734/// - a0 - gets the fitted value of a0 parameter
2735/// - a0Err - gets error value of a0 parameter
2736/// - a1 - gets the fitted value of a1 parameter
2737/// - a1Err - gets error value of a1 parameter
2738/// - a2 - gets the fitted value of a2 parameter
2739/// - a2Err - gets error value of a2 parameter
2740
2750
2751////////////////////////////////////////////////////////////////////////////////
2752/// This function gets the tail parameters and their errors
2753/// - t - gets the fitted value of t parameter
2754/// - tErr - gets error value of t parameter
2755/// - b - gets the fitted value of b parameter
2756/// - bErr - gets error value of b parameter
2757/// - s - gets the fitted value of s parameter
2758/// - sErr - gets error value of s parameter
2759///////////////////////////////////////////////////////////////////////////////
2760
2761void TSpectrumFit::GetTailParameters(Double_t &t, Double_t &tErr, Double_t &b, Double_t &bErr, Double_t &s, Double_t &sErr)
2762{
2763 t = fTCalc;
2764 tErr = fTErr;
2765 b = fBCalc;
2766 bErr = fBErr;
2767 s = fSCalc;
2768 sErr = fSErr;
2769}
2770
d
#define d(i)
Definition RSha256.hxx:102
b
#define b(i)
Definition RSha256.hxx:100
f
#define f(i)
Definition RSha256.hxx:104
c
#define c(i)
Definition RSha256.hxx:101
a
#define a(i)
Definition RSha256.hxx:99
e
#define e(i)
Definition RSha256.hxx:103
size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Double_t
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
TRangeDynCast
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Definition TCollection.h:358
p
winID h TVirtualViewer3D TVirtualGLPainter p
Definition TGWin32VirtualGLProxy.cxx:51
r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Definition TGWin32VirtualXProxy.cxx:168
xmin
float xmin
Definition THbookFile.cxx:95
q
float * q
Definition THbookFile.cxx:89
xmax
float xmax
Definition THbookFile.cxx:95
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
TSpectrumFit::fFixSigma
Bool_t fFixSigma
logical value of sigma parameter, which allows to fix the parameter (not to fit).
Definition TSpectrumFit.h:61
TSpectrumFit::Ders
Double_t Ders(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma)
This function calculates derivative of peaks shape function (see manual) according to relative amplit...
Definition TSpectrumFit.cxx:438
TSpectrumFit::Area
Double_t Area(Double_t a, Double_t sigma, Double_t t, Double_t b)
This function calculates area of a peak Function parameters:
Definition TSpectrumFit.cxx:570
TSpectrumFit::fAmpErr
Double_t * fAmpErr
[fNPeaks] array of amplitude errors
Definition TSpectrumFit.h:35
TSpectrumFit::SetPeakParameters
void SetPeakParameters(Double_t sigma, Bool_t fixSigma, const Double_t *positionInit, const Bool_t *fixPosition, const Double_t *ampInit, const Bool_t *fixAmp)
This function sets the following fitting parameters of peaks:
Definition TSpectrumFit.cxx:2657
TSpectrumFit::fAlpha
Double_t fAlpha
convergence coefficient, input parameter, it should be positive number and <=1, for details see refer...
Definition TSpectrumFit.h:28
TSpectrumFit::Deri0
Double_t Deri0(Double_t i, Double_t amp, Double_t i0, Double_t sigma, Double_t t, Double_t s, Double_t b)
This function calculates derivative of peak shape function (see manual) according to peak position.
Definition TSpectrumFit.cxx:267
TSpectrumFit::Ourpowl
Double_t Ourpowl(Double_t a, Int_t pw)
Power function.
Definition TSpectrumFit.cxx:675
TSpectrumFit::fAlphaOptim
Int_t fAlphaOptim
optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal
Definition TSpectrumFit.h:25
TSpectrumFit::fA1Init
Double_t fA1Init
initial value of background a1 parameter(backgroud is estimated as a0+a1*x+a2*x*x)
Definition TSpectrumFit.h:53
TSpectrumFit::FitStiefel
void FitStiefel(Double_t *source)
This function fits the source spectrum.
Definition TSpectrumFit.cxx:1858
TSpectrumFit::Derdersigma
Double_t Derdersigma(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma)
This function calculates second derivative of peaks shape function (see manual) according to sigma of...
Definition TSpectrumFit.cxx:377
TSpectrumFit::fPositionCalc
Double_t * fPositionCalc
[fNPeaks] array of calculated values of fitted positions, output parameters
Definition TSpectrumFit.h:31
TSpectrumFit::fSInit
Double_t fSInit
initial value of s parameter (relative amplitude of step), for details see html manual and references
Definition TSpectrumFit.h:47
TSpectrumFit::fA1Calc
Double_t fA1Calc
calculated value of background a1 parameter
Definition TSpectrumFit.h:54
TSpectrumFit::fSigmaErr
Double_t fSigmaErr
error value of sigma parameter
Definition TSpectrumFit.h:40
TSpectrumFit::fFixAmp
Bool_t * fFixAmp
[fNPeaks] array of logical values which allow to fix appropriate amplitudes (not fit)....
Definition TSpectrumFit.h:60
TSpectrumFit::fAmpCalc
Double_t * fAmpCalc
[fNPeaks] array of calculated values of fitted amplitudes, output parameters
Definition TSpectrumFit.h:34
TSpectrumFit::fFixB
Bool_t fFixB
logical value of b parameter, which allows to fix the parameter (not to fit).
Definition TSpectrumFit.h:63
TSpectrumFit::Dersigma
Double_t Dersigma(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma, Double_t t, Double_t s, Double_t b)
This function calculates derivative of peaks shape function (see manual) according to sigma of peaks.
Definition TSpectrumFit.cxx:332
TSpectrumFit::fBErr
Double_t fBErr
error value of b parameter
Definition TSpectrumFit.h:46
TSpectrumFit::FitAwmi
void FitAwmi(Double_t *source)
This function fits the source spectrum.
Definition TSpectrumFit.cxx:821
TSpectrumFit::fNumberIterations
Int_t fNumberIterations
number of iterations in fitting procedure, input parameter, it should be > 0
Definition TSpectrumFit.h:21
TSpectrumFit::fAmpInit
Double_t * fAmpInit
[fNPeaks] array of initial values of peaks amplitudes, input parameters
Definition TSpectrumFit.h:33
TSpectrumFit::fA1Err
Double_t fA1Err
error value of background a1 parameter
Definition TSpectrumFit.h:55
TSpectrumFit::fA2Err
Double_t fA2Err
error value of background a2 parameter
Definition TSpectrumFit.h:58
TSpectrumFit::fFixPosition
Bool_t * fFixPosition
[fNPeaks] array of logical values which allow to fix appropriate positions (not fit)....
Definition TSpectrumFit.h:59
TSpectrumFit::fA0Calc
Double_t fA0Calc
calculated value of background a0 parameter
Definition TSpectrumFit.h:51
TSpectrumFit::fNPeaks
Int_t fNPeaks
number of peaks present in fit, input parameter, it should be > 0
Definition TSpectrumFit.h:20
TSpectrumFit::fBCalc
Double_t fBCalc
calculated value of b parameter
Definition TSpectrumFit.h:45
TSpectrumFit::fFixA2
Bool_t fFixA2
logical value of a2 parameter, which allows to fix the parameter (not to fit).
Definition TSpectrumFit.h:67
TSpectrumFit::Erfc
Double_t Erfc(Double_t x)
Definition TSpectrumFit.cxx:178
TSpectrumFit::TSpectrumFit
TSpectrumFit(void)
Default constructor.
Definition TSpectrumFit.cxx:36
TSpectrumFit::~TSpectrumFit
~TSpectrumFit() override
Destructor.
Definition TSpectrumFit.cxx:161
TSpectrumFit::fA2Init
Double_t fA2Init
initial value of background a2 parameter(backgroud is estimated as a0+a1*x+a2*x*x)
Definition TSpectrumFit.h:56
TSpectrumFit::Derpt
Double_t Derpt(Double_t a, Double_t sigma, Double_t b)
This function calculates derivative of the area of peak according to t parameter.
Definition TSpectrumFit.cxx:637
TSpectrumFit::fFixA1
Bool_t fFixA1
logical value of a1 parameter, which allows to fix the parameter (not to fit).
Definition TSpectrumFit.h:66
TSpectrumFit::fA2Calc
Double_t fA2Calc
calculated value of background a2 parameter
Definition TSpectrumFit.h:57
TSpectrumFit::Derpb
Double_t Derpb(Double_t a, Double_t sigma, Double_t t, Double_t b)
This function calculates derivative of the area of peak according to b parameter.
Definition TSpectrumFit.cxx:659
TSpectrumFit::fPositionErr
Double_t * fPositionErr
[fNPeaks] array of position errors
Definition TSpectrumFit.h:32
TSpectrumFit::fTInit
Double_t fTInit
initial value of t parameter (relative amplitude of tail), for details see html manual and references
Definition TSpectrumFit.h:41
TSpectrumFit::fFixT
Bool_t fFixT
logical value of t parameter, which allows to fix the parameter (not to fit).
Definition TSpectrumFit.h:62
TSpectrumFit::GetSigma
void GetSigma(Double_t &sigma, Double_t &sigmaErr)
This function gets the sigma parameter and its error.
Definition TSpectrumFit.cxx:2726
TSpectrumFit::fFitTaylor
Int_t fFitTaylor
order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond....
Definition TSpectrumFit.h:27
TSpectrumFit::Dera1
Double_t Dera1(Double_t i)
Derivative of background according to a1.
Definition TSpectrumFit.cxx:493
TSpectrumFit::fSigmaCalc
Double_t fSigmaCalc
calculated value of sigma parameter
Definition TSpectrumFit.h:39
TSpectrumFit::Derpsigma
Double_t Derpsigma(Double_t a, Double_t t, Double_t b)
This function calculates derivative of the area of peak according to sigma of peaks.
Definition TSpectrumFit.cxx:615
TSpectrumFit::fSErr
Double_t fSErr
error value of s parameter
Definition TSpectrumFit.h:49
TSpectrumFit::fFixA0
Bool_t fFixA0
logical value of a0 parameter, which allows to fix the parameter (not to fit).
Definition TSpectrumFit.h:65
TSpectrumFit::fPower
Int_t fPower
possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting ...
Definition TSpectrumFit.h:26
TSpectrumFit::fChi
Double_t fChi
here the fitting functions return resulting chi square
Definition TSpectrumFit.h:29
TSpectrumFit::fArea
Double_t * fArea
[fNPeaks] array of calculated areas of peaks
Definition TSpectrumFit.h:36
TSpectrumFit::Dera2
Double_t Dera2(Double_t i)
Derivative of background according to a2.
Definition TSpectrumFit.cxx:501
TSpectrumFit::fBInit
Double_t fBInit
initial value of b parameter (slope), for details see html manual and references
Definition TSpectrumFit.h:44
TSpectrumFit::fPositionInit
Double_t * fPositionInit
[fNPeaks] array of initial values of peaks positions, input parameters
Definition TSpectrumFit.h:30
TSpectrumFit::GetTailParameters
void GetTailParameters(Double_t &t, Double_t &tErr, Double_t &b, Double_t &bErr, Double_t &s, Double_t &sErr)
This function gets the tail parameters and their errors.
Definition TSpectrumFit.cxx:2761
TSpectrumFit::fFixS
Bool_t fFixS
logical value of s parameter, which allows to fix the parameter (not to fit).
Definition TSpectrumFit.h:64
TSpectrumFit::fA0Init
Double_t fA0Init
initial value of background a0 parameter(backgroud is estimated as a0+a1*x+a2*x*x)
Definition TSpectrumFit.h:50
TSpectrumFit::fXmin
Int_t fXmin
first fitted channel
Definition TSpectrumFit.h:22
TSpectrumFit::SetFitParameters
void SetFitParameters(Int_t xmin, Int_t xmax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor)
This function sets the following fitting parameters:
Definition TSpectrumFit.cxx:2609
TSpectrumFit::fSCalc
Double_t fSCalc
calculated value of s parameter
Definition TSpectrumFit.h:48
TSpectrumFit::SetTailParameters
void SetTailParameters(Double_t tInit, Bool_t fixT, Double_t bInit, Bool_t fixB, Double_t sInit, Bool_t fixS)
This function sets the following fitting parameters of tails of peaks.
Definition TSpectrumFit.cxx:2711
TSpectrumFit::fXmax
Int_t fXmax
last fitted channel
Definition TSpectrumFit.h:23
TSpectrumFit::fStatisticType
Int_t fStatisticType
type of statistics, possible values kFitOptimChiCounts (chi square statistics with counts as weightin...
Definition TSpectrumFit.h:24
TSpectrumFit::fTErr
Double_t fTErr
error value of t parameter
Definition TSpectrumFit.h:43
TSpectrumFit::Dert
Double_t Dert(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma, Double_t b)
This function calculates derivative of peaks shape function (see manual) according to relative amplit...
Definition TSpectrumFit.cxx:410
TSpectrumFit::Derb
Double_t Derb(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma, Double_t t, Double_t b)
This function calculates derivative of peaks shape function (see manual) according to slope b.
Definition TSpectrumFit.cxx:464
TSpectrumFit::Derpa
Double_t Derpa(Double_t sigma, Double_t t, Double_t b)
This function calculates derivative of the area of peak according to its amplitude.
Definition TSpectrumFit.cxx:593
TSpectrumFit::fSigmaInit
Double_t fSigmaInit
initial value of sigma parameter
Definition TSpectrumFit.h:38
TSpectrumFit::Deramp
Double_t Deramp(Double_t i, Double_t i0, Double_t sigma, Double_t t, Double_t s, Double_t b)
This function calculates derivative of peak shape function (see manual) according to amplitude of pea...
Definition TSpectrumFit.cxx:230
TSpectrumFit::Derderi0
Double_t Derderi0(Double_t i, Double_t amp, Double_t i0, Double_t sigma)
This function calculates second derivative of peak shape function (see manual) according to peak posi...
Definition TSpectrumFit.cxx:304
TSpectrumFit::fA0Err
Double_t fA0Err
error value of background a0 parameter
Definition TSpectrumFit.h:52
TSpectrumFit::fTCalc
Double_t fTCalc
calculated value of t parameter
Definition TSpectrumFit.h:42
TSpectrumFit::Shape
Double_t Shape(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma, Double_t t, Double_t s, Double_t b, Double_t a0, Double_t a1, Double_t a2)
This function calculates peaks shape function (see manual) Function parameters:
Definition TSpectrumFit.cxx:517
TSpectrumFit::fAreaErr
Double_t * fAreaErr
[fNPeaks] array of errors of peak areas
Definition TSpectrumFit.h:37
TSpectrumFit::StiefelInversion
void StiefelInversion(Double_t **a, Int_t rozmer)
This function calculates solution of the system of linear equations.
Definition TSpectrumFit.cxx:1722
TSpectrumFit::Derfc
Double_t Derfc(Double_t x)
This function calculates derivative of error function of x.
Definition TSpectrumFit.cxx:201
TSpectrumFit::SetBackgroundParameters
void SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t a1Init, Bool_t fixA1, Double_t a2Init, Bool_t fixA2)
This function sets the following fitting parameters of background:
Definition TSpectrumFit.cxx:2692
TSpectrumFit::GetBackgroundParameters
void GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &a1, Double_t &a1Err, Double_t &a2, Double_t &a2Err)
This function gets the background parameters and their errors.
Definition TSpectrumFit.cxx:2741
sigma
const Double_t sigma
Definition h1analysisProxy.h:11
x
Double_t x[n]
Definition legend1.C:17
TMath::Log
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition TMath.h:767
TMath::Sqrt
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
TMath::Abs
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124