ROOT  6.06/09
Reference Guide
TSpectrumFit.cxx
Go to the documentation of this file.
1 // @(#)root/spectrum:$Id$
2 // Author: Miroslav Morhac 25/09/06
3 
4 //__________________________________________________________________________
5 // THIS CLASS CONTAINS ADVANCED SPECTRA FITTING FUNCTIONS. //
6 // //
7 // //
8 // These functions were written by: //
9 // Miroslav Morhac //
10 // Institute of Physics //
11 // Slovak Academy of Sciences //
12 // Dubravska cesta 9, 842 28 BRATISLAVA //
13 // SLOVAKIA //
14 // //
15 // email:fyzimiro@savba.sk, fax:+421 7 54772479 //
16 // //
17 // The original code in C has been repackaged as a C++ class by R.Brun //
18 // //
19 // The algorithms in this class have been published in the following //
20 // references: //
21 // [1] M. Morhac et al.: Efficient fitting algorithms applied to //
22 // analysis of coincidence gamma-ray spectra. Computer Physics //
23 // Communications, Vol 172/1 (2005) pp. 19-41. //
24 // //
25 // [2] M. Morhac et al.: Study of fitting algorithms applied to //
26 // simultaneous analysis of large number of peaks in gamma-ray spectra. //
27 // Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003. //
28 // //
29 // //
30 //____________________________________________________________________________
31 
32 /** \class TSpectrumFit
33  \ingroup Hist
34  \brief Advanced 1-dimentional spectra fitting functions
35  \author Miroslav Morhac
36 
37  The original code in C has been repackaged as a C++ class by R.Brun
38 
39  The algorithms in this class have been published in the following
40  references:
41  1. M. Morhac et al.: Efficient fitting algorithms applied to
42  analysis of coincidence gamma-ray spectra. Computer Physics
43  Communications, Vol 172/1 (2005) pp. 19-41.
44 
45  2. M. Morhac et al.: Study of fitting algorithms applied to
46  simultaneous analysis of large number of peaks in gamma-ray spectra.
47  Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003.
48 
49 */
50 
51 
52 #include "TSpectrumFit.h"
53 #include "TMath.h"
54 
55 ClassImp(TSpectrumFit)
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 ///default constructor
59 
60 TSpectrumFit::TSpectrumFit() :TNamed("SpectrumFit", "Miroslav Morhac peak fitter")
61 {
62  fNPeaks = 0;
63  fNumberIterations = 1;
64  fXmin = 0;
65  fXmax = 100;
66  fStatisticType = kFitOptimChiCounts;
67  fAlphaOptim = kFitAlphaHalving;
68  fPower = kFitPower2;
69  fFitTaylor = kFitTaylorOrderFirst;
70  fAlpha =1;
71  fChi = 0;
72  fPositionInit = 0;
73  fPositionCalc = 0;
74  fPositionErr = 0;
75  fFixPosition = 0;
76  fAmpInit = 0;
77  fAmpCalc = 0;
78  fAmpErr = 0;
79  fFixAmp = 0;
80  fArea = 0;
81  fAreaErr = 0;
82  fSigmaInit = 2;
83  fSigmaCalc = 1;
84  fSigmaErr = 0;
85  fTInit = 0;
86  fTCalc = 0;
87  fTErr = 0;
88  fBInit = 1;
89  fBCalc = 0;
90  fBErr = 0;
91  fSInit = 0;
92  fSCalc = 0;
93  fSErr = 0;
94  fA0Init = 0;
95  fA0Calc = 0;
96  fA0Err = 0;
97  fA1Init = 0;
98  fA1Calc = 0;
99  fA1Err = 0;
100  fA2Init = 0;
101  fA2Calc = 0;
102  fA2Err = 0;
103  fFixSigma = false;
104  fFixT = true;
105  fFixB = true;
106  fFixS = true;
107  fFixA0 = true;
108  fFixA1 = true;
109  fFixA2 = true;
110 }
111 
112 ////////////////////////////////////////////////////////////////////////////////
113 ///numberPeaks: number of fitted peaks (must be greater than zero)
114 ///the constructor allocates arrays for all fitted parameters (peak positions, amplitudes etc) and sets the member
115 ///variables to their default values. One can change these variables by member functions (setters) of TSpectrumFit class.
116 ///Begin_Html <!--
117 
118 TSpectrumFit::TSpectrumFit(Int_t numberPeaks) :TNamed("SpectrumFit", "Miroslav Morhac peak fitter")
119 {
120 /* -->
121 <div class=Section1>
122 
123 <p class=MsoNormal style='text-align:justify'>Shape function of the fitted
124 peaks is </p>
125 
126 <p class=MsoNormal style='text-align:justify'>
127 
128 <table cellpadding=0 cellspacing=0 align=left>
129  <tr>
130  <td width=68 height=6></td>
131  </tr>
132  <tr>
133  <td></td>
134  <td><img width=388 height=132 src="gif/spectrumfit_constructor_image001.gif"></td>
135  </tr>
136 </table>
137 
138 <span style='font-family:Arial'>&nbsp;</span></p>
139 
140 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
141 
142 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
143 
144 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
145 
146 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
147 
148 <p class=MsoNormal><i>&nbsp;</i></p>
149 
150 <p class=MsoNormal><i>&nbsp;</i></p>
151 
152 <p class=MsoNormal><i>&nbsp;</i></p>
153 
154 <br clear=ALL>
155 
156 <p class=MsoNormal style='text-align:justify'>where a represents vector of
157 fitted parameters (positions p(j), amplitudes A(j), sigma, relative amplitudes
158 T, S and slope B).</p>
159 
160 <p class=MsoNormal><span style='font-size:16.0pt'>&nbsp;</span></p>
161 
162 </div>
163 
164 <!-- */
165 // --> End_Html
166 
167  if (numberPeaks <= 0){
168  Error ("TSpectrumFit","Invalid number of peaks, must be > than 0");
169  return;
170  }
171  fNPeaks = numberPeaks;
172  fNumberIterations = 1;
173  fXmin = 0;
174  fXmax = 100;
175  fStatisticType = kFitOptimChiCounts;
176  fAlphaOptim = kFitAlphaHalving;
177  fPower = kFitPower2;
178  fFitTaylor = kFitTaylorOrderFirst;
179  fAlpha =1;
180  fChi = 0;
181  fPositionInit = new Double_t[numberPeaks];
182  fPositionCalc = new Double_t[numberPeaks];
183  fPositionErr = new Double_t[numberPeaks];
184  fFixPosition = new Bool_t[numberPeaks];
185  fAmpInit = new Double_t[numberPeaks];
186  fAmpCalc = new Double_t[numberPeaks];
187  fAmpErr = new Double_t[numberPeaks];
188  fFixAmp = new Bool_t[numberPeaks];
189  fArea = new Double_t[numberPeaks];
190  fAreaErr = new Double_t[numberPeaks];
191  fSigmaInit = 2;
192  fSigmaCalc = 1;
193  fSigmaErr = 0;
194  fTInit = 0;
195  fTCalc = 0;
196  fTErr = 0;
197  fBInit = 1;
198  fBCalc = 0;
199  fBErr = 0;
200  fSInit = 0;
201  fSCalc = 0;
202  fSErr = 0;
203  fA0Init = 0;
204  fA0Calc = 0;
205  fA0Err = 0;
206  fA1Init = 0;
207  fA1Calc = 0;
208  fA1Err = 0;
209  fA2Init = 0;
210  fA2Calc = 0;
211  fA2Err = 0;
212  fFixSigma = false;
213  fFixT = true;
214  fFixB = true;
215  fFixS = true;
216  fFixA0 = true;
217  fFixA1 = true;
218  fFixA2 = true;
219 }
220 
221 
222 
223 ////////////////////////////////////////////////////////////////////////////////
224 ///destructor
225 
226 TSpectrumFit::~TSpectrumFit()
227 {
228  delete [] fPositionInit;
229  delete [] fPositionCalc;
230  delete [] fPositionErr;
231  delete [] fFixPosition;
232  delete [] fAmpInit;
233  delete [] fAmpCalc;
234  delete [] fAmpErr;
235  delete [] fFixAmp;
236  delete [] fArea;
237  delete [] fAreaErr;
238 }
239 
240 //_____________________________________________________________________________
241 /////////////////BEGINNING OF AUXILIARY FUNCTIONS USED BY FITTING FUNCTION Fit1//////////////////////////
242 Double_t TSpectrumFit::Erfc(Double_t x)
243 {
244 //////////////////////////////////////////////////////////////////////////////
245 // AUXILIARY FUNCTION //
246 // //
247 // This function calculates error function of x. //
248 // //
249 //////////////////////////////////////////////////////////////////////////////
250  Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap = 0.47047;
251  Double_t a, t, c, w;
252  a = TMath::Abs(x);
253  w = 1. + dap * a;
254  t = 1. / w;
255  w = a * a;
256  if (w < 700)
257  c = exp(-w);
258 
259  else {
260  c = 0;
261  }
262  c = c * t * (da1 + t * (da2 + t * da3));
263  if (x < 0)
264  c = 1. - c;
265  return (c);
266 }
267 
268 ////////////////////////////////////////////////////////////////////////////////
269 ///////////////////////////////////////////////////////////////////////////////
270 /// AUXILIARY FUNCTION //
271 /// //
272 /// This function calculates derivative of error function of x. //
273 /// //
274 ///////////////////////////////////////////////////////////////////////////////
275 
276 Double_t TSpectrumFit::Derfc(Double_t x)
277 {
278  Double_t a, t, c, w;
279  Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap = 0.47047;
280  a = TMath::Abs(x);
281  w = 1. + dap * a;
282  t = 1. / w;
283  w = a * a;
284  if (w < 700)
285  c = exp(-w);
286 
287  else {
288  c = 0;
289  }
290  c = (-1.) * dap * c * t * t * (da1 + t * (2. * da2 + t * 3. * da3)) -
291  2. * a * Erfc(a);
292  return (c);
293 }
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 ///////////////////////////////////////////////////////////////////////////////
297 /// AUXILIARY FUNCTION //
298 /// //
299 /// This function calculates derivative of peak shape function (see manual) //
300 /// according to amplitude of peak. //
301 /// Function parameters: //
302 /// -i-channel //
303 /// -i0-position of peak //
304 /// -sigma-sigma of peak //
305 /// -t, s-relative amplitudes //
306 /// -b-slope //
307 /// //
308 ///////////////////////////////////////////////////////////////////////////////
309 
310 Double_t TSpectrumFit::Deramp(Double_t i, Double_t i0, Double_t sigma, Double_t t,
311  Double_t s, Double_t b)
312 {
313  Double_t p, q, r, a;
314  p = (i - i0) / sigma;
315  if ((p * p) < 700)
316  q = exp(-p * p);
317 
318  else {
319  q = 0;
320  }
321  r = 0;
322  if (t != 0) {
323  a = p / b;
324  if (a > 700)
325  a = 700;
326  r = t * exp(a) / 2.;
327  }
328  if (r != 0)
329  r = r * Erfc(p + 1. / (2. * b));
330  q = q + r;
331  if (s != 0)
332  q = q + s * Erfc(p) / 2.;
333  return (q);
334 }
335 
336 ////////////////////////////////////////////////////////////////////////////////
337 ///////////////////////////////////////////////////////////////////////////////
338 /// AUXILIARY FUNCTION //
339 /// //
340 /// This function calculates derivative of peak shape function (see manual) //
341 /// according to peak position. //
342 /// Function parameters: //
343 /// -i-channel //
344 /// -amp-amplitude of peak //
345 /// -i0-position of peak //
346 /// -sigma-sigma of peak //
347 /// -t, s-relative amplitudes //
348 /// -b-slope //
349 /// //
350 ///////////////////////////////////////////////////////////////////////////////
351 
352 Double_t TSpectrumFit::Deri0(Double_t i, Double_t amp, Double_t i0, Double_t sigma,
353  Double_t t, Double_t s, Double_t b)
354 {
355  Double_t p, r1, r2, r3, r4, c, d, e;
356  p = (i - i0) / sigma;
357  d = 2. * sigma;
358  if ((p * p) < 700)
359  r1 = 2. * p * exp(-p * p) / sigma;
360 
361  else {
362  r1 = 0;
363  }
364  r2 = 0, r3 = 0;
365  if (t != 0) {
366  c = p + 1. / (2. * b);
367  e = p / b;
368  if (e > 700)
369  e = 700;
370  r2 = -t * exp(e) * Erfc(c) / (d * b);
371  r3 = -t * exp(e) * Derfc(c) / d;
372  }
373  r4 = 0;
374  if (s != 0)
375  r4 = -s * Derfc(p) / d;
376  r1 = amp * (r1 + r2 + r3 + r4);
377  return (r1);
378 }
379 
380 ////////////////////////////////////////////////////////////////////////////////
381 ///////////////////////////////////////////////////////////////////////////////
382 /// AUXILIARY FUNCTION //
383 /// //
384 /// This function calculates second derivative of peak shape function //
385 /// (see manual) according to peak position. //
386 /// Function parameters: //
387 /// -i-channel //
388 /// -amp-amplitude of peak //
389 /// -i0-position of peak //
390 /// -sigma-width of peak //
391 /// //
392 ///////////////////////////////////////////////////////////////////////////////
393 
394 Double_t TSpectrumFit::Derderi0(Double_t i, Double_t amp, Double_t i0,
395  Double_t sigma)
396 {
397  Double_t p, r1, r2, r3, r4;
398  p = (i - i0) / sigma;
399  if ((p * p) < 700)
400  r1 = exp(-p * p);
401 
402  else {
403  r1 = 0;
404  }
405  r1 = r1 * (4 * p * p - 2) / (sigma * sigma);
406  r2 = 0, r3 = 0, r4 = 0;
407  r1 = amp * (r1 + r2 + r3 + r4);
408  return (r1);
409 }
410 
411 ////////////////////////////////////////////////////////////////////////////////
412 ///////////////////////////////////////////////////////////////////////////////////
413 /// AUXILIARY FUNCTION //
414 /// //
415 /// This function calculates derivative of peaks shape function (see manual) //
416 /// according to sigma of peaks. //
417 /// Function parameters: //
418 /// -num_of_fitted_peaks-number of fitted peaks //
419 /// -i-channel //
420 /// -parameter-array of peaks parameters (amplitudes and positions) //
421 /// -sigma-sigma of peak //
422 /// -t, s-relative amplitudes //
423 /// -b-slope //
424 /// //
425 ///////////////////////////////////////////////////////////////////////////////////
426 
427 Double_t TSpectrumFit::Dersigma(Int_t num_of_fitted_peaks, Double_t i,
428  const Double_t *parameter, Double_t sigma,
429  Double_t t, Double_t s, Double_t b)
430 {
431  Int_t j;
432  Double_t r, p, r1, r2, r3, r4, c, d, e;
433  r = 0;
434  d = 2. * sigma;
435  for (j = 0; j < num_of_fitted_peaks; j++) {
436  p = (i - parameter[2 * j + 1]) / sigma;
437  r1 = 0;
438  if (TMath::Abs(p) < 3) {
439  if ((p * p) < 700)
440  r1 = 2. * p * p * exp(-p * p) / sigma;
441 
442  else {
443  r1 = 0;
444  }
445  }
446  r2 = 0, r3 = 0;
447  if (t != 0) {
448  c = p + 1. / (2. * b);
449  e = p / b;
450  if (e > 700)
451  e = 700;
452  r2 = -t * p * exp(e) * Erfc(c) / (d * b);
453  r3 = -t * p * exp(e) * Derfc(c) / d;
454  }
455  r4 = 0;
456  if (s != 0)
457  r4 = -s * p * Derfc(p) / d;
458  r = r + parameter[2 * j] * (r1 + r2 + r3 + r4);
459  }
460  return (r);
461 }
462 
463 ////////////////////////////////////////////////////////////////////////////////
464 ///////////////////////////////////////////////////////////////////////////////////
465 /// AUXILIARY FUNCTION //
466 /// //
467 /// This function calculates second derivative of peaks shape function //
468 /// (see manual) according to sigma of peaks. //
469 /// Function parameters: //
470 /// -num_of_fitted_peaks-number of fitted peaks //
471 /// -i-channel //
472 /// -parameter-array of peaks parameters (amplitudes and positions) //
473 /// -sigma-sigma of peak //
474 /// //
475 ///////////////////////////////////////////////////////////////////////////////////
476 
477 Double_t TSpectrumFit::Derdersigma(Int_t num_of_fitted_peaks, Double_t i,
478  const Double_t *parameter, Double_t sigma)
479 {
480  Int_t j;
481  Double_t r, p, r1, r2, r3, r4;
482  r = 0;
483  for (j = 0; j < num_of_fitted_peaks; j++) {
484  p = (i - parameter[2 * j + 1]) / sigma;
485  r1 = 0;
486  if (TMath::Abs(p) < 3) {
487  if ((p * p) < 700)
488  r1 = exp(-p * p) * p * p * (4. * p * p - 6) / (sigma * sigma);
489 
490  else {
491  r1 = 0;
492  }
493  }
494  r2 = 0, r3 = 0, r4 = 0;
495  r = r + parameter[2 * j] * (r1 + r2 + r3 + r4);
496  }
497  return (r);
498 }
499 
500 ////////////////////////////////////////////////////////////////////////////////
501 ///////////////////////////////////////////////////////////////////////////////////
502 /// AUXILIARY FUNCTION //
503 /// //
504 /// This function calculates derivative of peaks shape function (see manual) //
505 /// according to relative amplitude t. //
506 /// Function parameters: //
507 /// -num_of_fitted_peaks-number of fitted peaks //
508 /// -i-channel //
509 /// -parameter-array of peaks parameters (amplitudes and positions) //
510 /// -sigma-sigma of peak //
511 /// -b-slope //
512 /// //
513 ///////////////////////////////////////////////////////////////////////////////////
514 
515 Double_t TSpectrumFit::Dert(Int_t num_of_fitted_peaks, Double_t i,
516  const Double_t *parameter, Double_t sigma, Double_t b)
517 {
518  Int_t j;
519  Double_t r, p, r1, c, e;
520  r = 0;
521  for (j = 0; j < num_of_fitted_peaks; j++) {
522  p = (i - parameter[2 * j + 1]) / sigma;
523  c = p + 1. / (2. * b);
524  e = p / b;
525  if (e > 700)
526  e = 700;
527  r1 = exp(e) * Erfc(c);
528  r = r + parameter[2 * j] * r1;
529  }
530  r = r / 2.;
531  return (r);
532 }
533 
534 ////////////////////////////////////////////////////////////////////////////////
535 ///////////////////////////////////////////////////////////////////////////////////
536 /// AUXILIARY FUNCTION //
537 /// //
538 /// This function calculates derivative of peaks shape function (see manual) //
539 /// according to relative amplitude s. //
540 /// Function parameters: //
541 /// -num_of_fitted_peaks-number of fitted peaks //
542 /// -i-channel //
543 /// -parameter-array of peaks parameters (amplitudes and positions) //
544 /// -sigma-sigma of peak //
545 /// //
546 ///////////////////////////////////////////////////////////////////////////////////
547 
548 Double_t TSpectrumFit::Ders(Int_t num_of_fitted_peaks, Double_t i,
549  const Double_t *parameter, Double_t sigma)
550 {
551  Int_t j;
552  Double_t r, p, r1;
553  r = 0;
554  for (j = 0; j < num_of_fitted_peaks; j++) {
555  p = (i - parameter[2 * j + 1]) / sigma;
556  r1 = Erfc(p);
557  r = r + parameter[2 * j] * r1;
558  }
559  r = r / 2.;
560  return (r);
561 }
562 
563 ////////////////////////////////////////////////////////////////////////////////
564 ///////////////////////////////////////////////////////////////////////////////////
565 /// AUXILIARY FUNCTION //
566 /// //
567 /// This function calculates derivative of peaks shape function (see manual) //
568 /// according to slope b. //
569 /// Function parameters: //
570 /// -num_of_fitted_peaks-number of fitted peaks //
571 /// -i-channel //
572 /// -parameter-array of peaks parameters (amplitudes and positions) //
573 /// -sigma-sigma of peak //
574 /// -t-relative amplitude //
575 /// -b-slope //
576 /// //
577 ///////////////////////////////////////////////////////////////////////////////////
578 
579 Double_t TSpectrumFit::Derb(Int_t num_of_fitted_peaks, Double_t i,
580  const Double_t *parameter, Double_t sigma, Double_t t,
581  Double_t b)
582 {
583  Int_t j;
584  Double_t r, p, r1, c, e;
585  r = 0;
586  for (j = 0; j < num_of_fitted_peaks && t != 0; j++) {
587  p = (i - parameter[2 * j + 1]) / sigma;
588  c = p + 1. / (2. * b);
589  e = p / b;
590  r1 = p * Erfc(c);
591  r1 = r1 + Derfc(c) / 2.;
592  if (e > 700)
593  e = 700;
594  if (e < -700)
595  r1 = 0;
596 
597  else
598  r1 = r1 * exp(e);
599  r = r + parameter[2 * j] * r1;
600  }
601  r = -r * t / (2. * b * b);
602  return (r);
603 }
604 
605 ////////////////////////////////////////////////////////////////////////////////
606 ///derivative of background according to a1
607 
608 Double_t TSpectrumFit::Dera1(Double_t i)
609 {
610  return (i);
611 }
612 
613 ////////////////////////////////////////////////////////////////////////////////
614 ///derivative of background according to a2
615 
616 Double_t TSpectrumFit::Dera2(Double_t i)
617 {
618  return (i * i);
619 }
620 
621 ////////////////////////////////////////////////////////////////////////////////
622 ///////////////////////////////////////////////////////////////////////////////////
623 /// AUXILIARY FUNCTION //
624 /// //
625 /// This function calculates peaks shape function (see manual) //
626 /// Function parameters: //
627 /// -num_of_fitted_peaks-number of fitted peaks //
628 /// -i-channel //
629 /// -parameter-array of peaks parameters (amplitudes and positions) //
630 /// -sigma-sigma of peak //
631 /// -t, s-relative amplitudes //
632 /// -b-slope //
633 /// -a0, a1, a2- background coefficients //
634 /// //
635 ///////////////////////////////////////////////////////////////////////////////////
636 
637 Double_t TSpectrumFit::Shape(Int_t num_of_fitted_peaks, Double_t i,
638  const Double_t *parameter, Double_t sigma, Double_t t,
639  Double_t s, Double_t b, Double_t a0, Double_t a1,
640  Double_t a2)
641 {
642  Int_t j;
643  Double_t r, p, r1, r2, r3, c, e;
644  r = 0;
645  for (j = 0; j < num_of_fitted_peaks; j++) {
646  if (sigma > 0.0001)
647  p = (i - parameter[2 * j + 1]) / sigma;
648 
649  else {
650  if (i == parameter[2 * j + 1])
651  p = 0;
652 
653  else
654  p = 10;
655  }
656  r1 = 0;
657  if (TMath::Abs(p) < 3) {
658  if ((p * p) < 700)
659  r1 = exp(-p * p);
660 
661  else {
662  r1 = 0;
663  }
664  }
665  r2 = 0;
666  if (t != 0) {
667  c = p + 1. / (2. * b);
668  e = p / b;
669  if (e > 700)
670  e = 700;
671  r2 = t * exp(e) * Erfc(c) / 2.;
672  }
673  r3 = 0;
674  if (s != 0)
675  r3 = s * Erfc(p) / 2.;
676  r = r + parameter[2 * j] * (r1 + r2 + r3);
677  }
678  r = r + a0 + a1 * i + a2 * i * i;
679  return (r);
680 }
681 
682 ////////////////////////////////////////////////////////////////////////////////
683 ///////////////////////////////////////////////////////////////////////////////////
684 /// AUXILIARY FUNCTION //
685 /// //
686 /// This function calculates area of a peak //
687 /// Function parameters: //
688 /// -a-amplitude of the peak //
689 /// -sigma-sigma of peak //
690 /// -t-relative amplitude //
691 /// -b-slope //
692 /// //
693 ///////////////////////////////////////////////////////////////////////////////////
694 
695 Double_t TSpectrumFit::Area(Double_t a, Double_t sigma, Double_t t, Double_t b)
696 {
697  Double_t odm_pi = 1.7724538, r = 0;
698  if (b != 0)
699  r = 0.5 / b;
700  r = (-1.) * r * r;
701  if (TMath::Abs(r) < 700)
702  r = a * sigma * (odm_pi + t * b * exp(r));
703 
704  else {
705  r = a * sigma * odm_pi;
706  }
707  return (r);
708 }
709 
710 ////////////////////////////////////////////////////////////////////////////////
711 ///////////////////////////////////////////////////////////////////////////////////
712 /// AUXILIARY FUNCTION //
713 /// //
714 /// This function calculates derivative of the area of peak //
715 /// according to its amplitude. //
716 /// Function parameters: //
717 /// -sigma-sigma of peak //
718 /// -t-relative amplitudes //
719 /// -b-slope //
720 /// //
721 ///////////////////////////////////////////////////////////////////////////////////
722 
723 Double_t TSpectrumFit::Derpa(Double_t sigma, Double_t t, Double_t b)
724 {
725  Double_t odm_pi = 1.7724538, r;
726  r = 0.5 / b;
727  r = (-1.) * r * r;
728  if (TMath::Abs(r) < 700)
729  r = sigma * (odm_pi + t * b * exp(r));
730 
731  else {
732  r = sigma * odm_pi;
733  }
734  return (r);
735 }
736 Double_t TSpectrumFit::Derpsigma(Double_t a, Double_t t, Double_t b)
737 {
738 //////////////////////////////////////////////////////////////////////////////////
739 // AUXILIARY FUNCTION //
740 // //
741 // This function calculates derivative of the area of peak //
742 // according to sigma of peaks. //
743 // Function parameters: //
744 // -a-amplitude of peak //
745 // -t-relative amplitudes //
746 // -b-slope //
747 // //
748 //////////////////////////////////////////////////////////////////////////////////
749  Double_t odm_pi = 1.7724538, r;
750  r = 0.5 / b;
751  r = (-1.) * r * r;
752  if (TMath::Abs(r) < 700)
753  r = a * (odm_pi + t * b * exp(r));
754 
755  else {
756  r = a * odm_pi;
757  }
758  return (r);
759 }
760 
761 ////////////////////////////////////////////////////////////////////////////////
762 ///////////////////////////////////////////////////////////////////////////////////
763 /// AUXILIARY FUNCTION //
764 /// //
765 /// This function calculates derivative of the area of peak //
766 /// according to t parameter. //
767 /// Function parameters: //
768 /// -sigma-sigma of peak //
769 /// -t-relative amplitudes //
770 /// -b-slope //
771 /// //
772 ///////////////////////////////////////////////////////////////////////////////////
773 
774 Double_t TSpectrumFit::Derpt(Double_t a, Double_t sigma, Double_t b)
775 {
776  Double_t r;
777  r = 0.5 / b;
778  r = (-1.) * r * r;
779  if (TMath::Abs(r) < 700)
780  r = a * sigma * b * exp(r);
781 
782  else {
783  r = 0;
784  }
785  return (r);
786 }
787 
788 ////////////////////////////////////////////////////////////////////////////////
789 ///////////////////////////////////////////////////////////////////////////////////
790 /// AUXILIARY FUNCTION //
791 /// //
792 /// This function calculates derivative of the area of peak //
793 /// according to b parameter. //
794 /// Function parameters: //
795 /// -sigma-sigma of peak //
796 /// -t-relative amplitudes //
797 /// -b-slope //
798 /// //
799 ///////////////////////////////////////////////////////////////////////////////////
800 
801 Double_t TSpectrumFit::Derpb(Double_t a, Double_t sigma, Double_t t, Double_t b)
802 {
803  Double_t r;
804  r = (-1) * 0.25 / (b * b);
805  if (TMath::Abs(r) < 700)
806  r = a * sigma * t * exp(r) * (1 - 2 * r);
807 
808  else {
809  r = 0;
810  }
811  return (r);
812 }
813 
814 ////////////////////////////////////////////////////////////////////////////////
815 ///power function
816 
817 Double_t TSpectrumFit::Ourpowl(Double_t a, Int_t pw)
818 {
819  Double_t c;
820  Double_t a2 = a*a;
821  c = 1;
822  if (pw > 0) c *= a2;
823  if (pw > 2) c *= a2;
824  if (pw > 4) c *= a2;
825  if (pw > 6) c *= a2;
826  if (pw > 8) c *= a2;
827  if (pw > 10) c *= a2;
828  if (pw > 12) c *= a2;
829  return (c);
830 }
831 
832 /////////////////END OF AUXILIARY FUNCTIONS USED BY FITTING FUNCTIONS FitAWMI, FitStiefel//////////////////////////
833 /////////////////FITTING FUNCTION WITHOUT MATRIX INVERSION///////////////////////////////////////
834 
835 ////////////////////////////////////////////////////////////////////////////////
836 //////////////////////////////////////////////////////////////////////////////
837 /// ONE-DIMENSIONAL FIT FUNCTION
838 /// ALGORITHM WITHOUT MATRIX INVERSION
839 /// This function fits the source spectrum. The calling program should
840 /// fill in input parameters of the TSpectrumFit class
841 /// The fitted parameters are written into
842 /// TSpectrumFit class output parameters and fitted data are written into
843 /// source spectrum.
844 ///
845 /// Function parameters:
846 /// source-pointer to the vector of source spectrum
847 ///
848 //////////////////////////////////////////////////////////////////////////////
849 ///
850 ///Begin_Html <!--
851 
852 void TSpectrumFit::FitAwmi(Double_t *source)
853 {
854 /* -->
855 <div class=Section2>
856 
857 <p class=MsoNormal><b><span style='font-size:14.0pt'>Fitting</span></b></p>
858 
859 <p class=MsoNormal style='text-align:justify'><i>Goal: to estimate
860 simultaneously peak shape parameters in spectra with large number of peaks</i></p>
861 
862 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
863 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
864 </span>peaks can be fitted separately, each peak (or multiplets) in a region or
865 together all peaks in a spectrum. To fit separately each peak one needs to
866 determine the fitted region. However it can happen that the regions of
867 neighboring peaks are overlapping. Then the results of fitting are very poor.
868 On the other hand, when fitting together all peaks found in a  spectrum, one
869 needs to have a method that is  stable (converges) and fast enough to carry out
870 fitting in reasonable time </p>
871 
872 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
873 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
874 </span>we have implemented the nonsymmetrical semiempirical peak shape function
875 [1]</p>
876 
877 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
878 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
879 </span>it contains the symmetrical Gaussian as well as nonsymmetrical terms.</p>
880 
881 <p class=MsoNormal style='text-align:justify'>
882 
883 <table cellpadding=0 cellspacing=0 align=left>
884  <tr>
885  <td width=84 height=18></td>
886  </tr>
887  <tr>
888  <td></td>
889  <td><img width=372 height=127 src="gif/spectrumfit_awni_image001.gif"></td>
890  </tr>
891 </table>
892 
893 <span style='font-size:16.0pt'>&nbsp;</span></p>
894 
895 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>&nbsp;</span></p>
896 
897 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>&nbsp;</span></p>
898 
899 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>&nbsp;</span></p>
900 
901 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>&nbsp;</span></p>
902 
903 <br clear=ALL>
904 
905 <p class=MsoNormal style='text-indent:34.2pt'>where T and S are relative amplitudes
906 and B is slope.</p>
907 
908 <p class=MsoNormal>&nbsp;</p>
909 
910 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
911 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
912 </span>algorithm without matrix inversion (AWMI) allows fitting tens, hundreds
913 of peaks simultaneously that represent sometimes thousands of parameters [2],
914 [5]. </p>
915 
916 <p class=MsoNormal><i>Function:</i></p>
917 
918 <p class=MsoNormal style='text-align:justify'>void <a
919 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrumFit::FitAwmi</b></a>(<a
920 href="http://root.cern.ch/root/html/ListOfTypes.html#double"><b>double</b></a> *fSource)
921 </p>
922 
923 <p class=MsoNormal style='text-align:justify'>This function fits the source
924 spectrum using AWMI algorithm. The calling program should fill in input fitting
925 parameters of the TSpectrumFit class using a set of TSpectrumFit setters. The
926 fitted parameters are written into the class and the fitted data are written
927 into source spectrum. </p>
928 
929 <p class=MsoNormal>&nbsp;</p>
930 
931 <p class=MsoNormal><i><span style='color:red'>Parameter:</span></i></p>
932 
933 <p class=MsoNormal style='text-align:justify'>        <b>fSource</b>-pointer to
934 the vector of source spectrum                  </p>
935 
936 <p class=MsoNormal style='text-align:justify'>        </p>
937 
938 <p class=MsoNormal><i><span style='color:red'>Member variables of the
939 TSpectrumFit class:</span></i></p>
940 
941 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
942 Int_t     fNPeaks;                    //number of peaks present in fit, input
943 parameter, it should be &gt; 0</span></p>
944 
945 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
946 Int_t     fNumberIterations;          //number of iterations in fitting
947 procedure, input parameter, it should be &gt; 0</span></p>
948 
949 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
950 Int_t     fXmin;                      //first fitted channel</span></p>
951 
952 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
953 Int_t     fXmax;                      //last fitted channel</span></p>
954 
955 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
956 Int_t     fStatisticType;             //type of statistics, possible values
957 kFitOptimChiCounts (chi square statistics with counts as weighting
958 coefficients), kFitOptimChiFuncValues (chi square statistics with function
959 values as weighting coefficients),kFitOptimMaxLikelihood</span></p>
960 
961 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
962 Int_t     fAlphaOptim;                //optimization of convergence algorithm, possible
963 values kFitAlphaHalving, kFitAlphaOptimal</span></p>
964 
965 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
966 Int_t     fPower;                     //possible values kFitPower2,4,6,8,10,12,
967 for details see references. It applies only for Awmi fitting function.</span></p>
968 
969 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
970 Int_t     fFitTaylor;                 //order of Taylor expansion, possible
971 values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi
972 fitting function.</span></p>
973 
974 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
975 Double_t  fAlpha;                     //convergence coefficient, input
976 parameter, it should be positive number and &lt;=1, for details see references</span></p>
977 
978 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
979 Double_t  fChi;                       //here the fitting functions return
980 resulting chi square   </span></p>
981 
982 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
983 Double_t *fPositionInit;              //[fNPeaks] array of initial values of
984 peaks positions, input parameters</span></p>
985 
986 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
987 Double_t *fPositionCalc;              //[fNPeaks] array of calculated values of
988 fitted positions, output parameters</span></p>
989 
990 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
991 Double_t *fPositionErr;               //[fNPeaks] array of position errors</span></p>
992 
993 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
994 Double_t *fAmpInit;                   //[fNPeaks] array of initial values of
995 peaks amplitudes, input parameters</span></p>
996 
997 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
998 Double_t *fAmpCalc;                   //[fNPeaks] array of calculated values of
999 fitted amplitudes, output parameters</span></p>
1000 
1001 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1002 Double_t *fAmpErr;                    //[fNPeaks] array of amplitude errors</span></p>
1003 
1004 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1005 Double_t *fArea;                      //[fNPeaks] array of calculated areas of
1006 peaks</span></p>
1007 
1008 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1009 Double_t *fAreaErr;                   //[fNPeaks] array of errors of peak areas</span></p>
1010 
1011 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1012 Double_t  fSigmaInit;                 //initial value of sigma parameter</span></p>
1013 
1014 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1015 Double_t  fSigmaCalc;                 //calculated value of sigma parameter</span></p>
1016 
1017 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1018 Double_t  fSigmaErr;                  //error value of sigma parameter</span></p>
1019 
1020 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1021 Double_t  fTInit;                     //initial value of t parameter (relative
1022 amplitude of tail), for details see html manual and references</span></p>
1023 
1024 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1025 Double_t  fTCalc;                     //calculated value of t parameter</span></p>
1026 
1027 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1028 Double_t  fTErr;                      //error value of t parameter</span></p>
1029 
1030 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1031 Double_t  fBInit;                     //initial value of b parameter (slope),
1032 for details see html manual and references</span></p>
1033 
1034 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1035 Double_t  fBCalc;                     //calculated value of b parameter</span></p>
1036 
1037 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1038 Double_t  fBErr;                      //error value of b parameter</span></p>
1039 
1040 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1041 Double_t  fSInit;                     //initial value of s parameter (relative
1042 amplitude of step), for details see html manual and references</span></p>
1043 
1044 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1045 Double_t  fSCalc;                     //calculated value of s parameter</span></p>
1046 
1047 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1048 Double_t  fSErr;                      //error value of s parameter</span></p>
1049 
1050 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1051 Double_t  fA0Init;                    //initial value of background a0
1052 parameter(backgroud is estimated as a0+a1*x+a2*x*x)</span></p>
1053 
1054 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1055 Double_t  fA0Calc;                    //calculated value of background a0
1056 parameter</span></p>
1057 
1058 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1059 Double_t  fA0Err;                     //error value of background a0 parameter</span></p>
1060 
1061 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1062 Double_t  fA1Init;                    //initial value of background a1
1063 parameter(backgroud is estimated as a0+a1*x+a2*x*x)</span></p>
1064 
1065 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1066 Double_t  fA1Calc;                    //calculated value of background a1
1067 parameter</span></p>
1068 
1069 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1070 Double_t  fA1Err;                     //error value of background a1 parameter</span></p>
1071 
1072 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1073 Double_t  fA2Init;                    //initial value of background a2
1074 parameter(backgroud is estimated as a0+a1*x+a2*x*x)</span></p>
1075 
1076 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1077 Double_t  fA2Calc;                    //calculated value of background a2
1078 parameter</span></p>
1079 
1080 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1081 Double_t  fA2Err;                     //error value of background a2 parameter</span></p>
1082 
1083 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1084 Bool_t   *fFixPosition;               //[fNPeaks] array of logical values which
1085 allow to fix appropriate positions (not fit). However they are present in the
1086 estimated functional   </span></p>
1087 
1088 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1089 Bool_t   *fFixAmp;                    //[fNPeaks] array of logical values which
1090 allow to fix appropriate amplitudes (not fit). However they are present in the
1091 estimated functional      </span></p>
1092 
1093 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1094 Bool_t    fFixSigma;                  //logical value of sigma parameter, which
1095 allows to fix the parameter (not to fit).   </span></p>
1096 
1097 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1098 Bool_t    fFixT;                      //logical value of t parameter, which
1099 allows to fix the parameter (not to fit).      </span></p>
1100 
1101 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1102 Bool_t    fFixB;                      //logical value of b parameter, which
1103 allows to fix the parameter (not to fit).   </span></p>
1104 
1105 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1106 Bool_t    fFixS;                      //logical value of s parameter, which
1107 allows to fix the parameter (not to fit).      </span></p>
1108 
1109 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1110 Bool_t    fFixA0;                     //logical value of a0 parameter, which
1111 allows to fix the parameter (not to fit).</span></p>
1112 
1113 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1114 Bool_t    fFixA1;                     //logical value of a1 parameter, which
1115 allows to fix the parameter (not to fit).   </span></p>
1116 
1117 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
1118 Bool_t    fFixA2;                     //logical value of a2 parameter, which
1119 allows to fix the parameter (not to fit).</span></p>
1120 
1121 <p class=MsoNormal style='text-align:justify'><b><i>&nbsp;</i></b></p>
1122 
1123 <p class=MsoNormal style='text-align:justify'><b><i>References:</i></b></p>
1124 
1125 <p class=MsoNormal style='text-align:justify'>[1] Phillps G.W., Marlow K.W.,
1126 NIM 137 (1976) 525.</p>
1127 
1128 <p class=MsoNormal style='text-align:justify'>[2] I. A. Slavic: Nonlinear
1129 least-squares fitting without matrix inversion applied to complex Gaussian
1130 spectra analysis. NIM 134 (1976) 285-289.</p>
1131 
1132 <p class=MsoNormal style='text-align:justify'>[3] T. Awaya: A new method for
1133 curve fitting to the data with low statistics not using chi-square method. NIM
1134 165 (1979) 317-323.</p>
1135 
1136 <p class=MsoNormal style='text-align:justify'>[4] T. Hauschild, M. Jentschel:
1137 Comparison of maximum likelihood estimation and chi-square statistics applied
1138 to counting experiments. NIM A 457 (2001) 384-401.</p>
1139 
1140 <p class=MsoNormal style='text-align:justify'> [5]  M. Morháč,  J.
1141 Kliman,  M. Jandel,  &#317;. Krupa, V. Matoušek: Study of fitting algorithms
1142 applied to simultaneous analysis of large number of peaks in -ray spectra. <span
1143 lang=EN-GB>Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003</span></p>
1144 
1145 <p class=MsoNormal style='text-align:justify'> </p>
1146 
1147 <p class=MsoNormal style='text-align:justify'><i>Example  – script FitAwmi.c:</i></p>
1148 
1149 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
1150 border=0 width=601 height=402 src="gif/spectrumfit_awni_image002.jpg"></span></i></p>
1151 
1152 <p class=MsoNormal style='text-align:justify'><b>Fig. 1 Original spectrum
1153 (black line) and fitted spectrum using AWMI algorithm (red line) and number of
1154 iteration steps = 1000. Positions of fitted peaks are denoted by markers</b></p>
1155 
1156 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
1157 
1158 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
1159 fitting function using AWMI algorithm.</span></p>
1160 
1161 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
1162 do</span></p>
1163 
1164 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x FitAwmi.C</span></p>
1165 
1166 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>&nbsp;</span></p>
1167 
1168 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>void FitAwmi() {</span></p>
1169 
1170 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t a;</span></p>
1171 
1172 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t
1173 i,nfound=0,bin;</span></p>
1174 
1175 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t nbins = 256;</span></p>
1176 
1177 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmin  = 0;</span></p>
1178 
1179 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmax  =
1180 nbins;</span></p>
1181 
1182 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t * source =
1183 new Double_t[nbins];</span></p>
1184 
1185 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t * dest =
1186 new Double_t[nbins];   </span></p>
1187 
1188 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TH1F *h = new
1189 TH1F(&quot;h&quot;,&quot;Fitting using AWMI algorithm&quot;,nbins,xmin,xmax);</span></p>
1190 
1191 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TH1F *d = new
1192 TH1F(&quot;d&quot;,&quot;&quot;,nbins,xmin,xmax);      </span></p>
1193 
1194 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TFile *f = new
1195 TFile(&quot;TSpectrum.root&quot;);</span></p>
1196 
1197 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   h=(TH1F*)
1198 f-&gt;Get(&quot;fit;1&quot;);   </span></p>
1199 
1200 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for (i = 0; i &lt;
1201 nbins; i++) source[i]=h-&gt;GetBinContent(i + 1);      </span></p>
1202 
1203 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TCanvas *Fit1 =
1204 gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Fit1&quot;);</span></p>
1205 
1206 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   if (!Fit1) Fit1 =
1207 new TCanvas(&quot;Fit1&quot;,&quot;Fit1&quot;,10,10,1000,700);</span></p>
1208 
1209 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1210 h-&gt;Draw(&quot;L&quot;);</span></p>
1211 
1212 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TSpectrum *s = new
1213 TSpectrum();</span></p>
1214 
1215 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   //searching for
1216 candidate peaks positions</span></p>
1217 
1218 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   nfound =
1219 s-&gt;SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);</span></p>
1220 
1221 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Bool_t *FixPos =
1222 new Bool_t[nfound];</span></p>
1223 
1224 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Bool_t *FixAmp =
1225 new Bool_t[nfound];      </span></p>
1226 
1227 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for(i = 0; i&lt;
1228 nfound ; i++){</span></p>
1229 
1230 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>      FixPos[i] =
1231 kFALSE;</span></p>
1232 
1233 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>      FixAmp[i] =
1234 kFALSE;    </span></p>
1235 
1236 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }</span></p>
1237 
1238 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   //filling in the
1239 initial estimates of the input parameters</span></p>
1240 
1241 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t *PosX =
1242 new Double_t[nfound];         </span></p>
1243 
1244 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t *PosY =
1245 new Double_t[nfound];</span></p>
1246 
1247 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   PosX =
1248 s-&gt;GetPositionX();</span></p>
1249 
1250 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for (i = 0; i &lt;
1251 nfound; i++) {</span></p>
1252 
1253 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>                                                a=PosX[i];</span></p>
1254 
1255 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        bin = 1 +
1256 Int_t(a + 0.5);</span></p>
1257 
1258 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        PosY[i] =
1259 h-&gt;GetBinContent(bin);</span></p>
1260 
1261 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }   </span></p>
1262 
1263 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TSpectrumFit
1264 *pfit=new TSpectrumFit(nfound);</span></p>
1265 
1266 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1267 pfit-&gt;SetFitParameters(xmin, xmax-1, 1000, 0.1, pfit-&gt;kFitOptimChiCounts,
1268 pfit-&gt;kFitAlphaHalving, pfit-&gt;kFitPower2,
1269 pfit-&gt;kFitTaylorOrderFirst);   </span></p>
1270 
1271 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   pfit-&gt;SetPeakParameters(2,
1272 kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);   </span></p>
1273 
1274 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1275 pfit-&gt;FitAwmi(source);</span></p>
1276 
1277 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t
1278 *CalcPositions = new Double_t[nfound];      </span></p>
1279 
1280 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t
1281 *CalcAmplitudes = new Double_t[nfound];         </span></p>
1282 
1283 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1284 CalcPositions=pfit-&gt;GetPositions();</span></p>
1285 
1286 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1287 CalcAmplitudes=pfit-&gt;GetAmplitudes();   </span></p>
1288 
1289 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for (i = 0; i &lt;
1290 nbins; i++) d-&gt;SetBinContent(i + 1,source[i]);</span></p>
1291 
1292 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1293 d-&gt;SetLineColor(kRed);   </span></p>
1294 
1295 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1296 d-&gt;Draw(&quot;SAME L&quot;);  </span></p>
1297 
1298 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for (i = 0; i &lt;
1299 nfound; i++) {</span></p>
1300 
1301 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>                                                a=CalcPositions[i];</span></p>
1302 
1303 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        bin = 1 +
1304 Int_t(a + 0.5);                </span></p>
1305 
1306 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        PosX[i] =
1307 d-&gt;GetBinCenter(bin);</span></p>
1308 
1309 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        PosY[i] =
1310 d-&gt;GetBinContent(bin);</span></p>
1311 
1312 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }</span></p>
1313 
1314 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TPolyMarker * pm =
1315 (TPolyMarker*)h-&gt;GetListOfFunctions()-&gt;FindObject(&quot;TPolyMarker&quot;);</span></p>
1316 
1317 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   if (pm) {</span></p>
1318 
1319 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>     
1320 h-&gt;GetListOfFunctions()-&gt;Remove(pm);</span></p>
1321 
1322 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>      delete pm;</span></p>
1323 
1324 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }</span></p>
1325 
1326 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   pm = new
1327 TPolyMarker(nfound, PosX, PosY);</span></p>
1328 
1329 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1330 h-&gt;GetListOfFunctions()-&gt;Add(pm);</span></p>
1331 
1332 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1333 pm-&gt;SetMarkerStyle(23);</span></p>
1334 
1335 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1336 pm-&gt;SetMarkerColor(kRed);</span></p>
1337 
1338 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
1339 pm-&gt;SetMarkerSize(1);   </span></p>
1340 
1341 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>}</span></p>
1342 
1343 </div>
1344 
1345 <!-- */
1346 // --> End_Html
1347 
1348  Int_t i, j, k, shift =
1349  2 * fNPeaks + 7, peak_vel, rozmer, iter, pw, regul_cycle,
1350  flag;
1351  Double_t a, b, c, d = 0, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
1352  0, pi, pmin = 0, chi_cel = 0, chi_er;
1353  Double_t *working_space = new Double_t[5 * (2 * fNPeaks + 7)];
1354  for (i = 0, j = 0; i < fNPeaks; i++) {
1355  working_space[2 * i] = fAmpInit[i]; //vector parameter
1356  if (fFixAmp[i] == false) {
1357  working_space[shift + j] = fAmpInit[i]; //vector xk
1358  j += 1;
1359  }
1360  working_space[2 * i + 1] = fPositionInit[i]; //vector parameter
1361  if (fFixPosition[i] == false) {
1362  working_space[shift + j] = fPositionInit[i]; //vector xk
1363  j += 1;
1364  }
1365  }
1366  peak_vel = 2 * i;
1367  working_space[2 * i] = fSigmaInit; //vector parameter
1368  if (fFixSigma == false) {
1369  working_space[shift + j] = fSigmaInit; //vector xk
1370  j += 1;
1371  }
1372  working_space[2 * i + 1] = fTInit; //vector parameter
1373  if (fFixT == false) {
1374  working_space[shift + j] = fTInit; //vector xk
1375  j += 1;
1376  }
1377  working_space[2 * i + 2] = fBInit; //vector parameter
1378  if (fFixB == false) {
1379  working_space[shift + j] = fBInit; //vector xk
1380  j += 1;
1381  }
1382  working_space[2 * i + 3] = fSInit; //vector parameter
1383  if (fFixS == false) {
1384  working_space[shift + j] = fSInit; //vector xk
1385  j += 1;
1386  }
1387  working_space[2 * i + 4] = fA0Init; //vector parameter
1388  if (fFixA0 == false) {
1389  working_space[shift + j] = fA0Init; //vector xk
1390  j += 1;
1391  }
1392  working_space[2 * i + 5] = fA1Init; //vector parameter
1393  if (fFixA1 == false) {
1394  working_space[shift + j] = fA1Init; //vector xk
1395  j += 1;
1396  }
1397  working_space[2 * i + 6] = fA2Init; //vector parameter
1398  if (fFixA2 == false) {
1399  working_space[shift + j] = fA2Init; //vector xk
1400  j += 1;
1401  }
1402  rozmer = j;
1403  if (rozmer == 0){
1404  delete [] working_space;
1405  Error ("FitAwmi","All parameters are fixed");
1406  return;
1407  }
1408  if (rozmer >= fXmax - fXmin + 1){
1409  delete [] working_space;
1410  Error ("FitAwmi","Number of fitted parameters is larger than # of fitted points");
1411  return;
1412  }
1413  for (iter = 0; iter < fNumberIterations; iter++) {
1414  for (j = 0; j < rozmer; j++) {
1415  working_space[2 * shift + j] = 0, working_space[3 * shift + j] = 0; //der,temp
1416  }
1417 
1418  //filling vectors
1419  alpha = fAlpha;
1420  chi_opt = 0, pw = fPower - 2;
1421  for (i = fXmin; i <= fXmax; i++) {
1422  yw = source[i];
1423  ywm = yw;
1424  f = Shape(fNPeaks, (Double_t) i, working_space,
1425  working_space[peak_vel], working_space[peak_vel + 1],
1426  working_space[peak_vel + 3],
1427  working_space[peak_vel + 2],
1428  working_space[peak_vel + 4],
1429  working_space[peak_vel + 5],
1430  working_space[peak_vel + 6]);
1431  if (fStatisticType == kFitOptimMaxLikelihood) {
1432  if (f > 0.00001)
1433  chi_opt += yw * TMath::Log(f) - f;
1434  }
1435 
1436  else {
1437  if (ywm != 0)
1438  chi_opt += (yw - f) * (yw - f) / ywm;
1439  }
1440  if (fStatisticType == kFitOptimChiFuncValues) {
1441  ywm = f;
1442  if (f < 0.00001)
1443  ywm = 0.00001;
1444  }
1445 
1446  else if (fStatisticType == kFitOptimMaxLikelihood) {
1447  ywm = f;
1448  if (f < 0.001)
1449  ywm = 0.001;
1450  }
1451 
1452  else {
1453  if (ywm == 0)
1454  ywm = 1;
1455  }
1456 
1457  //calculation of gradient vector
1458  for (j = 0, k = 0; j < fNPeaks; j++) {
1459  if (fFixAmp[j] == false) {
1460  a = Deramp((Double_t) i, working_space[2 * j + 1],
1461  working_space[peak_vel],
1462  working_space[peak_vel + 1],
1463  working_space[peak_vel + 3],
1464  working_space[peak_vel + 2]);
1465  if (ywm != 0) {
1466  c = Ourpowl(a, pw);
1467  if (fStatisticType == kFitOptimChiFuncValues) {
1468  b = a * (yw * yw - f * f) / (ywm * ywm);
1469  working_space[2 * shift + k] += b * c; //der
1470  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1471  working_space[3 * shift + k] += b * c; //temp
1472  }
1473 
1474  else {
1475  b = a * (yw - f) / ywm;
1476  working_space[2 * shift + k] += b * c; //der
1477  b = a * a / ywm;
1478  working_space[3 * shift + k] += b * c; //temp
1479  }
1480  }
1481  k += 1;
1482  }
1483  if (fFixPosition[j] == false) {
1484  a = Deri0((Double_t) i, working_space[2 * j],
1485  working_space[2 * j + 1],
1486  working_space[peak_vel],
1487  working_space[peak_vel + 1],
1488  working_space[peak_vel + 3],
1489  working_space[peak_vel + 2]);
1490  if (fFitTaylor == kFitTaylorOrderSecond)
1491  d = Derderi0((Double_t) i, working_space[2 * j],
1492  working_space[2 * j + 1],
1493  working_space[peak_vel]);
1494  if (ywm != 0) {
1495  c = Ourpowl(a, pw);
1496  if (TMath::Abs(a) > 0.00000001
1497  && fFitTaylor == kFitTaylorOrderSecond) {
1498  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
1499  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
1500  d = 0;
1501  }
1502 
1503  else
1504  d = 0;
1505  a = a + d;
1506  if (fStatisticType == kFitOptimChiFuncValues) {
1507  b = a * (yw * yw - f * f) / (ywm * ywm);
1508  working_space[2 * shift + k] += b * c; //der
1509  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1510  working_space[3 * shift + k] += b * c; //temp
1511  }
1512 
1513  else {
1514  b = a * (yw - f) / ywm;
1515  working_space[2 * shift + k] += b * c; //der
1516  b = a * a / ywm;
1517  working_space[3 * shift + k] += b * c; //temp
1518  }
1519  }
1520  k += 1;
1521  }
1522  }
1523  if (fFixSigma == false) {
1524  a = Dersigma(fNPeaks, (Double_t) i, working_space,
1525  working_space[peak_vel],
1526  working_space[peak_vel + 1],
1527  working_space[peak_vel + 3],
1528  working_space[peak_vel + 2]);
1529  if (fFitTaylor == kFitTaylorOrderSecond)
1530  d = Derdersigma(fNPeaks, (Double_t) i,
1531  working_space, working_space[peak_vel]);
1532  if (ywm != 0) {
1533  c = Ourpowl(a, pw);
1534  if (TMath::Abs(a) > 0.00000001
1535  && fFitTaylor == kFitTaylorOrderSecond) {
1536  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
1537  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
1538  d = 0;
1539  }
1540 
1541  else
1542  d = 0;
1543  a = a + d;
1544  if (fStatisticType == kFitOptimChiFuncValues) {
1545  b = a * (yw * yw - f * f) / (ywm * ywm);
1546  working_space[2 * shift + k] += b * c; //der
1547  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1548  working_space[3 * shift + k] += b * c; //temp
1549  }
1550 
1551  else {
1552  b = a * (yw - f) / ywm;
1553  working_space[2 * shift + k] += b * c; //der
1554  b = a * a / ywm;
1555  working_space[3 * shift + k] += b * c; //temp
1556  }
1557  }
1558  k += 1;
1559  }
1560  if (fFixT == false) {
1561  a = Dert(fNPeaks, (Double_t) i, working_space,
1562  working_space[peak_vel],
1563  working_space[peak_vel + 2]);
1564  if (ywm != 0) {
1565  c = Ourpowl(a, pw);
1566  if (fStatisticType == kFitOptimChiFuncValues) {
1567  b = a * (yw * yw - f * f) / (ywm * ywm);
1568  working_space[2 * shift + k] += b * c; //der
1569  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1570  working_space[3 * shift + k] += b * c; //temp
1571  }
1572 
1573  else {
1574  b = a * (yw - f) / ywm;
1575  working_space[2 * shift + k] += b * c; //der
1576  b = a * a / ywm;
1577  working_space[3 * shift + k] += b * c; //temp
1578  }
1579  }
1580  k += 1;
1581  }
1582  if (fFixB == false) {
1583  a = Derb(fNPeaks, (Double_t) i, working_space,
1584  working_space[peak_vel], working_space[peak_vel + 1],
1585  working_space[peak_vel + 2]);
1586  if (ywm != 0) {
1587  c = Ourpowl(a, pw);
1588  if (fStatisticType == kFitOptimChiFuncValues) {
1589  b = a * (yw * yw - f * f) / (ywm * ywm);
1590  working_space[2 * shift + k] += b * c; //der
1591  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1592  working_space[3 * shift + k] += b * c; //temp
1593  }
1594 
1595  else {
1596  b = a * (yw - f) / ywm;
1597  working_space[2 * shift + k] += b * c; //der
1598  b = a * a / ywm;
1599  working_space[3 * shift + k] += b * c; //temp
1600  }
1601  }
1602  k += 1;
1603  }
1604  if (fFixS == false) {
1605  a = Ders(fNPeaks, (Double_t) i, working_space,
1606  working_space[peak_vel]);
1607  if (ywm != 0) {
1608  c = Ourpowl(a, pw);
1609  if (fStatisticType == kFitOptimChiFuncValues) {
1610  b = a * (yw * yw - f * f) / (ywm * ywm);
1611  working_space[2 * shift + k] += b * c; //der
1612  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1613  working_space[3 * shift + k] += b * c; //temp
1614  }
1615 
1616  else {
1617  b = a * (yw - f) / ywm;
1618  working_space[2 * shift + k] += b * c; //der
1619  b = a * a / ywm;
1620  working_space[3 * shift + k] += b * c; //temp
1621  }
1622  }
1623  k += 1;
1624  }
1625  if (fFixA0 == false) {
1626  a = 1.;
1627  if (ywm != 0) {
1628  c = Ourpowl(a, pw);
1629  if (fStatisticType == kFitOptimChiFuncValues) {
1630  b = a * (yw * yw - f * f) / (ywm * ywm);
1631  working_space[2 * shift + k] += b * c; //der
1632  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1633  working_space[3 * shift + k] += b * c; //temp
1634  }
1635 
1636  else {
1637  b = a * (yw - f) / ywm;
1638  working_space[2 * shift + k] += b * c; //der
1639  b = a * a / ywm;
1640  working_space[3 * shift + k] += b * c; //temp
1641  }
1642  }
1643  k += 1;
1644  }
1645  if (fFixA1 == false) {
1646  a = Dera1((Double_t) i);
1647  if (ywm != 0) {
1648  c = Ourpowl(a, pw);
1649  if (fStatisticType == kFitOptimChiFuncValues) {
1650  b = a * (yw * yw - f * f) / (ywm * ywm);
1651  working_space[2 * shift + k] += b * c; //der
1652  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1653  working_space[3 * shift + k] += b * c; //temp
1654  }
1655 
1656  else {
1657  b = a * (yw - f) / ywm;
1658  working_space[2 * shift + k] += b * c; //der
1659  b = a * a / ywm;
1660  working_space[3 * shift + k] += b * c; //temp
1661  }
1662  }
1663  k += 1;
1664  }
1665  if (fFixA2 == false) {
1666  a = Dera2((Double_t) i);
1667  if (ywm != 0) {
1668  c = Ourpowl(a, pw);
1669  if (fStatisticType == kFitOptimChiFuncValues) {
1670  b = a * (yw * yw - f * f) / (ywm * ywm);
1671  working_space[2 * shift + k] += b * c; //der
1672  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
1673  working_space[3 * shift + k] += b * c; //temp
1674  }
1675 
1676  else {
1677  b = a * (yw - f) / ywm;
1678  working_space[2 * shift + k] += b * c; //der
1679  b = a * a / ywm;
1680  working_space[3 * shift + k] += b * c; //temp
1681  }
1682  }
1683  k += 1;
1684  }
1685  }
1686  for (j = 0; j < rozmer; j++) {
1687  if (TMath::Abs(working_space[3 * shift + j]) > 0.000001)
1688  working_space[2 * shift + j] = working_space[2 * shift + j] / TMath::Abs(working_space[3 * shift + j]); //der[j]=der[j]/temp[j]
1689  else
1690  working_space[2 * shift + j] = 0; //der[j]
1691  }
1692 
1693  //calculate chi_opt
1694  chi2 = chi_opt;
1695  chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
1696 
1697  //calculate new parameters
1698  regul_cycle = 0;
1699  for (j = 0; j < rozmer; j++) {
1700  working_space[4 * shift + j] = working_space[shift + j]; //temp_xk[j]=xk[j]
1701  }
1702 
1703  do {
1704  if (fAlphaOptim == kFitAlphaOptimal) {
1705  if (fStatisticType != kFitOptimMaxLikelihood)
1706  chi_min = 10000 * chi2;
1707 
1708  else
1709  chi_min = 0.1 * chi2;
1710  flag = 0;
1711  for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
1712  for (j = 0; j < rozmer; j++) {
1713  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]
1714  }
1715  for (i = 0, j = 0; i < fNPeaks; i++) {
1716  if (fFixAmp[i] == false) {
1717  if (working_space[shift + j] < 0) //xk[j]
1718  working_space[shift + j] = 0; //xk[j]
1719  working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
1720  j += 1;
1721  }
1722  if (fFixPosition[i] == false) {
1723  if (working_space[shift + j] < fXmin) //xk[j]
1724  working_space[shift + j] = fXmin; //xk[j]
1725  if (working_space[shift + j] > fXmax) //xk[j]
1726  working_space[shift + j] = fXmax; //xk[j]
1727  working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
1728  j += 1;
1729  }
1730  }
1731  if (fFixSigma == false) {
1732  if (working_space[shift + j] < 0.001) { //xk[j]
1733  working_space[shift + j] = 0.001; //xk[j]
1734  }
1735  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
1736  j += 1;
1737  }
1738  if (fFixT == false) {
1739  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
1740  j += 1;
1741  }
1742  if (fFixB == false) {
1743  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
1744  if (working_space[shift + j] < 0) //xk[j]
1745  working_space[shift + j] = -0.001; //xk[j]
1746  else
1747  working_space[shift + j] = 0.001; //xk[j]
1748  }
1749  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
1750  j += 1;
1751  }
1752  if (fFixS == false) {
1753  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
1754  j += 1;
1755  }
1756  if (fFixA0 == false) {
1757  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
1758  j += 1;
1759  }
1760  if (fFixA1 == false) {
1761  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
1762  j += 1;
1763  }
1764  if (fFixA2 == false) {
1765  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
1766  j += 1;
1767  }
1768  chi2 = 0;
1769  for (i = fXmin; i <= fXmax; i++) {
1770  yw = source[i];
1771  ywm = yw;
1772  f = Shape(fNPeaks, (Double_t) i, working_space,
1773  working_space[peak_vel],
1774  working_space[peak_vel + 1],
1775  working_space[peak_vel + 3],
1776  working_space[peak_vel + 2],
1777  working_space[peak_vel + 4],
1778  working_space[peak_vel + 5],
1779  working_space[peak_vel + 6]);
1780  if (fStatisticType == kFitOptimChiFuncValues) {
1781  ywm = f;
1782  if (f < 0.00001)
1783  ywm = 0.00001;
1784  }
1785  if (fStatisticType == kFitOptimMaxLikelihood) {
1786  if (f > 0.00001)
1787  chi2 += yw * TMath::Log(f) - f;
1788  }
1789 
1790  else {
1791  if (ywm != 0)
1792  chi2 += (yw - f) * (yw - f) / ywm;
1793  }
1794  }
1795  if ((chi2 < chi_min
1796  && fStatisticType != kFitOptimMaxLikelihood)
1797  || (chi2 > chi_min
1798  && fStatisticType == kFitOptimMaxLikelihood)) {
1799  pmin = pi, chi_min = chi2;
1800  }
1801 
1802  else
1803  flag = 1;
1804  if (pi == 0.1)
1805  chi_min = chi2;
1806  chi = chi_min;
1807  }
1808  if (pmin != 0.1) {
1809  for (j = 0; j < rozmer; j++) {
1810  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]
1811  }
1812  for (i = 0, j = 0; i < fNPeaks; i++) {
1813  if (fFixAmp[i] == false) {
1814  if (working_space[shift + j] < 0) //xk[j]
1815  working_space[shift + j] = 0; //xk[j]
1816  working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
1817  j += 1;
1818  }
1819  if (fFixPosition[i] == false) {
1820  if (working_space[shift + j] < fXmin) //xk[j]
1821  working_space[shift + j] = fXmin; //xk[j]
1822  if (working_space[shift + j] > fXmax) //xk[j]
1823  working_space[shift + j] = fXmax; //xk[j]
1824  working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
1825  j += 1;
1826  }
1827  }
1828  if (fFixSigma == false) {
1829  if (working_space[shift + j] < 0.001) { //xk[j]
1830  working_space[shift + j] = 0.001; //xk[j]
1831  }
1832  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
1833  j += 1;
1834  }
1835  if (fFixT == false) {
1836  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
1837  j += 1;
1838  }
1839  if (fFixB == false) {
1840  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
1841  if (working_space[shift + j] < 0) //xk[j]
1842  working_space[shift + j] = -0.001; //xk[j]
1843  else
1844  working_space[shift + j] = 0.001; //xk[j]
1845  }
1846  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
1847  j += 1;
1848  }
1849  if (fFixS == false) {
1850  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
1851  j += 1;
1852  }
1853  if (fFixA0 == false) {
1854  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
1855  j += 1;
1856  }
1857  if (fFixA1 == false) {
1858  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
1859  j += 1;
1860  }
1861  if (fFixA2 == false) {
1862  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
1863  j += 1;
1864  }
1865  chi = chi_min;
1866  }
1867  }
1868 
1869  else {
1870  for (j = 0; j < rozmer; j++) {
1871  working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
1872  }
1873  for (i = 0, j = 0; i < fNPeaks; i++) {
1874  if (fFixAmp[i] == false) {
1875  if (working_space[shift + j] < 0) //xk[j]
1876  working_space[shift + j] = 0; //xk[j]
1877  working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
1878  j += 1;
1879  }
1880  if (fFixPosition[i] == false) {
1881  if (working_space[shift + j] < fXmin) //xk[j]
1882  working_space[shift + j] = fXmin; //xk[j]
1883  if (working_space[shift + j] > fXmax) //xk[j]
1884  working_space[shift + j] = fXmax; //xk[j]
1885  working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
1886  j += 1;
1887  }
1888  }
1889  if (fFixSigma == false) {
1890  if (working_space[shift + j] < 0.001) { //xk[j]
1891  working_space[shift + j] = 0.001; //xk[j]
1892  }
1893  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
1894  j += 1;
1895  }
1896  if (fFixT == false) {
1897  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
1898  j += 1;
1899  }
1900  if (fFixB == false) {
1901  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
1902  if (working_space[shift + j] < 0) //xk[j]
1903  working_space[shift + j] = -0.001; //xk[j]
1904  else
1905  working_space[shift + j] = 0.001; //xk[j]
1906  }
1907  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
1908  j += 1;
1909  }
1910  if (fFixS == false) {
1911  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
1912  j += 1;
1913  }
1914  if (fFixA0 == false) {
1915  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
1916  j += 1;
1917  }
1918  if (fFixA1 == false) {
1919  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
1920  j += 1;
1921  }
1922  if (fFixA2 == false) {
1923  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
1924  j += 1;
1925  }
1926  chi = 0;
1927  for (i = fXmin; i <= fXmax; i++) {
1928  yw = source[i];
1929  ywm = yw;
1930  f = Shape(fNPeaks, (Double_t) i, working_space,
1931  working_space[peak_vel],
1932  working_space[peak_vel + 1],
1933  working_space[peak_vel + 3],
1934  working_space[peak_vel + 2],
1935  working_space[peak_vel + 4],
1936  working_space[peak_vel + 5],
1937  working_space[peak_vel + 6]);
1938  if (fStatisticType == kFitOptimChiFuncValues) {
1939  ywm = f;
1940  if (f < 0.00001)
1941  ywm = 0.00001;
1942  }
1943  if (fStatisticType == kFitOptimMaxLikelihood) {
1944  if (f > 0.00001)
1945  chi += yw * TMath::Log(f) - f;
1946  }
1947 
1948  else {
1949  if (ywm != 0)
1950  chi += (yw - f) * (yw - f) / ywm;
1951  }
1952  }
1953  }
1954  chi2 = chi;
1955  chi = TMath::Sqrt(TMath::Abs(chi));
1956  if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
1957  alpha = alpha * chi_opt / (2 * chi);
1958 
1959  else if (fAlphaOptim == kFitAlphaOptimal)
1960  alpha = alpha / 10.0;
1961  iter += 1;
1962  regul_cycle += 1;
1963  } while (((chi > chi_opt
1964  && fStatisticType != kFitOptimMaxLikelihood)
1965  || (chi < chi_opt
1966  && fStatisticType == kFitOptimMaxLikelihood))
1967  && regul_cycle < kFitNumRegulCycles);
1968  for (j = 0; j < rozmer; j++) {
1969  working_space[4 * shift + j] = 0; //temp_xk[j]
1970  working_space[2 * shift + j] = 0; //der[j]
1971  }
1972  for (i = fXmin, chi_cel = 0; i <= fXmax; i++) {
1973  yw = source[i];
1974  if (yw == 0)
1975  yw = 1;
1976  f = Shape(fNPeaks, (Double_t) i, working_space,
1977  working_space[peak_vel], working_space[peak_vel + 1],
1978  working_space[peak_vel + 3],
1979  working_space[peak_vel + 2],
1980  working_space[peak_vel + 4],
1981  working_space[peak_vel + 5],
1982  working_space[peak_vel + 6]);
1983  chi_opt = (yw - f) * (yw - f) / yw;
1984  chi_cel += (yw - f) * (yw - f) / yw;
1985 
1986  //calculate gradient vector
1987  for (j = 0, k = 0; j < fNPeaks; j++) {
1988  if (fFixAmp[j] == false) {
1989  a = Deramp((Double_t) i, working_space[2 * j + 1],
1990  working_space[peak_vel],
1991  working_space[peak_vel + 1],
1992  working_space[peak_vel + 3],
1993  working_space[peak_vel + 2]);
1994  if (yw != 0) {
1995  c = Ourpowl(a, pw);
1996  working_space[2 * shift + k] += chi_opt * c; //der[k]
1997  b = a * a / yw;
1998  working_space[4 * shift + k] += b * c; //temp_xk[k]
1999  }
2000  k += 1;
2001  }
2002  if (fFixPosition[j] == false) {
2003  a = Deri0((Double_t) i, working_space[2 * j],
2004  working_space[2 * j + 1],
2005  working_space[peak_vel],
2006  working_space[peak_vel + 1],
2007  working_space[peak_vel + 3],
2008  working_space[peak_vel + 2]);
2009  if (yw != 0) {
2010  c = Ourpowl(a, pw);
2011  working_space[2 * shift + k] += chi_opt * c; //der[k]
2012  b = a * a / yw;
2013  working_space[4 * shift + k] += b * c; //temp_xk[k]
2014  }
2015  k += 1;
2016  }
2017  }
2018  if (fFixSigma == false) {
2019  a = Dersigma(fNPeaks, (Double_t) i, working_space,
2020  working_space[peak_vel],
2021  working_space[peak_vel + 1],
2022  working_space[peak_vel + 3],
2023  working_space[peak_vel + 2]);
2024  if (yw != 0) {
2025  c = Ourpowl(a, pw);
2026  working_space[2 * shift + k] += chi_opt * c; //der[k]
2027  b = a * a / yw;
2028  working_space[4 * shift + k] += b * c; //temp_xk[k]
2029  }
2030  k += 1;
2031  }
2032  if (fFixT == false) {
2033  a = Dert(fNPeaks, (Double_t) i, working_space,
2034  working_space[peak_vel],
2035  working_space[peak_vel + 2]);
2036  if (yw != 0) {
2037  c = Ourpowl(a, pw);
2038  working_space[2 * shift + k] += chi_opt * c; //der[k]
2039  b = a * a / yw;
2040  working_space[4 * shift + k] += b * c; //temp_xk[k]
2041  }
2042  k += 1;
2043  }
2044  if (fFixB == false) {
2045  a = Derb(fNPeaks, (Double_t) i, working_space,
2046  working_space[peak_vel], working_space[peak_vel + 1],
2047  working_space[peak_vel + 2]);
2048  if (yw != 0) {
2049  c = Ourpowl(a, pw);
2050  working_space[2 * shift + k] += chi_opt * c; //der[k]
2051  b = a * a / yw;
2052  working_space[4 * shift + k] += b * c; //temp_xk[k]
2053  }
2054  k += 1;
2055  }
2056  if (fFixS == false) {
2057  a = Ders(fNPeaks, (Double_t) i, working_space,
2058  working_space[peak_vel]);
2059  if (yw != 0) {
2060  c = Ourpowl(a, pw);
2061  working_space[2 * shift + k] += chi_opt * c; //der[k]
2062  b = a * a / yw;
2063  working_space[4 * shift + k] += b * c; //tem_xk[k]
2064  }
2065  k += 1;
2066  }
2067  if (fFixA0 == false) {
2068  a = 1.0;
2069  if (yw != 0) {
2070  c = Ourpowl(a, pw);
2071  working_space[2 * shift + k] += chi_opt * c; //der[k]
2072  b = a * a / yw;
2073  working_space[4 * shift + k] += b * c; //temp_xk[k]
2074  }
2075  k += 1;
2076  }
2077  if (fFixA1 == false) {
2078  a = Dera1((Double_t) i);
2079  if (yw != 0) {
2080  c = Ourpowl(a, pw);
2081  working_space[2 * shift + k] += chi_opt * c; //der[k]
2082  b = a * a / yw;
2083  working_space[4 * shift + k] += b * c; //temp_xk[k]
2084  }
2085  k += 1;
2086  }
2087  if (fFixA2 == false) {
2088  a = Dera2((Double_t) i);
2089  if (yw != 0) {
2090  c = Ourpowl(a, pw);
2091  working_space[2 * shift + k] += chi_opt * c; //der[k]
2092  b = a * a / yw;
2093  working_space[4 * shift + k] += b * c; //temp_xk[k]
2094  }
2095  k += 1;
2096  }
2097  }
2098  }
2099  b = fXmax - fXmin + 1 - rozmer;
2100  chi_er = chi_cel / b;
2101  for (i = 0, j = 0; i < fNPeaks; i++) {
2102  fArea[i] =
2103  Area(working_space[2 * i], working_space[peak_vel],
2104  working_space[peak_vel + 1], working_space[peak_vel + 2]);
2105  if (fFixAmp[i] == false) {
2106  fAmpCalc[i] = working_space[shift + j]; //xk[j]
2107  if (working_space[3 * shift + j] != 0)
2108  fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2109  if (fArea[i] > 0) {
2110  a = Derpa(working_space[peak_vel],
2111  working_space[peak_vel + 1],
2112  working_space[peak_vel + 2]);
2113  b = working_space[4 * shift + j]; //temp_xk[j]
2114  if (b == 0)
2115  b = 1;
2116 
2117  else
2118  b = 1 / b;
2119  fAreaErr[i] = TMath::Sqrt(TMath::Abs(a * a * b * chi_er));
2120  }
2121 
2122  else
2123  fAreaErr[i] = 0;
2124  j += 1;
2125  }
2126 
2127  else {
2128  fAmpCalc[i] = fAmpInit[i];
2129  fAmpErr[i] = 0;
2130  fAreaErr[i] = 0;
2131  }
2132  if (fFixPosition[i] == false) {
2133  fPositionCalc[i] = working_space[shift + j]; //xk[j]
2134  if (working_space[3 * shift + j] != 0) //temp[j]
2135  fPositionErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2136  j += 1;
2137  }
2138 
2139  else {
2140  fPositionCalc[i] = fPositionInit[i];
2141  fPositionErr[i] = 0;
2142  }
2143  }
2144  if (fFixSigma == false) {
2145  fSigmaCalc = working_space[shift + j]; //xk[j]
2146  if (working_space[3 * shift + j] != 0) //temp[j]
2147  fSigmaErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2148  j += 1;
2149  }
2150 
2151  else {
2152  fSigmaCalc = fSigmaInit;
2153  fSigmaErr = 0;
2154  }
2155  if (fFixT == false) {
2156  fTCalc = working_space[shift + j]; //xk[j]
2157  if (working_space[3 * shift + j] != 0) //temp[j]
2158  fTErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2159  j += 1;
2160  }
2161 
2162  else {
2163  fTCalc = fTInit;
2164  fTErr = 0;
2165  }
2166  if (fFixB == false) {
2167  fBCalc = working_space[shift + j]; //xk[j]
2168  if (working_space[3 * shift + j] != 0) //temp[j]
2169  fBErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2170  j += 1;
2171  }
2172 
2173  else {
2174  fBCalc = fBInit;
2175  fBErr = 0;
2176  }
2177  if (fFixS == false) {
2178  fSCalc = working_space[shift + j]; //xk[j]
2179  if (working_space[3 * shift + j] != 0) //temp[j]
2180  fSErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2181  j += 1;
2182  }
2183 
2184  else {
2185  fSCalc = fSInit;
2186  fSErr = 0;
2187  }
2188  if (fFixA0 == false) {
2189  fA0Calc = working_space[shift + j]; //xk[j]
2190  if (working_space[3 * shift + j] != 0) //temp[j]
2191  fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2192  j += 1;
2193  }
2194 
2195  else {
2196  fA0Calc = fA0Init;
2197  fA0Err = 0;
2198  }
2199  if (fFixA1 == false) {
2200  fA1Calc = working_space[shift + j]; //xk[j]
2201  if (working_space[3 * shift + j] != 0) //temp[j]
2202  fA1Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2203  j += 1;
2204  }
2205 
2206  else {
2207  fA1Calc = fA1Init;
2208  fA1Err = 0;
2209  }
2210  if (fFixA2 == false) {
2211  fA2Calc = working_space[shift + j]; //xk[j]
2212  if (working_space[3 * shift + j] != 0) //temp[j]
2213  fA2Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
2214  j += 1;
2215  }
2216 
2217  else {
2218  fA2Calc = fA2Init;
2219  fA2Err = 0;
2220  }
2221  b = fXmax - fXmin + 1 - rozmer;
2222  fChi = chi_cel / b;
2223  for (i = fXmin; i <= fXmax; i++) {
2224  f = Shape(fNPeaks, (Double_t) i, working_space,
2225  working_space[peak_vel], working_space[peak_vel + 1],
2226  working_space[peak_vel + 3], working_space[peak_vel + 2],
2227  working_space[peak_vel + 4], working_space[peak_vel + 5],
2228  working_space[peak_vel + 6]);
2229  source[i] = f;
2230  }
2231  delete[]working_space;
2232  return;
2233 }
2234 
2235 /////////////////FITTING FUNCTION WITH MATRIX INVERSION///////////////////////////////////////
2236 
2237 ////////////////////////////////////////////////////////////////////////////////
2238 ///////////////////////////////////////////////////////////////////////////////////
2239 /// AUXILIARY FUNCTION //
2240 /// //
2241 /// This function calculates solution of the system of linear equations. //
2242 /// The matrix a should have a dimension size*(size+4) //
2243 /// The calling function should fill in the matrix, the column size should //
2244 /// contain vector y (right side of the system of equations). The result is //
2245 /// placed into size+1 column of the matrix. //
2246 /// according to sigma of peaks. //
2247 /// Function parameters: //
2248 /// -a-matrix with dimension size*(size+4) // //
2249 /// -size-number of rows of the matrix //
2250 /// //
2251 ///////////////////////////////////////////////////////////////////////////////////
2252 
2253 void TSpectrumFit::StiefelInversion(Double_t **a, Int_t size)
2254 {
2255  Int_t i, j, k = 0;
2256  Double_t sk = 0, b, lambdak, normk, normk_old = 0;
2257 
2258  do {
2259  normk = 0;
2260 
2261  //calculation of rk and norm
2262  for (i = 0; i < size; i++) {
2263  a[i][size + 2] = -a[i][size]; //rk=-C
2264  for (j = 0; j < size; j++) {
2265  a[i][size + 2] += a[i][j] * a[j][size + 1]; //A*xk-C
2266  }
2267  normk += a[i][size + 2] * a[i][size + 2]; //calculation normk
2268  }
2269 
2270  //calculation of sk
2271  if (k != 0) {
2272  sk = normk / normk_old;
2273  }
2274 
2275  //calculation of uk
2276  for (i = 0; i < size; i++) {
2277  a[i][size + 3] = -a[i][size + 2] + sk * a[i][size + 3]; //uk=-rk+sk*uk-1
2278  }
2279 
2280  //calculation of lambdak
2281  lambdak = 0;
2282  for (i = 0; i < size; i++) {
2283  for (j = 0, b = 0; j < size; j++) {
2284  b += a[i][j] * a[j][size + 3]; //A*uk
2285  }
2286  lambdak += b * a[i][size + 3];
2287  }
2288  if (TMath::Abs(lambdak) > 1e-50) //computer zero
2289  lambdak = normk / lambdak;
2290 
2291  else
2292  lambdak = 0;
2293  for (i = 0; i < size; i++)
2294  a[i][size + 1] += lambdak * a[i][size + 3]; //xk+1=xk+lambdak*uk
2295  normk_old = normk;
2296  k += 1;
2297  } while (k < size && TMath::Abs(normk) > 1e-50); //computer zero
2298  return;
2299 }
2300 
2301 ////////////////////////////////////////////////////////////////////////////////
2302 //////////////////////////////////////////////////////////////////////////////
2303 /// ONE-DIMENSIONAL FIT FUNCTION
2304 /// ALGORITHM WITH MATRIX INVERSION (STIEFEL-HESTENS METHOD)
2305 /// This function fits the source spectrum. The calling program should
2306 /// fill in input parameters
2307 /// The fitted parameters are written into
2308 /// output parameters and fitted data are written into
2309 /// source spectrum.
2310 ///
2311 /// Function parameters:
2312 /// source-pointer to the vector of source spectrum
2313 ///
2314 //////////////////////////////////////////////////////////////////////////////
2315 ///Begin_Html <!--
2316 
2317 void TSpectrumFit::FitStiefel(Double_t *source)
2318 {
2319 /* -->
2320 <div class=Section3>
2321 
2322 <p class=MsoNormal><b><span style='font-size:14.0pt'>Stiefel fitting algorithm</span></b></p>
2323 
2324 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>&nbsp;</span></i></p>
2325 
2326 <p class=MsoNormal><i>Function:</i></p>
2327 
2328 <p class=MsoNormal style='text-align:justify'>void <a
2329 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrumFit::</b></a>FitStiefel(<a
2330 href="http://root.cern.ch/root/html/ListOfTypes.html#double"><b>double</b></a> *fSource)
2331 </p>
2332 
2333 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
2334 
2335 <p class=MsoNormal style='text-align:justify'>This function fits the source
2336 spectrum using Stiefel-Hestens method [1] (see Awmi function).  The calling
2337 program should fill in input fitting parameters of the TSpectrumFit class using
2338 a set of TSpectrumFit setters. The fitted parameters are written into the class
2339 and the fitted data are written into source spectrum. It converges faster than
2340 Awmi method.</p>
2341 
2342 <p class=MsoNormal>&nbsp;</p>
2343 
2344 <p class=MsoNormal><i><span style='color:red'>Parameter:</span></i></p>
2345 
2346 <p class=MsoNormal style='text-align:justify'>        <b>fSource</b>-pointer to
2347 the vector of source spectrum                  </p>
2348 
2349 <p class=MsoNormal style='text-align:justify'>        </p>
2350 
2351 <p class=MsoNormal style='text-align:justify'><i>Example – script FitStiefel.c:</i></p>
2352 
2353 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
2354 border=0 width=601 height=402 src="gif/spectrumfit_stiefel_image001.jpg"></span></i></p>
2355 
2356 <p class=MsoNormal style='text-align:justify'><b>Fig. 2 Original spectrum
2357 (black line) and fitted spectrum using Stiefel-Hestens method (red line) and
2358 number of iteration steps = 100. Positions of fitted peaks are denoted by
2359 markers</b></p>
2360 
2361 <p class=MsoNormal><b><span style='color:#339966'>&nbsp;</span></b></p>
2362 
2363 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
2364 
2365 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
2366 fitting function using Stiefel-Hestens method.</span></p>
2367 
2368 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example, do</span></p>
2369 
2370 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x FitStiefel.C</span></p>
2371 
2372 <p class=MsoNormal><span style='font-size:10.0pt'>&nbsp;</span></p>
2373 
2374 <p class=MsoNormal><span style='font-size:10.0pt'>void FitStiefel() {</span></p>
2375 
2376 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t a;</span></p>
2377 
2378 <p class=MsoNormal><span style='font-size:10.0pt'>   Int_t i,nfound=0,bin;</span></p>
2379 
2380 <p class=MsoNormal><span style='font-size:10.0pt'>   Int_t nbins = 256;</span></p>
2381 
2382 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
2383 style='font-size:10.0pt'>Int_t xmin  = 0;</span></p>
2384 
2385 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmax  =
2386 nbins;</span></p>
2387 
2388 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
2389 style='font-size:10.0pt'>Double_t * source = new Double_t[nbins];</span></p>
2390 
2391 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t * dest = new
2392 Double_t[nbins];   </span></p>
2393 
2394 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *h = new TH1F(&quot;h&quot;,&quot;Fitting
2395 using AWMI algorithm&quot;,nbins,xmin,xmax);</span></p>
2396 
2397 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *d = new
2398 TH1F(&quot;d&quot;,&quot;&quot;,nbins,xmin,xmax);      </span></p>
2399 
2400 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
2401 TFile(&quot;TSpectrum.root&quot;);</span></p>
2402 
2403 <p class=MsoNormal><span style='font-size:10.0pt'>   h=(TH1F*)
2404 f-&gt;Get(&quot;fit;1&quot;);   </span></p>
2405 
2406 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
2407 i++) source[i]=h-&gt;GetBinContent(i + 1);      </span></p>
2408 
2409 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Fit1 =
2410 gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Fit1&quot;);</span></p>
2411 
2412 <p class=MsoNormal><span style='font-size:10.0pt'>   if (!Fit1) Fit1 = new
2413 TCanvas(&quot;Fit1&quot;,&quot;Fit1&quot;,10,10,1000,700);</span></p>
2414 
2415 <p class=MsoNormal><span style='font-size:10.0pt'>   h-&gt;Draw(&quot;L&quot;);</span></p>
2416 
2417 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrum *s = new
2418 TSpectrum();</span></p>
2419 
2420 <p class=MsoNormal><span style='font-size:10.0pt'>   //searching for candidate
2421 peaks positions</span></p>
2422 
2423 <p class=MsoNormal><span style='font-size:10.0pt'>   nfound =
2424 s-&gt;SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);</span></p>
2425 
2426 <p class=MsoNormal><span style='font-size:10.0pt'>   Bool_t *FixPos = new
2427 Bool_t[nfound];</span></p>
2428 
2429 <p class=MsoNormal><span style='font-size:10.0pt'>   Bool_t *FixAmp = new
2430 Bool_t[nfound];      </span></p>
2431 
2432 <p class=MsoNormal><span style='font-size:10.0pt'>   for(i = 0; i&lt; nfound ;
2433 i++){</span></p>
2434 
2435 <p class=MsoNormal><span style='font-size:10.0pt'>      FixPos[i] = kFALSE;</span></p>
2436 
2437 <p class=MsoNormal><span style='font-size:10.0pt'>      FixAmp[i] = kFALSE;    </span></p>
2438 
2439 <p class=MsoNormal><span style='font-size:10.0pt'>   }</span></p>
2440 
2441 <p class=MsoNormal><span style='font-size:10.0pt'>   //filling in the initial
2442 estimates of the input parameters</span></p>
2443 
2444 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t *PosX = new
2445 Double_t[nfound];         </span></p>
2446 
2447 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t *PosY = new
2448 Double_t[nfound];</span></p>
2449 
2450 <p class=MsoNormal><span style='font-size:10.0pt'>   PosX =
2451 s-&gt;GetPositionX();</span></p>
2452 
2453 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nfound;
2454 i++) {</span></p>
2455 
2456 <p class=MsoNormal><span style='font-size:10.0pt'>                                                a=PosX[i];</span></p>
2457 
2458 <p class=MsoNormal><span style='font-size:10.0pt'>        bin = 1 + Int_t(a +
2459 0.5);</span></p>
2460 
2461 <p class=MsoNormal><span style='font-size:10.0pt'>        PosY[i] =
2462 h-&gt;GetBinContent(bin);</span></p>
2463 
2464 <p class=MsoNormal><span style='font-size:10.0pt'>   }   </span></p>
2465 
2466 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrumFit *pfit=new
2467 TSpectrumFit(nfound);</span></p>
2468 
2469 <p class=MsoNormal><span style='font-size:10.0pt'>  
2470 pfit-&gt;SetFitParameters(xmin, xmax-1, 1000, 0.1, pfit-&gt;kFitOptimChiCounts,
2471 pfit-&gt;kFitAlphaHalving, pfit-&gt;kFitPower2,
2472 pfit-&gt;kFitTaylorOrderFirst);   </span></p>
2473 
2474 <p class=MsoNormal><span style='font-size:10.0pt'>   pfit-&gt;SetPeakParameters(2,
2475 kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);   </span></p>
2476 
2477 <p class=MsoNormal><span style='font-size:10.0pt'>  
2478 pfit-&gt;FitStiefel(source);</span></p>
2479 
2480 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t *CalcPositions =
2481 new Double_t[nfound];      </span></p>
2482 
2483 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
2484 style='font-size:10.0pt'>Double_t *CalcAmplitudes = new
2485 Double_t[nfound];         </span></p>
2486 
2487 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
2488 style='font-size:10.0pt'>CalcPositions=pfit-&gt;GetPositions();</span></p>
2489 
2490 <p class=MsoNormal><span style='font-size:10.0pt'>  
2491 CalcAmplitudes=pfit-&gt;GetAmplitudes();   </span></p>
2492 
2493 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
2494 i++) d-&gt;SetBinContent(i + 1,source[i]);</span></p>
2495 
2496 <p class=MsoNormal><span style='font-size:10.0pt'>  
2497 d-&gt;SetLineColor(kRed);   </span></p>
2498 
2499 <p class=MsoNormal><span style='font-size:10.0pt'>   d-&gt;Draw(&quot;SAME
2500 L&quot;);  </span></p>
2501 
2502 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nfound;
2503 i++) {</span></p>
2504 
2505 <p class=MsoNormal><span style='font-size:10.0pt'>                                                a=CalcPositions[i];</span></p>
2506 
2507 <p class=MsoNormal><span style='font-size:10.0pt'>        bin = 1 + Int_t(a +
2508 0.5);                </span></p>
2509 
2510 <p class=MsoNormal><span style='font-size:10.0pt'>        PosX[i] =
2511 d-&gt;GetBinCenter(bin);</span></p>
2512 
2513 <p class=MsoNormal><span style='font-size:10.0pt'>        PosY[i] =
2514 d-&gt;GetBinContent(bin);</span></p>
2515 
2516 <p class=MsoNormal><span style='font-size:10.0pt'>   }</span></p>
2517 
2518 <p class=MsoNormal><span style='font-size:10.0pt'>   TPolyMarker * pm =
2519 (TPolyMarker*)h-&gt;GetListOfFunctions()-&gt;FindObject(&quot;TPolyMarker&quot;);</span></p>
2520 
2521 <p class=MsoNormal><span style='font-size:10.0pt'>   if (pm) {</span></p>
2522 
2523 <p class=MsoNormal><span style='font-size:10.0pt'>     
2524 h-&gt;GetListOfFunctions()-&gt;Remove(pm);</span></p>
2525 
2526 <p class=MsoNormal><span style='font-size:10.0pt'>      delete pm;</span></p>
2527 
2528 <p class=MsoNormal><span style='font-size:10.0pt'>   }</span></p>
2529 
2530 <p class=MsoNormal><span style='font-size:10.0pt'>   pm = new
2531 TPolyMarker(nfound, PosX, PosY);</span></p>
2532 
2533 <p class=MsoNormal><span style='font-size:10.0pt'>  
2534 h-&gt;GetListOfFunctions()-&gt;Add(pm);</span></p>
2535 
2536 <p class=MsoNormal><span style='font-size:10.0pt'>   pm-&gt;SetMarkerStyle(23);</span></p>
2537 
2538 <p class=MsoNormal><span style='font-size:10.0pt'>  
2539 pm-&gt;SetMarkerColor(kRed);</span></p>
2540 
2541 <p class=MsoNormal><span style='font-size:10.0pt'>   pm-&gt;SetMarkerSize(1);  
2542 </span></p>
2543 
2544 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
2545 
2546 </div>
2547 
2548 <!-- */
2549 // --> End_Html
2550 
2551  Int_t i, j, k, shift =
2552  2 * fNPeaks + 7, peak_vel, rozmer, iter, regul_cycle,
2553  flag;
2554  Double_t a, b, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
2555  0, pi, pmin = 0, chi_cel = 0, chi_er;
2556  Double_t *working_space = new Double_t[5 * (2 * fNPeaks + 7)];
2557  for (i = 0, j = 0; i < fNPeaks; i++) {
2558  working_space[2 * i] = fAmpInit[i]; //vector parameter
2559  if (fFixAmp[i] == false) {
2560  working_space[shift + j] = fAmpInit[i]; //vector xk
2561  j += 1;
2562  }
2563  working_space[2 * i + 1] = fPositionInit[i]; //vector parameter
2564  if (fFixPosition[i] == false) {
2565  working_space[shift + j] = fPositionInit[i]; //vector xk
2566  j += 1;
2567  }
2568  }
2569  peak_vel = 2 * i;
2570  working_space[2 * i] = fSigmaInit; //vector parameter
2571  if (fFixSigma == false) {
2572  working_space[shift + j] = fSigmaInit; //vector xk
2573  j += 1;
2574  }
2575  working_space[2 * i + 1] = fTInit; //vector parameter
2576  if (fFixT == false) {
2577  working_space[shift + j] = fTInit; //vector xk
2578  j += 1;
2579  }
2580  working_space[2 * i + 2] = fBInit; //vector parameter
2581  if (fFixB == false) {
2582  working_space[shift + j] = fBInit; //vector xk
2583  j += 1;
2584  }
2585  working_space[2 * i + 3] = fSInit; //vector parameter
2586  if (fFixS == false) {
2587  working_space[shift + j] = fSInit; //vector xk
2588  j += 1;
2589  }
2590  working_space[2 * i + 4] = fA0Init; //vector parameter
2591  if (fFixA0 == false) {
2592  working_space[shift + j] = fA0Init; //vector xk
2593  j += 1;
2594  }
2595  working_space[2 * i + 5] = fA1Init; //vector parameter
2596  if (fFixA1 == false) {
2597  working_space[shift + j] = fA1Init; //vector xk
2598  j += 1;
2599  }
2600  working_space[2 * i + 6] = fA2Init; //vector parameter
2601  if (fFixA2 == false) {
2602  working_space[shift + j] = fA2Init; //vector xk
2603  j += 1;
2604  }
2605  rozmer = j;
2606  if (rozmer == 0){
2607  Error ("FitAwmi","All parameters are fixed");
2608  delete [] working_space;
2609  return;
2610  }
2611  if (rozmer >= fXmax - fXmin + 1){
2612  Error ("FitAwmi","Number of fitted parameters is larger than # of fitted points");
2613  delete [] working_space;
2614  return;
2615  }
2616  Double_t **working_matrix = new Double_t *[rozmer];
2617  for (i = 0; i < rozmer; i++)
2618  working_matrix[i] = new Double_t[rozmer + 4];
2619  for (iter = 0; iter < fNumberIterations; iter++) {
2620  for (j = 0; j < rozmer; j++) {
2621  working_space[3 * shift + j] = 0; //temp
2622  for (k = 0; k <= rozmer; k++) {
2623  working_matrix[j][k] = 0;
2624  }
2625  }
2626 
2627  //filling working matrix
2628  alpha = fAlpha;
2629  chi_opt = 0;
2630  for (i = fXmin; i <= fXmax; i++) {
2631 
2632  //calculation of gradient vector
2633  for (j = 0, k = 0; j < fNPeaks; j++) {
2634  if (fFixAmp[j] == false) {
2635  working_space[2 * shift + k] =
2636  Deramp((Double_t) i, working_space[2 * j + 1],
2637  working_space[peak_vel],
2638  working_space[peak_vel + 1],
2639  working_space[peak_vel + 3],
2640  working_space[peak_vel + 2]);
2641  k += 1;
2642  }
2643  if (fFixPosition[j] == false) {
2644  working_space[2 * shift + k] =
2645  Deri0((Double_t) i, working_space[2 * j],
2646  working_space[2 * j + 1], working_space[peak_vel],
2647  working_space[peak_vel + 1],
2648  working_space[peak_vel + 3],
2649  working_space[peak_vel + 2]);
2650  k += 1;
2651  }
2652  } if (fFixSigma == false) {
2653  working_space[2 * shift + k] =
2654  Dersigma(fNPeaks, (Double_t) i, working_space,
2655  working_space[peak_vel],
2656  working_space[peak_vel + 1],
2657  working_space[peak_vel + 3],
2658  working_space[peak_vel + 2]);
2659  k += 1;
2660  }
2661  if (fFixT == false) {
2662  working_space[2 * shift + k] =
2663  Dert(fNPeaks, (Double_t) i, working_space,
2664  working_space[peak_vel], working_space[peak_vel + 2]);
2665  k += 1;
2666  }
2667  if (fFixB == false) {
2668  working_space[2 * shift + k] =
2669  Derb(fNPeaks, (Double_t) i, working_space,
2670  working_space[peak_vel], working_space[peak_vel + 1],
2671  working_space[peak_vel + 2]);
2672  k += 1;
2673  }
2674  if (fFixS == false) {
2675  working_space[2 * shift + k] =
2676  Ders(fNPeaks, (Double_t) i, working_space,
2677  working_space[peak_vel]);
2678  k += 1;
2679  }
2680  if (fFixA0 == false) {
2681  working_space[2 * shift + k] = 1.;
2682  k += 1;
2683  }
2684  if (fFixA1 == false) {
2685  working_space[2 * shift + k] = Dera1((Double_t) i);
2686  k += 1;
2687  }
2688  if (fFixA2 == false) {
2689  working_space[2 * shift + k] = Dera2((Double_t) i);
2690  k += 1;
2691  }
2692  yw = source[i];
2693  ywm = yw;
2694  f = Shape(fNPeaks, (Double_t) i, working_space,
2695  working_space[peak_vel], working_space[peak_vel + 1],
2696  working_space[peak_vel + 3],
2697  working_space[peak_vel + 2],
2698  working_space[peak_vel + 4],
2699  working_space[peak_vel + 5],
2700  working_space[peak_vel + 6]);
2701  if (fStatisticType == kFitOptimMaxLikelihood) {
2702  if (f > 0.00001)
2703  chi_opt += yw * TMath::Log(f) - f;
2704  }
2705 
2706  else {
2707  if (ywm != 0)
2708  chi_opt += (yw - f) * (yw - f) / ywm;
2709  }
2710  if (fStatisticType == kFitOptimChiFuncValues) {
2711  ywm = f;
2712  if (f < 0.00001)
2713  ywm = 0.00001;
2714  }
2715 
2716  else if (fStatisticType == kFitOptimMaxLikelihood) {
2717  ywm = f;
2718  if (f < 0.00001)
2719  ywm = 0.00001;
2720  }
2721 
2722  else {
2723  if (ywm == 0)
2724  ywm = 1;
2725  }
2726  for (j = 0; j < rozmer; j++) {
2727  for (k = 0; k < rozmer; k++) {
2728  b = working_space[2 * shift +
2729  j] * working_space[2 * shift + k] / ywm;
2730  if (fStatisticType == kFitOptimChiFuncValues)
2731  b = b * (4 * yw - 2 * f) / ywm;
2732  working_matrix[j][k] += b;
2733  if (j == k)
2734  working_space[3 * shift + j] += b;
2735  }
2736  }
2737  if (fStatisticType == kFitOptimChiFuncValues)
2738  b = (f * f - yw * yw) / (ywm * ywm);
2739 
2740  else
2741  b = (f - yw) / ywm;
2742  for (j = 0; j < rozmer; j++) {
2743  working_matrix[j][rozmer] -= b * working_space[2 * shift + j];
2744  }
2745  }
2746  for (i = 0; i < rozmer; i++) {
2747  working_matrix[i][rozmer + 1] = 0; //xk
2748  }
2749  StiefelInversion(working_matrix, rozmer);
2750  for (i = 0; i < rozmer; i++) {
2751  working_space[2 * shift + i] = working_matrix[i][rozmer + 1]; //der
2752  }
2753 
2754  //calculate chi_opt
2755  chi2 = chi_opt;
2756  chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
2757 
2758  //calculate new parameters
2759  regul_cycle = 0;
2760  for (j = 0; j < rozmer; j++) {
2761  working_space[4 * shift + j] = working_space[shift + j]; //temp_xk[j]=xk[j]
2762  }
2763 
2764  do {
2765  if (fAlphaOptim == kFitAlphaOptimal) {
2766  if (fStatisticType != kFitOptimMaxLikelihood)
2767  chi_min = 10000 * chi2;
2768 
2769  else
2770  chi_min = 0.1 * chi2;
2771  flag = 0;
2772  for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
2773  for (j = 0; j < rozmer; j++) {
2774  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]
2775  }
2776  for (i = 0, j = 0; i < fNPeaks; i++) {
2777  if (fFixAmp[i] == false) {
2778  if (working_space[shift + j] < 0) //xk[j]
2779  working_space[shift + j] = 0; //xk[j]
2780  working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
2781  j += 1;
2782  }
2783  if (fFixPosition[i] == false) {
2784  if (working_space[shift + j] < fXmin) //xk[j]
2785  working_space[shift + j] = fXmin; //xk[j]
2786  if (working_space[shift + j] > fXmax) //xk[j]
2787  working_space[shift + j] = fXmax; //xk[j]
2788  working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
2789  j += 1;
2790  }
2791  }
2792  if (fFixSigma == false) {
2793  if (working_space[shift + j] < 0.001) { //xk[j]
2794  working_space[shift + j] = 0.001; //xk[j]
2795  }
2796  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
2797  j += 1;
2798  }
2799  if (fFixT == false) {
2800  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
2801  j += 1;
2802  }
2803  if (fFixB == false) {
2804  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2805  if (working_space[shift + j] < 0) //xk[j]
2806  working_space[shift + j] = -0.001; //xk[j]
2807  else
2808  working_space[shift + j] = 0.001; //xk[j]
2809  }
2810  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
2811  j += 1;
2812  }
2813  if (fFixS == false) {
2814  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
2815  j += 1;
2816  }
2817  if (fFixA0 == false) {
2818  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
2819  j += 1;
2820  }
2821  if (fFixA1 == false) {
2822  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
2823  j += 1;
2824  }
2825  if (fFixA2 == false) {
2826  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
2827  j += 1;
2828  }
2829  chi2 = 0;
2830  for (i = fXmin; i <= fXmax; i++) {
2831  yw = source[i];
2832  ywm = yw;
2833  f = Shape(fNPeaks, (Double_t) i, working_space,
2834  working_space[peak_vel],
2835  working_space[peak_vel + 1],
2836  working_space[peak_vel + 3],
2837  working_space[peak_vel + 2],
2838  working_space[peak_vel + 4],
2839  working_space[peak_vel + 5],
2840  working_space[peak_vel + 6]);
2841  if (fStatisticType == kFitOptimChiFuncValues) {
2842  ywm = f;
2843  if (f < 0.00001)
2844  ywm = 0.00001;
2845  }
2846  if (fStatisticType == kFitOptimMaxLikelihood) {
2847  if (f > 0.00001)
2848  chi2 += yw * TMath::Log(f) - f;
2849  }
2850 
2851  else {
2852  if (ywm != 0)
2853  chi2 += (yw - f) * (yw - f) / ywm;
2854  }
2855  }
2856  if ((chi2 < chi_min
2857  && fStatisticType != kFitOptimMaxLikelihood)
2858  || (chi2 > chi_min
2859  && fStatisticType == kFitOptimMaxLikelihood)) {
2860  pmin = pi, chi_min = chi2;
2861  }
2862 
2863  else
2864  flag = 1;
2865  if (pi == 0.1)
2866  chi_min = chi2;
2867  chi = chi_min;
2868  }
2869  if (pmin != 0.1) {
2870  for (j = 0; j < rozmer; j++) {
2871  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]
2872  }
2873  for (i = 0, j = 0; i < fNPeaks; i++) {
2874  if (fFixAmp[i] == false) {
2875  if (working_space[shift + j] < 0) //xk[j]
2876  working_space[shift + j] = 0; //xk[j]
2877  working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
2878  j += 1;
2879  }
2880  if (fFixPosition[i] == false) {
2881  if (working_space[shift + j] < fXmin) //xk[j]
2882  working_space[shift + j] = fXmin; //xk[j]
2883  if (working_space[shift + j] > fXmax) //xk[j]
2884  working_space[shift + j] = fXmax; //xk[j]
2885  working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
2886  j += 1;
2887  }
2888  }
2889  if (fFixSigma == false) {
2890  if (working_space[shift + j] < 0.001) { //xk[j]
2891  working_space[shift + j] = 0.001; //xk[j]
2892  }
2893  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
2894  j += 1;
2895  }
2896  if (fFixT == false) {
2897  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
2898  j += 1;
2899  }
2900  if (fFixB == false) {
2901  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2902  if (working_space[shift + j] < 0) //xk[j]
2903  working_space[shift + j] = -0.001; //xk[j]
2904  else
2905  working_space[shift + j] = 0.001; //xk[j]
2906  }
2907  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
2908  j += 1;
2909  }
2910  if (fFixS == false) {
2911  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
2912  j += 1;
2913  }
2914  if (fFixA0 == false) {
2915  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
2916  j += 1;
2917  }
2918  if (fFixA1 == false) {
2919  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
2920  j += 1;
2921  }
2922  if (fFixA2 == false) {
2923  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
2924  j += 1;
2925  }
2926  chi = chi_min;
2927  }
2928  }
2929 
2930  else {
2931  for (j = 0; j < rozmer; j++) {
2932  working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+alpha*der[j]
2933  }
2934  for (i = 0, j = 0; i < fNPeaks; i++) {
2935  if (fFixAmp[i] == false) {
2936  if (working_space[shift + j] < 0) //xk[j]
2937  working_space[shift + j] = 0; //xk[j]
2938  working_space[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
2939  j += 1;
2940  }
2941  if (fFixPosition[i] == false) {
2942  if (working_space[shift + j] < fXmin) //xk[j]
2943  working_space[shift + j] = fXmin; //xk[j]
2944  if (working_space[shift + j] > fXmax) //xk[j]
2945  working_space[shift + j] = fXmax; //xk[j]
2946  working_space[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
2947  j += 1;
2948  }
2949  }
2950  if (fFixSigma == false) {
2951  if (working_space[shift + j] < 0.001) { //xk[j]
2952  working_space[shift + j] = 0.001; //xk[j]
2953  }
2954  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
2955  j += 1;
2956  }
2957  if (fFixT == false) {
2958  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
2959  j += 1;
2960  }
2961  if (fFixB == false) {
2962  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2963  if (working_space[shift + j] < 0) //xk[j]
2964  working_space[shift + j] = -0.001; //xk[j]
2965  else
2966  working_space[shift + j] = 0.001; //xk[j]
2967  }
2968  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
2969  j += 1;
2970  }
2971  if (fFixS == false) {
2972  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
2973  j += 1;
2974  }
2975  if (fFixA0 == false) {
2976  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
2977  j += 1;
2978  }
2979  if (fFixA1 == false) {
2980  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
2981  j += 1;
2982  }
2983  if (fFixA2 == false) {
2984  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
2985  j += 1;
2986  }
2987  chi = 0;
2988  for (i = fXmin; i <= fXmax; i++) {
2989  yw = source[i];
2990  ywm = yw;
2991  f = Shape(fNPeaks, (Double_t) i, working_space,
2992  working_space[peak_vel],
2993  working_space[peak_vel + 1],
2994  working_space[peak_vel + 3],
2995  working_space[peak_vel + 2],
2996  working_space[peak_vel + 4],
2997  working_space[peak_vel + 5],
2998  working_space[peak_vel + 6]);
2999  if (fStatisticType == kFitOptimChiFuncValues) {
3000  ywm = f;
3001  if (f < 0.00001)
3002  ywm = 0.00001;
3003  }
3004  if (fStatisticType == kFitOptimMaxLikelihood) {
3005  if (f > 0.00001)
3006  chi += yw * TMath::Log(f) - f;
3007  }
3008 
3009  else {
3010  if (ywm != 0)
3011  chi += (yw - f) * (yw - f) / ywm;
3012  }
3013  }
3014  }
3015  chi2 = chi;
3016  chi = TMath::Sqrt(TMath::Abs(chi));
3017  if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
3018  alpha = alpha * chi_opt / (2 * chi);
3019 
3020  else if (fAlphaOptim == kFitAlphaOptimal)
3021  alpha = alpha / 10.0;
3022  iter += 1;
3023  regul_cycle += 1;
3024  } while (((chi > chi_opt
3025  && fStatisticType != kFitOptimMaxLikelihood)
3026  || (chi < chi_opt
3027  && fStatisticType == kFitOptimMaxLikelihood))
3028  && regul_cycle < kFitNumRegulCycles);
3029  for (j = 0; j < rozmer; j++) {
3030  working_space[4 * shift + j] = 0; //temp_xk[j]
3031  working_space[2 * shift + j] = 0; //der[j]
3032  }
3033  for (i = fXmin, chi_cel = 0; i <= fXmax; i++) {
3034  yw = source[i];
3035  if (yw == 0)
3036  yw = 1;
3037  f = Shape(fNPeaks, (Double_t) i, working_space,
3038  working_space[peak_vel], working_space[peak_vel + 1],
3039  working_space[peak_vel + 3],
3040  working_space[peak_vel + 2],
3041  working_space[peak_vel + 4],
3042  working_space[peak_vel + 5],
3043  working_space[peak_vel + 6]);
3044  chi_opt = (yw - f) * (yw - f) / yw;
3045  chi_cel += (yw - f) * (yw - f) / yw;
3046 
3047  //calculate gradient vector
3048  for (j = 0, k = 0; j < fNPeaks; j++) {
3049  if (fFixAmp[j] == false) {
3050  a = Deramp((Double_t) i, working_space[2 * j + 1],
3051  working_space[peak_vel],
3052  working_space[peak_vel + 1],
3053  working_space[peak_vel + 3],
3054  working_space[peak_vel + 2]);
3055  if (yw != 0) {
3056  working_space[2 * shift + k] += chi_opt; //der[k]
3057  b = a * a / yw;
3058  working_space[4 * shift + k] += b; //temp_xk[k]
3059  }
3060  k += 1;
3061  }
3062  if (fFixPosition[j] == false) {
3063  a = Deri0((Double_t) i, working_space[2 * j],
3064  working_space[2 * j + 1],
3065  working_space[peak_vel],
3066  working_space[peak_vel + 1],
3067  working_space[peak_vel + 3],
3068  working_space[peak_vel + 2]);
3069  if (yw != 0) {
3070  working_space[2 * shift + k] += chi_opt; //der[k]
3071  b = a * a / yw;
3072  working_space[4 * shift + k] += b; //temp_xk[k]
3073  }
3074  k += 1;
3075  }
3076  }
3077  if (fFixSigma == false) {
3078  a = Dersigma(fNPeaks, (Double_t) i, working_space,
3079  working_space[peak_vel],
3080  working_space[peak_vel + 1],
3081  working_space[peak_vel + 3],
3082  working_space[peak_vel + 2]);
3083  if (yw != 0) {
3084  working_space[2 * shift + k] += chi_opt; //der[k]
3085  b = a * a / yw;
3086  working_space[4 * shift + k] += b; //temp_xk[k]
3087  }
3088  k += 1;
3089  }
3090  if (fFixT == false) {
3091  a = Dert(fNPeaks, (Double_t) i, working_space,
3092  working_space[peak_vel],
3093  working_space[peak_vel + 2]);
3094  if (yw != 0) {
3095  working_space[2 * shift + k] += chi_opt; //der[k]
3096  b = a * a / yw;
3097  working_space[4 * shift + k] += b; //temp_xk[k]
3098  }
3099  k += 1;
3100  }
3101  if (fFixB == false) {
3102  a = Derb(fNPeaks, (Double_t) i, working_space,
3103  working_space[peak_vel], working_space[peak_vel + 1],
3104  working_space[peak_vel + 2]);
3105  if (yw != 0) {
3106  working_space[2 * shift + k] += chi_opt; //der[k]
3107  b = a * a / yw;
3108  working_space[4 * shift + k] += b; //temp_xk[k]
3109  }
3110  k += 1;
3111  }
3112  if (fFixS == false) {
3113  a = Ders(fNPeaks, (Double_t) i, working_space,
3114  working_space[peak_vel]);
3115  if (yw != 0) {
3116  working_space[2 * shift + k] += chi_opt; //der[k]
3117  b = a * a / yw;
3118  working_space[4 * shift + k] += b; //tem_xk[k]
3119  }
3120  k += 1;
3121  }
3122  if (fFixA0 == false) {
3123  a = 1.0;
3124  if (yw != 0) {
3125  working_space[2 * shift + k] += chi_opt; //der[k]
3126  b = a * a / yw;
3127  working_space[4 * shift + k] += b; //temp_xk[k]
3128  }
3129  k += 1;
3130  }
3131  if (fFixA1 == false) {
3132  a = Dera1((Double_t) i);
3133  if (yw != 0) {
3134  working_space[2 * shift + k] += chi_opt; //der[k]
3135  b = a * a / yw;
3136  working_space[4 * shift + k] += b; //temp_xk[k]
3137  }
3138  k += 1;
3139  }
3140  if (fFixA2 == false) {
3141  a = Dera2((Double_t) i);
3142  if (yw != 0) {
3143  working_space[2 * shift + k] += chi_opt; //der[k]
3144  b = a * a / yw;
3145  working_space[4 * shift + k] += b; //temp_xk[k]
3146  }
3147  k += 1;
3148  }
3149  }
3150  }
3151  b = fXmax - fXmin + 1 - rozmer;
3152  chi_er = chi_cel / b;
3153  for (i = 0, j = 0; i < fNPeaks; i++) {
3154  fArea[i] =
3155  Area(working_space[2 * i], working_space[peak_vel],
3156  working_space[peak_vel + 1], working_space[peak_vel + 2]);
3157  if (fFixAmp[i] == false) {
3158  fAmpCalc[i] = working_space[shift + j]; //xk[j]
3159  if (working_space[3 * shift + j] != 0)
3160  fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3161  if (fArea[i] > 0) {
3162  a = Derpa(working_space[peak_vel],
3163  working_space[peak_vel + 1],
3164  working_space[peak_vel + 2]);
3165  b = working_space[4 * shift + j]; //temp_xk[j]
3166  if (b == 0)
3167  b = 1;
3168 
3169  else
3170  b = 1 / b;
3171  fAreaErr[i] = TMath::Sqrt(TMath::Abs(a * a * b * chi_er));
3172  }
3173 
3174  else
3175  fAreaErr[i] = 0;
3176  j += 1;
3177  }
3178 
3179  else {
3180  fAmpCalc[i] = fAmpInit[i];
3181  fAmpErr[i] = 0;
3182  fAreaErr[i] = 0;
3183  }
3184  if (fFixPosition[i] == false) {
3185  fPositionCalc[i] = working_space[shift + j]; //xk[j]
3186  if (working_space[3 * shift + j] != 0) //temp[j]
3187  fPositionErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //Der[j]/temp[j]
3188  j += 1;
3189  }
3190 
3191  else {
3192  fPositionCalc[i] = fPositionInit[i];
3193  fPositionErr[i] = 0;
3194  }
3195  }
3196  if (fFixSigma == false) {
3197  fSigmaCalc = working_space[shift + j]; //xk[j]
3198  if (working_space[3 * shift + j] != 0) //temp[j]
3199  fSigmaErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3200  j += 1;
3201  }
3202 
3203  else {
3204  fSigmaCalc = fSigmaInit;
3205  fSigmaErr = 0;
3206  }
3207  if (fFixT == false) {
3208  fTCalc = working_space[shift + j]; //xk[j]
3209  if (working_space[3 * shift + j] != 0) //temp[j]
3210  fTErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3211  j += 1;
3212  }
3213 
3214  else {
3215  fTCalc = fTInit;
3216  fTErr = 0;
3217  }
3218  if (fFixB == false) {
3219  fBCalc = working_space[shift + j]; //xk[j]
3220  if (working_space[3 * shift + j] != 0) //temp[j]
3221  fBErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3222  j += 1;
3223  }
3224 
3225  else {
3226  fBCalc = fBInit;
3227  fBErr = 0;
3228  }
3229  if (fFixS == false) {
3230  fSCalc = working_space[shift + j]; //xk[j]
3231  if (working_space[3 * shift + j] != 0) //temp[j]
3232  fSErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3233  j += 1;
3234  }
3235 
3236  else {
3237  fSCalc = fSInit;
3238  fSErr = 0;
3239  }
3240  if (fFixA0 == false) {
3241  fA0Calc = working_space[shift + j]; //xk[j]
3242  if (working_space[3 * shift + j] != 0) //temp[j]
3243  fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3244  j += 1;
3245  }
3246 
3247  else {
3248  fA0Calc = fA0Init;
3249  fA0Err = 0;
3250  }
3251  if (fFixA1 == false) {
3252  fA1Calc = working_space[shift + j]; //xk[j]
3253  if (working_space[3 * shift + j] != 0) //temp[j]
3254  fA1Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3255  j += 1;
3256  }
3257 
3258  else {
3259  fA1Calc = fA1Init;
3260  fA1Err = 0;
3261  }
3262  if (fFixA2 == false) {
3263  fA2Calc = working_space[shift + j]; //xk[j]
3264  if (working_space[3 * shift + j] != 0) //temp[j]
3265  fA2Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3266  j += 1;
3267  }
3268 
3269  else {
3270  fA2Calc = fA2Init;
3271  fA2Err = 0;
3272  }
3273  b = fXmax - fXmin + 1 - rozmer;
3274  fChi = chi_cel / b;
3275  for (i = fXmin; i <= fXmax; i++) {
3276  f = Shape(fNPeaks, (Double_t) i, working_space,
3277  working_space[peak_vel], working_space[peak_vel + 1],
3278  working_space[peak_vel + 3], working_space[peak_vel + 2],
3279  working_space[peak_vel + 4], working_space[peak_vel + 5],
3280  working_space[peak_vel + 6]);
3281  source[i] = f;
3282  }
3283  for (i = 0; i < rozmer; i++)
3284  delete [] working_matrix[i];
3285  delete [] working_matrix;
3286  delete [] working_space;
3287  return;
3288 }
3289 
3290 ////////////////////////////////////////////////////////////////////////////////
3291 ///////////////////////////////////////////////////////////////////////////////
3292 /// SETTER FUNCTION
3293 ///
3294 /// This function sets the following fitting parameters:
3295 /// -xmin, xmax - fitting region
3296 /// -numberIterations - # of desired iterations in the fit
3297 /// -alpha - convergence coefficient, it should be positive number and <=1, for details see references
3298 /// -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
3299 /// -alphaOptim - optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal
3300 /// -power - possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting function.
3301 /// -fitTaylor - order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi fitting function.
3302 ///////////////////////////////////////////////////////////////////////////////
3303 
3304 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)
3305 {
3306  if(xmin<0 || xmax <= xmin){
3307  Error("SetFitParameters", "Wrong range");
3308  return;
3309  }
3310  if (numberIterations <= 0){
3311  Error("SetFitParameters","Invalid number of iterations, must be positive");
3312  return;
3313  }
3314  if (alpha <= 0 || alpha > 1){
3315  Error ("SetFitParameters","Invalid step coefficient alpha, must be > than 0 and <=1");
3316  return;
3317  }
3318  if (statisticType != kFitOptimChiCounts
3319  && statisticType != kFitOptimChiFuncValues
3320  && statisticType != kFitOptimMaxLikelihood){
3321  Error("SetFitParameters","Wrong type of statistic");
3322  return;
3323  }
3324  if (alphaOptim != kFitAlphaHalving
3325  && alphaOptim != kFitAlphaOptimal){
3326  Error("SetFitParameters","Wrong optimization algorithm");
3327  return;
3328  }
3329  if (power != kFitPower2 && power != kFitPower4
3330  && power != kFitPower6 && power != kFitPower8
3331  && power != kFitPower10 && power != kFitPower12){
3332  Error("SetFitParameters","Wrong power");
3333  return;
3334  }
3335  if (fitTaylor != kFitTaylorOrderFirst
3336  && fitTaylor != kFitTaylorOrderSecond){
3337  Error("SetFitParameters","Wrong order of Taylor development");
3338  return;
3339  }
3340  fXmin=xmin,fXmax=xmax,fNumberIterations=numberIterations,fAlpha=alpha,fStatisticType=statisticType,fAlphaOptim=alphaOptim,fPower=power,fFitTaylor=fitTaylor;
3341 }
3342 
3343 ////////////////////////////////////////////////////////////////////////////////
3344 ///////////////////////////////////////////////////////////////////////////////
3345 /// SETTER FUNCTION
3346 ///
3347 /// This function sets the following fitting parameters of peaks:
3348 /// -sigma - initial value of sigma parameter
3349 /// -fixSigma - logical value of sigma parameter, which allows to fix the parameter (not to fit)
3350 /// -positionInit - aray of initial values of peaks positions
3351 /// -fixPosition - array of logical values which allow to fix appropriate positions (not fit). However they are present in the estimated functional.
3352 /// -ampInit - aray of initial values of peaks amplitudes
3353 /// -fixAmp - aray of logical values which allow to fix appropriate amplitudes (not fit). However they are present in the estimated functional
3354 ///////////////////////////////////////////////////////////////////////////////
3355 
3356 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)
3357 {
3358  Int_t i;
3359  if (sigma <= 0){
3360  Error ("SetPeakParameters","Invalid sigma, must be > than 0");
3361  return;
3362  }
3363  for(i=0; i < fNPeaks; i++){
3364  if((Int_t) positionInit[i] < fXmin || (Int_t) positionInit[i] > fXmax){
3365  Error ("SetPeakParameters","Invalid peak position, must be in the range fXmin, fXmax");
3366  return;
3367  }
3368  if(ampInit[i] < 0){
3369  Error ("SetPeakParameters","Invalid peak amplitude, must be > than 0");
3370  return;
3371  }
3372  }
3373  fSigmaInit = sigma,fFixSigma = fixSigma;
3374  for(i=0; i < fNPeaks; i++){
3375  fPositionInit[i] = (Double_t) positionInit[i];
3376  fFixPosition[i] = fixPosition[i];
3377  fAmpInit[i] = (Double_t) ampInit[i];
3378  fFixAmp[i] = fixAmp[i];
3379  }
3380 }
3381 
3382 ////////////////////////////////////////////////////////////////////////////////
3383 ///////////////////////////////////////////////////////////////////////////////
3384 /// SETTER FUNCTION
3385 ///
3386 /// This function sets the following fitting parameters of background:
3387 /// -a0Init - initial value of a0 parameter (backgroud is estimated as a0+a1*x+a2*x*x)
3388 /// -fixA0 - logical value of a0 parameter, which allows to fix the parameter (not to fit)
3389 /// -a1Init - initial value of a1 parameter
3390 /// -fixA1 - logical value of a1 parameter, which allows to fix the parameter (not to fit)
3391 /// -a2Init - initial value of a2 parameter
3392 /// -fixA2 - logical value of a2 parameter, which allows to fix the parameter (not to fit)
3393 ///////////////////////////////////////////////////////////////////////////////
3394 
3395 void TSpectrumFit::SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t a1Init, Bool_t fixA1, Double_t a2Init, Bool_t fixA2)
3396 {
3397  fA0Init = a0Init;
3398  fFixA0 = fixA0;
3399  fA1Init = a1Init;
3400  fFixA1 = fixA1;
3401  fA2Init = a2Init;
3402  fFixA2 = fixA2;
3403 }
3404 
3405 ////////////////////////////////////////////////////////////////////////////////
3406 ///////////////////////////////////////////////////////////////////////////////
3407 /// SETTER FUNCTION
3408 ///
3409 /// This function sets the following fitting parameters of tails of peaks
3410 /// -tInit - initial value of t parameter
3411 /// -fixT - logical value of t parameter, which allows to fix the parameter (not to fit)
3412 /// -bInit - initial value of b parameter
3413 /// -fixB - logical value of b parameter, which allows to fix the parameter (not to fit)
3414 /// -sInit - initial value of s parameter
3415 /// -fixS - logical value of s parameter, which allows to fix the parameter (not to fit)
3416 ///////////////////////////////////////////////////////////////////////////////
3417 
3418 void TSpectrumFit::SetTailParameters(Double_t tInit, Bool_t fixT, Double_t bInit, Bool_t fixB, Double_t sInit, Bool_t fixS)
3419 {
3420  fTInit = tInit;
3421  fFixT = fixT;
3422  fBInit = bInit;
3423  fFixB = fixB;
3424  fSInit = sInit;
3425  fFixS = fixS;
3426 }
3427 
3428 ////////////////////////////////////////////////////////////////////////////////
3429 ///////////////////////////////////////////////////////////////////////////////
3430 /// GETTER FUNCTION
3431 ///
3432 /// This function gets the sigma parameter and its error
3433 /// -sigma - gets the fitted value of sigma parameter
3434 /// -sigmaErr - gets error value of sigma parameter
3435 ///////////////////////////////////////////////////////////////////////////////
3436 
3437 void TSpectrumFit::GetSigma(Double_t &sigma, Double_t &sigmaErr)
3438 {
3439  sigma=fSigmaCalc;
3440  sigmaErr=fSigmaErr;
3441 }
3442 
3443 ////////////////////////////////////////////////////////////////////////////////
3444 ///////////////////////////////////////////////////////////////////////////////
3445 /// GETTER FUNCTION
3446 ///
3447 /// This function gets the background parameters and their errors
3448 /// -a0 - gets the fitted value of a0 parameter
3449 /// -a0Err - gets error value of a0 parameter
3450 /// -a1 - gets the fitted value of a1 parameter
3451 /// -a1Err - gets error value of a1 parameter
3452 /// -a2 - gets the fitted value of a2 parameter
3453 /// -a2Err - gets error value of a2 parameter
3454 ///////////////////////////////////////////////////////////////////////////////
3455 
3456 void TSpectrumFit::GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &a1, Double_t &a1Err, Double_t &a2, Double_t &a2Err)
3457 {
3458  a0 = fA0Calc;
3459  a0Err = fA0Err;
3460  a1 = fA1Calc;
3461  a1Err = fA1Err;
3462  a2 = fA2Calc;
3463  a2Err = fA2Err;
3464 }
3465 
3466 ////////////////////////////////////////////////////////////////////////////////
3467 ///////////////////////////////////////////////////////////////////////////////
3468 /// GETTER FUNCTION
3469 ///
3470 /// This function gets the tail parameters and their errors
3471 /// -t - gets the fitted value of t parameter
3472 /// -tErr - gets error value of t parameter
3473 /// -b - gets the fitted value of b parameter
3474 /// -bErr - gets error value of b parameter
3475 /// -s - gets the fitted value of s parameter
3476 /// -sErr - gets error value of s parameter
3477 ///////////////////////////////////////////////////////////////////////////////
3478 
3479 void TSpectrumFit::GetTailParameters(Double_t &t, Double_t &tErr, Double_t &b, Double_t &bErr, Double_t &s, Double_t &sErr)
3480 {
3481  t = fTCalc;
3482  tErr = fTErr;
3483  b = fBCalc;
3484  bErr = fBErr;
3485  s = fSCalc;
3486  sErr = fSErr;
3487 }
3488 
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)
AUXILIARY FUNCTION // // This function calculates peaks shape function (see manual) // Function param...
Definition: TSpectrumFit.cxx:637
TSpectrumFit::fAlpha
Double_t fAlpha
Definition: TSpectrumFit.h:42
xmin
float xmin
Definition: THbookFile.cxx:93
TSpectrumFit::fAreaErr
Double_t * fAreaErr
Definition: TSpectrumFit.h:51
TSpectrumFit::~TSpectrumFit
virtual ~TSpectrumFit()
destructor
Definition: TSpectrumFit.cxx:226
TSpectrumFit::fFixSigma
Bool_t fFixSigma
Definition: TSpectrumFit.h:75
TSpectrumFit::Derpb
Double_t Derpb(Double_t a, Double_t sigma, Double_t t, Double_t b)
AUXILIARY FUNCTION // // This function calculates derivative of the area of peak // according to b pa...
Definition: TSpectrumFit.cxx:801
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)
SETTER FUNCTION.
Definition: TSpectrumFit.cxx:3356
TSpectrumFit::fAmpInit
Double_t * fAmpInit
Definition: TSpectrumFit.h:47
TMath::Log
Double_t Log(Double_t x)
Definition: TMath.h:526
TSpectrumFit::fPositionCalc
Double_t * fPositionCalc
Definition: TSpectrumFit.h:45
pi
const double pi
Definition: rotationApplication.cxx:282
TSpectrumFit::Area
Double_t Area(Double_t a, Double_t sigma, Double_t t, Double_t b)
AUXILIARY FUNCTION // // This function calculates area of a peak // Function parameters: // -a-amplit...
Definition: TSpectrumFit.cxx:695
TSpectrumFit::StiefelInversion
void StiefelInversion(Double_t **a, Int_t rozmer)
AUXILIARY FUNCTION // // This function calculates solution of the system of linear equations...
Definition: TSpectrumFit.cxx:2253
TSpectrumFit::fFitTaylor
Int_t fFitTaylor
Definition: TSpectrumFit.h:41
TSpectrumFit::fStatisticType
Int_t fStatisticType
Definition: TSpectrumFit.h:38
Int_t
int Int_t
Definition: RtypesCore.h:41
TSpectrumFit::Dert
Double_t Dert(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma, Double_t b)
AUXILIARY FUNCTION // // This function calculates derivative of peaks shape function (see manual) // ...
Definition: TSpectrumFit.cxx:515
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
a
TArc * a
Definition: textangle.C:12
TSpectrumFit::Derpsigma
Double_t Derpsigma(Double_t a, Double_t t, Double_t b)
Definition: TSpectrumFit.cxx:736
TSpectrumFit::fPositionErr
Double_t * fPositionErr
Definition: TSpectrumFit.h:46
TSpectrumFit::fTErr
Double_t fTErr
Definition: TSpectrumFit.h:57
TSpectrumFit::GetTailParameters
void GetTailParameters(Double_t &t, Double_t &tErr, Double_t &b, Double_t &bErr, Double_t &s, Double_t &sErr)
GETTER FUNCTION.
Definition: TSpectrumFit.cxx:3479
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
TSpectrumFit::fA0Init
Double_t fA0Init
Definition: TSpectrumFit.h:64
TSpectrumFit::fFixS
Bool_t fFixS
Definition: TSpectrumFit.h:78
TSpectrumFit::fAlphaOptim
Int_t fAlphaOptim
Definition: TSpectrumFit.h:39
x
Double_t x[n]
Definition: legend1.C:17
TNamed
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:33
r3
unsigned int r3[N_CITIES]
Definition: simanTSP.cxx:323
TSpectrumFit::fA1Calc
Double_t fA1Calc
Definition: TSpectrumFit.h:68
TSpectrumFit::fA2Init
Double_t fA2Init
Definition: TSpectrumFit.h:70
TSpectrumFit::SetBackgroundParameters
void SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t a1Init, Bool_t fixA1, Double_t a2Init, Bool_t fixA2)
SETTER FUNCTION.
Definition: TSpectrumFit.cxx:3395
iter
std::map< std::string, std::string >::const_iterator iter
Definition: TAlienJob.cxx:54
TSpectrumFit::fPositionInit
Double_t * fPositionInit
Definition: TSpectrumFit.h:44
TSpectrumFit::Dera2
Double_t Dera2(Double_t i)
derivative of background according to a2
Definition: TSpectrumFit.cxx:616
TSpectrumFit::fSigmaCalc
Double_t fSigmaCalc
Definition: TSpectrumFit.h:53
TSpectrumFit::fNumberIterations
Int_t fNumberIterations
Definition: TSpectrumFit.h:35
TObject::Error
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
TSpectrumFit::Derfc
Double_t Derfc(Double_t x)
AUXILIARY FUNCTION // // This function calculates derivative of error function of x...
Definition: TSpectrumFit.cxx:276
TSpectrumFit::fSCalc
Double_t fSCalc
Definition: TSpectrumFit.h:62
TSpectrumFit::fChi
Double_t fChi
Definition: TSpectrumFit.h:43
TSpectrumFit::fA2Err
Double_t fA2Err
Definition: TSpectrumFit.h:72
TSpectrumFit::fFixB
Bool_t fFixB
Definition: TSpectrumFit.h:77
TSpectrumFit::fA1Init
Double_t fA1Init
Definition: TSpectrumFit.h:67
TSpectrumFit::fTInit
Double_t fTInit
Definition: TSpectrumFit.h:55
TSpectrumFit::fSErr
Double_t fSErr
Definition: TSpectrumFit.h:63
TSpectrumFit::fBCalc
Double_t fBCalc
Definition: TSpectrumFit.h:59
TSpectrumFit::fAmpErr
Double_t * fAmpErr
Definition: TSpectrumFit.h:49
r
ROOT::R::TRInterface & r
Definition: Object.C:4
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)
AUXILIARY FUNCTION // // This function calculates derivative of peaks shape function (see manual) // ...
Definition: TSpectrumFit.cxx:427
r1
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
TSpectrumFit::fBInit
Double_t fBInit
Definition: TSpectrumFit.h:58
TMath::E
Double_t E()
Definition: TMath.h:54
TSpectrumFit::Erfc
Double_t Erfc(Double_t x)
Definition: TSpectrumFit.cxx:242
TSpectrumFit::fBErr
Double_t fBErr
Definition: TSpectrumFit.h:60
TSpectrumFit::Deramp
Double_t Deramp(Double_t i, Double_t i0, Double_t sigma, Double_t t, Double_t s, Double_t b)
AUXILIARY FUNCTION // // This function calculates derivative of peak shape function (see manual) // a...
Definition: TSpectrumFit.cxx:310
xmax
float xmax
Definition: THbookFile.cxx:93
TSpectrumFit::fA1Err
Double_t fA1Err
Definition: TSpectrumFit.h:69
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)
SETTER FUNCTION.
Definition: TSpectrumFit.cxx:3304
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)
AUXILIARY FUNCTION // // This function calculates derivative of peak shape function (see manual) // a...
Definition: TSpectrumFit.cxx:352
TSpectrumFit::fFixA0
Bool_t fFixA0
Definition: TSpectrumFit.h:79
ClassImp
ClassImp(TSpectrumFit) TSpectrumFit
default constructor
Definition: TSpectrumFit.cxx:55
TSpectrumFit::Derpa
Double_t Derpa(Double_t sigma, Double_t t, Double_t b)
AUXILIARY FUNCTION // // This function calculates derivative of the area of peak // according to its ...
Definition: TSpectrumFit.cxx:723
TSpectrumFit::Ders
Double_t Ders(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma)
AUXILIARY FUNCTION // // This function calculates derivative of peaks shape function (see manual) // ...
Definition: TSpectrumFit.cxx:548
f
double f(double x)
Definition: testIntegration.cxx:12
TSpectrumFit::fA0Err
Double_t fA0Err
Definition: TSpectrumFit.h:66
Double_t
double Double_t
Definition: RtypesCore.h:55
TSpectrumFit::Dera1
Double_t Dera1(Double_t i)
derivative of background according to a1
Definition: TSpectrumFit.cxx:608
TSpectrumFit::Derderi0
Double_t Derderi0(Double_t i, Double_t amp, Double_t i0, Double_t sigma)
AUXILIARY FUNCTION // // This function calculates second derivative of peak shape function // (see ma...
Definition: TSpectrumFit.cxx:394
TSpectrumFit::Derdersigma
Double_t Derdersigma(Int_t num_of_fitted_peaks, Double_t i, const Double_t *parameter, Double_t sigma)
AUXILIARY FUNCTION // // This function calculates second derivative of peaks shape function // (see m...
Definition: TSpectrumFit.cxx:477
TSpectrumFit::Derpt
Double_t Derpt(Double_t a, Double_t sigma, Double_t b)
AUXILIARY FUNCTION // // This function calculates derivative of the area of peak // according to t pa...
Definition: TSpectrumFit.cxx:774
TSpectrumFit::fFixPosition
Bool_t * fFixPosition
Definition: TSpectrumFit.h:73
TSpectrumFit::fFixAmp
Bool_t * fFixAmp
Definition: TSpectrumFit.h:74
TSpectrumFit::fSInit
Double_t fSInit
Definition: TSpectrumFit.h:61
TSpectrumFit::TSpectrumFit
TSpectrumFit(void)
TSpectrumFit::SetTailParameters
void SetTailParameters(Double_t tInit, Bool_t fixT, Double_t bInit, Bool_t fixB, Double_t sInit, Bool_t fixS)
SETTER FUNCTION.
Definition: TSpectrumFit.cxx:3418
TSpectrumFit::Ourpowl
Double_t Ourpowl(Double_t a, Int_t pw)
power function
Definition: TSpectrumFit.cxx:817
sample_config.c
tuple c
Definition: sample_config.py:3
TSpectrumFit
Advanced 1-dimentional spectra fitting functions.
Definition: TSpectrumFit.h:32
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)
AUXILIARY FUNCTION // // This function calculates derivative of peaks shape function (see manual) // ...
Definition: TSpectrumFit.cxx:579
TSpectrumFit::GetBackgroundParameters
void GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &a1, Double_t &a1Err, Double_t &a2, Double_t &a2Err)
GETTER FUNCTION.
Definition: TSpectrumFit.cxx:3456
TSpectrumFit::fArea
Double_t * fArea
Definition: TSpectrumFit.h:50
TSpectrumFit::fSigmaInit
Double_t fSigmaInit
Definition: TSpectrumFit.h:52
TSpectrumFit::fA2Calc
Double_t fA2Calc
Definition: TSpectrumFit.h:71
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
TSpectrumFit::fA0Calc
Double_t fA0Calc
Definition: TSpectrumFit.h:65
TSpectrumFit::fTCalc
Double_t fTCalc
Definition: TSpectrumFit.h:56
exp
double exp(double)
q
float * q
Definition: THbookFile.cxx:87
TSpectrumFit::GetSigma
void GetSigma(Double_t &sigma, Double_t &sigmaErr)
GETTER FUNCTION.
Definition: TSpectrumFit.cxx:3437
TSpectrumFit::FitAwmi
void FitAwmi(Double_t *source)
ONE-DIMENSIONAL FIT FUNCTION ALGORITHM WITHOUT MATRIX INVERSION This function fits the source spectru...
Definition: TSpectrumFit.cxx:852
TSpectrumFit::FitStiefel
void FitStiefel(Double_t *source)
ONE-DIMENSIONAL FIT FUNCTION ALGORITHM WITH MATRIX INVERSION (STIEFEL-HESTENS METHOD) This function f...
Definition: TSpectrumFit.cxx:2317
TSpectrumFit::fAmpCalc
Double_t * fAmpCalc
Definition: TSpectrumFit.h:48
r2
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
TSpectrumFit::fSigmaErr
Double_t fSigmaErr
Definition: TSpectrumFit.h:54
TSpectrumFit::fFixT
Bool_t fFixT
Definition: TSpectrumFit.h:76
TSpectrumFit::fFixA1
Bool_t fFixA1
Definition: TSpectrumFit.h:80
TSpectrumFit::fFixA2
Bool_t fFixA2
Definition: TSpectrumFit.h:81