Logo ROOT  

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