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