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