Logo ROOT  

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