Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TSpectrum2Fit.cxx
Go to the documentation of this file.
1 // @(#)root/spectrum:$Id$
2 // Author: Miroslav Morhac 25/09/2006
3 
4 /** \class TSpectrum2Fit
5  \ingroup Spectrum
6  \brief Advanced 2-dimensional spectra fitting functions
7  \author Miroslav Morhac
8 
9 Class for fitting 2D spectra using AWMI (algorithm without matrix
10 inversion) and conjugate gradient algorithms for symmetrical
11 matrices (Stiefel-Hestens method). AWMI method allows to fit
12 simultaneously 100s up to 1000s peaks. Stiefel method is very stable,
13 it converges faster, but is more time consuming.
14 
15 
16  The algorithms in this class have been published in the following references:
17 
18  1. M. Morhac et al.: Efficient fitting algorithms applied to
19  analysis of coincidence gamma-ray spectra. Computer Physics
20  Communications, Vol 172/1 (2005) pp. 19-41.
21 
22  2. M. Morhac et al.: Study of fitting algorithms applied to
23  simultaneous analysis of large number of peaks in gamma-ray spectra.
24  Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003.
25 */
26 
27 #include "TSpectrum2Fit.h"
28 #include "TMath.h"
29 
30 ClassImp(TSpectrum2Fit);
31 
32 ////////////////////////////////////////////////////////////////////////////////
33 /// Default constructor
34 
35 TSpectrum2Fit::TSpectrum2Fit() :TNamed("Spectrum2Fit", "Miroslav Morhac peak fitter")
36 {
37  fNPeaks = 0;
38  fNumberIterations = 1;
39  fXmin = 0;
40  fXmax = 100;
41  fYmin = 0;
42  fYmax = 100;
43  fStatisticType = kFitOptimChiCounts;
44  fAlphaOptim = kFitAlphaHalving;
45  fPower = kFitPower2;
46  fFitTaylor = kFitTaylorOrderFirst;
47  fAlpha = 1;
48  fChi = 0;
49  fPositionInitX = 0;
50  fPositionCalcX = 0;
51  fPositionErrX = 0;
52  fPositionInitY = 0;
53  fPositionCalcY = 0;
54  fPositionErrY = 0;
55  fPositionInitX1 = 0;
56  fPositionCalcX1 = 0;
57  fPositionErrX1 = 0;
58  fPositionInitY1 = 0;
59  fPositionCalcY1 = 0;
60  fPositionErrY1 = 0;
61  fAmpInit = 0;
62  fAmpCalc = 0;
63  fAmpErr = 0;
64  fAmpInitX1 = 0;
65  fAmpCalcX1 = 0;
66  fAmpErrX1 = 0;
67  fAmpInitY1 = 0;
68  fAmpCalcY1 = 0;
69  fAmpErrY1 = 0;
70  fVolume = 0;
71  fVolumeErr = 0;
72  fSigmaInitX = 2;
73  fSigmaCalcX = 0;
74  fSigmaErrX = 0;
75  fSigmaInitY = 2;
76  fSigmaCalcY = 0;
77  fSigmaErrY = 0;
78  fRoInit = 0;
79  fRoCalc = 0;
80  fRoErr = 0;
81  fTxyInit = 0;
82  fTxyCalc = 0;
83  fTxyErr = 0;
84  fTxInit = 0;
85  fTxCalc = 0;
86  fTxErr = 0;
87  fTyInit = 0;
88  fTyCalc = 0;
89  fTyErr = 0;
90  fBxInit = 1;
91  fBxCalc = 0;
92  fBxErr = 0;
93  fByInit = 1;
94  fByCalc = 0;
95  fByErr = 0;
96  fSxyInit = 0;
97  fSxyCalc = 0;
98  fSxyErr = 0;
99  fSxInit = 0;
100  fSxCalc = 0;
101  fSxErr = 0;
102  fSyInit = 0;
103  fSyCalc = 0;
104  fSyErr = 0;
105  fA0Init = 0;
106  fA0Calc = 0;
107  fA0Err = 0;
108  fAxInit = 0;
109  fAxCalc = 0;
110  fAxErr = 0;
111  fAyInit = 0;
112  fAyCalc = 0;
113  fAyErr = 0;
114  fFixPositionX = 0;
115  fFixPositionY = 0;
116  fFixPositionX1 = 0;
117  fFixPositionY1 = 0;
118  fFixAmp = 0;
119  fFixAmpX1 = 0;
120  fFixAmpY1 = 0;
121  fFixSigmaX = false;
122  fFixSigmaY = false;
123  fFixRo = true;
124  fFixTxy = true;
125  fFixTx = true;
126  fFixTy = true;
127  fFixBx = true;
128  fFixBy = true;
129  fFixSxy = true;
130  fFixSx = true;
131  fFixSy = true;
132  fFixA0 = true;
133  fFixAx = true;
134  fFixAy = true;
135 
136 }
137 
138 ////////////////////////////////////////////////////////////////////////////////
139 /// numberPeaks: number of fitted peaks (must be greater than zero)
140 /// the constructor allocates arrays for all fitted parameters (peak positions,
141 /// amplitudes etc) and sets the member variables to their default values. One
142 /// can change these variables by member functions (setters) of TSpectrumFit class.
143 ///
144 /// Shape function of the fitted
145 /// peaks contains the two-dimensional symmetrical Gaussian two one-dimensional
146 /// symmetrical Gaussian ridges as well as non-symmetrical terms and background.
147 ///
148 /// \image html spectrum2fit_constructor_image001.gif
149 
150 TSpectrum2Fit::TSpectrum2Fit(Int_t numberPeaks) :TNamed("Spectrum2Fit", "Miroslav Morhac peak fitter")
151 {
152  if (numberPeaks <= 0){
153  Error ("TSpectrum2Fit","Invalid number of peaks, must be > than 0");
154  return;
155  }
156  fNPeaks = numberPeaks;
157  fNumberIterations = 1;
158  fXmin = 0;
159  fXmax = 100;
160  fYmin = 0;
161  fYmax = 100;
162  fStatisticType = kFitOptimChiCounts;
163  fAlphaOptim = kFitAlphaHalving;
164  fPower = kFitPower2;
165  fFitTaylor = kFitTaylorOrderFirst;
166  fAlpha = 1;
167  fChi = 0;
168  fPositionInitX = new Double_t[numberPeaks];
169  fPositionCalcX = new Double_t[numberPeaks];
170  fPositionErrX = new Double_t[numberPeaks];
171  fPositionInitY = new Double_t[numberPeaks];
172  fPositionCalcY = new Double_t[numberPeaks];
173  fPositionErrY = new Double_t[numberPeaks];
174  fPositionInitX1 = new Double_t[numberPeaks];
175  fPositionCalcX1 = new Double_t[numberPeaks];
176  fPositionErrX1 = new Double_t[numberPeaks];
177  fPositionInitY1 = new Double_t[numberPeaks];
178  fPositionCalcY1 = new Double_t[numberPeaks];
179  fPositionErrY1 = new Double_t[numberPeaks];
180  fAmpInit = new Double_t[numberPeaks];
181  fAmpCalc = new Double_t[numberPeaks];
182  fAmpErr = new Double_t[numberPeaks];
183  fAmpInitX1 = new Double_t[numberPeaks];
184  fAmpCalcX1 = new Double_t[numberPeaks];
185  fAmpErrX1 = new Double_t[numberPeaks];
186  fAmpInitY1 = new Double_t[numberPeaks];
187  fAmpCalcY1 = new Double_t[numberPeaks];
188  fAmpErrY1 = new Double_t[numberPeaks];
189  fVolume = new Double_t[numberPeaks];
190  fVolumeErr = new Double_t[numberPeaks];
191  fSigmaInitX = 2;
192  fSigmaCalcX = 0;
193  fSigmaErrX = 0;
194  fSigmaInitY = 2;
195  fSigmaCalcY = 0;
196  fSigmaErrY = 0;
197  fRoInit = 0;
198  fRoCalc = 0;
199  fRoErr = 0;
200  fTxyInit = 0;
201  fTxyCalc = 0;
202  fTxyErr = 0;
203  fTxInit = 0;
204  fTxCalc = 0;
205  fTxErr = 0;
206  fTyInit = 0;
207  fTyCalc = 0;
208  fTyErr = 0;
209  fBxInit = 1;
210  fBxCalc = 0;
211  fBxErr = 0;
212  fByInit = 1;
213  fByCalc = 0;
214  fByErr = 0;
215  fSxyInit = 0;
216  fSxyCalc = 0;
217  fSxyErr = 0;
218  fSxInit = 0;
219  fSxCalc = 0;
220  fSxErr = 0;
221  fSyInit = 0;
222  fSyCalc = 0;
223  fSyErr = 0;
224  fA0Init = 0;
225  fA0Calc = 0;
226  fA0Err = 0;
227  fAxInit = 0;
228  fAxCalc = 0;
229  fAxErr = 0;
230  fAyInit = 0;
231  fAyCalc = 0;
232  fAyErr = 0;
233  fFixPositionX = new Bool_t[numberPeaks];
234  fFixPositionY = new Bool_t[numberPeaks];
235  fFixPositionX1 = new Bool_t[numberPeaks];
236  fFixPositionY1 = new Bool_t[numberPeaks];
237  fFixAmp = new Bool_t[numberPeaks];
238  fFixAmpX1 = new Bool_t[numberPeaks];
239  fFixAmpY1 = new Bool_t[numberPeaks];
240  fFixSigmaX = false;
241  fFixSigmaY = false;
242  fFixRo = true;
243  fFixTxy = true;
244  fFixTx = true;
245  fFixTy = true;
246  fFixBx = true;
247  fFixBy = true;
248  fFixSxy = true;
249  fFixSx = true;
250  fFixSy = true;
251  fFixA0 = true;
252  fFixAx = true;
253  fFixAy = true;
254 }
255 
256 ////////////////////////////////////////////////////////////////////////////////
257 /// Destructor
258 
259 TSpectrum2Fit::~TSpectrum2Fit()
260 {
261  delete [] fPositionInitX;
262  delete [] fPositionCalcX;
263  delete [] fPositionErrX;
264  delete [] fFixPositionX;
265  delete [] fPositionInitY;
266  delete [] fPositionCalcY;
267  delete [] fPositionErrY;
268  delete [] fFixPositionY;
269  delete [] fPositionInitX1;
270  delete [] fPositionCalcX1;
271  delete [] fPositionErrX1;
272  delete [] fFixPositionX1;
273  delete [] fPositionInitY1;
274  delete [] fPositionCalcY1;
275  delete [] fPositionErrY1;
276  delete [] fFixPositionY1;
277  delete [] fAmpInit;
278  delete [] fAmpCalc;
279  delete [] fAmpErr;
280  delete [] fFixAmp;
281  delete [] fAmpInitX1;
282  delete [] fAmpCalcX1;
283  delete [] fAmpErrX1;
284  delete [] fFixAmpX1;
285  delete [] fAmpInitY1;
286  delete [] fAmpCalcY1;
287  delete [] fAmpErrY1;
288  delete [] fFixAmpY1;
289  delete [] fVolume;
290  delete [] fVolumeErr;
291 }
292 
293 
294 ////////////////////////////////////////////////////////////////////////////////
295 /// This function calculates error function of x.
296 
297 
298 Double_t TSpectrum2Fit::Erfc(Double_t x)
299 {
300  Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap =
301  0.47047;
302  Double_t a, t, c, w;
303  a = TMath::Abs(x);
304  w = 1. + dap * a;
305  t = 1. / w;
306  w = a * a;
307  if (w < 700)
308  c = exp(-w);
309 
310  else {
311  c = 0;
312  }
313  c = c * t * (da1 + t * (da2 + t * da3));
314  if (x < 0)
315  c = 1. - c;
316  return (c);
317 }
318 
319 ////////////////////////////////////////////////////////////////////////////////
320 /// This function calculates derivative of error function of x.
321 
322 Double_t TSpectrum2Fit::Derfc(Double_t x)
323 {
324  Double_t a, t, c, w;
325  Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap =
326  0.47047;
327  a = TMath::Abs(x);
328  w = 1. + dap * a;
329  t = 1. / w;
330  w = a * a;
331  if (w < 700)
332  c = exp(-w);
333 
334  else {
335  c = 0;
336  }
337  c = (-1.) * dap * c * t * t * (da1 + t * (2. * da2 + t * 3. * da3)) -
338  2. * a * Erfc(a);
339  return (c);
340 }
341 
342 ////////////////////////////////////////////////////////////////////////////////
343 /// power function
344 
345 Double_t TSpectrum2Fit::Ourpowl(Double_t a, Int_t pw)
346 {
347  Double_t c;
348  Double_t a2 = a*a;
349  c = 1;
350  if (pw > 0) c *= a2;
351  if (pw > 2) c *= a2;
352  if (pw > 4) c *= a2;
353  if (pw > 6) c *= a2;
354  if (pw > 8) c *= a2;
355  if (pw > 10) c *= a2;
356  if (pw > 12) c *= a2;
357  return (c);
358 }
359 
360 ////////////////////////////////////////////////////////////////////////////////
361 /// This function calculates solution of the system of linear equations.
362 /// The matrix a should have a dimension size*(size+4)
363 /// The calling function should fill in the matrix, the column size should
364 /// contain vector y (right side of the system of equations). The result is
365 /// placed into size+1 column of the matrix.
366 /// according to sigma of peaks.
367 ///
368 /// Function parameters:
369 /// - a-matrix with dimension size*(size+4)
370 /// - size-number of rows of the matrix
371 
372 void TSpectrum2Fit::StiefelInversion(Double_t **a, Int_t size)
373 {
374  Int_t i, j, k = 0;
375  Double_t sk = 0, b, lambdak, normk, normk_old = 0;
376 
377  do {
378  normk = 0;
379 
380  //calculation of rk and norm
381  for (i = 0; i < size; i++) {
382  a[i][size + 2] = -a[i][size]; //rk=-C
383  for (j = 0; j < size; j++) {
384  a[i][size + 2] += a[i][j] * a[j][size + 1]; //A*xk-C
385  }
386  normk += a[i][size + 2] * a[i][size + 2]; //calculation normk
387  }
388 
389  //calculation of sk
390  if (k != 0) {
391  sk = normk / normk_old;
392  }
393 
394  //calculation of uk
395  for (i = 0; i < size; i++) {
396  a[i][size + 3] = -a[i][size + 2] + sk * a[i][size + 3]; //uk=-rk+sk*uk-1
397  }
398 
399  //calculation of lambdak
400  lambdak = 0;
401  for (i = 0; i < size; i++) {
402  for (j = 0, b = 0; j < size; j++) {
403  b += a[i][j] * a[j][size + 3]; //A*uk
404  }
405  lambdak += b * a[i][size + 3];
406  }
407  if (TMath::Abs(lambdak) > 1e-50) //computer zero
408  lambdak = normk / lambdak;
409 
410  else
411  lambdak = 0;
412  for (i = 0; i < size; i++)
413  a[i][size + 1] += lambdak * a[i][size + 3]; //xk+1=xk+lambdak*uk
414  normk_old = normk;
415  k += 1;
416  } while (k < size && TMath::Abs(normk) > 1e-50); //computer zero
417  return;
418 }
419 
420 ////////////////////////////////////////////////////////////////////////////////
421 /// This function calculates 2D peaks shape function (see manual)
422 ///
423 /// Function parameters:
424 /// - numOfFittedPeaks-number of fitted peaks
425 /// - x-channel in x-dimension
426 /// - y-channel in y-dimension
427 /// - parameter-array of peaks parameters (amplitudes and positions)
428 /// - sigmax-sigmax of peaks
429 /// - sigmay-sigmay of peaks
430 /// - ro-correlation coefficient
431 /// - a0,ax,ay-bac kground coefficients
432 /// - txy,tx,ty, sxy,sy,sx-relative amplitudes
433 /// - bx, by-slopes
434 
435 Double_t TSpectrum2Fit::Shape2(Int_t numOfFittedPeaks, Double_t x, Double_t y,
436  const Double_t *parameter, Double_t sigmax,
437  Double_t sigmay, Double_t ro, Double_t a0, Double_t ax,
438  Double_t ay, Double_t txy, Double_t sxy, Double_t tx,
439  Double_t ty, Double_t sx, Double_t sy, Double_t bx,
440  Double_t by)
441 {
442  Int_t j;
443  Double_t r, p, r1, e, ex, ey, vx, s2, px, py, rx, ry, erx, ery;
444  vx = 0;
445  s2 = TMath::Sqrt(2.0);
446  for (j = 0; j < numOfFittedPeaks; j++) {
447  p = (x - parameter[7 * j + 1]) / sigmax;
448  r = (y - parameter[7 * j + 2]) / sigmay;
449  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
450  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
451  if (e < 700)
452  r1 = exp(-e);
453 
454  else {
455  r1 = 0;
456  }
457  if (txy != 0) {
458  px = 0, py = 0;
459  erx = Erfc(p / s2 + 1 / (2 * bx)), ery =
460  Erfc(r / s2 + 1 / (2 * by));
461  ex = p / (s2 * bx), ey = r / (s2 * by);
462  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
463  px = exp(ex) * erx, py = exp(ey) * ery;
464  }
465  r1 += 0.5 * txy * px * py;
466  }
467  if (sxy != 0) {
468  rx = Erfc(p / s2), ry = Erfc(r / s2);
469  r1 += 0.5 * sxy * rx * ry;
470  }
471  vx = vx + parameter[7 * j] * r1;
472  }
473  p = (x - parameter[7 * j + 5]) / sigmax;
474  if (TMath::Abs(p) < 3) {
475  e = p * p / 2;
476  if (e < 700)
477  r1 = exp(-e);
478 
479  else {
480  r1 = 0;
481  }
482  if (tx != 0) {
483  px = 0;
484  erx = Erfc(p / s2 + 1 / (2 * bx));
485  ex = p / (s2 * bx);
486  if (TMath::Abs(ex) < 9) {
487  px = exp(ex) * erx;
488  }
489  r1 += 0.5 * tx * px;
490  }
491  if (sx != 0) {
492  rx = Erfc(p / s2);
493  r1 += 0.5 * sx * rx;
494  }
495  vx = vx + parameter[7 * j + 3] * r1;
496  }
497  r = (y - parameter[7 * j + 6]) / sigmay;
498  if (TMath::Abs(r) < 3) {
499  e = r * r / 2;
500  if (e < 700)
501  r1 = exp(-e);
502 
503  else {
504  r1 = 0;
505  }
506  if (ty != 0) {
507  py = 0;
508  ery = Erfc(r / s2 + 1 / (2 * by));
509  ey = r / (s2 * by);
510  if (TMath::Abs(ey) < 9) {
511  py = exp(ey) * ery;
512  }
513  r1 += 0.5 * ty * py;
514  }
515  if (sy != 0) {
516  ry = Erfc(r / s2);
517  r1 += 0.5 * sy * ry;
518  }
519  vx = vx + parameter[7 * j + 4] * r1;
520  }
521  }
522  vx = vx + a0 + ax * x + ay * y;
523  return (vx);
524 }
525 
526 ////////////////////////////////////////////////////////////////////////////////
527 /// This function calculates derivative of 2D peaks shape function (see manual)
528 /// according to amplitude of 2D peak
529 ///
530 /// Function parameters:
531 /// - x-channel in x-dimension
532 /// - y-channel in y-dimension
533 /// - x0-position of peak in x-dimension
534 /// - y0-position of peak in y-dimension
535 /// - sigmax-sigmax of peaks
536 /// - sigmay-sigmay of peaks
537 /// - ro-correlation coefficient
538 /// - txy, sxy-relative amplitudes
539 /// - bx, by-slopes
540 
541 Double_t TSpectrum2Fit::Deramp2(Double_t x, Double_t y, Double_t x0, Double_t y0,
542  Double_t sigmax, Double_t sigmay, Double_t ro,
543  Double_t txy, Double_t sxy, Double_t bx, Double_t by)
544 {
545  Double_t p, r, r1 = 0, e, ex, ey, px, py, rx, ry, erx, ery, s2;
546  p = (x - x0) / sigmax;
547  r = (y - y0) / sigmay;
548  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
549  s2 = TMath::Sqrt(2.0);
550  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
551  if (e < 700)
552  r1 = exp(-e);
553 
554  else {
555  r1 = 0;
556  }
557  if (txy != 0) {
558  px = 0, py = 0;
559  erx = Erfc(p / s2 + 1 / (2 * bx)), ery =
560  Erfc(r / s2 + 1 / (2 * by));
561  ex = p / (s2 * bx), ey = r / (s2 * by);
562  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
563  px = exp(ex) * erx, py = exp(ey) * ery;
564  }
565  r1 += 0.5 * txy * px * py;
566  }
567  if (sxy != 0) {
568  rx = Erfc(p / s2), ry = Erfc(r / s2);
569  r1 += 0.5 * sxy * rx * ry;
570  }
571  }
572  return (r1);
573 }
574 
575 ////////////////////////////////////////////////////////////////////////////////
576 /// This function calculates derivative of 2D peaks shape function (see manual)
577 /// according to amplitude of the ridge
578 ///
579 /// Function parameters:
580 /// - x-channel in x-dimension
581 /// - x0-position of peak in x-dimension
582 /// - y0-position of peak in y-dimension
583 /// - sigmax-sigmax of peaks
584 /// - ro-correlation coefficient
585 /// - tx, sx-relative amplitudes
586 /// - bx-slope
587 
588 Double_t TSpectrum2Fit::Derampx(Double_t x, Double_t x0, Double_t sigmax, Double_t tx,
589  Double_t sx, Double_t bx)
590 {
591  Double_t p, r1 = 0, px, erx, rx, ex, s2;
592  p = (x - x0) / sigmax;
593  if (TMath::Abs(p) < 3) {
594  s2 = TMath::Sqrt(2.0);
595  p = p * p / 2;
596  if (p < 700)
597  r1 = exp(-p);
598 
599  else {
600  r1 = 0;
601  }
602  if (tx != 0) {
603  px = 0;
604  erx = Erfc(p / s2 + 1 / (2 * bx));
605  ex = p / (s2 * bx);
606  if (TMath::Abs(ex) < 9) {
607  px = exp(ex) * erx;
608  }
609  r1 += 0.5 * tx * px;
610  }
611  if (sx != 0) {
612  rx = Erfc(p / s2);
613  r1 += 0.5 * sx * rx;
614  }
615  }
616  return (r1);
617 }
618 
619 ////////////////////////////////////////////////////////////////////////////////
620 /// This function calculates derivative of 2D peaks shape function (see manual)
621 /// according to x position of 2D peak
622 ///
623 /// Function parameters:
624 /// - x-channel in x-dimension
625 /// - y-channel in y-dimension
626 /// - a-amplitude
627 /// - x0-position of peak in x-dimension
628 /// - y0-position of peak in y-dimension
629 /// - sigmax-sigmax of peaks
630 /// - sigmay-sigmay of peaks
631 /// - ro-correlation coefficient
632 /// - txy, sxy-relative amplitudes
633 /// - bx, by-slopes
634 
635 Double_t TSpectrum2Fit::Deri02(Double_t x, Double_t y, Double_t a, Double_t x0,
636  Double_t y0, Double_t sigmax, Double_t sigmay,
637  Double_t ro, Double_t txy, Double_t sxy, Double_t bx,
638  Double_t by)
639 {
640  Double_t p, r, r1 = 0, e, ex, ey, px, py, rx, ry, erx, ery, s2;
641  p = (x - x0) / sigmax;
642  r = (y - y0) / sigmay;
643  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
644  s2 = TMath::Sqrt(2.0);
645  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
646  if (e < 700)
647  r1 = exp(-e);
648 
649  else {
650  r1 = 0;
651  }
652  e = -(ro * r - p) / sigmax;
653  e = e / (1 - ro * ro);
654  r1 = r1 * e;
655  if (txy != 0) {
656  px = 0, py = 0;
657  erx =
658  (-Erfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * sigmax) -
659  Derfc(p / s2 + 1 / (2 * bx)) / (s2 * sigmax)), ery =
660  Erfc(r / s2 + 1 / (2 * by));
661  ex = p / (s2 * bx), ey = r / (s2 * by);
662  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
663  px = exp(ex) * erx, py = exp(ey) * ery;
664  }
665  r1 += 0.5 * txy * px * py;
666  }
667  if (sxy != 0) {
668  rx = -Derfc(p / s2) / (s2 * sigmax), ry = Erfc(r / s2);
669  r1 += 0.5 * sxy * rx * ry;
670  }
671  r1 = a * r1;
672  }
673  return (r1);
674 }
675 
676 ////////////////////////////////////////////////////////////////////////////////
677 /// This function calculates second derivative of 2D peaks shape function
678 /// (see manual) according to x position of 2D peak
679 ///
680 /// Function parameters:
681 /// - x-channel in x-dimension
682 /// - y-channel in y-dimension
683 /// - a-amplitude
684 /// - x0-position of peak in x-dimension
685 /// - y0-position of peak in y-dimension
686 /// - sigmax-sigmax of peaks
687 /// - sigmay-sigmay of peaks
688 /// - ro-correlation coefficient
689 
690 Double_t TSpectrum2Fit::Derderi02(Double_t x, Double_t y, Double_t a, Double_t x0,
691  Double_t y0, Double_t sigmax, Double_t sigmay,
692  Double_t ro)
693 {
694  Double_t p, r, r1 = 0, e;
695  p = (x - x0) / sigmax;
696  r = (y - y0) / sigmay;
697  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
698  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
699  if (e < 700)
700  r1 = exp(-e);
701 
702  else {
703  r1 = 0;
704  }
705  e = -(ro * r - p) / sigmax;
706  e = e / (1 - ro * ro);
707  r1 = r1 * (e * e - 1 / ((1 - ro * ro) * sigmax * sigmax));
708  r1 = a * r1;
709  }
710  return (r1);
711 }
712 
713 ////////////////////////////////////////////////////////////////////////////////
714 /// This function calculates derivative of 2D peaks shape function (see manual)
715 /// according to y position of 2D peak
716 /// Function parameters:
717 /// - x-channel in x-dimension
718 /// - y-channel in y-dimension
719 /// - a-amplitude
720 /// - x0-position of peak in x-dimension
721 /// - y0-position of peak in y-dimension
722 /// - sigmax-sigmax of peaks
723 /// - sigmay-sigmay of peaks
724 /// - ro-correlation coefficient
725 /// - txy, sxy-relative amplitudes
726 /// - bx, by-slopes
727 
728 Double_t TSpectrum2Fit::Derj02(Double_t x, Double_t y, Double_t a, Double_t x0,
729  Double_t y0, Double_t sigmax, Double_t sigmay,
730  Double_t ro, Double_t txy, Double_t sxy, Double_t bx,
731  Double_t by)
732 {
733 
734 
735  Double_t p, r, r1 = 0, e, ex, ey, px, py, rx, ry, erx, ery, s2;
736  p = (x - x0) / sigmax;
737  r = (y - y0) / sigmay;
738  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
739  s2 = TMath::Sqrt(2.0);
740  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
741  if (e < 700)
742  r1 = exp(-e);
743 
744  else {
745  r1 = 0;
746  }
747  e = -(ro * p - r) / sigmay;
748  e = e / (1 - ro * ro);
749  r1 = r1 * e;
750  if (txy != 0) {
751  px = 0, py = 0;
752  ery =
753  (-Erfc(r / s2 + 1 / (2 * by)) / (s2 * by * sigmay) -
754  Derfc(r / s2 + 1 / (2 * by)) / (s2 * sigmay)), erx =
755  Erfc(p / s2 + 1 / (2 * bx));
756  ex = p / (s2 * bx), ey = r / (s2 * by);
757  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
758  px = exp(ex) * erx, py = exp(ey) * ery;
759  }
760  r1 += 0.5 * txy * px * py;
761  }
762  if (sxy != 0) {
763  ry = -Derfc(r / s2) / (s2 * sigmay), rx = Erfc(p / s2);
764  r1 += 0.5 * sxy * rx * ry;
765  }
766  r1 = a * r1;
767  }
768  return (r1);
769 }
770 
771 ////////////////////////////////////////////////////////////////////////////////
772 /// This function calculates second derivative of 2D peaks shape function
773 /// (see manual) according to y position of 2D peak
774 ///
775 /// Function parameters:
776 /// - x-channel in x-dimension
777 /// - y-channel in y-dimension
778 /// - a-amplitude
779 /// - x0-position of peak in x-dimension
780 /// - y0-position of peak in y-dimension
781 /// - sigmax-sigmax of peaks
782 /// - sigmay-sigmay of peaks
783 /// - ro-correlation coefficient
784 
785 Double_t TSpectrum2Fit::Derderj02(Double_t x, Double_t y, Double_t a, Double_t x0,
786  Double_t y0, Double_t sigmax, Double_t sigmay,
787  Double_t ro)
788 {
789  Double_t p, r, r1 = 0, e;
790  p = (x - x0) / sigmax;
791  r = (y - y0) / sigmay;
792  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
793  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
794  if (e < 700)
795  r1 = exp(-e);
796 
797  else {
798  r1 = 0;
799  }
800  e = -(ro * p - r) / sigmay;
801  e = e / (1 - ro * ro);
802  r1 = r1 * (e * e - 1 / ((1 - ro * ro) * sigmay * sigmay));
803  r1 = a * r1;
804  }
805  return (r1);
806 }
807 
808 ////////////////////////////////////////////////////////////////////////////////
809 /// This function calculates derivative of 2D peaks shape function (see manual)
810 /// according to x position of 1D ridge
811 ///
812 /// Function parameters:
813 /// - x-channel in x-dimension
814 /// - ax-amplitude of ridge
815 /// - x0-position of peak in x-dimension
816 /// - sigmax-sigmax of peaks
817 /// - ro-correlation coefficient
818 /// - tx, sx-relative amplitudes
819 /// - bx-slope
820 
821 Double_t TSpectrum2Fit::Deri01(Double_t x, Double_t ax, Double_t x0, Double_t sigmax,
822  Double_t tx, Double_t sx, Double_t bx)
823 {
824  Double_t p, e, r1 = 0, px, rx, erx, ex, s2;
825  p = (x - x0) / sigmax;
826  if (TMath::Abs(p) < 3) {
827  s2 = TMath::Sqrt(2.0);
828  e = p * p / 2;
829  if (e < 700)
830  r1 = exp(-e);
831 
832  else {
833  r1 = 0;
834  }
835  r1 = r1 * p / sigmax;
836  if (tx != 0) {
837  px = 0;
838  erx =
839  (-Erfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * sigmax) -
840  Derfc(p / s2 + 1 / (2 * bx)) / (s2 * sigmax));
841  ex = p / (s2 * bx);
842  if (TMath::Abs(ex) < 9)
843  px = exp(ex) * erx;
844  r1 += 0.5 * tx * px;
845  }
846  if (sx != 0) {
847  rx = -Derfc(p / s2) / (s2 * sigmax);
848  r1 += 0.5 * sx * rx;
849  }
850  r1 = ax * r1;
851  }
852  return (r1);
853 }
854 
855 ////////////////////////////////////////////////////////////////////////////////
856 /// This function calculates second derivative of 2D peaks shape function
857 /// (see manual) according to x position of 1D ridge
858 ///
859 /// Function parameters:
860 /// - x-channel in x-dimension
861 /// - ax-amplitude of ridge
862 /// - x0-position of peak in x-dimension
863 /// - sigmax-sigmax of peaks
864 
865 Double_t TSpectrum2Fit::Derderi01(Double_t x, Double_t ax, Double_t x0,
866  Double_t sigmax)
867 {
868  Double_t p, e, r1 = 0;
869  p = (x - x0) / sigmax;
870  if (TMath::Abs(p) < 3) {
871  e = p * p / 2;
872  if (e < 700)
873  r1 = exp(-e);
874 
875  else {
876  r1 = 0;
877  }
878  r1 = r1 * (p * p / (sigmax * sigmax) - 1 / (sigmax * sigmax));
879  r1 = ax * r1;
880  }
881  return (r1);
882 }
883 
884 ////////////////////////////////////////////////////////////////////////////////
885 /// This function calculates derivative of peaks shape function (see manual)
886 /// according to sigmax of peaks.
887 ///
888 /// Function parameters:
889 /// - numOfFittedPeaks-number of fitted peaks
890 /// - x,y-position of channel
891 /// - parameter-array of peaks parameters (amplitudes and positions)
892 /// - sigmax-sigmax of peaks
893 /// - sigmay-sigmay of peaks
894 /// - ro-correlation coefficient
895 /// - txy, sxy, tx, sx-relative amplitudes
896 /// - bx, by-slopes
897 
898 Double_t TSpectrum2Fit::Dersigmax(Int_t numOfFittedPeaks, Double_t x, Double_t y,
899  const Double_t *parameter, Double_t sigmax,
900  Double_t sigmay, Double_t ro, Double_t txy,
901  Double_t sxy, Double_t tx, Double_t sx, Double_t bx,
902  Double_t by)
903 {
904  Double_t p, r, r1 =
905  0, e, a, b, x0, y0, s2, px, py, rx, ry, erx, ery, ex, ey;
906  Int_t j;
907  s2 = TMath::Sqrt(2.0);
908  for (j = 0; j < numOfFittedPeaks; j++) {
909  a = parameter[7 * j];
910  x0 = parameter[7 * j + 1];
911  y0 = parameter[7 * j + 2];
912  p = (x - x0) / sigmax;
913  r = (y - y0) / sigmay;
914  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
915  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
916  if (e < 700)
917  e = exp(-e);
918 
919  else {
920  e = 0;
921  }
922  b = -(ro * p * r - p * p) / sigmax;
923  e = e * b / (1 - ro * ro);
924  if (txy != 0) {
925  px = 0, py = 0;
926  erx =
927  -Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * sigmax) -
928  Derfc(p / s2 + 1 / (2 * bx)) * p / (s2 * sigmax), ery =
929  Erfc(r / s2 + 1 / (2 * by));
930  ex = p / (s2 * bx), ey = r / (s2 * by);
931  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
932  px = exp(ex) * erx, py = exp(ey) * ery;
933  }
934  e += 0.5 * txy * px * py;
935  }
936  if (sxy != 0) {
937  rx = -Derfc(p / s2) * p / (s2 * sigmax), ry = Erfc(r / s2);
938  e += 0.5 * sxy * rx * ry;
939  }
940  r1 = r1 + a * e;
941  }
942  if (TMath::Abs(p) < 3) {
943  x0 = parameter[7 * j + 5];
944  p = (x - x0) / sigmax;
945  b = p * p / 2;
946  if (b < 700)
947  e = exp(-b);
948 
949  else {
950  e = 0;
951  }
952  e = 2 * b * e / sigmax;
953  if (tx != 0) {
954  px = 0;
955  erx =
956  (-Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * sigmax) -
957  Derfc(p / s2 + 1 / (2 * bx)) * p / (s2 * sigmax));
958  ex = p / (s2 * bx);
959  if (TMath::Abs(ex) < 9)
960  px = exp(ex) * erx;
961  e += 0.5 * tx * px;
962  }
963  if (sx != 0) {
964  rx = -Derfc(p / s2) * p / (s2 * sigmax);
965  e += 0.5 * sx * rx;
966  }
967  r1 += parameter[7 * j + 3] * e;
968  }
969  }
970  return (r1);
971 }
972 
973 ////////////////////////////////////////////////////////////////////////////////
974 /// This function calculates second derivative of peaks shape function
975 /// (see manual) according to sigmax of peaks.
976 ///
977 /// Function parameters:
978 /// - numOfFittedPeaks-number of fitted peaks
979 /// - x,y-position of channel
980 /// - parameter-array of peaks parameters (amplitudes and positions)
981 /// - sigmax-sigmax of peaks
982 /// - sigmay-sigmay of peaks
983 /// - ro-correlation coefficient
984 
985 Double_t TSpectrum2Fit::Derdersigmax(Int_t numOfFittedPeaks, Double_t x,
986  Double_t y, const Double_t *parameter,
987  Double_t sigmax, Double_t sigmay,
988  Double_t ro)
989 {
990  Double_t p, r, r1 = 0, e, a, b, x0, y0;
991  Int_t j;
992  for (j = 0; j < numOfFittedPeaks; j++) {
993  a = parameter[7 * j];
994  x0 = parameter[7 * j + 1];
995  y0 = parameter[7 * j + 2];
996  p = (x - x0) / sigmax;
997  r = (y - y0) / sigmay;
998  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
999  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
1000  if (e < 700)
1001  e = exp(-e);
1002 
1003  else {
1004  e = 0;
1005  }
1006  b = -(ro * p * r - p * p) / sigmax;
1007  e = e * (b * b / (1 - ro * ro) -
1008  (3 * p * p - 2 * ro * p * r) / (sigmax * sigmax)) / (1 -
1009  ro
1010  *
1011  ro);
1012  r1 = r1 + a * e;
1013  }
1014  if (TMath::Abs(p) < 3) {
1015  x0 = parameter[7 * j + 5];
1016  p = (x - x0) / sigmax;
1017  b = p * p / 2;
1018  if (b < 700)
1019  e = exp(-b);
1020 
1021  else {
1022  e = 0;
1023  }
1024  e = e * (4 * b * b - 6 * b) / (sigmax * sigmax);
1025  r1 += parameter[7 * j + 3] * e;
1026  }
1027  }
1028  return (r1);
1029 }
1030 
1031 ////////////////////////////////////////////////////////////////////////////////
1032 /// This function calculates derivative of peaks shape function (see manual)
1033 /// according to sigmax of peaks.
1034 ///
1035 /// Function parameters:
1036 /// - numOfFittedPeaks-number of fitted peaks
1037 /// - x,y-position of channel
1038 /// - parameter-array of peaks parameters (amplitudes and positions)
1039 /// - sigmax-sigmax of peaks
1040 /// - sigmay-sigmay of peaks
1041 /// - ro-correlation coefficient
1042 /// - txy, sxy, ty, sy-relative amplitudes
1043 /// - bx, by-slopes
1044 
1045 Double_t TSpectrum2Fit::Dersigmay(Int_t numOfFittedPeaks, Double_t x, Double_t y,
1046  const Double_t *parameter, Double_t sigmax,
1047  Double_t sigmay, Double_t ro, Double_t txy,
1048  Double_t sxy, Double_t ty, Double_t sy, Double_t bx,
1049  Double_t by)
1050 {
1051  Double_t p, r, r1 =
1052  0, e, a, b, x0, y0, s2, px, py, rx, ry, erx, ery, ex, ey;
1053  Int_t j;
1054  s2 = TMath::Sqrt(2.0);
1055  for (j = 0; j < numOfFittedPeaks; j++) {
1056  a = parameter[7 * j];
1057  x0 = parameter[7 * j + 1];
1058  y0 = parameter[7 * j + 2];
1059  p = (x - x0) / sigmax;
1060  r = (y - y0) / sigmay;
1061  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
1062  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
1063  if (e < 700)
1064  e = exp(-e);
1065 
1066  else {
1067  e = 0;
1068  }
1069  b = -(ro * p * r - r * r) / sigmay;
1070  e = e * b / (1 - ro * ro);
1071  if (txy != 0) {
1072  px = 0, py = 0;
1073  ery =
1074  -Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * sigmay) -
1075  Derfc(r / s2 + 1 / (2 * by)) * r / (s2 * sigmay), erx =
1076  Erfc(p / s2 + 1 / (2 * bx));
1077  ex = p / (s2 * bx), ey = r / (s2 * by);
1078  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
1079  px = exp(ex) * erx, py = exp(ey) * ery;
1080  }
1081  e += 0.5 * txy * px * py;
1082  }
1083  if (sxy != 0) {
1084  ry = -Derfc(r / s2) * r / (s2 * sigmay), rx = Erfc(p / s2);
1085  e += 0.5 * sxy * rx * ry;
1086  }
1087  r1 = r1 + a * e;
1088  }
1089  if (TMath::Abs(r) < 3) {
1090  y0 = parameter[7 * j + 6];
1091  r = (y - y0) / sigmay;
1092  b = r * r / 2;
1093  if (b < 700)
1094  e = exp(-b);
1095 
1096  else {
1097  e = 0;
1098  }
1099  e = 2 * b * e / sigmay;
1100  if (ty != 0) {
1101  py = 0;
1102  ery =
1103  (-Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * sigmay) -
1104  Derfc(r / s2 + 1 / (2 * by)) * r / (s2 * sigmay));
1105  ey = r / (s2 * by);
1106  if (TMath::Abs(ey) < 9)
1107  py = exp(ey) * ery;
1108  e += 0.5 * ty * py;
1109  }
1110  if (sy != 0) {
1111  ry = -Derfc(r / s2) * r / (s2 * sigmay);
1112  e += 0.5 * sy * ry;
1113  }
1114  r1 += parameter[7 * j + 4] * e;
1115  }
1116  }
1117  return (r1);
1118 }
1119 
1120 ////////////////////////////////////////////////////////////////////////////////
1121 /// This function calculates second derivative of peaks shape function
1122 /// (see manual) according to sigmay of peaks.
1123 ///
1124 /// Function parameters:
1125 /// - numOfFittedPeaks-number of fitted peaks
1126 /// - x,y-position of channel
1127 /// - parameter-array of peaks parameters (amplitudes and positions)
1128 /// - sigmax-sigmax of peaks
1129 /// - sigmay-sigmay of peaks
1130 /// - ro-correlation coefficient
1131 
1132 Double_t TSpectrum2Fit::Derdersigmay(Int_t numOfFittedPeaks, Double_t x,
1133  Double_t y, const Double_t *parameter,
1134  Double_t sigmax, Double_t sigmay,
1135  Double_t ro)
1136 {
1137  Double_t p, r, r1 = 0, e, a, b, x0, y0;
1138  Int_t j;
1139  for (j = 0; j < numOfFittedPeaks; j++) {
1140  a = parameter[7 * j];
1141  x0 = parameter[7 * j + 1];
1142  y0 = parameter[7 * j + 2];
1143  p = (x - x0) / sigmax;
1144  r = (y - y0) / sigmay;
1145  if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
1146  e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
1147  if (e < 700)
1148  e = exp(-e);
1149 
1150  else {
1151  e = 0;
1152  }
1153  b = -(ro * p * r - r * r) / sigmay;
1154  e = e * (b * b / (1 - ro * ro) -
1155  (3 * r * r - 2 * ro * r * p) / (sigmay * sigmay)) / (1 -
1156  ro
1157  *
1158  ro);
1159  r1 = r1 + a * e;
1160  }
1161  if (TMath::Abs(r) < 3) {
1162  y0 = parameter[7 * j + 6];
1163  r = (y - y0) / sigmay;
1164  b = r * r / 2;
1165  if (b < 700)
1166  e = exp(-b);
1167 
1168  else {
1169  e = 0;
1170  }
1171  e = e * (4 * b * b - 6 * b) / (sigmay * sigmay);
1172  r1 += parameter[7 * j + 4] * e;
1173  }
1174  }
1175  return (r1);
1176 }
1177 
1178 ////////////////////////////////////////////////////////////////////////////////
1179 /// This function calculates derivative of peaks shape function (see manual)
1180 /// according to correlation coefficient ro.
1181 ///
1182 /// Function parameters:
1183 /// - numOfFittedPeaks-number of fitted peaks
1184 /// - x,y-position of channel
1185 /// - parameter-array of peaks parameters (amplitudes and positions)
1186 /// - sx-sigmax of peaks
1187 /// - sy-sigmay of peaks
1188 /// - r-correlation coefficient ro
1189 
1190 Double_t TSpectrum2Fit::Derro(Int_t numOfFittedPeaks, Double_t x, Double_t y,
1191  const Double_t *parameter, Double_t sx, Double_t sy,
1192  Double_t r)
1193 {
1194  Double_t px, qx, rx, vx, x0, y0, a, ex, tx;
1195  Int_t j;
1196  vx = 0;
1197  for (j = 0; j < numOfFittedPeaks; j++) {
1198  a = parameter[7 * j];
1199  x0 = parameter[7 * j + 1];
1200  y0 = parameter[7 * j + 2];
1201  px = (x - x0) / sx;
1202  qx = (y - y0) / sy;
1203  if (TMath::Abs(px) < 3 && TMath::Abs(qx) < 3) {
1204  rx = (px * px - 2 * r * px * qx + qx * qx);
1205  ex = rx / (2 * (1 - r * r));
1206  if ((ex) < 700)
1207  ex = exp(-ex);
1208 
1209  else {
1210  ex = 0;
1211  }
1212  tx = px * qx / (1 - r * r);
1213  tx = tx - r * rx / ((1 - r * r) * (1 - r * r));
1214  vx = vx + a * ex * tx;
1215  }
1216  }
1217  return (vx);
1218 }
1219 
1220 ////////////////////////////////////////////////////////////////////////////////
1221 /// This function calculates derivative of peaks shape function (see manual)
1222 /// according to relative amplitude txy.
1223 ///
1224 /// Function parameters:
1225 /// - numOfFittedPeaks-number of fitted peaks
1226 /// - x,y-position of channel
1227 /// - parameter-array of peaks parameters (amplitudes and positions)
1228 /// - sigmax-sigmax of peaks
1229 /// - sigmay-sigmay of peaks
1230 /// - bx, by-slopes
1231 
1232 Double_t TSpectrum2Fit::Dertxy(Int_t numOfFittedPeaks, Double_t x, Double_t y,
1233  const Double_t *parameter, Double_t sigmax,
1234  Double_t sigmay, Double_t bx, Double_t by)
1235 {
1236  Double_t p, r, r1 = 0, ex, ey, px, py, erx, ery, s2, x0, y0, a;
1237  Int_t j;
1238  s2 = TMath::Sqrt(2.0);
1239  for (j = 0; j < numOfFittedPeaks; j++) {
1240  a = parameter[7 * j];
1241  x0 = parameter[7 * j + 1];
1242  y0 = parameter[7 * j + 2];
1243  p = (x - x0) / sigmax;
1244  r = (y - y0) / sigmay;
1245  px = 0, py = 0;
1246  erx = Erfc(p / s2 + 1 / (2 * bx)), ery =
1247  Erfc(r / s2 + 1 / (2 * by));
1248  ex = p / (s2 * bx), ey = r / (s2 * by);
1249  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
1250  px = exp(ex) * erx, py = exp(ey) * ery;
1251  }
1252  r1 += 0.5 * a * px * py;
1253  }
1254  return (r1);
1255 }
1256 
1257 ////////////////////////////////////////////////////////////////////////////////
1258 /// This function calculates derivative of peaks shape function (see manual)
1259 /// according to relative amplitude sxy.
1260 ///
1261 /// Function parameters:
1262 /// - numOfFittedPeaks-number of fitted peaks
1263 /// - x,y-position of channel
1264 /// - parameter-array of peaks parameters (amplitudes and positions)
1265 /// - sigmax-sigmax of peaks
1266 /// - sigmay-sigmay of peaks
1267 
1268 Double_t TSpectrum2Fit::Dersxy(Int_t numOfFittedPeaks, Double_t x, Double_t y,
1269  const Double_t *parameter, Double_t sigmax,
1270  Double_t sigmay)
1271 {
1272  Double_t p, r, r1 = 0, rx, ry, x0, y0, a, s2;
1273  Int_t j;
1274  s2 = TMath::Sqrt(2.0);
1275  for (j = 0; j < numOfFittedPeaks; j++) {
1276  a = parameter[7 * j];
1277  x0 = parameter[7 * j + 1];
1278  y0 = parameter[7 * j + 2];
1279  p = (x - x0) / sigmax;
1280  r = (y - y0) / sigmay;
1281  rx = Erfc(p / s2), ry = Erfc(r / s2);
1282  r1 += 0.5 * a * rx * ry;
1283  }
1284  return (r1);
1285 }
1286 
1287 ////////////////////////////////////////////////////////////////////////////////
1288 /// This function calculates derivative of peaks shape function (see manual)
1289 /// according to relative amplitude tx.
1290 ///
1291 /// Function parameters:
1292 /// - numOfFittedPeaks-number of fitted peaks
1293 /// - x-position of channel
1294 /// - parameter-array of peaks parameters (amplitudes and positions)
1295 /// - sigmax-sigma of 1D ridge
1296 /// - bx-slope
1297 
1298 Double_t TSpectrum2Fit::Dertx(Int_t numOfFittedPeaks, Double_t x,
1299  const Double_t *parameter, Double_t sigmax,
1300  Double_t bx)
1301 {
1302  Double_t p, r1 = 0, ex, px, erx, s2, ax, x0;
1303  Int_t j;
1304  s2 = TMath::Sqrt(2.0);
1305  for (j = 0; j < numOfFittedPeaks; j++) {
1306  ax = parameter[7 * j + 3];
1307  x0 = parameter[7 * j + 5];
1308  p = (x - x0) / sigmax;
1309  px = 0;
1310  erx = Erfc(p / s2 + 1 / (2 * bx));
1311  ex = p / (s2 * bx);
1312  if (TMath::Abs(ex) < 9) {
1313  px = exp(ex) * erx;
1314  }
1315  r1 += 0.5 * ax * px;
1316  }
1317  return (r1);
1318 }
1319 
1320 ////////////////////////////////////////////////////////////////////////////////
1321 /// This function calculates derivative of peaks shape function (see manual)
1322 /// according to relative amplitude ty.
1323 ///
1324 /// Function parameters:
1325 /// - numOfFittedPeaks-number of fitted peaks
1326 /// - x-position of channel
1327 /// - parameter-array of peaks parameters (amplitudes and positions)
1328 /// - sigmax-sigma of 1D ridge
1329 /// - bx-slope
1330 
1331 Double_t TSpectrum2Fit::Derty(Int_t numOfFittedPeaks, Double_t x,
1332  const Double_t *parameter, Double_t sigmax,
1333  Double_t bx)
1334 {
1335  Double_t p, r1 = 0, ex, px, erx, s2, ax, x0;
1336  Int_t j;
1337  s2 = TMath::Sqrt(2.0);
1338  for (j = 0; j < numOfFittedPeaks; j++) {
1339  ax = parameter[7 * j + 4];
1340  x0 = parameter[7 * j + 6];
1341  p = (x - x0) / sigmax;
1342  px = 0;
1343  erx = Erfc(p / s2 + 1 / (2 * bx));
1344  ex = p / (s2 * bx);
1345  if (TMath::Abs(ex) < 9) {
1346  px = exp(ex) * erx;
1347  }
1348  r1 += 0.5 * ax * px;
1349  }
1350  return (r1);
1351 }
1352 
1353 ////////////////////////////////////////////////////////////////////////////////
1354 /// This function calculates derivative of peaks shape function (see manual)
1355 /// according to relative amplitude sx.
1356 ///
1357 /// Function parameters:
1358 /// - numOfFittedPeaks-number of fitted peaks
1359 /// - x-position of channel
1360 /// - parameter-array of peaks parameters (amplitudes and positions)
1361 /// - sigmax-sigma of 1D ridge
1362 
1363 Double_t TSpectrum2Fit::Dersx(Int_t numOfFittedPeaks, Double_t x,
1364  const Double_t *parameter, Double_t sigmax)
1365 {
1366  Double_t p, r1 = 0, rx, ax, x0, s2;
1367  Int_t j;
1368  s2 = TMath::Sqrt(2.0);
1369  for (j = 0; j < numOfFittedPeaks; j++) {
1370  ax = parameter[7 * j + 3];
1371  x0 = parameter[7 * j + 5];
1372  p = (x - x0) / sigmax;
1373  s2 = TMath::Sqrt(2.0);
1374  rx = Erfc(p / s2);
1375  r1 += 0.5 * ax * rx;
1376  }
1377  return (r1);
1378 }
1379 
1380 ////////////////////////////////////////////////////////////////////////////////
1381 /// This function calculates derivative of peaks shape function (see manual)
1382 /// according to relative amplitude sy.
1383 ///
1384 /// Function parameters:
1385 /// - numOfFittedPeaks-number of fitted peaks
1386 /// - x-position of channel
1387 /// - parameter-array of peaks parameters (amplitudes and positions)
1388 /// - sigmax-sigma of 1D ridge
1389 
1390 Double_t TSpectrum2Fit::Dersy(Int_t numOfFittedPeaks, Double_t x,
1391  const Double_t *parameter, Double_t sigmax)
1392 {
1393  Double_t p, r1 = 0, rx, ax, x0, s2;
1394  Int_t j;
1395  s2 = TMath::Sqrt(2.0);
1396  for (j = 0; j < numOfFittedPeaks; j++) {
1397  ax = parameter[7 * j + 4];
1398  x0 = parameter[7 * j + 6];
1399  p = (x - x0) / sigmax;
1400  s2 = TMath::Sqrt(2.0);
1401  rx = Erfc(p / s2);
1402  r1 += 0.5 * ax * rx;
1403  }
1404  return (r1);
1405 }
1406 
1407 ////////////////////////////////////////////////////////////////////////////////
1408 /// This function calculates derivative of peaks shape function (see manual)
1409 /// according to slope bx.
1410 ///
1411 /// Function parameters:
1412 /// - numOfFittedPeaks-number of fitted peaks
1413 /// - x,y-position of channel
1414 /// - parameter-array of peaks parameters (amplitudes and positions)
1415 /// - sigmax-sigmax of peaks
1416 /// - sigmay-sigmay of peaks
1417 /// - txy, tx-relative amplitudes
1418 /// - bx, by-slopes
1419 
1420 Double_t TSpectrum2Fit::Derbx(Int_t numOfFittedPeaks, Double_t x, Double_t y,
1421  const Double_t *parameter, Double_t sigmax,
1422  Double_t sigmay, Double_t txy, Double_t tx, Double_t bx,
1423  Double_t by)
1424 {
1425  Double_t p, r, r1 = 0, a, x0, y0, s2, px, py, erx, ery, ex, ey;
1426  Int_t j;
1427  s2 = TMath::Sqrt(2.0);
1428  for (j = 0; j < numOfFittedPeaks; j++) {
1429  a = parameter[7 * j];
1430  x0 = parameter[7 * j + 1];
1431  y0 = parameter[7 * j + 2];
1432  p = (x - x0) / sigmax;
1433  r = (y - y0) / sigmay;
1434  if (txy != 0) {
1435  px = 0, py = 0;
1436  erx =
1437  -Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * bx) -
1438  Derfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * bx), ery =
1439  Erfc(r / s2 + 1 / (2 * by));
1440  ex = p / (s2 * bx), ey = r / (s2 * by);
1441  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
1442  px = exp(ex) * erx, py = exp(ey) * ery;
1443  }
1444  r1 += 0.5 * a * txy * px * py;
1445  }
1446  a = parameter[7 * j + 3];
1447  x0 = parameter[7 * j + 5];
1448  p = (x - x0) / sigmax;
1449  if (tx != 0) {
1450  px = 0;
1451  erx =
1452  (-Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * bx) -
1453  Derfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * bx));
1454  ex = p / (s2 * bx);
1455  if (TMath::Abs(ex) < 9)
1456  px = exp(ex) * erx;
1457  r1 += 0.5 * a * tx * px;
1458  }
1459  }
1460  return (r1);
1461 }
1462 
1463 ////////////////////////////////////////////////////////////////////////////////
1464 /// This function calculates derivative of peaks shape function (see manual)
1465 /// according to slope by.
1466 ///
1467 /// Function parameters:
1468 /// - numOfFittedPeaks-number of fitted peaks
1469 /// - x,y-position of channel
1470 /// - parameter-array of peaks parameters (amplitudes and positions)
1471 /// - sigmax-sigmax of peaks
1472 /// - sigmay-sigmay of peaks
1473 /// - txy, ty-relative amplitudes
1474 /// - bx, by-slopes
1475 
1476 Double_t TSpectrum2Fit::Derby(Int_t numOfFittedPeaks, Double_t x, Double_t y,
1477  const Double_t *parameter, Double_t sigmax,
1478  Double_t sigmay, Double_t txy, Double_t ty, Double_t bx,
1479  Double_t by)
1480 {
1481  Double_t p, r, r1 = 0, a, x0, y0, s2, px, py, erx, ery, ex, ey;
1482  Int_t j;
1483  s2 = TMath::Sqrt(2.0);
1484  for (j = 0; j < numOfFittedPeaks; j++) {
1485  a = parameter[7 * j];
1486  x0 = parameter[7 * j + 1];
1487  y0 = parameter[7 * j + 2];
1488  p = (x - x0) / sigmax;
1489  r = (y - y0) / sigmay;
1490  if (txy != 0) {
1491  px = 0, py = 0;
1492  ery =
1493  -Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * by) -
1494  Derfc(r / s2 + 1 / (2 * by)) / (s2 * by * by), erx =
1495  Erfc(p / s2 + 1 / (2 * bx));
1496  ex = p / (s2 * bx), ey = r / (s2 * by);
1497  if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
1498  px = exp(ex) * erx, py = exp(ey) * ery;
1499  }
1500  r1 += 0.5 * a * txy * px * py;
1501  }
1502  a = parameter[7 * j + 4];
1503  y0 = parameter[7 * j + 6];
1504  r = (y - y0) / sigmay;
1505  if (ty != 0) {
1506  py = 0;
1507  ery =
1508  (-Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * by) -
1509  Derfc(r / s2 + 1 / (2 * by)) / (s2 * by * by));
1510  ey = r / (s2 * by);
1511  if (TMath::Abs(ey) < 9)
1512  py = exp(ey) * ery;
1513  r1 += 0.5 * a * ty * py;
1514  }
1515  }
1516  return (r1);
1517 }
1518 
1519 ////////////////////////////////////////////////////////////////////////////////
1520 /// This function calculates volume of a peak
1521 ///
1522 /// Function parameters:
1523 /// - a-amplitude of the peak
1524 /// - sx,sy-sigmas of peak
1525 /// - ro-correlation coefficient
1526 
1527 Double_t TSpectrum2Fit::Volume(Double_t a, Double_t sx, Double_t sy, Double_t ro)
1528 {
1529  Double_t pi = 3.1415926535, r;
1530  r = 1 - ro * ro;
1531  if (r > 0)
1532  r = TMath::Sqrt(r);
1533 
1534  else {
1535  return (0);
1536  }
1537  r = 2 * a * pi * sx * sy * r;
1538  return (r);
1539 }
1540 
1541 ////////////////////////////////////////////////////////////////////////////////
1542 /// This function calculates derivative of the volume of a peak
1543 /// according to amplitude
1544 ///
1545 /// Function parameters:
1546 /// - sx,sy-sigmas of peak
1547 /// - ro-correlation coefficient
1548 
1549 Double_t TSpectrum2Fit::Derpa2(Double_t sx, Double_t sy, Double_t ro)
1550 {
1551  Double_t pi = 3.1415926535, r;
1552  r = 1 - ro * ro;
1553  if (r > 0)
1554  r = TMath::Sqrt(r);
1555 
1556  else {
1557  return (0);
1558  }
1559  r = 2 * pi * sx * sy * r;
1560  return (r);
1561 }
1562 
1563 ////////////////////////////////////////////////////////////////////////////////
1564 /// This function calculates derivative of the volume of a peak
1565 /// according to sigmax
1566 ///
1567 /// Function parameters:
1568 /// - a-amplitude of peak
1569 /// - sy-sigma of peak
1570 /// - ro-correlation coefficient
1571 
1572 Double_t TSpectrum2Fit::Derpsigmax(Double_t a, Double_t sy, Double_t ro)
1573 {
1574  Double_t pi = 3.1415926535, r;
1575  r = 1 - ro * ro;
1576  if (r > 0)
1577  r = TMath::Sqrt(r);
1578 
1579  else {
1580  return (0);
1581  }
1582  r = a * 2 * pi * sy * r;
1583  return (r);
1584 }
1585 
1586 ////////////////////////////////////////////////////////////////////////////////
1587 /// This function calculates derivative of the volume of a peak
1588 /// according to sigmay
1589 ///
1590 /// Function parameters:
1591 /// - a-amplitude of peak
1592 /// - sx-sigma of peak
1593 /// - ro-correlation coefficient
1594 
1595 Double_t TSpectrum2Fit::Derpsigmay(Double_t a, Double_t sx, Double_t ro)
1596 {
1597  Double_t pi = 3.1415926535, r;
1598  r = 1 - ro * ro;
1599  if (r > 0)
1600  r = TMath::Sqrt(r);
1601 
1602  else {
1603  return (0);
1604  }
1605  r = a * 2 * pi * sx * r;
1606  return (r);
1607 }
1608 
1609 ////////////////////////////////////////////////////////////////////////////////
1610 /// This function calculates derivative of the volume of a peak
1611 /// according to ro
1612 ///
1613 /// Function parameters:
1614 /// - a-amplitude of peak
1615 /// - sx,sy-sigmas of peak
1616 /// - ro-correlation coefficient
1617 
1618 Double_t TSpectrum2Fit::Derpro(Double_t a, Double_t sx, Double_t sy, Double_t ro)
1619 {
1620  Double_t pi = 3.1415926535, r;
1621  r = 1 - ro * ro;
1622  if (r > 0)
1623  r = TMath::Sqrt(r);
1624 
1625  else {
1626  return (0);
1627  }
1628  r = -a * 2 * pi * sx * sy * ro / r;
1629  return (r);
1630 }
1631 
1632 
1633 ////////////////////////////////////////////////////////////////////////////////
1634 /// This function fits the source spectrum. The calling program should
1635 /// fill in input parameters of the TSpectrum2Fit class.
1636 /// The fitted parameters are written into
1637 /// TSpectrum2Fit class output parameters and fitted data are written into
1638 /// source spectrum.
1639 ///
1640 /// Function parameters:
1641 /// - source-pointer to the matrix of source spectrum
1642 ///
1643 /// ### Fitting
1644 ///
1645 /// Goal: to estimate simultaneously peak shape parameters in spectra with large
1646 /// number of peaks
1647 ///
1648 /// - peaks can be fitted separately, each peak (or multiplets) in a region or
1649 /// together all peaks in a spectrum. To fit separately each peak one needs to
1650 /// determine the fitted region. However it can happen that the regions of
1651 /// neighbouring peaks are overlapping. Then the results of fitting are very poor.
1652 /// On the other hand, when fitting together all peaks found in a spectrum, one
1653 /// needs to have a method that is stable (converges) and fast enough to carry out
1654 /// fitting in reasonable time
1655 ///
1656 /// - we have implemented the non-symmetrical semiempirical peak shape function
1657 ///
1658 /// - it contains the two-dimensional symmetrical Gaussian two one-dimensional
1659 /// symmetrical Gaussian ridges as well as non-symmetrical terms and background.
1660 ///
1661 /// \image html spectrum2fit_awmi_image001.gif
1662 ///
1663 /// where Txy, Tx, Ty, Sxy, Sx, Sy are relative amplitudes and Bx, By are slopes.
1664 ///
1665 /// - algorithm without matrix inversion (AWMI) allows fitting tens, hundreds
1666 /// of peaks simultaneously that represent sometimes thousands of parameters [2],[5].
1667 ///
1668 /// #### References:
1669 ///
1670 /// [1] Phillps G.W., Marlow K.W.,
1671 /// NIM 137 (1976) 525.
1672 ///
1673 /// [2] I. A. Slavic: Nonlinear
1674 /// least-squares fitting without matrix inversion applied to complex Gaussian
1675 /// spectra analysis. NIM 134 (1976) 285-289.
1676 ///
1677 /// [3] T. Awaya: A new method for
1678 /// curve fitting to the data with low statistics not using chi-square method. NIM
1679 /// 165 (1979) 317-323.
1680 ///
1681 /// [4] T. Hauschild, M. Jentschel:
1682 /// Comparison of maximum likelihood estimation and chi-square statistics applied
1683 /// to counting experiments. NIM A 457 (2001) 384-401.
1684 ///
1685 /// [5] M. Morh, J.
1686 /// Kliman, M. Jandel, Krupa, V. Matouoek: Study of fitting algorithms
1687 /// applied to simultaneous analysis of large number of peaks in -ray spectra.
1688 /// Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003
1689 ///
1690 /// ### Example 1 - script FitAwmi2.c:
1691 ///
1692 /// \image html spectrum2fit_awmi_image002.jpg Original Fig. 1 two-dimensional spectrum with found peaks (using TSpectrum2 peak searching function). The positions of peaks were used as initial estimates in fitting procedure.
1693 ///
1694 /// \image html spectrum2fit_awmi_image003.jpg Fig. 2 Fitted (generated from fitted parameters) spectrum of the data from Fig. 1 using Algorithm Without Matrix Inversion. Each peak was represented by 7 parameters, which together with Sigmax, Sigmay and a0 resulted in 38 fitted parameters. The chi-squareafter 1000 iterations was 0.642342.
1695 ///
1696 /// #### Script:
1697 ///
1698 /// Example to illustrate fitting function, algorithm without matrix inversion (AWMI) (class TSpectrumFit2).
1699 /// To execute this example, do
1700 ///
1701 /// `root > .x FitAwmi2.C`
1702 ///
1703 /// ~~~ {.cpp}
1704 /// void FitAwmi2() {
1705 /// Int_t i, j, nfound;
1706 /// Int_t nbinsx = 64;
1707 /// Int_t nbinsy = 64;
1708 /// Int_t xmin = 0;
1709 /// Int_t xmax = nbinsx;
1710 /// Int_t ymin = 0;
1711 /// Int_t ymax = nbinsy;
1712 /// Double_t ** source = new float *[nbinsx];
1713 /// Double_t ** dest = new float *[nbinsx];
1714 /// for (i=0;i<nbinsx;i++)
1715 /// source[i]=new float[nbinsy];
1716 /// for (i=0;i<nbinsx;i++)
1717 /// dest[i]=new float[nbinsy];
1718 /// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1719 /// TFile *f = new TFile("TSpectrum2.root");
1720 /// search=(TH2F*) f->Get("search4;1");
1721 /// TCanvas *Searching = new TCanvas("Searching","Two-dimensional fitting using Algorithm Without Matrix Inversion",10,10,1000,700);
1722 /// TSpectrum2 *s = new TSpectrum2();
1723 /// for (i = 0; i < nbinsx; i++){
1724 /// for (j = 0; j < nbinsy; j++){
1725 /// source[i][j] = search->GetBinContent(i + 1,j + 1);
1726 /// }
1727 /// }
1728 /// //searching for candidate peaks positions
1729 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);
1730 /// Bool_t *FixPosX = new Bool_t[nfound];
1731 /// Bool_t *FixPosY = new Bool_t[nfound];
1732 /// Bool_t *FixAmp = new Bool_t[nfound];
1733 /// Double_t *PosX = new Double_t[nfound];
1734 /// Double_t *PosY = new Double_t[nfound];
1735 /// Double_t *Amp = new Double_t[nfound];
1736 /// Double_t *AmpXY = new Double_t[nfound];
1737 /// PosX = s->GetPositionX();
1738 /// PosY = s->GetPositionY();
1739 /// printf("Found %d candidate peaks\n",nfound);
1740 /// for(i = 0; i< nfound ; i++){
1741 /// FixPosX[i] = kFALSE;
1742 /// FixPosY[i] = kFALSE;
1743 /// FixAmp[i] = kFALSE;
1744 /// Amp[i] = source[(Int_t)(PosX[i]+0.5)][(Int_t)(PosY[i]+0.5)]; //initial values of peaks amplitudes, input parameters
1745 /// AmpXY[i] = 0;
1746 /// }
1747 /// //filling in the initial estimates of the input parameters
1748 /// TSpectrumFit2 *pfit=new TSpectrumFit2(nfound);
1749 /// pfit->SetFitParameters(xmin, xmax-1, ymin, ymax-1, 1000, 0.1, pfit->kFitOptimChiCounts,
1750 /// pfit->kFitAlphaHalving, pfit->kFitPower2,
1751 /// pfit->kFitTaylorOrderFirst);
1752 /// pfit->SetPeakParameters(2, kFALSE, 2, kFALSE, 0, kTRUE, PosX, (Bool_t *)
1753 /// FixPosX, PosY, (Bool_t *) FixPosY, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *)
1754 /// FixPosY, Amp, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp, AmpXY, (Bool_t *)
1755 /// FixAmp);
1756 /// pfit->SetBackgroundParameters(0, kFALSE, 0, kTRUE, 0, kTRUE);
1757 /// pfit->FitAwmi(source);
1758 /// for (i = 0; i < nbinsx; i++){
1759 /// for (j = 0; j < nbinsy; j++){
1760 /// search->SetBinContent(i + 1, j + 1,source[i][j]);
1761 /// }
1762 /// }
1763 /// search->Draw("SURF");
1764 /// }
1765 /// ~~~
1766 ///
1767 /// ### Example 2 - script FitAwmi2.c:
1768 ///
1769 /// \image html spectrum2fit_awmi_image004.jpg Fig. 3 Original two-dimensional gamma-gamma-ray spectrum with found peaks (using TSpectrum2 peak searching function).
1770 ///
1771 /// \image html spectrum2fit_awmi_image005.jpg Fig. 4 Fitted (generated from fitted parameters) spectrum of the data from Fig. 3 using Algorithm Without Matrix Inversion. 152 peaks were identified. Each peak was represented by 7 parameters, which together with Sigmax, Sigmay and a0 resulted in 1067 fitted parameters. The chi-square after 1000 iterations was 0.728675. One can observe good correspondence with the original data.
1772 ///
1773 /// #### Script:
1774 ///
1775 ////
1776 /// Example to illustrate fitting function, algorithm without matrix inversion
1777 /// (AWMI) (class TSpectrumFit2). To execute this example, do:
1778 ///
1779 /// `root > .x FitA2.C`
1780 ///
1781 /// ~~~ {.cpp}
1782 /// void FitA2() {
1783 /// Int_t i, j, nfound;
1784 /// Int_t nbinsx = 256;
1785 /// Int_t nbinsy = 256;
1786 /// Int_t xmin = 0;
1787 /// Int_t xmax = nbinsx;
1788 /// Int_t ymin = 0;
1789 /// Int_t ymax = nbinsy;
1790 /// Double_t ** source = new float *[nbinsx];
1791 /// Double_t ** dest = new float *[nbinsx];
1792 /// for (i=0;i<nbinsx;i++)
1793 /// source[i]=new
1794 /// float[nbinsy];
1795 /// for (i=0;i<nbinsx;i++)
1796 /// dest[i]=new
1797 /// float[nbinsy];
1798 /// TH2F *search = new TH2F("search","High resolution peak
1799 /// searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
1800 /// TFile *f = new TFile("TSpectrum2.root");
1801 /// search=(TH2F*) f->Get("fit1;1");
1802 /// TCanvas *Searching = new TCanvas("Searching","Two-dimensional fitting using Algorithm Without Matrix Inversion",10,10,1000,700);
1803 /// TSpectrum2 *s = new TSpectrum2(1000,1);
1804 /// for (i = 0; i < nbinsx; i++){
1805 /// for (j = 0; j < nbinsy; j++){
1806 /// source[i][j] = search->GetBinContent(i + 1,j + 1);
1807 /// }
1808 /// }
1809 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 2, kTRUE, 100, kFALSE, 3);
1810 /// printf("Found %d candidate peaks\n",nfound);
1811 /// Bool_t *FixPosX = new Bool_t[nfound];
1812 /// Bool_t *FixPosY = new Bool_t[nfound];
1813 /// Bool_t *FixAmp = new Bool_t[nfound];
1814 /// Double_t *PosX = new Double_t[nfound];
1815 /// Double_t *PosY = new Double_t[nfound];
1816 /// Double_t *Amp = new Double_t[nfound];
1817 /// Double_t *AmpXY = new Double_t[nfound];
1818 /// PosX = s->GetPositionX();
1819 /// PosY = s->GetPositionY();
1820 /// for(i = 0; i< nfound ; i++){
1821 /// FixPosX[i] = kFALSE;
1822 /// FixPosY[i] = kFALSE;
1823 /// FixAmp[i] = kFALSE;
1824 /// Amp[i] = source[(Int_t)(PosX[i]+0.5)][(Int_t)(PosY[i]+0.5)]; //initial values of peaks amplitudes, input parameters
1825 /// AmpXY[i] = 0;
1826 /// }
1827 /// //filling in the initial estimates of the input parameters
1828 /// TSpectrumFit2 *pfit=new TSpectrumFit2(nfound);
1829 /// pfit->SetFitParameters(xmin, xmax-1, ymin, ymax-1, 1000, 0.1,
1830 /// pfit->kFitOptimChiCounts, pfit->kFitAlphaHalving, pfit->kFitPower2,
1831 /// pfit->kFitTaylorOrderFirst);
1832 /// pfit->SetPeakParameters(2, kFALSE, 2, kFALSE, 0, kTRUE, PosX, (Bool_t *)
1833 /// FixPosX, PosY, (Bool_t *) FixPosY, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *)
1834 /// FixPosY, Amp, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp, AmpXY, (Bool_t *)
1835 /// FixAmp);
1836 /// pfit->SetBackgroundParameters(0, kFALSE, 0, kTRUE, 0, kTRUE);
1837 /// pfit->FitAwmi(source);
1838 /// for (i = 0; i < nbinsx; i++){
1839 /// for (j = 0; j < nbinsy; j++){
1840 /// search->SetBinContent(i + 1, j + 1,source[i][j]);
1841 /// }
1842 /// }
1843 /// search->Draw("SURF");
1844 /// }
1845 /// ~~~
1846 
1847 void TSpectrum2Fit::FitAwmi(Double_t **source)
1848 {
1849 
1850 
1851  Int_t i, i1, i2, j, k, shift =
1852  7 * fNPeaks + 14, peak_vel, size, iter, pw,
1853  regul_cycle, flag;
1854  Double_t a, b, c, d = 0, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
1855  0, pi, pmin = 0, chi_cel = 0, chi_er;
1856  Double_t *working_space = new Double_t[5 * (7 * fNPeaks + 14)];
1857  for (i = 0, j = 0; i < fNPeaks; i++) {
1858  working_space[7 * i] = fAmpInit[i]; //vector parameter
1859  if (fFixAmp[i] == false) {
1860  working_space[shift + j] = fAmpInit[i]; //vector xk
1861  j += 1;
1862  }
1863  working_space[7 * i + 1] = fPositionInitX[i]; //vector parameter
1864  if (fFixPositionX[i] == false) {
1865  working_space[shift + j] = fPositionInitX[i]; //vector xk
1866  j += 1;
1867  }
1868  working_space[7 * i + 2] = fPositionInitY[i]; //vector parameter
1869  if (fFixPositionY[i] == false) {
1870  working_space[shift + j] = fPositionInitY[i]; //vector xk
1871  j += 1;
1872  }
1873  working_space[7 * i + 3] = fAmpInitX1[i]; //vector parameter
1874  if (fFixAmpX1[i] == false) {
1875  working_space[shift + j] = fAmpInitX1[i]; //vector xk
1876  j += 1;
1877  }
1878  working_space[7 * i + 4] = fAmpInitY1[i]; //vector parameter
1879  if (fFixAmpY1[i] == false) {
1880  working_space[shift + j] = fAmpInitY1[i]; //vector xk
1881  j += 1;
1882  }
1883  working_space[7 * i + 5] = fPositionInitX1[i]; //vector parameter
1884  if (fFixPositionX1[i] == false) {
1885  working_space[shift + j] = fPositionInitX1[i]; //vector xk
1886  j += 1;
1887  }
1888  working_space[7 * i + 6] = fPositionInitY1[i]; //vector parameter
1889  if (fFixPositionY1[i] == false) {
1890  working_space[shift + j] = fPositionInitY1[i]; //vector xk
1891  j += 1;
1892  }
1893  }
1894  peak_vel = 7 * i;
1895  working_space[7 * i] = fSigmaInitX; //vector parameter
1896  if (fFixSigmaX == false) {
1897  working_space[shift + j] = fSigmaInitX; //vector xk
1898  j += 1;
1899  }
1900  working_space[7 * i + 1] = fSigmaInitY; //vector parameter
1901  if (fFixSigmaY == false) {
1902  working_space[shift + j] = fSigmaInitY; //vector xk
1903  j += 1;
1904  }
1905  working_space[7 * i + 2] = fRoInit; //vector parameter
1906  if (fFixRo == false) {
1907  working_space[shift + j] = fRoInit; //vector xk
1908  j += 1;
1909  }
1910  working_space[7 * i + 3] = fA0Init; //vector parameter
1911  if (fFixA0 == false) {
1912  working_space[shift + j] = fA0Init; //vector xk
1913  j += 1;
1914  }
1915  working_space[7 * i + 4] = fAxInit; //vector parameter
1916  if (fFixAx == false) {
1917  working_space[shift + j] = fAxInit; //vector xk
1918  j += 1;
1919  }
1920  working_space[7 * i + 5] = fAyInit; //vector parameter
1921  if (fFixAy == false) {
1922  working_space[shift + j] = fAyInit; //vector xk
1923  j += 1;
1924  }
1925  working_space[7 * i + 6] = fTxyInit; //vector parameter
1926  if (fFixTxy == false) {
1927  working_space[shift + j] = fTxyInit; //vector xk
1928  j += 1;
1929  }
1930  working_space[7 * i + 7] = fSxyInit; //vector parameter
1931  if (fFixSxy == false) {
1932  working_space[shift + j] = fSxyInit; //vector xk
1933  j += 1;
1934  }
1935  working_space[7 * i + 8] = fTxInit; //vector parameter
1936  if (fFixTx == false) {
1937  working_space[shift + j] = fTxInit; //vector xk
1938  j += 1;
1939  }
1940  working_space[7 * i + 9] = fTyInit; //vector parameter
1941  if (fFixTy == false) {
1942  working_space[shift + j] = fTyInit; //vector xk
1943  j += 1;
1944  }
1945  working_space[7 * i + 10] = fSxyInit; //vector parameter
1946  if (fFixSx == false) {
1947  working_space[shift + j] = fSxInit; //vector xk
1948  j += 1;
1949  }
1950  working_space[7 * i + 11] = fSyInit; //vector parameter
1951  if (fFixSy == false) {
1952  working_space[shift + j] = fSyInit; //vector xk
1953  j += 1;
1954  }
1955  working_space[7 * i + 12] = fBxInit; //vector parameter
1956  if (fFixBx == false) {
1957  working_space[shift + j] = fBxInit; //vector xk
1958  j += 1;
1959  }
1960  working_space[7 * i + 13] = fByInit; //vector parameter
1961  if (fFixBy == false) {
1962  working_space[shift + j] = fByInit; //vector xk
1963  j += 1;
1964  }
1965  size = j;
1966  for (iter = 0; iter < fNumberIterations; iter++) {
1967  for (j = 0; j < size; j++) {
1968  working_space[2 * shift + j] = 0, working_space[3 * shift + j] = 0; //der,temp
1969  }
1970 
1971  //filling vectors
1972  alpha = fAlpha;
1973  chi_opt = 0, pw = fPower - 2;
1974  for (i1 = fXmin; i1 <= fXmax; i1++) {
1975  for (i2 = fYmin; i2 <= fYmax; i2++) {
1976  yw = source[i1][i2];
1977  ywm = yw;
1978  f = Shape2(fNPeaks, i1, i2,
1979  working_space, working_space[peak_vel],
1980  working_space[peak_vel + 1],
1981  working_space[peak_vel + 2],
1982  working_space[peak_vel + 3],
1983  working_space[peak_vel + 4],
1984  working_space[peak_vel + 5],
1985  working_space[peak_vel + 6],
1986  working_space[peak_vel + 7],
1987  working_space[peak_vel + 8],
1988  working_space[peak_vel + 9],
1989  working_space[peak_vel + 10],
1990  working_space[peak_vel + 11],
1991  working_space[peak_vel + 12],
1992  working_space[peak_vel + 13]);
1993  if (fStatisticType == kFitOptimMaxLikelihood) {
1994  if (f > 0.00001)
1995  chi_opt += yw * TMath::Log(f) - f;
1996  }
1997 
1998  else {
1999  if (ywm != 0)
2000  chi_opt += (yw - f) * (yw - f) / ywm;
2001  }
2002  if (fStatisticType == kFitOptimChiFuncValues) {
2003  ywm = f;
2004  if (f < 0.00001)
2005  ywm = 0.00001;
2006  }
2007 
2008  else if (fStatisticType == kFitOptimMaxLikelihood) {
2009  ywm = f;
2010  if (f < 0.00001)
2011  ywm = 0.00001;
2012  }
2013 
2014  else {
2015  if (ywm == 0)
2016  ywm = 1;
2017  }
2018 
2019  //calculation of gradient vector
2020  for (j = 0, k = 0; j < fNPeaks; j++) {
2021  if (fFixAmp[j] == false) {
2022  a = Deramp2(i1, i2,
2023  working_space[7 * j + 1],
2024  working_space[7 * j + 2],
2025  working_space[peak_vel],
2026  working_space[peak_vel + 1],
2027  working_space[peak_vel + 2],
2028  working_space[peak_vel + 6],
2029  working_space[peak_vel + 7],
2030  working_space[peak_vel + 12],
2031  working_space[peak_vel + 13]);
2032  if (ywm != 0) {
2033  c = Ourpowl(a, pw);
2034  if (fStatisticType == kFitOptimChiFuncValues) {
2035  b = a * (yw * yw - f * f) / (ywm * ywm);
2036  working_space[2 * shift + k] += b * c; //der
2037  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2038  working_space[3 * shift + k] += b * c; //temp
2039  }
2040 
2041  else {
2042  b = a * (yw - f) / ywm;
2043  working_space[2 * shift + k] += b * c; //der
2044  b = a * a / ywm;
2045  working_space[3 * shift + k] += b * c; //temp
2046  }
2047  }
2048  k += 1;
2049  }
2050  if (fFixPositionX[j] == false) {
2051  a = Deri02(i1, i2,
2052  working_space[7 * j],
2053  working_space[7 * j + 1],
2054  working_space[7 * j + 2],
2055  working_space[peak_vel],
2056  working_space[peak_vel + 1],
2057  working_space[peak_vel + 2],
2058  working_space[peak_vel + 6],
2059  working_space[peak_vel + 7],
2060  working_space[peak_vel + 12],
2061  working_space[peak_vel + 13]);
2062  if (fFitTaylor == kFitTaylorOrderSecond)
2063  d = Derderi02(i1, i2,
2064  working_space[7 * j],
2065  working_space[7 * j + 1],
2066  working_space[7 * j + 2],
2067  working_space[peak_vel],
2068  working_space[peak_vel + 1],
2069  working_space[peak_vel + 2]);
2070  if (ywm != 0) {
2071  c = Ourpowl(a, pw);
2072  if (TMath::Abs(a) > 0.00000001
2073  && fFitTaylor == kFitTaylorOrderSecond) {
2074  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
2075  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
2076  && a <= 0))
2077  d = 0;
2078  }
2079 
2080  else
2081  d = 0;
2082  a = a + d;
2083  if (fStatisticType == kFitOptimChiFuncValues) {
2084  b = a * (yw * yw - f * f) / (ywm * ywm);
2085  working_space[2 * shift + k] += b * c; //der
2086  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2087  working_space[3 * shift + k] += b * c; //temp
2088  }
2089 
2090  else {
2091  b = a * (yw - f) / ywm;
2092  working_space[2 * shift + k] += b * c; //der
2093  b = a * a / ywm;
2094  working_space[3 * shift + k] += b * c; //temp
2095  }
2096  }
2097  k += 1;
2098  }
2099  if (fFixPositionY[j] == false) {
2100  a = Derj02(i1, i2,
2101  working_space[7 * j],
2102  working_space[7 * j + 1],
2103  working_space[7 * j + 2],
2104  working_space[peak_vel],
2105  working_space[peak_vel + 1],
2106  working_space[peak_vel + 2],
2107  working_space[peak_vel + 6],
2108  working_space[peak_vel + 7],
2109  working_space[peak_vel + 12],
2110  working_space[peak_vel + 13]);
2111  if (fFitTaylor == kFitTaylorOrderSecond)
2112  d = Derderj02(i1, i2,
2113  working_space[7 * j],
2114  working_space[7 * j + 1],
2115  working_space[7 * j + 2],
2116  working_space[peak_vel],
2117  working_space[peak_vel + 1],
2118  working_space[peak_vel + 2]);
2119  if (ywm != 0) {
2120  c = Ourpowl(a, pw);
2121  if (TMath::Abs(a) > 0.00000001
2122  && fFitTaylor == kFitTaylorOrderSecond) {
2123  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
2124  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
2125  && a <= 0))
2126  d = 0;
2127  }
2128 
2129  else
2130  d = 0;
2131  a = a + d;
2132  if (fStatisticType == kFitOptimChiFuncValues) {
2133  b = a * (yw * yw - f * f) / (ywm * ywm);
2134  working_space[2 * shift + k] += b * c; //der
2135  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2136  working_space[3 * shift + k] += b * c; //temp
2137  }
2138 
2139  else {
2140  b = a * (yw - f) / ywm;
2141  working_space[2 * shift + k] += b * c; //der
2142  b = a * a / ywm;
2143  working_space[3 * shift + k] += b * c; //temp
2144  }
2145  }
2146  k += 1;
2147  }
2148  if (fFixAmpX1[j] == false) {
2149  a = Derampx(i1, working_space[7 * j + 5],
2150  working_space[peak_vel],
2151  working_space[peak_vel + 8],
2152  working_space[peak_vel + 10],
2153  working_space[peak_vel + 12]);
2154  if (ywm != 0) {
2155  c = Ourpowl(a, pw);
2156  if (fStatisticType == kFitOptimChiFuncValues) {
2157  b = a * (yw * yw - f * f) / (ywm * ywm);
2158  working_space[2 * shift + k] += b * c; //der
2159  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2160  working_space[3 * shift + k] += b * c; //temp
2161  }
2162 
2163  else {
2164  b = a * (yw - f) / ywm;
2165  working_space[2 * shift + k] += b * c; //der
2166  b = a * a / ywm;
2167  working_space[3 * shift + k] += b * c; //temp
2168  }
2169  }
2170  k += 1;
2171  }
2172  if (fFixAmpY1[j] == false) {
2173  a = Derampx(i2, working_space[7 * j + 6],
2174  working_space[peak_vel + 1],
2175  working_space[peak_vel + 9],
2176  working_space[peak_vel + 11],
2177  working_space[peak_vel + 13]);
2178  if (ywm != 0) {
2179  c = Ourpowl(a, pw);
2180  if (fStatisticType == kFitOptimChiFuncValues) {
2181  b = a * (yw * yw - f * f) / (ywm * ywm);
2182  working_space[2 * shift + k] += b * c; //der
2183  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2184  working_space[3 * shift + k] += b * c; //temp
2185  }
2186 
2187  else {
2188  b = a * (yw - f) / ywm;
2189  working_space[2 * shift + k] += b * c; //der
2190  b = a * a / ywm;
2191  working_space[3 * shift + k] += b * c; //temp
2192  }
2193  }
2194  k += 1;
2195  }
2196  if (fFixPositionX1[j] == false) {
2197  a = Deri01(i1, working_space[7 * j + 3],
2198  working_space[7 * j + 5],
2199  working_space[peak_vel],
2200  working_space[peak_vel + 8],
2201  working_space[peak_vel + 10],
2202  working_space[peak_vel + 12]);
2203  if (fFitTaylor == kFitTaylorOrderSecond)
2204  d = Derderi01(i1, working_space[7 * j + 3],
2205  working_space[7 * j + 5],
2206  working_space[peak_vel]);
2207  if (ywm != 0) {
2208  c = Ourpowl(a, pw);
2209  if (TMath::Abs(a) > 0.00000001
2210  && fFitTaylor == kFitTaylorOrderSecond) {
2211  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
2212  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
2213  && a <= 0))
2214  d = 0;
2215  }
2216 
2217  else
2218  d = 0;
2219  a = a + d;
2220  if (fStatisticType == kFitOptimChiFuncValues) {
2221  b = a * (yw * yw - f * f) / (ywm * ywm);
2222  working_space[2 * shift + k] += b * c; //Der
2223  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2224  working_space[3 * shift + k] += b * c; //temp
2225  }
2226 
2227  else {
2228  b = a * (yw - f) / ywm;
2229  working_space[2 * shift + k] += b * c; //Der
2230  b = a * a / ywm;
2231  working_space[3 * shift + k] += b * c; //temp
2232  }
2233  }
2234  k += 1;
2235  }
2236  if (fFixPositionY1[j] == false) {
2237  a = Deri01(i2, working_space[7 * j + 4],
2238  working_space[7 * j + 6],
2239  working_space[peak_vel + 1],
2240  working_space[peak_vel + 9],
2241  working_space[peak_vel + 11],
2242  working_space[peak_vel + 13]);
2243  if (fFitTaylor == kFitTaylorOrderSecond)
2244  d = Derderi01(i2, working_space[7 * j + 4],
2245  working_space[7 * j + 6],
2246  working_space[peak_vel + 1]);
2247  if (ywm != 0) {
2248  c = Ourpowl(a, pw);
2249  if (TMath::Abs(a) > 0.00000001
2250  && fFitTaylor == kFitTaylorOrderSecond) {
2251  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
2252  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
2253  && a <= 0))
2254  d = 0;
2255  }
2256 
2257  else
2258  d = 0;
2259  a = a + d;
2260  if (fStatisticType == kFitOptimChiFuncValues) {
2261  b = a * (yw * yw - f * f) / (ywm * ywm);
2262  working_space[2 * shift + k] += b * c; //der
2263  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2264  working_space[3 * shift + k] += b * c; //temp
2265  }
2266 
2267  else {
2268  b = a * (yw - f) / ywm;
2269  working_space[2 * shift + k] += b * c; //der
2270  b = a * a / ywm;
2271  working_space[3 * shift + k] += b * c; //temp
2272  }
2273  }
2274  k += 1;
2275  }
2276  }
2277  if (fFixSigmaX == false) {
2278  a = Dersigmax(fNPeaks, i1, i2,
2279  working_space, working_space[peak_vel],
2280  working_space[peak_vel + 1],
2281  working_space[peak_vel + 2],
2282  working_space[peak_vel + 6],
2283  working_space[peak_vel + 7],
2284  working_space[peak_vel + 8],
2285  working_space[peak_vel + 10],
2286  working_space[peak_vel + 12],
2287  working_space[peak_vel + 13]);
2288  if (fFitTaylor == kFitTaylorOrderSecond)
2289  d = Derdersigmax(fNPeaks, i1,
2290  i2, working_space,
2291  working_space[peak_vel],
2292  working_space[peak_vel + 1],
2293  working_space[peak_vel + 2]);
2294  if (ywm != 0) {
2295  c = Ourpowl(a, pw);
2296  if (TMath::Abs(a) > 0.00000001
2297  && fFitTaylor == kFitTaylorOrderSecond) {
2298  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
2299  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
2300  d = 0;
2301  }
2302 
2303  else
2304  d = 0;
2305  a = a + d;
2306  if (fStatisticType == kFitOptimChiFuncValues) {
2307  b = a * (yw * yw - f * f) / (ywm * ywm);
2308  working_space[2 * shift + k] += b * c; //der
2309  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2310  working_space[3 * shift + k] += b * c; //temp
2311  }
2312 
2313  else {
2314  b = a * (yw - f) / ywm;
2315  working_space[2 * shift + k] += b * c; //der
2316  b = a * a / ywm;
2317  working_space[3 * shift + k] += b * c; //temp
2318  }
2319  }
2320  k += 1;
2321  }
2322  if (fFixSigmaY == false) {
2323  a = Dersigmay(fNPeaks, i1, i2,
2324  working_space, working_space[peak_vel],
2325  working_space[peak_vel + 1],
2326  working_space[peak_vel + 2],
2327  working_space[peak_vel + 6],
2328  working_space[peak_vel + 7],
2329  working_space[peak_vel + 9],
2330  working_space[peak_vel + 11],
2331  working_space[peak_vel + 12],
2332  working_space[peak_vel + 13]);
2333  if (fFitTaylor == kFitTaylorOrderSecond)
2334  d = Derdersigmay(fNPeaks, i1,
2335  i2, working_space,
2336  working_space[peak_vel],
2337  working_space[peak_vel + 1],
2338  working_space[peak_vel + 2]);
2339  if (ywm != 0) {
2340  c = Ourpowl(a, pw);
2341  if (TMath::Abs(a) > 0.00000001
2342  && fFitTaylor == kFitTaylorOrderSecond) {
2343  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
2344  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
2345  d = 0;
2346  }
2347 
2348  else
2349  d = 0;
2350  a = a + d;
2351  if (fStatisticType == kFitOptimChiFuncValues) {
2352  b = a * (yw * yw - f * f) / (ywm * ywm);
2353  working_space[2 * shift + k] += b * c; //der
2354  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2355  working_space[3 * shift + k] += b * c; //temp
2356  }
2357 
2358  else {
2359  b = a * (yw - f) / ywm;
2360  working_space[2 * shift + k] += b * c; //der
2361  b = a * a / ywm;
2362  working_space[3 * shift + k] += b * c; //temp
2363  }
2364  }
2365  k += 1;
2366  }
2367  if (fFixRo == false) {
2368  a = Derro(fNPeaks, i1, i2,
2369  working_space, working_space[peak_vel],
2370  working_space[peak_vel + 1],
2371  working_space[peak_vel + 2]);
2372  if (ywm != 0) {
2373  c = Ourpowl(a, pw);
2374  if (TMath::Abs(a) > 0.00000001
2375  && fFitTaylor == kFitTaylorOrderSecond) {
2376  d = d * TMath::Abs(yw - f) / (2 * a * ywm);
2377  if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
2378  d = 0;
2379  }
2380 
2381  else
2382  d = 0;
2383  a = a + d;
2384  if (fStatisticType == kFitOptimChiFuncValues) {
2385  b = a * (yw * yw - f * f) / (ywm * ywm);
2386  working_space[2 * shift + k] += b * c; //der
2387  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2388  working_space[3 * shift + k] += b * c; //temp
2389  }
2390 
2391  else {
2392  b = a * (yw - f) / ywm;
2393  working_space[2 * shift + k] += b * c; //der
2394  b = a * a / ywm;
2395  working_space[3 * shift + k] += b * c; //temp
2396  }
2397  }
2398  k += 1;
2399  }
2400  if (fFixA0 == false) {
2401  a = 1.;
2402  if (ywm != 0) {
2403  c = Ourpowl(a, pw);
2404  if (fStatisticType == kFitOptimChiFuncValues) {
2405  b = a * (yw * yw - f * f) / (ywm * ywm);
2406  working_space[2 * shift + k] += b * c; //der
2407  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2408  working_space[3 * shift + k] += b * c; //temp
2409  }
2410 
2411  else {
2412  b = a * (yw - f) / ywm;
2413  working_space[2 * shift + k] += b * c; //der
2414  b = a * a / ywm;
2415  working_space[3 * shift + k] += b * c; //temp
2416  }
2417  }
2418  k += 1;
2419  }
2420  if (fFixAx == false) {
2421  a = i1;
2422  if (ywm != 0) {
2423  c = Ourpowl(a, pw);
2424  if (fStatisticType == kFitOptimChiFuncValues) {
2425  b = a * (yw * yw - f * f) / (ywm * ywm);
2426  working_space[2 * shift + k] += b * c; //der
2427  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2428  working_space[3 * shift + k] += b * c; //temp
2429  }
2430 
2431  else {
2432  b = a * (yw - f) / ywm;
2433  working_space[2 * shift + k] += b * c; //der
2434  b = a * a / ywm;
2435  working_space[3 * shift + k] += b * c; //temp
2436  }
2437  }
2438  k += 1;
2439  }
2440  if (fFixAy == false) {
2441  a = i2;
2442  if (ywm != 0) {
2443  c = Ourpowl(a, pw);
2444  if (fStatisticType == kFitOptimChiFuncValues) {
2445  b = a * (yw * yw - f * f) / (ywm * ywm);
2446  working_space[2 * shift + k] += b * c; //der
2447  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2448  working_space[3 * shift + k] += b * c; //temp
2449  }
2450 
2451  else {
2452  b = a * (yw - f) / ywm;
2453  working_space[2 * shift + k] += b * c; //der
2454  b = a * a / ywm;
2455  working_space[3 * shift + k] += b * c; //temp
2456  }
2457  }
2458  k += 1;
2459  }
2460  if (fFixTxy == false) {
2461  a = Dertxy(fNPeaks, i1, i2,
2462  working_space, working_space[peak_vel],
2463  working_space[peak_vel + 1],
2464  working_space[peak_vel + 12],
2465  working_space[peak_vel + 13]);
2466  if (ywm != 0) {
2467  c = Ourpowl(a, pw);
2468  if (fStatisticType == kFitOptimChiFuncValues) {
2469  b = a * (yw * yw - f * f) / (ywm * ywm);
2470  working_space[2 * shift + k] += b * c; //der
2471  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2472  working_space[3 * shift + k] += b * c; //temp
2473  }
2474 
2475  else {
2476  b = a * (yw - f) / ywm;
2477  working_space[2 * shift + k] += b * c; //der
2478  b = a * a / ywm;
2479  working_space[3 * shift + k] += b * c; //temp
2480  }
2481  }
2482  k += 1;
2483  }
2484  if (fFixSxy == false) {
2485  a = Dersxy(fNPeaks, i1, i2,
2486  working_space, working_space[peak_vel],
2487  working_space[peak_vel + 1]);
2488  if (ywm != 0) {
2489  c = Ourpowl(a, pw);
2490  if (fStatisticType == kFitOptimChiFuncValues) {
2491  b = a * (yw * yw - f * f) / (ywm * ywm);
2492  working_space[2 * shift + k] += b * c; //der
2493  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2494  working_space[3 * shift + k] += b * c; //temp
2495  }
2496 
2497  else {
2498  b = a * (yw - f) / ywm;
2499  working_space[2 * shift + k] += b * c; //der
2500  b = a * a / ywm;
2501  working_space[3 * shift + k] += b * c; //temp
2502  }
2503  }
2504  k += 1;
2505  }
2506  if (fFixTx == false) {
2507  a = Dertx(fNPeaks, i1, working_space,
2508  working_space[peak_vel],
2509  working_space[peak_vel + 12]);
2510  if (ywm != 0) {
2511  c = Ourpowl(a, pw);
2512  if (fStatisticType == kFitOptimChiFuncValues) {
2513  b = a * (yw * yw - f * f) / (ywm * ywm);
2514  working_space[2 * shift + k] += b * c; //der
2515  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2516  working_space[3 * shift + k] += b * c; //temp
2517  }
2518 
2519  else {
2520  b = a * (yw - f) / ywm;
2521  working_space[2 * shift + k] += b * c; //der
2522  b = a * a / ywm;
2523  working_space[3 * shift + k] += b * c; //temp
2524  }
2525  }
2526  k += 1;
2527  }
2528  if (fFixTy == false) {
2529  a = Derty(fNPeaks, i2, working_space,
2530  working_space[peak_vel + 1],
2531  working_space[peak_vel + 13]);
2532  if (ywm != 0) {
2533  c = Ourpowl(a, pw);
2534  if (fStatisticType == kFitOptimChiFuncValues) {
2535  b = a * (yw * yw - f * f) / (ywm * ywm);
2536  working_space[2 * shift + k] += b * c; //der
2537  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2538  working_space[3 * shift + k] += b * c; //temp
2539  }
2540 
2541  else {
2542  b = a * (yw - f) / ywm;
2543  working_space[2 * shift + k] += b * c; //der
2544  b = a * a / ywm;
2545  working_space[3 * shift + k] += b * c; //temp
2546  }
2547  }
2548  k += 1;
2549  }
2550  if (fFixSx == false) {
2551  a = Dersx(fNPeaks, i1, working_space,
2552  working_space[peak_vel]);
2553  if (ywm != 0) {
2554  c = Ourpowl(a, pw);
2555  if (fStatisticType == kFitOptimChiFuncValues) {
2556  b = a * (yw * yw - f * f) / (ywm * ywm);
2557  working_space[2 * shift + k] += b * c; //der
2558  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2559  working_space[3 * shift + k] += b * c; //temp
2560  }
2561 
2562  else {
2563  b = a * (yw - f) / ywm;
2564  working_space[2 * shift + k] += b * c; //der
2565  b = a * a / ywm;
2566  working_space[3 * shift + k] += b * c; //temp
2567  }
2568  }
2569  k += 1;
2570  }
2571  if (fFixSy == false) {
2572  a = Dersy(fNPeaks, i2, working_space,
2573  working_space[peak_vel + 1]);
2574  if (ywm != 0) {
2575  c = Ourpowl(a, pw);
2576  if (fStatisticType == kFitOptimChiFuncValues) {
2577  b = a * (yw * yw - f * f) / (ywm * ywm);
2578  working_space[2 * shift + k] += b * c; //der
2579  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2580  working_space[3 * shift + k] += b * c; //temp
2581  }
2582 
2583  else {
2584  b = a * (yw - f) / ywm;
2585  working_space[2 * shift + k] += b * c; //der
2586  b = a * a / ywm;
2587  working_space[3 * shift + k] += b * c; //temp
2588  }
2589  }
2590  k += 1;
2591  }
2592  if (fFixBx == false) {
2593  a = Derbx(fNPeaks, i1, i2,
2594  working_space, working_space[peak_vel],
2595  working_space[peak_vel + 1],
2596  working_space[peak_vel + 6],
2597  working_space[peak_vel + 8],
2598  working_space[peak_vel + 12],
2599  working_space[peak_vel + 13]);
2600  if (ywm != 0) {
2601  c = Ourpowl(a, pw);
2602  if (fStatisticType == kFitOptimChiFuncValues) {
2603  b = a * (yw * yw - f * f) / (ywm * ywm);
2604  working_space[2 * shift + k] += b * c; //der
2605  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2606  working_space[3 * shift + k] += b * c; //temp
2607  }
2608 
2609  else {
2610  b = a * (yw - f) / ywm;
2611  working_space[2 * shift + k] += b * c; //der
2612  b = a * a / ywm;
2613  working_space[3 * shift + k] += b * c; //temp
2614  }
2615  }
2616  k += 1;
2617  }
2618  if (fFixBy == false) {
2619  a = Derby(fNPeaks, i1, i2,
2620  working_space, working_space[peak_vel],
2621  working_space[peak_vel + 1],
2622  working_space[peak_vel + 6],
2623  working_space[peak_vel + 8],
2624  working_space[peak_vel + 12],
2625  working_space[peak_vel + 13]);
2626  if (ywm != 0) {
2627  c = Ourpowl(a, pw);
2628  if (fStatisticType == kFitOptimChiFuncValues) {
2629  b = a * (yw * yw - f * f) / (ywm * ywm);
2630  working_space[2 * shift + k] += b * c; //der
2631  b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
2632  working_space[3 * shift + k] += b * c; //temp
2633  }
2634 
2635  else {
2636  b = a * (yw - f) / ywm;
2637  working_space[2 * shift + k] += b * c; //der
2638  b = a * a / ywm;
2639  working_space[3 * shift + k] += b * c; //temp
2640  }
2641  }
2642  k += 1;
2643  }
2644  }
2645  }
2646  for (j = 0; j < size; j++) {
2647  if (TMath::Abs(working_space[3 * shift + j]) > 0.000001)
2648  working_space[2 * shift + j] = working_space[2 * shift + j] / TMath::Abs(working_space[3 * shift + j]); //der[j]=der[j]/temp[j]
2649  else
2650  working_space[2 * shift + j] = 0; //der[j]
2651  }
2652 
2653  //calculate chi_opt
2654  chi2 = chi_opt;
2655  chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
2656 
2657  //calculate new parameters
2658  regul_cycle = 0;
2659  for (j = 0; j < size; j++) {
2660  working_space[4 * shift + j] = working_space[shift + j]; //temp_xk[j]=xk[j]
2661  }
2662 
2663  do {
2664  if (fAlphaOptim == kFitAlphaOptimal) {
2665  if (fStatisticType != kFitOptimMaxLikelihood)
2666  chi_min = 10000 * chi2;
2667 
2668  else
2669  chi_min = 0.1 * chi2;
2670  flag = 0;
2671  for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
2672  for (j = 0; j < size; j++) {
2673  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]
2674  }
2675  for (i = 0, j = 0; i < fNPeaks; i++) {
2676  if (fFixAmp[i] == false) {
2677  if (working_space[shift + j] < 0) //xk[j]
2678  working_space[shift + j] = 0; //xk[j]
2679  working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
2680  j += 1;
2681  }
2682  if (fFixPositionX[i] == false) {
2683  if (working_space[shift + j] < fXmin) //xk[j]
2684  working_space[shift + j] = fXmin; //xk[j]
2685  if (working_space[shift + j] > fXmax) //xk[j]
2686  working_space[shift + j] = fXmax; //xk[j]
2687  working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
2688  j += 1;
2689  }
2690  if (fFixPositionY[i] == false) {
2691  if (working_space[shift + j] < fYmin) //xk[j]
2692  working_space[shift + j] = fYmin; //xk[j]
2693  if (working_space[shift + j] > fYmax) //xk[j]
2694  working_space[shift + j] = fYmax; //xk[j]
2695  working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
2696  j += 1;
2697  }
2698  if (fFixAmpX1[i] == false) {
2699  if (working_space[shift + j] < 0) //xk[j]
2700  working_space[shift + j] = 0; //xk[j]
2701  working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
2702  j += 1;
2703  }
2704  if (fFixAmpY1[i] == false) {
2705  if (working_space[shift + j] < 0) //xk[j]
2706  working_space[shift + j] = 0; //xk[j]
2707  working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
2708  j += 1;
2709  }
2710  if (fFixPositionX1[i] == false) {
2711  if (working_space[shift + j] < fXmin) //xk[j]
2712  working_space[shift + j] = fXmin; //xk[j]
2713  if (working_space[shift + j] > fXmax) //xk[j]
2714  working_space[shift + j] = fXmax; //xk[j]
2715  working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
2716  j += 1;
2717  }
2718  if (fFixPositionY1[i] == false) {
2719  if (working_space[shift + j] < fYmin) //xk[j]
2720  working_space[shift + j] = fYmin; //xk[j]
2721  if (working_space[shift + j] > fYmax) //xk[j]
2722  working_space[shift + j] = fYmax; //xk[j]
2723  working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
2724  j += 1;
2725  }
2726  }
2727  if (fFixSigmaX == false) {
2728  if (working_space[shift + j] < 0.001) { //xk[j]
2729  working_space[shift + j] = 0.001; //xk[j]
2730  }
2731  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
2732  j += 1;
2733  }
2734  if (fFixSigmaY == false) {
2735  if (working_space[shift + j] < 0.001) { //xk[j]
2736  working_space[shift + j] = 0.001; //xk[j]
2737  }
2738  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
2739  j += 1;
2740  }
2741  if (fFixRo == false) {
2742  if (working_space[shift + j] < -1) { //xk[j]
2743  working_space[shift + j] = -1; //xk[j]
2744  }
2745  if (working_space[shift + j] > 1) { //xk[j]
2746  working_space[shift + j] = 1; //xk[j]
2747  }
2748  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
2749  j += 1;
2750  }
2751  if (fFixA0 == false) {
2752  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
2753  j += 1;
2754  }
2755  if (fFixAx == false) {
2756  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
2757  j += 1;
2758  }
2759  if (fFixAy == false) {
2760  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
2761  j += 1;
2762  }
2763  if (fFixTxy == false) {
2764  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
2765  j += 1;
2766  }
2767  if (fFixSxy == false) {
2768  working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
2769  j += 1;
2770  }
2771  if (fFixTx == false) {
2772  working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
2773  j += 1;
2774  }
2775  if (fFixTy == false) {
2776  working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
2777  j += 1;
2778  }
2779  if (fFixSx == false) {
2780  working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
2781  j += 1;
2782  }
2783  if (fFixSy == false) {
2784  working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
2785  j += 1;
2786  }
2787  if (fFixBx == false) {
2788  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2789  if (working_space[shift + j] < 0) //xk[j]
2790  working_space[shift + j] = -0.001; //xk[j]
2791  else
2792  working_space[shift + j] = 0.001; //xk[j]
2793  }
2794  working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
2795  j += 1;
2796  }
2797  if (fFixBy == 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 + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
2805  j += 1;
2806  }
2807  chi2 = 0;
2808  for (i1 = fXmin; i1 <= fXmax; i1++) {
2809  for (i2 = fYmin; i2 <= fYmax; i2++) {
2810  yw = source[i1][i2];
2811  ywm = yw;
2812  f = Shape2(fNPeaks, i1,
2813  i2, working_space,
2814  working_space[peak_vel],
2815  working_space[peak_vel + 1],
2816  working_space[peak_vel + 2],
2817  working_space[peak_vel + 3],
2818  working_space[peak_vel + 4],
2819  working_space[peak_vel + 5],
2820  working_space[peak_vel + 6],
2821  working_space[peak_vel + 7],
2822  working_space[peak_vel + 8],
2823  working_space[peak_vel + 9],
2824  working_space[peak_vel + 10],
2825  working_space[peak_vel + 11],
2826  working_space[peak_vel + 12],
2827  working_space[peak_vel + 13]);
2828  if (fStatisticType == kFitOptimChiFuncValues) {
2829  ywm = f;
2830  if (f < 0.00001)
2831  ywm = 0.00001;
2832  }
2833  if (fStatisticType == kFitOptimMaxLikelihood) {
2834  if (f > 0.00001)
2835  chi2 += yw * TMath::Log(f) - f;
2836  }
2837 
2838  else {
2839  if (ywm != 0)
2840  chi2 += (yw - f) * (yw - f) / ywm;
2841  }
2842  }
2843  }
2844  if ((chi2 < chi_min
2845  && fStatisticType != kFitOptimMaxLikelihood)
2846  || (chi2 > chi_min
2847  && fStatisticType == kFitOptimMaxLikelihood)) {
2848  pmin = pi, chi_min = chi2;
2849  }
2850 
2851  else
2852  flag = 1;
2853  if (pi == 0.1)
2854  chi_min = chi2;
2855  chi = chi_min;
2856  }
2857  if (pmin != 0.1) {
2858  for (j = 0; j < size; j++) {
2859  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]
2860  }
2861  for (i = 0, j = 0; i < fNPeaks; i++) {
2862  if (fFixAmp[i] == false) {
2863  if (working_space[shift + j] < 0) //xk[j]
2864  working_space[shift + j] = 0; //xk[j]
2865  working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
2866  j += 1;
2867  }
2868  if (fFixPositionX[i] == false) {
2869  if (working_space[shift + j] < fXmin) //xk[j]
2870  working_space[shift + j] = fXmin; //xk[j]
2871  if (working_space[shift + j] > fXmax) //xk[j]
2872  working_space[shift + j] = fXmax; //xk[j]
2873  working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
2874  j += 1;
2875  }
2876  if (fFixPositionY[i] == false) {
2877  if (working_space[shift + j] < fYmin) //xk[j]
2878  working_space[shift + j] = fYmin; //xk[j]
2879  if (working_space[shift + j] > fYmax) //xk[j]
2880  working_space[shift + j] = fYmax; //xk[j]
2881  working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
2882  j += 1;
2883  }
2884  if (fFixAmpX1[i] == false) {
2885  if (working_space[shift + j] < 0) //xk[j]
2886  working_space[shift + j] = 0; //xk[j]
2887  working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
2888  j += 1;
2889  }
2890  if (fFixAmpY1[i] == false) {
2891  if (working_space[shift + j] < 0) //xk[j]
2892  working_space[shift + j] = 0; //xk[j]
2893  working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
2894  j += 1;
2895  }
2896  if (fFixPositionX1[i] == false) {
2897  if (working_space[shift + j] < fXmin) //xk[j]
2898  working_space[shift + j] = fXmin; //xk[j]
2899  if (working_space[shift + j] > fXmax) //xk[j]
2900  working_space[shift + j] = fXmax; //xk[j]
2901  working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
2902  j += 1;
2903  }
2904  if (fFixPositionY1[i] == false) {
2905  if (working_space[shift + j] < fYmin) //xk[j]
2906  working_space[shift + j] = fYmin; //xk[j]
2907  if (working_space[shift + j] > fYmax) //xk[j]
2908  working_space[shift + j] = fYmax; //xk[j]
2909  working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
2910  j += 1;
2911  }
2912  }
2913  if (fFixSigmaX == false) {
2914  if (working_space[shift + j] < 0.001) { //xk[j]
2915  working_space[shift + j] = 0.001; //xk[j]
2916  }
2917  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
2918  j += 1;
2919  }
2920  if (fFixSigmaY == false) {
2921  if (working_space[shift + j] < 0.001) { //xk[j]
2922  working_space[shift + j] = 0.001; //xk[j]
2923  }
2924  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
2925  j += 1;
2926  }
2927  if (fFixRo == false) {
2928  if (working_space[shift + j] < -1) { //xk[j]
2929  working_space[shift + j] = -1; //xk[j]
2930  }
2931  if (working_space[shift + j] > 1) { //xk[j]
2932  working_space[shift + j] = 1; //xk[j]
2933  }
2934  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
2935  j += 1;
2936  }
2937  if (fFixA0 == false) {
2938  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
2939  j += 1;
2940  }
2941  if (fFixAx == false) {
2942  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
2943  j += 1;
2944  }
2945  if (fFixAy == false) {
2946  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
2947  j += 1;
2948  }
2949  if (fFixTxy == false) {
2950  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
2951  j += 1;
2952  }
2953  if (fFixSxy == false) {
2954  working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
2955  j += 1;
2956  }
2957  if (fFixTx == false) {
2958  working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
2959  j += 1;
2960  }
2961  if (fFixTy == false) {
2962  working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
2963  j += 1;
2964  }
2965  if (fFixSx == false) {
2966  working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
2967  j += 1;
2968  }
2969  if (fFixSy == false) {
2970  working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
2971  j += 1;
2972  }
2973  if (fFixBx == false) {
2974  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2975  if (working_space[shift + j] < 0) //xk[j]
2976  working_space[shift + j] = -0.001; //xk[j]
2977  else
2978  working_space[shift + j] = 0.001; //xk[j]
2979  }
2980  working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
2981  j += 1;
2982  }
2983  if (fFixBy == false) {
2984  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
2985  if (working_space[shift + j] < 0) //xk[j]
2986  working_space[shift + j] = -0.001; //xk[j]
2987  else
2988  working_space[shift + j] = 0.001; //xk[j]
2989  }
2990  working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
2991  j += 1;
2992  }
2993  chi = chi_min;
2994  }
2995  }
2996 
2997  else {
2998  for (j = 0; j < size; j++) {
2999  working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
3000  }
3001  for (i = 0, j = 0; i < fNPeaks; i++) {
3002  if (fFixAmp[i] == false) {
3003  if (working_space[shift + j] < 0) //xk[j]
3004  working_space[shift + j] = 0; //xk[j]
3005  working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
3006  j += 1;
3007  }
3008  if (fFixPositionX[i] == false) {
3009  if (working_space[shift + j] < fXmin) //xk[j]
3010  working_space[shift + j] = fXmin; //xk[j]
3011  if (working_space[shift + j] > fXmax) //xk[j]
3012  working_space[shift + j] = fXmax; //xk[j]
3013  working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
3014  j += 1;
3015  }
3016  if (fFixPositionY[i] == false) {
3017  if (working_space[shift + j] < fYmin) //xk[j]
3018  working_space[shift + j] = fYmin; //xk[j]
3019  if (working_space[shift + j] > fYmax) //xk[j]
3020  working_space[shift + j] = fYmax; //xk[j]
3021  working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
3022  j += 1;
3023  }
3024  if (fFixAmpX1[i] == false) {
3025  if (working_space[shift + j] < 0) //xk[j]
3026  working_space[shift + j] = 0; //xk[j]
3027  working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
3028  j += 1;
3029  }
3030  if (fFixAmpY1[i] == false) {
3031  if (working_space[shift + j] < 0) //xk[j]
3032  working_space[shift + j] = 0; //xk[j]
3033  working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
3034  j += 1;
3035  }
3036  if (fFixPositionX1[i] == false) {
3037  if (working_space[shift + j] < fXmin) //xk[j]
3038  working_space[shift + j] = fXmin; //xk[j]
3039  if (working_space[shift + j] > fXmax) //xk[j]
3040  working_space[shift + j] = fXmax; //xk[j]
3041  working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
3042  j += 1;
3043  }
3044  if (fFixPositionY1[i] == false) {
3045  if (working_space[shift + j] < fYmin) //xk[j]
3046  working_space[shift + j] = fYmin; //xk[j]
3047  if (working_space[shift + j] > fYmax) //xk[j]
3048  working_space[shift + j] = fYmax; //xk[j]
3049  working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
3050  j += 1;
3051  }
3052  }
3053  if (fFixSigmaX == false) {
3054  if (working_space[shift + j] < 0.001) { //xk[j]
3055  working_space[shift + j] = 0.001; //xk[j]
3056  }
3057  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
3058  j += 1;
3059  }
3060  if (fFixSigmaY == false) {
3061  if (working_space[shift + j] < 0.001) { //xk[j]
3062  working_space[shift + j] = 0.001; //xk[j]
3063  }
3064  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
3065  j += 1;
3066  }
3067  if (fFixRo == false) {
3068  if (working_space[shift + j] < -1) { //xk[j]
3069  working_space[shift + j] = -1; //xk[j]
3070  }
3071  if (working_space[shift + j] > 1) { //xk[j]
3072  working_space[shift + j] = 1; //xk[j]
3073  }
3074  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
3075  j += 1;
3076  }
3077  if (fFixA0 == false) {
3078  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
3079  j += 1;
3080  }
3081  if (fFixAx == false) {
3082  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
3083  j += 1;
3084  }
3085  if (fFixAy == false) {
3086  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
3087  j += 1;
3088  }
3089  if (fFixTxy == false) {
3090  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
3091  j += 1;
3092  }
3093  if (fFixSxy == false) {
3094  working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
3095  j += 1;
3096  }
3097  if (fFixTx == false) {
3098  working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
3099  j += 1;
3100  }
3101  if (fFixTy == false) {
3102  working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
3103  j += 1;
3104  }
3105  if (fFixSx == false) {
3106  working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
3107  j += 1;
3108  }
3109  if (fFixSy == false) {
3110  working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
3111  j += 1;
3112  }
3113  if (fFixBx == false) {
3114  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
3115  if (working_space[shift + j] < 0) //xk[j]
3116  working_space[shift + j] = -0.001; //xk[j]
3117  else
3118  working_space[shift + j] = 0.001; //xk[j]
3119  }
3120  working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
3121  j += 1;
3122  }
3123  if (fFixBy == false) {
3124  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
3125  if (working_space[shift + j] < 0) //xk[j]
3126  working_space[shift + j] = -0.001; //xk[j]
3127  else
3128  working_space[shift + j] = 0.001; //xk[j]
3129  }
3130  working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
3131  j += 1;
3132  }
3133  chi = 0;
3134  for (i1 = fXmin; i1 <= fXmax; i1++) {
3135  for (i2 = fYmin; i2 <= fYmax; i2++) {
3136  yw = source[i1][i2];
3137  ywm = yw;
3138  f = Shape2(fNPeaks, i1, i2,
3139  working_space, working_space[peak_vel],
3140  working_space[peak_vel + 1],
3141  working_space[peak_vel + 2],
3142  working_space[peak_vel + 3],
3143  working_space[peak_vel + 4],
3144  working_space[peak_vel + 5],
3145  working_space[peak_vel + 6],
3146  working_space[peak_vel + 7],
3147  working_space[peak_vel + 8],
3148  working_space[peak_vel + 9],
3149  working_space[peak_vel + 10],
3150  working_space[peak_vel + 11],
3151  working_space[peak_vel + 12],
3152  working_space[peak_vel + 13]);
3153  if (fStatisticType == kFitOptimChiFuncValues) {
3154  ywm = f;
3155  if (f < 0.00001)
3156  ywm = 0.00001;
3157  }
3158  if (fStatisticType == kFitOptimMaxLikelihood) {
3159  if (f > 0.00001)
3160  chi += yw * TMath::Log(f) - f;
3161  }
3162 
3163  else {
3164  if (ywm != 0)
3165  chi += (yw - f) * (yw - f) / ywm;
3166  }
3167  }
3168  }
3169  }
3170  chi2 = chi;
3171  chi = TMath::Sqrt(TMath::Abs(chi));
3172  if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
3173  alpha = alpha * chi_opt / (2 * chi);
3174 
3175  else if (fAlphaOptim == kFitAlphaOptimal)
3176  alpha = alpha / 10.0;
3177  iter += 1;
3178  regul_cycle += 1;
3179  } while (((chi > chi_opt
3180  && fStatisticType != kFitOptimMaxLikelihood)
3181  || (chi < chi_opt
3182  && fStatisticType == kFitOptimMaxLikelihood))
3183  && regul_cycle < kFitNumRegulCycles);
3184  for (j = 0; j < size; j++) {
3185  working_space[4 * shift + j] = 0; //temp_xk[j]
3186  working_space[2 * shift + j] = 0; //der[j]
3187  }
3188  for (i1 = fXmin, chi_cel = 0; i1 <= fXmax; i1++) {
3189  for (i2 = fYmin; i2 <= fYmax; i2++) {
3190  yw = source[i1][i2];
3191  if (yw == 0)
3192  yw = 1;
3193  f = Shape2(fNPeaks, i1, i2,
3194  working_space, working_space[peak_vel],
3195  working_space[peak_vel + 1],
3196  working_space[peak_vel + 2],
3197  working_space[peak_vel + 3],
3198  working_space[peak_vel + 4],
3199  working_space[peak_vel + 5],
3200  working_space[peak_vel + 6],
3201  working_space[peak_vel + 7],
3202  working_space[peak_vel + 8],
3203  working_space[peak_vel + 9],
3204  working_space[peak_vel + 10],
3205  working_space[peak_vel + 11],
3206  working_space[peak_vel + 12],
3207  working_space[peak_vel + 13]);
3208  chi_opt = (yw - f) * (yw - f) / yw;
3209  chi_cel += (yw - f) * (yw - f) / yw;
3210 
3211  //calculate gradient vector
3212  for (j = 0, k = 0; j < fNPeaks; j++) {
3213  if (fFixAmp[j] == false) {
3214  a = Deramp2(i1, i2,
3215  working_space[7 * j + 1],
3216  working_space[7 * j + 2],
3217  working_space[peak_vel],
3218  working_space[peak_vel + 1],
3219  working_space[peak_vel + 2],
3220  working_space[peak_vel + 6],
3221  working_space[peak_vel + 7],
3222  working_space[peak_vel + 12],
3223  working_space[peak_vel + 13]);
3224  if (yw != 0) {
3225  c = Ourpowl(a, pw);
3226  working_space[2 * shift + k] += chi_opt * c; //der[k]
3227  b = a * a / yw;
3228  working_space[4 * shift + k] += b * c; //temp_xk[k]
3229  }
3230  k += 1;
3231  }
3232  if (fFixPositionX[j] == false) {
3233  a = Deri02(i1, i2,
3234  working_space[7 * j],
3235  working_space[7 * j + 1],
3236  working_space[7 * j + 2],
3237  working_space[peak_vel],
3238  working_space[peak_vel + 1],
3239  working_space[peak_vel + 2],
3240  working_space[peak_vel + 6],
3241  working_space[peak_vel + 7],
3242  working_space[peak_vel + 12],
3243  working_space[peak_vel + 13]);
3244  if (yw != 0) {
3245  c = Ourpowl(a, pw);
3246  working_space[2 * shift + k] += chi_opt * c; //der[k]
3247  b = a * a / yw;
3248  working_space[4 * shift + k] += b * c; //temp_xk[k]
3249  }
3250  k += 1;
3251  }
3252  if (fFixPositionY[j] == false) {
3253  a = Derj02(i1, i2,
3254  working_space[7 * j],
3255  working_space[7 * j + 1],
3256  working_space[7 * j + 2],
3257  working_space[peak_vel],
3258  working_space[peak_vel + 1],
3259  working_space[peak_vel + 2],
3260  working_space[peak_vel + 6],
3261  working_space[peak_vel + 7],
3262  working_space[peak_vel + 12],
3263  working_space[peak_vel + 13]);
3264  if (yw != 0) {
3265  c = Ourpowl(a, pw);
3266  working_space[2 * shift + k] += chi_opt * c; //der[k]
3267  b = a * a / yw;
3268  working_space[4 * shift + k] += b * c; //temp_xk[k]
3269  }
3270  k += 1;
3271  }
3272  if (fFixAmpX1[j] == false) {
3273  a = Derampx(i1, working_space[7 * j + 5],
3274  working_space[peak_vel],
3275  working_space[peak_vel + 8],
3276  working_space[peak_vel + 10],
3277  working_space[peak_vel + 12]);
3278  if (yw != 0) {
3279  c = Ourpowl(a, pw);
3280  working_space[2 * shift + k] += chi_opt * c; //der[k]
3281  b = a * a / yw;
3282  working_space[4 * shift + k] += b * c; //temp_xk[k]
3283  }
3284  k += 1;
3285  }
3286  if (fFixAmpY1[j] == false) {
3287  a = Derampx(i2, working_space[7 * j + 6],
3288  working_space[peak_vel + 1],
3289  working_space[peak_vel + 9],
3290  working_space[peak_vel + 11],
3291  working_space[peak_vel + 13]);
3292  if (yw != 0) {
3293  c = Ourpowl(a, pw);
3294  working_space[2 * shift + k] += chi_opt * c; //der[k]
3295  b = a * a / yw;
3296  working_space[4 * shift + k] += b * c; //temp_xk[k]
3297  }
3298  k += 1;
3299  }
3300  if (fFixPositionX1[j] == false) {
3301  a = Deri01(i1, working_space[7 * j + 3],
3302  working_space[7 * j + 5],
3303  working_space[peak_vel],
3304  working_space[peak_vel + 8],
3305  working_space[peak_vel + 10],
3306  working_space[peak_vel + 12]);
3307  if (yw != 0) {
3308  c = Ourpowl(a, pw);
3309  working_space[2 * shift + k] += chi_opt * c; //der[k]
3310  b = a * a / yw;
3311  working_space[4 * shift + k] += b * c; //temp_xk[k]
3312  }
3313  k += 1;
3314  }
3315  if (fFixPositionY1[j] == false) {
3316  a = Deri01(i2, working_space[7 * j + 4],
3317  working_space[7 * j + 6],
3318  working_space[peak_vel + 1],
3319  working_space[peak_vel + 9],
3320  working_space[peak_vel + 11],
3321  working_space[peak_vel + 13]);
3322  if (yw != 0) {
3323  c = Ourpowl(a, pw);
3324  working_space[2 * shift + k] += chi_opt * c; //der[k]
3325  b = a * a / yw;
3326  working_space[4 * shift + k] += b * c; //temp_xk[k]
3327  }
3328  k += 1;
3329  }
3330  }
3331  if (fFixSigmaX == false) {
3332  a = Dersigmax(fNPeaks, i1, i2,
3333  working_space, working_space[peak_vel],
3334  working_space[peak_vel + 1],
3335  working_space[peak_vel + 2],
3336  working_space[peak_vel + 6],
3337  working_space[peak_vel + 7],
3338  working_space[peak_vel + 8],
3339  working_space[peak_vel + 10],
3340  working_space[peak_vel + 12],
3341  working_space[peak_vel + 13]);
3342  if (yw != 0) {
3343  c = Ourpowl(a, pw);
3344  working_space[2 * shift + k] += chi_opt * c; //der[k]
3345  b = a * a / yw;
3346  working_space[4 * shift + k] += b * c; //temp_xk[k]
3347  }
3348  k += 1;
3349  }
3350  if (fFixSigmaY == false) {
3351  a = Dersigmay(fNPeaks, i1, i2,
3352  working_space, working_space[peak_vel],
3353  working_space[peak_vel + 1],
3354  working_space[peak_vel + 2],
3355  working_space[peak_vel + 6],
3356  working_space[peak_vel + 7],
3357  working_space[peak_vel + 9],
3358  working_space[peak_vel + 11],
3359  working_space[peak_vel + 12],
3360  working_space[peak_vel + 13]);
3361  if (yw != 0) {
3362  c = Ourpowl(a, pw);
3363  working_space[2 * shift + k] += chi_opt * c; //der[k]
3364  b = a * a / yw;
3365  working_space[4 * shift + k] += b * c; //temp_xk[k]
3366  }
3367  k += 1;
3368  }
3369  if (fFixRo == false) {
3370  a = Derro(fNPeaks, i1, i2,
3371  working_space, working_space[peak_vel],
3372  working_space[peak_vel + 1],
3373  working_space[peak_vel + 2]);
3374  if (yw != 0) {
3375  c = Ourpowl(a, pw);
3376  working_space[2 * shift + k] += chi_opt * c; //der[k]
3377  b = a * a / yw;
3378  working_space[4 * shift + k] += b * c; //temp_xk[k]
3379  }
3380  k += 1;
3381  }
3382  if (fFixA0 == false) {
3383  a = 1.;
3384  if (yw != 0) {
3385  c = Ourpowl(a, pw);
3386  working_space[2 * shift + k] += chi_opt * c; //der[k]
3387  b = a * a / yw;
3388  working_space[4 * shift + k] += b * c; //temp_xk[k]
3389  }
3390  k += 1;
3391  }
3392  if (fFixAx == false) {
3393  a = i1;
3394  if (yw != 0) {
3395  c = Ourpowl(a, pw);
3396  working_space[2 * shift + k] += chi_opt * c; //der[k]
3397  b = a * a / yw;
3398  working_space[4 * shift + k] += b * c; //temp_xk[k]
3399  }
3400  k += 1;
3401  }
3402  if (fFixAy == false) {
3403  a = i2;
3404  if (yw != 0) {
3405  c = Ourpowl(a, pw);
3406  working_space[2 * shift + k] += chi_opt * c; //der[k]
3407  b = a * a / yw;
3408  working_space[4 * shift + k] += b * c; //temp_xk[k]
3409  }
3410  k += 1;
3411  }
3412  if (fFixTxy == false) {
3413  a = Dertxy(fNPeaks, i1, i2,
3414  working_space, working_space[peak_vel],
3415  working_space[peak_vel + 1],
3416  working_space[peak_vel + 12],
3417  working_space[peak_vel + 13]);
3418  if (yw != 0) {
3419  c = Ourpowl(a, pw);
3420  working_space[2 * shift + k] += chi_opt * c; //der[k]
3421  b = a * a / yw;
3422  working_space[4 * shift + k] += b * c; //temp_xk[k]
3423  }
3424  k += 1;
3425  }
3426  if (fFixSxy == false) {
3427  a = Dersxy(fNPeaks, i1, i2,
3428  working_space, working_space[peak_vel],
3429  working_space[peak_vel + 1]);
3430  if (yw != 0) {
3431  c = Ourpowl(a, pw);
3432  working_space[2 * shift + k] += chi_opt * c; //der[k]
3433  b = a * a / yw;
3434  working_space[4 * shift + k] += b * c; //temp_xk[k]
3435  }
3436  k += 1;
3437  }
3438  if (fFixTx == false) {
3439  a = Dertx(fNPeaks, i1, working_space,
3440  working_space[peak_vel],
3441  working_space[peak_vel + 12]);
3442  if (yw != 0) {
3443  c = Ourpowl(a, pw);
3444  working_space[2 * shift + k] += chi_opt * c; //der[k]
3445  b = a * a / yw;
3446  working_space[4 * shift + k] += b * c; //temp_xk[k]
3447  }
3448  k += 1;
3449  }
3450  if (fFixTy == false) {
3451  a = Derty(fNPeaks, i2, working_space,
3452  working_space[peak_vel + 1],
3453  working_space[peak_vel + 13]);
3454  if (yw != 0) {
3455  c = Ourpowl(a, pw);
3456  working_space[2 * shift + k] += chi_opt * c; //der[k]
3457  b = a * a / yw;
3458  working_space[4 * shift + k] += b * c; //temp_xk[k]
3459  }
3460  k += 1;
3461  }
3462  if (fFixSx == false) {
3463  a = Dersx(fNPeaks, i1, working_space,
3464  working_space[peak_vel]);
3465  if (yw != 0) {
3466  c = Ourpowl(a, pw);
3467  working_space[2 * shift + k] += chi_opt * c; //der[k]
3468  b = a * a / yw;
3469  working_space[4 * shift + k] += b * c; //temp_xk[k]
3470  }
3471  k += 1;
3472  }
3473  if (fFixSy == false) {
3474  a = Dersy(fNPeaks, i2, working_space,
3475  working_space[peak_vel + 1]);
3476  if (yw != 0) {
3477  c = Ourpowl(a, pw);
3478  working_space[2 * shift + k] += chi_opt * c; //der[k]
3479  b = a * a / yw;
3480  working_space[4 * shift + k] += b * c; //temp_xk[k]
3481  }
3482  k += 1;
3483  }
3484  if (fFixBx == false) {
3485  a = Derbx(fNPeaks, i1, i2,
3486  working_space, working_space[peak_vel],
3487  working_space[peak_vel + 1],
3488  working_space[peak_vel + 6],
3489  working_space[peak_vel + 8],
3490  working_space[peak_vel + 12],
3491  working_space[peak_vel + 13]);
3492  if (yw != 0) {
3493  c = Ourpowl(a, pw);
3494  working_space[2 * shift + k] += chi_opt * c; //der[k]
3495  b = a * a / yw;
3496  working_space[4 * shift + k] += b * c; //temp_xk[k]
3497  }
3498  k += 1;
3499  }
3500  if (fFixBy == false) {
3501  a = Derby(fNPeaks, i1, i2,
3502  working_space, working_space[peak_vel],
3503  working_space[peak_vel + 1],
3504  working_space[peak_vel + 6],
3505  working_space[peak_vel + 8],
3506  working_space[peak_vel + 12],
3507  working_space[peak_vel + 13]);
3508  if (yw != 0) {
3509  c = Ourpowl(a, pw);
3510  working_space[2 * shift + k] += chi_opt * c; //der[k]
3511  b = a * a / yw;
3512  working_space[4 * shift + k] += b * c; //temp_xk[k]
3513  }
3514  k += 1;
3515  }
3516  }
3517  }
3518  }
3519  b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
3520  chi_er = chi_cel / b;
3521  for (i = 0, j = 0; i < fNPeaks; i++) {
3522  fVolume[i] =
3523  Volume(working_space[7 * i], working_space[peak_vel],
3524  working_space[peak_vel + 1], working_space[peak_vel + 2]);
3525  if (fVolume[i] > 0) {
3526  c = 0;
3527  if (fFixAmp[i] == false) {
3528  a = Derpa2(working_space[peak_vel],
3529  working_space[peak_vel + 1],
3530  working_space[peak_vel + 2]);
3531  b = working_space[4 * shift + j]; //temp_xk[j]
3532  if (b == 0)
3533  b = 1;
3534 
3535  else
3536  b = 1 / b;
3537  c = c + a * a * b;
3538  }
3539  if (fFixSigmaX == false) {
3540  a = Derpsigmax(working_space[shift + j],
3541  working_space[peak_vel + 1],
3542  working_space[peak_vel + 2]);
3543  b = working_space[4 * shift + peak_vel]; //temp_xk[j]
3544  if (b == 0)
3545  b = 1;
3546 
3547  else
3548  b = 1 / b;
3549  c = c + a * a * b;
3550  }
3551  if (fFixSigmaY == false) {
3552  a = Derpsigmay(working_space[shift + j],
3553  working_space[peak_vel],
3554  working_space[peak_vel + 2]);
3555  b = working_space[4 * shift + peak_vel + 1]; //temp_xk[j]
3556  if (b == 0)
3557  b = 1;
3558 
3559  else
3560  b = 1 / b;
3561  c = c + a * a * b;
3562  }
3563  if (fFixRo == false) {
3564  a = Derpro(working_space[shift + j], working_space[peak_vel],
3565  working_space[peak_vel + 1],
3566  working_space[peak_vel + 2]);
3567  b = working_space[4 * shift + peak_vel + 2]; //temp_xk[j]
3568  if (b == 0)
3569  b = 1;
3570 
3571  else
3572  b = 1 / b;
3573  c = c + a * a * b;
3574  }
3575  fVolumeErr[i] = TMath::Sqrt(TMath::Abs(chi_er * c));
3576  }
3577 
3578  else {
3579  fVolumeErr[i] = 0;
3580  }
3581  if (fFixAmp[i] == false) {
3582  fAmpCalc[i] = working_space[shift + j]; //xk[j]
3583  if (working_space[3 * shift + j] != 0)
3584  fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3585  j += 1;
3586  }
3587 
3588  else {
3589  fAmpCalc[i] = fAmpInit[i];
3590  fAmpErr[i] = 0;
3591  }
3592  if (fFixPositionX[i] == false) {
3593  fPositionCalcX[i] = working_space[shift + j]; //xk[j]
3594  if (working_space[3 * shift + j] != 0)
3595  fPositionErrX[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3596  j += 1;
3597  }
3598 
3599  else {
3600  fPositionCalcX[i] = fPositionInitX[i];
3601  fPositionErrX[i] = 0;
3602  }
3603  if (fFixPositionY[i] == false) {
3604  fPositionCalcY[i] = working_space[shift + j]; //xk[j]
3605  if (working_space[3 * shift + j] != 0)
3606  fPositionErrY[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3607  j += 1;
3608  }
3609 
3610  else {
3611  fPositionCalcY[i] = fPositionInitY[i];
3612  fPositionErrY[i] = 0;
3613  }
3614  if (fFixAmpX1[i] == false) {
3615  fAmpCalcX1[i] = working_space[shift + j]; //xk[j]
3616  if (working_space[3 * shift + j] != 0)
3617  fAmpErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3618  j += 1;
3619  }
3620 
3621  else {
3622  fAmpCalcX1[i] = fAmpInitX1[i];
3623  fAmpErrX1[i] = 0;
3624  }
3625  if (fFixAmpY1[i] == false) {
3626  fAmpCalcY1[i] = working_space[shift + j]; //xk[j]
3627  if (working_space[3 * shift + j] != 0)
3628  fAmpErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3629  j += 1;
3630  }
3631 
3632  else {
3633  fAmpCalcY1[i] = fAmpInitY1[i];
3634  fAmpErrY1[i] = 0;
3635  }
3636  if (fFixPositionX1[i] == false) {
3637  fPositionCalcX1[i] = working_space[shift + j]; //xk[j]
3638  if (working_space[3 * shift + j] != 0)
3639  fPositionErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3640  j += 1;
3641  }
3642 
3643  else {
3644  fPositionCalcX1[i] = fPositionInitX1[i];
3645  fPositionErrX1[i] = 0;
3646  }
3647  if (fFixPositionY1[i] == false) {
3648  fPositionCalcY1[i] = working_space[shift + j]; //xk[j]
3649  if (working_space[3 * shift + j] != 0)
3650  fPositionErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3651  j += 1;
3652  }
3653 
3654  else {
3655  fPositionCalcY1[i] = fPositionInitY1[i];
3656  fPositionErrY1[i] = 0;
3657  }
3658  }
3659  if (fFixSigmaX == false) {
3660  fSigmaCalcX = working_space[shift + j]; //xk[j]
3661  if (working_space[3 * shift + j] != 0) //temp[j]
3662  fSigmaErrX = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3663  j += 1;
3664  }
3665 
3666  else {
3667  fSigmaCalcX = fSigmaInitX;
3668  fSigmaErrX = 0;
3669  }
3670  if (fFixSigmaY == false) {
3671  fSigmaCalcY = working_space[shift + j]; //xk[j]
3672  if (working_space[3 * shift + j] != 0) //temp[j]
3673  fSigmaErrY = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3674  j += 1;
3675  }
3676 
3677  else {
3678  fSigmaCalcY = fSigmaInitY;
3679  fSigmaErrY = 0;
3680  }
3681  if (fFixRo == false) {
3682  fRoCalc = working_space[shift + j]; //xk[j]
3683  if (working_space[3 * shift + j] != 0) //temp[j]
3684  fRoErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3685  j += 1;
3686  }
3687 
3688  else {
3689  fRoCalc = fRoInit;
3690  fRoErr = 0;
3691  }
3692  if (fFixA0 == false) {
3693  fA0Calc = working_space[shift + j]; //xk[j]
3694  if (working_space[3 * shift + j] != 0) //temp[j]
3695  fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3696  j += 1;
3697  }
3698 
3699  else {
3700  fA0Calc = fA0Init;
3701  fA0Err = 0;
3702  }
3703  if (fFixAx == false) {
3704  fAxCalc = working_space[shift + j]; //xk[j]
3705  if (working_space[3 * shift + j] != 0) //temp[j]
3706  fAxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3707  j += 1;
3708  }
3709 
3710  else {
3711  fAxCalc = fAxInit;
3712  fAxErr = 0;
3713  }
3714  if (fFixAy == false) {
3715  fAyCalc = working_space[shift + j]; //xk[j]
3716  if (working_space[3 * shift + j] != 0) //temp[j]
3717  fAyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3718  j += 1;
3719  }
3720 
3721  else {
3722  fAyCalc = fAyInit;
3723  fAyErr = 0;
3724  }
3725  if (fFixTxy == false) {
3726  fTxyCalc = working_space[shift + j]; //xk[j]
3727  if (working_space[3 * shift + j] != 0) //temp[j]
3728  fTxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3729  j += 1;
3730  }
3731 
3732  else {
3733  fTxyCalc = fTxyInit;
3734  fTxyErr = 0;
3735  }
3736  if (fFixSxy == false) {
3737  fSxyCalc = working_space[shift + j]; //xk[j]
3738  if (working_space[3 * shift + j] != 0) //temp[j]
3739  fSxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3740  j += 1;
3741  }
3742 
3743  else {
3744  fSxyCalc = fSxyInit;
3745  fSxyErr = 0;
3746  }
3747  if (fFixTx == false) {
3748  fTxCalc = working_space[shift + j]; //xk[j]
3749  if (working_space[3 * shift + j] != 0) //temp[j]
3750  fTxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3751  j += 1;
3752  }
3753 
3754  else {
3755  fTxCalc = fTxInit;
3756  fTxErr = 0;
3757  }
3758  if (fFixTy == false) {
3759  fTyCalc = working_space[shift + j]; //xk[j]
3760  if (working_space[3 * shift + j] != 0) //temp[j]
3761  fTyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3762  j += 1;
3763  }
3764 
3765  else {
3766  fTyCalc = fTyInit;
3767  fTyErr = 0;
3768  }
3769  if (fFixSx == false) {
3770  fSxCalc = working_space[shift + j]; //xk[j]
3771  if (working_space[3 * shift + j] != 0) //temp[j]
3772  fSxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3773  j += 1;
3774  }
3775 
3776  else {
3777  fSxCalc = fSxInit;
3778  fSxErr = 0;
3779  }
3780  if (fFixSy == false) {
3781  fSyCalc = working_space[shift + j]; //xk[j]
3782  if (working_space[3 * shift + j] != 0) //temp[j]
3783  fSyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3784  j += 1;
3785  }
3786 
3787  else {
3788  fSyCalc = fSyInit;
3789  fSyErr = 0;
3790  }
3791  if (fFixBx == false) {
3792  fBxCalc = working_space[shift + j]; //xk[j]
3793  if (working_space[3 * shift + j] != 0) //temp[j]
3794  fBxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3795  j += 1;
3796  }
3797 
3798  else {
3799  fBxCalc = fBxInit;
3800  fBxErr = 0;
3801  }
3802  if (fFixBy == false) {
3803  fByCalc = working_space[shift + j]; //xk[j]
3804  if (working_space[3 * shift + j] != 0) //temp[j]
3805  fByErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
3806  j += 1;
3807  }
3808 
3809  else {
3810  fByCalc = fByInit;
3811  fByErr = 0;
3812  }
3813  b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
3814  fChi = chi_cel / b;
3815  for (i1 = fXmin; i1 <= fXmax; i1++) {
3816  for (i2 = fYmin; i2 <= fYmax; i2++) {
3817  f = Shape2(fNPeaks, i1, i2,
3818  working_space, working_space[peak_vel],
3819  working_space[peak_vel + 1],
3820  working_space[peak_vel + 2],
3821  working_space[peak_vel + 3],
3822  working_space[peak_vel + 4],
3823  working_space[peak_vel + 5],
3824  working_space[peak_vel + 6],
3825  working_space[peak_vel + 7],
3826  working_space[peak_vel + 8],
3827  working_space[peak_vel + 9],
3828  working_space[peak_vel + 10],
3829  working_space[peak_vel + 11],
3830  working_space[peak_vel + 12],
3831  working_space[peak_vel + 13]);
3832  source[i1][i2] = f;
3833  }
3834  }
3835  delete [] working_space;
3836  return;
3837 }
3838 
3839 
3840 
3841 ////////////////////////////////////////////////////////////////////////////////
3842 /// This function fits the source spectrum. The calling program should
3843 /// fill in input parameters of the TSpectrum2Fit class.
3844 /// The fitted parameters are written into
3845 /// TSpectrum2Fit class output parameters and fitted data are written into
3846 /// source spectrum.
3847 ///
3848 /// Function parameters:
3849 /// - source-pointer to the matrix of source spectrum
3850 ///
3851 /// ### Stiefel fitting algorithm
3852 ///
3853 /// This function fits the source
3854 /// spectrum using Stiefel-Hestens method [1]. The calling program should fill in
3855 /// input fitting parameters of the TSpectrumFit2 class using a set of
3856 /// TSpectrumFit2 setters. The fitted parameters are written into the class and the
3857 /// fitted data are written into source spectrum. It converges faster than Awmi
3858 /// method.
3859 ///
3860 /// #### Reference:
3861 ///
3862 /// [1] B. Mihaila: Analysis of
3863 /// complex gamma spectra, Rom. Jorn. Phys., Vol. 39, No. 2, (1994), 139-148.
3864 ///
3865 /// Example 1 - script FitS.c:
3866 ///
3867 /// \image html spectrum2fit_stiefel_image001.jpg Fig. 1 Original two-dimensional spectrum with found peaks (using TSpectrum2 peak searching function). The positions of peaks were used as initial estimates in fitting procedure.
3868 ///
3869 /// \image html spectrum2fit_stiefel_image002.jpg Fig. 2 Fitted (generated from fitted parameters) spectrum of the data from Fig. 1 using Stiefel-Hestens method. Each peak was represented by 7 parameters, which together with Sigmax, Sigmay and a0 resulted in 38 fitted parameters. The chi-square after 1000 iterations was 0.642157.
3870 ///
3871 /// #### Script:
3872 ///
3873 /// Example to illustrate fitting function, algorithm without matrix inversion (AWMI) (class
3874 /// TSpectrumFit2). To execute this example, do
3875 ///
3876 /// `root > .x FitStiefel2.C`
3877 ///
3878 /// ~~~ {.cpp}
3879 /// void FitStiefel2() {
3880 /// Int_t i, j, nfound;
3881 /// Int_t nbinsx = 64;
3882 /// Int_t nbinsy = 64;
3883 /// Int_t xmin = 0;
3884 /// Int_t xmax = nbinsx;
3885 /// Int_t ymin = 0;
3886 /// Int_t ymax = nbinsy;
3887 /// Double_t ** source = new float *[nbinsx];
3888 /// Double_t ** dest = new float *[nbinsx];
3889 /// for (i=0;i<nbinsx;i++)
3890 /// source[i]=new float[nbinsy];
3891 /// for (i=0;i<nbinsx;i++)
3892 /// dest[i]= new float[nbinsy];
3893 /// TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
3894 /// TFile *f = new TFile("TSpectrum2.root");
3895 /// search=(TH2F*)f->Get("search4;1");
3896 /// TCanvas *Searching = new TCanvas("Searching","Two-dimensional fitting using Stiefel-Hestens method",10,10,1000,700);
3897 /// TSpectrum2 *s = new TSpectrum2();
3898 /// for (i = 0; i < nbinsx; i++){
3899 /// for (j = 0; j < nbinsy; j++){
3900 /// source[i][j] = search->GetBinContent(i + 1,j + 1);
3901 /// }
3902 /// }
3903 /// nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);
3904 /// printf("Found %d candidate peaks\n",nfound);
3905 /// Bool_t *FixPosX = new Bool_t[nfound];
3906 /// Bool_t *FixPosY = new Bool_t[nfound];
3907 /// Bool_t *FixAmp = new Bool_t[nfound];
3908 /// Double_t *PosX = new Double_t[nfound];
3909 /// Double_t *PosY = new Double_t[nfound];
3910 /// Double_t *Amp = new Double_t[nfound];
3911 /// Double_t *AmpXY = new Double_t[nfound];
3912 /// PosX = s->GetPositionX();
3913 /// PosY = s->GetPositionY();
3914 /// for(i = 0; i< nfound ; i++){
3915 /// FixPosX[i] = kFALSE;
3916 /// FixPosY[i] = kFALSE;
3917 /// FixAmp[i] = kFALSE;
3918 /// Amp[i] = source[(Int_t)(PosX[i]+0.5)][(Int_t)(PosY[i]+0.5)]; //initial values of peaks amplitudes, input parameters
3919 /// AmpXY[i] = 0;
3920 /// }
3921 /// //filling in the initial estimates of the input parameters
3922 /// TSpectrumFit2 *pfit=new
3923 /// TSpectrumFit2(nfound);
3924 /// pfit->SetFitParameters(xmin, xmax-1, ymin, ymax-1, 1000, 0.1,
3925 /// pfit->kFitOptimChiCounts, pfit->kFitAlphaHalving, pfit->kFitPower2,
3926 /// pfit->kFitTaylorOrderFirst);
3927 /// pfit->SetPeakParameters(2, kFALSE, 2, kFALSE, 0, kTRUE, PosX, (Bool_t *)
3928 /// FixPosX, PosY, (Bool_t *) FixPosY, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *)
3929 /// FixPosY, Amp, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp, AmpXY, (Bool_t *)
3930 /// FixAmp);
3931 /// pfit->SetBackgroundParameters(0, kFALSE, 0, kTRUE, 0, kTRUE);
3932 /// pfit->FitStiefel(source);
3933 /// for (i = 0; i < nbinsx; i++){
3934 /// for (j = 0; j < nbinsy; j++){
3935 /// search->SetBinContent(i + 1, j + 1,source[i][j]);
3936 /// }
3937 /// }
3938 /// search->Draw("SURF");
3939 /// }
3940 /// ~~~
3941 
3942 void TSpectrum2Fit::FitStiefel(Double_t **source)
3943 {
3944 
3945  Int_t i, i1, i2, j, k, shift =
3946  7 * fNPeaks + 14, peak_vel, size, iter, regul_cycle,
3947  flag;
3948  Double_t a, b, c, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi = 0
3949  , pi, pmin = 0, chi_cel = 0, chi_er;
3950  Double_t *working_space = new Double_t[5 * (7 * fNPeaks + 14)];
3951  for (i = 0, j = 0; i < fNPeaks; i++) {
3952  working_space[7 * i] = fAmpInit[i]; //vector parameter
3953  if (fFixAmp[i] == false) {
3954  working_space[shift + j] = fAmpInit[i]; //vector xk
3955  j += 1;
3956  }
3957  working_space[7 * i + 1] = fPositionInitX[i]; //vector parameter
3958  if (fFixPositionX[i] == false) {
3959  working_space[shift + j] = fPositionInitX[i]; //vector xk
3960  j += 1;
3961  }
3962  working_space[7 * i + 2] = fPositionInitY[i]; //vector parameter
3963  if (fFixPositionY[i] == false) {
3964  working_space[shift + j] = fPositionInitY[i]; //vector xk
3965  j += 1;
3966  }
3967  working_space[7 * i + 3] = fAmpInitX1[i]; //vector parameter
3968  if (fFixAmpX1[i] == false) {
3969  working_space[shift + j] = fAmpInitX1[i]; //vector xk
3970  j += 1;
3971  }
3972  working_space[7 * i + 4] = fAmpInitY1[i]; //vector parameter
3973  if (fFixAmpY1[i] == false) {
3974  working_space[shift + j] = fAmpInitY1[i]; //vector xk
3975  j += 1;
3976  }
3977  working_space[7 * i + 5] = fPositionInitX1[i]; //vector parameter
3978  if (fFixPositionX1[i] == false) {
3979  working_space[shift + j] = fPositionInitX1[i]; //vector xk
3980  j += 1;
3981  }
3982  working_space[7 * i + 6] = fPositionInitY1[i]; //vector parameter
3983  if (fFixPositionY1[i] == false) {
3984  working_space[shift + j] = fPositionInitY1[i]; //vector xk
3985  j += 1;
3986  }
3987  }
3988  peak_vel = 7 * i;
3989  working_space[7 * i] = fSigmaInitX; //vector parameter
3990  if (fFixSigmaX == false) {
3991  working_space[shift + j] = fSigmaInitX; //vector xk
3992  j += 1;
3993  }
3994  working_space[7 * i + 1] = fSigmaInitY; //vector parameter
3995  if (fFixSigmaY == false) {
3996  working_space[shift + j] = fSigmaInitY; //vector xk
3997  j += 1;
3998  }
3999  working_space[7 * i + 2] = fRoInit; //vector parameter
4000  if (fFixRo == false) {
4001  working_space[shift + j] = fRoInit; //vector xk
4002  j += 1;
4003  }
4004  working_space[7 * i + 3] = fA0Init; //vector parameter
4005  if (fFixA0 == false) {
4006  working_space[shift + j] = fA0Init; //vector xk
4007  j += 1;
4008  }
4009  working_space[7 * i + 4] = fAxInit; //vector parameter
4010  if (fFixAx == false) {
4011  working_space[shift + j] = fAxInit; //vector xk
4012  j += 1;
4013  }
4014  working_space[7 * i + 5] = fAyInit; //vector parameter
4015  if (fFixAy == false) {
4016  working_space[shift + j] = fAyInit; //vector xk
4017  j += 1;
4018  }
4019  working_space[7 * i + 6] = fTxyInit; //vector parameter
4020  if (fFixTxy == false) {
4021  working_space[shift + j] = fTxyInit; //vector xk
4022  j += 1;
4023  }
4024  working_space[7 * i + 7] = fSxyInit; //vector parameter
4025  if (fFixSxy == false) {
4026  working_space[shift + j] = fSxyInit; //vector xk
4027  j += 1;
4028  }
4029  working_space[7 * i + 8] = fTxInit; //vector parameter
4030  if (fFixTx == false) {
4031  working_space[shift + j] = fTxInit; //vector xk
4032  j += 1;
4033  }
4034  working_space[7 * i + 9] = fTyInit; //vector parameter
4035  if (fFixTy == false) {
4036  working_space[shift + j] = fTyInit; //vector xk
4037  j += 1;
4038  }
4039  working_space[7 * i + 10] = fSxyInit; //vector parameter
4040  if (fFixSx == false) {
4041  working_space[shift + j] = fSxInit; //vector xk
4042  j += 1;
4043  }
4044  working_space[7 * i + 11] = fSyInit; //vector parameter
4045  if (fFixSy == false) {
4046  working_space[shift + j] = fSyInit; //vector xk
4047  j += 1;
4048  }
4049  working_space[7 * i + 12] = fBxInit; //vector parameter
4050  if (fFixBx == false) {
4051  working_space[shift + j] = fBxInit; //vector xk
4052  j += 1;
4053  }
4054  working_space[7 * i + 13] = fByInit; //vector parameter
4055  if (fFixBy == false) {
4056  working_space[shift + j] = fByInit; //vector xk
4057  j += 1;
4058  }
4059  size = j;
4060  Double_t **working_matrix = new Double_t *[size];
4061  for (i = 0; i < size; i++)
4062  working_matrix[i] = new Double_t[size + 4];
4063  for (iter = 0; iter < fNumberIterations; iter++) {
4064  for (j = 0; j < size; j++) {
4065  working_space[3 * shift + j] = 0; //temp
4066  for (k = 0; k < (size + 4); k++) {
4067  working_matrix[j][k] = 0;
4068  }
4069  }
4070 
4071  //filling working matrix
4072  alpha = fAlpha;
4073  chi_opt = 0;
4074  for (i1 = fXmin; i1 <= fXmax; i1++) {
4075  for (i2 = fYmin; i2 <= fYmax; i2++) {
4076  //calculation of gradient vector
4077  for (j = 0, k = 0; j < fNPeaks; j++) {
4078  if (fFixAmp[j] == false) {
4079  working_space[2 * shift + k] =
4080  Deramp2(i1, i2,
4081  working_space[7 * j + 1],
4082  working_space[7 * j + 2],
4083  working_space[peak_vel],
4084  working_space[peak_vel + 1],
4085  working_space[peak_vel + 2],
4086  working_space[peak_vel + 6],
4087  working_space[peak_vel + 7],
4088  working_space[peak_vel + 12],
4089  working_space[peak_vel + 13]);
4090  k += 1;
4091  }
4092  if (fFixPositionX[j] == false) {
4093  working_space[2 * shift + k] =
4094  Deri02(i1, i2,
4095  working_space[7 * j],
4096  working_space[7 * j + 1],
4097  working_space[7 * j + 2],
4098  working_space[peak_vel],
4099  working_space[peak_vel + 1],
4100  working_space[peak_vel + 2],
4101  working_space[peak_vel + 6],
4102  working_space[peak_vel + 7],
4103  working_space[peak_vel + 12],
4104  working_space[peak_vel + 13]);
4105  k += 1;
4106  }
4107  if (fFixPositionY[j] == false) {
4108  working_space[2 * shift + k] =
4109  Derj02(i1, i2,
4110  working_space[7 * j],
4111  working_space[7 * j + 1],
4112  working_space[7 * j + 2],
4113  working_space[peak_vel],
4114  working_space[peak_vel + 1],
4115  working_space[peak_vel + 2],
4116  working_space[peak_vel + 6],
4117  working_space[peak_vel + 7],
4118  working_space[peak_vel + 12],
4119  working_space[peak_vel + 13]);
4120  k += 1;
4121  }
4122  if (fFixAmpX1[j] == false) {
4123  working_space[2 * shift + k] =
4124  Derampx(i1, working_space[7 * j + 5],
4125  working_space[peak_vel],
4126  working_space[peak_vel + 8],
4127  working_space[peak_vel + 10],
4128  working_space[peak_vel + 12]);
4129  k += 1;
4130  }
4131  if (fFixAmpY1[j] == false) {
4132  working_space[2 * shift + k] =
4133  Derampx(i2, working_space[7 * j + 6],
4134  working_space[peak_vel + 1],
4135  working_space[peak_vel + 9],
4136  working_space[peak_vel + 11],
4137  working_space[peak_vel + 13]);
4138  k += 1;
4139  }
4140  if (fFixPositionX1[j] == false) {
4141  working_space[2 * shift + k] =
4142  Deri01(i1, working_space[7 * j + 3],
4143  working_space[7 * j + 5],
4144  working_space[peak_vel],
4145  working_space[peak_vel + 8],
4146  working_space[peak_vel + 10],
4147  working_space[peak_vel + 12]);
4148  k += 1;
4149  }
4150  if (fFixPositionY1[j] == false) {
4151  working_space[2 * shift + k] =
4152  Deri01(i2, working_space[7 * j + 4],
4153  working_space[7 * j + 6],
4154  working_space[peak_vel + 1],
4155  working_space[peak_vel + 9],
4156  working_space[peak_vel + 11],
4157  working_space[peak_vel + 13]);
4158  k += 1;
4159  }
4160  } if (fFixSigmaX == false) {
4161  working_space[2 * shift + k] =
4162  Dersigmax(fNPeaks, i1, i2,
4163  working_space, working_space[peak_vel],
4164  working_space[peak_vel + 1],
4165  working_space[peak_vel + 2],
4166  working_space[peak_vel + 6],
4167  working_space[peak_vel + 7],
4168  working_space[peak_vel + 8],
4169  working_space[peak_vel + 10],
4170  working_space[peak_vel + 12],
4171  working_space[peak_vel + 13]);
4172  k += 1;
4173  }
4174  if (fFixSigmaY == false) {
4175  working_space[2 * shift + k] =
4176  Dersigmay(fNPeaks, i1, i2,
4177  working_space, working_space[peak_vel],
4178  working_space[peak_vel + 1],
4179  working_space[peak_vel + 2],
4180  working_space[peak_vel + 6],
4181  working_space[peak_vel + 7],
4182  working_space[peak_vel + 9],
4183  working_space[peak_vel + 11],
4184  working_space[peak_vel + 12],
4185  working_space[peak_vel + 13]);
4186  k += 1;
4187  }
4188  if (fFixRo == false) {
4189  working_space[2 * shift + k] =
4190  Derro(fNPeaks, i1, i2,
4191  working_space, working_space[peak_vel],
4192  working_space[peak_vel + 1],
4193  working_space[peak_vel + 2]);
4194  k += 1;
4195  }
4196  if (fFixA0 == false) {
4197  working_space[2 * shift + k] = 1.;
4198  k += 1;
4199  }
4200  if (fFixAx == false) {
4201  working_space[2 * shift + k] = i1;
4202  k += 1;
4203  }
4204  if (fFixAy == false) {
4205  working_space[2 * shift + k] = i2;
4206  k += 1;
4207  }
4208  if (fFixTxy == false) {
4209  working_space[2 * shift + k] =
4210  Dertxy(fNPeaks, i1, i2,
4211  working_space, working_space[peak_vel],
4212  working_space[peak_vel + 1],
4213  working_space[peak_vel + 12],
4214  working_space[peak_vel + 13]);
4215  k += 1;
4216  }
4217  if (fFixSxy == false) {
4218  working_space[2 * shift + k] =
4219  Dersxy(fNPeaks, i1, i2,
4220  working_space, working_space[peak_vel],
4221  working_space[peak_vel + 1]);
4222  k += 1;
4223  }
4224  if (fFixTx == false) {
4225  working_space[2 * shift + k] =
4226  Dertx(fNPeaks, i1, working_space,
4227  working_space[peak_vel],
4228  working_space[peak_vel + 12]);
4229  k += 1;
4230  }
4231  if (fFixTy == false) {
4232  working_space[2 * shift + k] =
4233  Derty(fNPeaks, i2, working_space,
4234  working_space[peak_vel + 1],
4235  working_space[peak_vel + 13]);
4236  k += 1;
4237  }
4238  if (fFixSx == false) {
4239  working_space[2 * shift + k] =
4240  Dersx(fNPeaks, i1, working_space,
4241  working_space[peak_vel]);
4242  k += 1;
4243  }
4244  if (fFixSy == false) {
4245  working_space[2 * shift + k] =
4246  Dersy(fNPeaks, i2, working_space,
4247  working_space[peak_vel + 1]);
4248  k += 1;
4249  }
4250  if (fFixBx == false) {
4251  working_space[2 * shift + k] =
4252  Derbx(fNPeaks, i1, i2,
4253  working_space, working_space[peak_vel],
4254  working_space[peak_vel + 1],
4255  working_space[peak_vel + 6],
4256  working_space[peak_vel + 8],
4257  working_space[peak_vel + 12],
4258  working_space[peak_vel + 13]);
4259  k += 1;
4260  }
4261  if (fFixBy == false) {
4262  working_space[2 * shift + k] =
4263  Derby(fNPeaks, i1, i2,
4264  working_space, working_space[peak_vel],
4265  working_space[peak_vel + 1],
4266  working_space[peak_vel + 6],
4267  working_space[peak_vel + 8],
4268  working_space[peak_vel + 12],
4269  working_space[peak_vel + 13]);
4270  k += 1;
4271  }
4272  yw = source[i1][i2];
4273  ywm = yw;
4274  f = Shape2(fNPeaks, i1, i2,
4275  working_space, working_space[peak_vel],
4276  working_space[peak_vel + 1],
4277  working_space[peak_vel + 2],
4278  working_space[peak_vel + 3],
4279  working_space[peak_vel + 4],
4280  working_space[peak_vel + 5],
4281  working_space[peak_vel + 6],
4282  working_space[peak_vel + 7],
4283  working_space[peak_vel + 8],
4284  working_space[peak_vel + 9],
4285  working_space[peak_vel + 10],
4286  working_space[peak_vel + 11],
4287  working_space[peak_vel + 12],
4288  working_space[peak_vel + 13]);
4289  if (fStatisticType == kFitOptimMaxLikelihood) {
4290  if (f > 0.00001)
4291  chi_opt += yw * TMath::Log(f) - f;
4292  }
4293 
4294  else {
4295  if (ywm != 0)
4296  chi_opt += (yw - f) * (yw - f) / ywm;
4297  }
4298  if (fStatisticType == kFitOptimChiFuncValues) {
4299  ywm = f;
4300  if (f < 0.00001)
4301  ywm = 0.00001;
4302  }
4303 
4304  else if (fStatisticType == kFitOptimMaxLikelihood) {
4305  ywm = f;
4306  if (f < 0.00001)
4307  ywm = 0.00001;
4308  }
4309 
4310  else {
4311  if (ywm == 0)
4312  ywm = 1;
4313  }
4314  for (j = 0; j < size; j++) {
4315  for (k = 0; k < size; k++) {
4316  b = working_space[2 * shift +
4317  j] * working_space[2 * shift +
4318  k] / ywm;
4319  if (fStatisticType == kFitOptimChiFuncValues)
4320  b = b * (4 * yw - 2 * f) / ywm;
4321  working_matrix[j][k] += b;
4322  if (j == k)
4323  working_space[3 * shift + j] += b;
4324  }
4325  }
4326  if (fStatisticType == kFitOptimChiFuncValues)
4327  b = (f * f - yw * yw) / (ywm * ywm);
4328 
4329  else
4330  b = (f - yw) / ywm;
4331  for (j = 0; j < size; j++) {
4332  working_matrix[j][size] -=
4333  b * working_space[2 * shift + j];
4334  }
4335  }
4336  }
4337  for (i = 0; i < size; i++) {
4338  working_matrix[i][size + 1] = 0; //xk
4339  }
4340  StiefelInversion(working_matrix, size);
4341  for (i = 0; i < size; i++) {
4342  working_space[2 * shift + i] = working_matrix[i][size + 1]; //der
4343  }
4344 
4345  //calculate chi_opt
4346  chi2 = chi_opt;
4347  chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
4348 
4349  //calculate new parameters
4350  regul_cycle = 0;
4351  for (j = 0; j < size; j++) {
4352  working_space[4 * shift + j] = working_space[shift + j]; //temp_xk[j]=xk[j]
4353  }
4354 
4355  do {
4356  if (fAlphaOptim == kFitAlphaOptimal) {
4357  if (fStatisticType != kFitOptimMaxLikelihood)
4358  chi_min = 10000 * chi2;
4359 
4360  else
4361  chi_min = 0.1 * chi2;
4362  flag = 0;
4363  for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
4364  for (j = 0; j < size; j++) {
4365  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]
4366  }
4367  for (i = 0, j = 0; i < fNPeaks; i++) {
4368  if (fFixAmp[i] == false) {
4369  if (working_space[shift + j] < 0) //xk[j]
4370  working_space[shift + j] = 0; //xk[j]
4371  working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
4372  j += 1;
4373  }
4374  if (fFixPositionX[i] == false) {
4375  if (working_space[shift + j] < fXmin) //xk[j]
4376  working_space[shift + j] = fXmin; //xk[j]
4377  if (working_space[shift + j] > fXmax) //xk[j]
4378  working_space[shift + j] = fXmax; //xk[j]
4379  working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
4380  j += 1;
4381  }
4382  if (fFixPositionY[i] == false) {
4383  if (working_space[shift + j] < fYmin) //xk[j]
4384  working_space[shift + j] = fYmin; //xk[j]
4385  if (working_space[shift + j] > fYmax) //xk[j]
4386  working_space[shift + j] = fYmax; //xk[j]
4387  working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
4388  j += 1;
4389  }
4390  if (fFixAmpX1[i] == false) {
4391  if (working_space[shift + j] < 0) //xk[j]
4392  working_space[shift + j] = 0; //xk[j]
4393  working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
4394  j += 1;
4395  }
4396  if (fFixAmpY1[i] == false) {
4397  if (working_space[shift + j] < 0) //xk[j]
4398  working_space[shift + j] = 0; //xk[j]
4399  working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
4400  j += 1;
4401  }
4402  if (fFixPositionX1[i] == false) {
4403  if (working_space[shift + j] < fXmin) //xk[j]
4404  working_space[shift + j] = fXmin; //xk[j]
4405  if (working_space[shift + j] > fXmax) //xk[j]
4406  working_space[shift + j] = fXmax; //xk[j]
4407  working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
4408  j += 1;
4409  }
4410  if (fFixPositionY1[i] == false) {
4411  if (working_space[shift + j] < fYmin) //xk[j]
4412  working_space[shift + j] = fYmin; //xk[j]
4413  if (working_space[shift + j] > fYmax) //xk[j]
4414  working_space[shift + j] = fYmax; //xk[j]
4415  working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
4416  j += 1;
4417  }
4418  }
4419  if (fFixSigmaX == false) {
4420  if (working_space[shift + j] < 0.001) { //xk[j]
4421  working_space[shift + j] = 0.001; //xk[j]
4422  }
4423  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
4424  j += 1;
4425  }
4426  if (fFixSigmaY == false) {
4427  if (working_space[shift + j] < 0.001) { //xk[j]
4428  working_space[shift + j] = 0.001; //xk[j]
4429  }
4430  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
4431  j += 1;
4432  }
4433  if (fFixRo == false) {
4434  if (working_space[shift + j] < -1) { //xk[j]
4435  working_space[shift + j] = -1; //xk[j]
4436  }
4437  if (working_space[shift + j] > 1) { //xk[j]
4438  working_space[shift + j] = 1; //xk[j]
4439  }
4440  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
4441  j += 1;
4442  }
4443  if (fFixA0 == false) {
4444  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
4445  j += 1;
4446  }
4447  if (fFixAx == false) {
4448  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
4449  j += 1;
4450  }
4451  if (fFixAy == false) {
4452  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
4453  j += 1;
4454  }
4455  if (fFixTxy == false) {
4456  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
4457  j += 1;
4458  }
4459  if (fFixSxy == false) {
4460  working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
4461  j += 1;
4462  }
4463  if (fFixTx == false) {
4464  working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
4465  j += 1;
4466  }
4467  if (fFixTy == false) {
4468  working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
4469  j += 1;
4470  }
4471  if (fFixSx == false) {
4472  working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
4473  j += 1;
4474  }
4475  if (fFixSy == false) {
4476  working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
4477  j += 1;
4478  }
4479  if (fFixBx == false) {
4480  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
4481  if (working_space[shift + j] < 0) //xk[j]
4482  working_space[shift + j] = -0.001; //xk[j]
4483  else
4484  working_space[shift + j] = 0.001; //xk[j]
4485  }
4486  working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
4487  j += 1;
4488  }
4489  if (fFixBy == false) {
4490  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
4491  if (working_space[shift + j] < 0) //xk[j]
4492  working_space[shift + j] = -0.001; //xk[j]
4493  else
4494  working_space[shift + j] = 0.001; //xk[j]
4495  }
4496  working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
4497  j += 1;
4498  }
4499  chi2 = 0;
4500  for (i1 = fXmin; i1 <= fXmax; i1++) {
4501  for (i2 = fYmin; i2 <= fYmax; i2++) {
4502  yw = source[i1][i2];
4503  ywm = yw;
4504  f = Shape2(fNPeaks, i1,
4505  i2, working_space,
4506  working_space[peak_vel],
4507  working_space[peak_vel + 1],
4508  working_space[peak_vel + 2],
4509  working_space[peak_vel + 3],
4510  working_space[peak_vel + 4],
4511  working_space[peak_vel + 5],
4512  working_space[peak_vel + 6],
4513  working_space[peak_vel + 7],
4514  working_space[peak_vel + 8],
4515  working_space[peak_vel + 9],
4516  working_space[peak_vel + 10],
4517  working_space[peak_vel + 11],
4518  working_space[peak_vel + 12],
4519  working_space[peak_vel + 13]);
4520  if (fStatisticType == kFitOptimChiFuncValues) {
4521  ywm = f;
4522  if (f < 0.00001)
4523  ywm = 0.00001;
4524  }
4525  if (fStatisticType == kFitOptimMaxLikelihood) {
4526  if (f > 0.00001)
4527  chi2 += yw * TMath::Log(f) - f;
4528  }
4529 
4530  else {
4531  if (ywm != 0)
4532  chi2 += (yw - f) * (yw - f) / ywm;
4533  }
4534  }
4535  }
4536  if ((chi2 < chi_min
4537  && fStatisticType != kFitOptimMaxLikelihood)
4538  || (chi2 > chi_min
4539  && fStatisticType == kFitOptimMaxLikelihood)) {
4540  pmin = pi, chi_min = chi2;
4541  }
4542 
4543  else
4544  flag = 1;
4545  if (pi == 0.1)
4546  chi_min = chi2;
4547  chi = chi_min;
4548  }
4549  if (pmin != 0.1) {
4550  for (j = 0; j < size; j++) {
4551  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]
4552  }
4553  for (i = 0, j = 0; i < fNPeaks; i++) {
4554  if (fFixAmp[i] == false) {
4555  if (working_space[shift + j] < 0) //xk[j]
4556  working_space[shift + j] = 0; //xk[j]
4557  working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
4558  j += 1;
4559  }
4560  if (fFixPositionX[i] == false) {
4561  if (working_space[shift + j] < fXmin) //xk[j]
4562  working_space[shift + j] = fXmin; //xk[j]
4563  if (working_space[shift + j] > fXmax) //xk[j]
4564  working_space[shift + j] = fXmax; //xk[j]
4565  working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
4566  j += 1;
4567  }
4568  if (fFixPositionY[i] == false) {
4569  if (working_space[shift + j] < fYmin) //xk[j]
4570  working_space[shift + j] = fYmin; //xk[j]
4571  if (working_space[shift + j] > fYmax) //xk[j]
4572  working_space[shift + j] = fYmax; //xk[j]
4573  working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
4574  j += 1;
4575  }
4576  if (fFixAmpX1[i] == false) {
4577  if (working_space[shift + j] < 0) //xk[j]
4578  working_space[shift + j] = 0; //xk[j]
4579  working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
4580  j += 1;
4581  }
4582  if (fFixAmpY1[i] == false) {
4583  if (working_space[shift + j] < 0) //xk[j]
4584  working_space[shift + j] = 0; //xk[j]
4585  working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
4586  j += 1;
4587  }
4588  if (fFixPositionX1[i] == false) {
4589  if (working_space[shift + j] < fXmin) //xk[j]
4590  working_space[shift + j] = fXmin; //xk[j]
4591  if (working_space[shift + j] > fXmax) //xk[j]
4592  working_space[shift + j] = fXmax; //xk[j]
4593  working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
4594  j += 1;
4595  }
4596  if (fFixPositionY1[i] == false) {
4597  if (working_space[shift + j] < fYmin) //xk[j]
4598  working_space[shift + j] = fYmin; //xk[j]
4599  if (working_space[shift + j] > fYmax) //xk[j]
4600  working_space[shift + j] = fYmax; //xk[j]
4601  working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
4602  j += 1;
4603  }
4604  }
4605  if (fFixSigmaX == false) {
4606  if (working_space[shift + j] < 0.001) { //xk[j]
4607  working_space[shift + j] = 0.001; //xk[j]
4608  }
4609  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
4610  j += 1;
4611  }
4612  if (fFixSigmaY == false) {
4613  if (working_space[shift + j] < 0.001) { //xk[j]
4614  working_space[shift + j] = 0.001; //xk[j]
4615  }
4616  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
4617  j += 1;
4618  }
4619  if (fFixRo == false) {
4620  if (working_space[shift + j] < -1) { //xk[j]
4621  working_space[shift + j] = -1; //xk[j]
4622  }
4623  if (working_space[shift + j] > 1) { //xk[j]
4624  working_space[shift + j] = 1; //xk[j]
4625  }
4626  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
4627  j += 1;
4628  }
4629  if (fFixA0 == false) {
4630  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
4631  j += 1;
4632  }
4633  if (fFixAx == false) {
4634  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
4635  j += 1;
4636  }
4637  if (fFixAy == false) {
4638  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
4639  j += 1;
4640  }
4641  if (fFixTxy == false) {
4642  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
4643  j += 1;
4644  }
4645  if (fFixSxy == false) {
4646  working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
4647  j += 1;
4648  }
4649  if (fFixTx == false) {
4650  working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
4651  j += 1;
4652  }
4653  if (fFixTy == false) {
4654  working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
4655  j += 1;
4656  }
4657  if (fFixSx == false) {
4658  working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
4659  j += 1;
4660  }
4661  if (fFixSy == false) {
4662  working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
4663  j += 1;
4664  }
4665  if (fFixBx == false) {
4666  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
4667  if (working_space[shift + j] < 0) //xk[j]
4668  working_space[shift + j] = -0.001; //xk[j]
4669  else
4670  working_space[shift + j] = 0.001; //xk[j]
4671  }
4672  working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
4673  j += 1;
4674  }
4675  if (fFixBy == false) {
4676  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
4677  if (working_space[shift + j] < 0) //xk[j]
4678  working_space[shift + j] = -0.001; //xk[j]
4679  else
4680  working_space[shift + j] = 0.001; //xk[j]
4681  }
4682  working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
4683  j += 1;
4684  }
4685  chi = chi_min;
4686  }
4687  }
4688 
4689  else {
4690  for (j = 0; j < size; j++) {
4691  working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
4692  }
4693  for (i = 0, j = 0; i < fNPeaks; i++) {
4694  if (fFixAmp[i] == false) {
4695  if (working_space[shift + j] < 0) //xk[j]
4696  working_space[shift + j] = 0; //xk[j]
4697  working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
4698  j += 1;
4699  }
4700  if (fFixPositionX[i] == false) {
4701  if (working_space[shift + j] < fXmin) //xk[j]
4702  working_space[shift + j] = fXmin; //xk[j]
4703  if (working_space[shift + j] > fXmax) //xk[j]
4704  working_space[shift + j] = fXmax; //xk[j]
4705  working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
4706  j += 1;
4707  }
4708  if (fFixPositionY[i] == false) {
4709  if (working_space[shift + j] < fYmin) //xk[j]
4710  working_space[shift + j] = fYmin; //xk[j]
4711  if (working_space[shift + j] > fYmax) //xk[j]
4712  working_space[shift + j] = fYmax; //xk[j]
4713  working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
4714  j += 1;
4715  }
4716  if (fFixAmpX1[i] == false) {
4717  if (working_space[shift + j] < 0) //xk[j]
4718  working_space[shift + j] = 0; //xk[j]
4719  working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
4720  j += 1;
4721  }
4722  if (fFixAmpY1[i] == false) {
4723  if (working_space[shift + j] < 0) //xk[j]
4724  working_space[shift + j] = 0; //xk[j]
4725  working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
4726  j += 1;
4727  }
4728  if (fFixPositionX1[i] == false) {
4729  if (working_space[shift + j] < fXmin) //xk[j]
4730  working_space[shift + j] = fXmin; //xk[j]
4731  if (working_space[shift + j] > fXmax) //xk[j]
4732  working_space[shift + j] = fXmax; //xk[j]
4733  working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
4734  j += 1;
4735  }
4736  if (fFixPositionY1[i] == false) {
4737  if (working_space[shift + j] < fYmin) //xk[j]
4738  working_space[shift + j] = fYmin; //xk[j]
4739  if (working_space[shift + j] > fYmax) //xk[j]
4740  working_space[shift + j] = fYmax; //xk[j]
4741  working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
4742  j += 1;
4743  }
4744  }
4745  if (fFixSigmaX == false) {
4746  if (working_space[shift + j] < 0.001) { //xk[j]
4747  working_space[shift + j] = 0.001; //xk[j]
4748  }
4749  working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
4750  j += 1;
4751  }
4752  if (fFixSigmaY == false) {
4753  if (working_space[shift + j] < 0.001) { //xk[j]
4754  working_space[shift + j] = 0.001; //xk[j]
4755  }
4756  working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
4757  j += 1;
4758  }
4759  if (fFixRo == false) {
4760  if (working_space[shift + j] < -1) { //xk[j]
4761  working_space[shift + j] = -1; //xk[j]
4762  }
4763  if (working_space[shift + j] > 1) { //xk[j]
4764  working_space[shift + j] = 1; //xk[j]
4765  }
4766  working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
4767  j += 1;
4768  }
4769  if (fFixA0 == false) {
4770  working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
4771  j += 1;
4772  }
4773  if (fFixAx == false) {
4774  working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
4775  j += 1;
4776  }
4777  if (fFixAy == false) {
4778  working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
4779  j += 1;
4780  }
4781  if (fFixTxy == false) {
4782  working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
4783  j += 1;
4784  }
4785  if (fFixSxy == false) {
4786  working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
4787  j += 1;
4788  }
4789  if (fFixTx == false) {
4790  working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
4791  j += 1;
4792  }
4793  if (fFixTy == false) {
4794  working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
4795  j += 1;
4796  }
4797  if (fFixSx == false) {
4798  working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
4799  j += 1;
4800  }
4801  if (fFixSy == false) {
4802  working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
4803  j += 1;
4804  }
4805  if (fFixBx == false) {
4806  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
4807  if (working_space[shift + j] < 0) //xk[j]
4808  working_space[shift + j] = -0.001; //xk[j]
4809  else
4810  working_space[shift + j] = 0.001; //xk[j]
4811  }
4812  working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
4813  j += 1;
4814  }
4815  if (fFixBy == false) {
4816  if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
4817  if (working_space[shift + j] < 0) //xk[j]
4818  working_space[shift + j] = -0.001; //xk[j]
4819  else
4820  working_space[shift + j] = 0.001; //xk[j]
4821  }
4822  working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
4823  j += 1;
4824  }
4825  chi = 0;
4826  for (i1 = fXmin; i1 <= fXmax; i1++) {
4827  for (i2 = fYmin; i2 <= fYmax; i2++) {
4828  yw = source[i1][i2];
4829  ywm = yw;
4830  f = Shape2(fNPeaks, i1, i2,
4831  working_space, working_space[peak_vel],
4832  working_space[peak_vel + 1],
4833  working_space[peak_vel + 2],
4834  working_space[peak_vel + 3],
4835  working_space[peak_vel + 4],
4836  working_space[peak_vel + 5],
4837  working_space[peak_vel + 6],
4838  working_space[peak_vel + 7],
4839  working_space[peak_vel + 8],
4840  working_space[peak_vel + 9],
4841  working_space[peak_vel + 10],
4842  working_space[peak_vel + 11],
4843  working_space[peak_vel + 12],
4844  working_space[peak_vel + 13]);
4845  if (fStatisticType == kFitOptimChiFuncValues) {
4846  ywm = f;
4847  if (f < 0.00001)
4848  ywm = 0.00001;
4849  }
4850  if (fStatisticType == kFitOptimMaxLikelihood) {
4851  if (f > 0.00001)
4852  chi += yw * TMath::Log(f) - f;
4853  }
4854 
4855  else {
4856  if (ywm != 0)
4857  chi += (yw - f) * (yw - f) / ywm;
4858  }
4859  }
4860  }
4861  }
4862  chi2 = chi;
4863  chi = TMath::Sqrt(TMath::Abs(chi));
4864  if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
4865  alpha = alpha * chi_opt / (2 * chi);
4866 
4867  else if (fAlphaOptim == kFitAlphaOptimal)
4868  alpha = alpha / 10.0;
4869  iter += 1;
4870  regul_cycle += 1;
4871  } while (((chi > chi_opt
4872  && fStatisticType != kFitOptimMaxLikelihood)
4873  || (chi < chi_opt
4874  && fStatisticType == kFitOptimMaxLikelihood))
4875  && regul_cycle < kFitNumRegulCycles);
4876  for (j = 0; j < size; j++) {
4877  working_space[4 * shift + j] = 0; //temp_xk[j]
4878  working_space[2 * shift + j] = 0; //der[j]
4879  }
4880  for (i1 = fXmin, chi_cel = 0; i1 <= fXmax; i1++) {
4881  for (i2 = fYmin; i2 <= fYmax; i2++) {
4882  yw = source[i1][i2];
4883  if (yw == 0)
4884  yw = 1;
4885  f = Shape2(fNPeaks, i1, i2,
4886  working_space, working_space[peak_vel],
4887  working_space[peak_vel + 1],
4888  working_space[peak_vel + 2],
4889  working_space[peak_vel + 3],
4890  working_space[peak_vel + 4],
4891  working_space[peak_vel + 5],
4892  working_space[peak_vel + 6],
4893  working_space[peak_vel + 7],
4894  working_space[peak_vel + 8],
4895  working_space[peak_vel + 9],
4896  working_space[peak_vel + 10],
4897  working_space[peak_vel + 11],
4898  working_space[peak_vel + 12],
4899  working_space[peak_vel + 13]);
4900  chi_opt = (yw - f) * (yw - f) / yw;
4901  chi_cel += (yw - f) * (yw - f) / yw;
4902 
4903  //calculate gradient vector
4904  for (j = 0, k = 0; j < fNPeaks; j++) {
4905  if (fFixAmp[j] == false) {
4906  a = Deramp2(i1, i2,
4907  working_space[7 * j + 1],
4908  working_space[7 * j + 2],
4909  working_space[peak_vel],
4910  working_space[peak_vel + 1],
4911  working_space[peak_vel + 2],
4912  working_space[peak_vel + 6],
4913  working_space[peak_vel + 7],
4914  working_space[peak_vel + 12],
4915  working_space[peak_vel + 13]);
4916  if (yw != 0) {
4917  working_space[2 * shift + k] += chi_opt; //der[k]
4918  b = a * a / yw;
4919  working_space[4 * shift + k] += b; //temp_xk[k]
4920  }
4921  k += 1;
4922  }
4923  if (fFixPositionX[j] == false) {
4924  a = Deri02(i1, i2,
4925  working_space[7 * j],
4926  working_space[7 * j + 1],
4927  working_space[7 * j + 2],
4928  working_space[peak_vel],
4929  working_space[peak_vel + 1],
4930  working_space[peak_vel + 2],
4931  working_space[peak_vel + 6],
4932  working_space[peak_vel + 7],
4933  working_space[peak_vel + 12],
4934  working_space[peak_vel + 13]);
4935  if (yw != 0) {
4936  working_space[2 * shift + k] += chi_opt; //der[k]
4937  b = a * a / yw;
4938  working_space[4 * shift + k] += b; //temp_xk[k]
4939  }
4940  k += 1;
4941  }
4942  if (fFixPositionY[j] == false) {
4943  a = Derj02(i1, i2,
4944  working_space[7 * j],
4945  working_space[7 * j + 1],
4946  working_space[7 * j + 2],
4947  working_space[peak_vel],
4948  working_space[peak_vel + 1],
4949  working_space[peak_vel + 2],
4950  working_space[peak_vel + 6],
4951  working_space[peak_vel + 7],
4952  working_space[peak_vel + 12],
4953  working_space[peak_vel + 13]);
4954  if (yw != 0) {
4955  working_space[2 * shift + k] += chi_opt; //der[k]
4956  b = a * a / yw;
4957  working_space[4 * shift + k] += b; //temp_xk[k]
4958  }
4959  k += 1;
4960  }
4961  if (fFixAmpX1[j] == false) {
4962  a = Derampx(i1, working_space[7 * j + 5],
4963  working_space[peak_vel],
4964  working_space[peak_vel + 8],
4965  working_space[peak_vel + 10],
4966  working_space[peak_vel + 12]);
4967  if (yw != 0) {
4968  working_space[2 * shift + k] += chi_opt; //der[k]
4969  b = a * a / yw;
4970  working_space[4 * shift + k] += b; //temp_xk[k]
4971  }
4972  k += 1;
4973  }
4974  if (fFixAmpY1[j] == false) {
4975  a = Derampx(i2, working_space[7 * j + 6],
4976  working_space[peak_vel + 1],
4977  working_space[peak_vel + 9],
4978  working_space[peak_vel + 11],
4979  working_space[peak_vel + 13]);
4980  if (yw != 0) {
4981  working_space[2 * shift + k] += chi_opt; //der[k]
4982  b = a * a / yw;
4983  working_space[4 * shift + k] += b; //temp_xk[k]
4984  }
4985  k += 1;
4986  }
4987  if (fFixPositionX1[j] == false) {
4988  a = Deri01(i1, working_space[7 * j + 3],
4989  working_space[7 * j + 5],
4990  working_space[peak_vel],
4991  working_space[peak_vel + 8],
4992  working_space[peak_vel + 10],
4993  working_space[peak_vel + 12]);
4994  if (yw != 0) {
4995  working_space[2 * shift + k] += chi_opt; //der[k]
4996  b = a * a / yw;
4997  working_space[4 * shift + k] += b; //temp_xk[k]
4998  }
4999  k += 1;
5000  }
5001  if (fFixPositionY1[j] == false) {
5002  a = Deri01(i2, working_space[7 * j + 4],
5003  working_space[7 * j + 6],
5004  working_space[peak_vel + 1],
5005  working_space[peak_vel + 9],
5006  working_space[peak_vel + 11],
5007  working_space[peak_vel + 13]);
5008  if (yw != 0) {
5009  working_space[2 * shift + k] += chi_opt; //der[k]
5010  b = a * a / yw;
5011  working_space[4 * shift + k] += b; //temp_xk[k]
5012  }
5013  k += 1;
5014  }
5015  }
5016  if (fFixSigmaX == false) {
5017  a = Dersigmax(fNPeaks, i1, i2,
5018  working_space, working_space[peak_vel],
5019  working_space[peak_vel + 1],
5020  working_space[peak_vel + 2],
5021  working_space[peak_vel + 6],
5022  working_space[peak_vel + 7],
5023  working_space[peak_vel + 8],
5024  working_space[peak_vel + 10],
5025  working_space[peak_vel + 12],
5026  working_space[peak_vel + 13]);
5027  if (yw != 0) {
5028  working_space[2 * shift + k] += chi_opt; //der[k]
5029  b = a * a / yw;
5030  working_space[4 * shift + k] += b; //temp_xk[k]
5031  }
5032  k += 1;
5033  }
5034  if (fFixSigmaY == false) {
5035  a = Dersigmay(fNPeaks, i1, i2,
5036  working_space, working_space[peak_vel],
5037  working_space[peak_vel + 1],
5038  working_space[peak_vel + 2],
5039  working_space[peak_vel + 6],
5040  working_space[peak_vel + 7],
5041  working_space[peak_vel + 9],
5042  working_space[peak_vel + 11],
5043  working_space[peak_vel + 12],
5044  working_space[peak_vel + 13]);
5045  if (yw != 0) {
5046  working_space[2 * shift + k] += chi_opt; //der[k]
5047  b = a * a / yw;
5048  working_space[4 * shift + k] += b; //temp_xk[k]
5049  }
5050  k += 1;
5051  }
5052  if (fFixRo == false) {
5053  a = Derro(fNPeaks, i1, i2,
5054  working_space, working_space[peak_vel],
5055  working_space[peak_vel + 1],
5056  working_space[peak_vel + 2]);
5057  if (yw != 0) {
5058  working_space[2 * shift + k] += chi_opt; //der[k]
5059  b = a * a / yw;
5060  working_space[4 * shift + k] += b; //temp_xk[k]
5061  }
5062  k += 1;
5063  }
5064  if (fFixA0 == false) {
5065  a = 1.;
5066  if (yw != 0) {
5067  working_space[2 * shift + k] += chi_opt; //der[k]
5068  b = a * a / yw;
5069  working_space[4 * shift + k] += b; //temp_xk[k]
5070  }
5071  k += 1;
5072  }
5073  if (fFixAx == false) {
5074  a = i1;
5075  if (yw != 0) {
5076  working_space[2 * shift + k] += chi_opt; //der[k]
5077  b = a * a / yw;
5078  working_space[4 * shift + k] += b; //temp_xk[k]
5079  }
5080  k += 1;
5081  }
5082  if (fFixAy == false) {
5083  a = i2;
5084  if (yw != 0) {
5085  working_space[2 * shift + k] += chi_opt; //der[k]
5086  b = a * a / yw;
5087  working_space[4 * shift + k] += b; //temp_xk[k]
5088  }
5089  k += 1;
5090  }
5091  if (fFixTxy == false) {
5092  a = Dertxy(fNPeaks, i1, i2,
5093  working_space, working_space[peak_vel],
5094  working_space[peak_vel + 1],
5095  working_space[peak_vel + 12],
5096  working_space[peak_vel + 13]);
5097  if (yw != 0) {
5098  working_space[2 * shift + k] += chi_opt; //der[k]
5099  b = a * a / yw;
5100  working_space[4 * shift + k] += b; //temp_xk[k]
5101  }
5102  k += 1;
5103  }
5104  if (fFixSxy == false) {
5105  a = Dersxy(fNPeaks, i1, i2,
5106  working_space, working_space[peak_vel],
5107  working_space[peak_vel + 1]);
5108  if (yw != 0) {
5109  working_space[2 * shift + k] += chi_opt; //der[k]
5110  b = a * a / yw;
5111  working_space[4 * shift + k] += b; //temp_xk[k]
5112  }
5113  k += 1;
5114  }
5115  if (fFixTx == false) {
5116  a = Dertx(fNPeaks, i1, working_space,
5117  working_space[peak_vel],
5118  working_space[peak_vel + 12]);
5119  if (yw != 0) {
5120  working_space[2 * shift + k] += chi_opt; //der[k]
5121  b = a * a / yw;
5122  working_space[4 * shift + k] += b; //temp_xk[k]
5123  }
5124  k += 1;
5125  }
5126  if (fFixTy == false) {
5127  a = Derty(fNPeaks, i2, working_space,
5128  working_space[peak_vel + 1],
5129  working_space[peak_vel + 13]);
5130  if (yw != 0) {
5131  working_space[2 * shift + k] += chi_opt; //der[k]
5132  b = a * a / yw;
5133  working_space[4 * shift + k] += b; //temp_xk[k]
5134  }
5135  k += 1;
5136  }
5137  if (fFixSx == false) {
5138  a = Dersx(fNPeaks, i1, working_space,
5139  working_space[peak_vel]);
5140  if (yw != 0) {
5141  working_space[2 * shift + k] += chi_opt; //der[k]
5142  b = a * a / yw;
5143  working_space[4 * shift + k] += b; //temp_xk[k]
5144  }
5145  k += 1;
5146  }
5147  if (fFixSy == false) {
5148  a = Dersy(fNPeaks, i2, working_space,
5149  working_space[peak_vel + 1]);
5150  if (yw != 0) {
5151  working_space[2 * shift + k] += chi_opt; //der[k]
5152  b = a * a / yw;
5153  working_space[4 * shift + k] += b; //temp_xk[k]
5154  }
5155  k += 1;
5156  }
5157  if (fFixBx == false) {
5158  a = Derbx(fNPeaks, i1, i2,
5159  working_space, working_space[peak_vel],
5160  working_space[peak_vel + 1],
5161  working_space[peak_vel + 6],
5162  working_space[peak_vel + 8],
5163  working_space[peak_vel + 12],
5164  working_space[peak_vel + 13]);
5165  if (yw != 0) {
5166  working_space[2 * shift + k] += chi_opt; //der[k]
5167  b = a * a / yw;
5168  working_space[4 * shift + k] += b; //temp_xk[k]
5169  }
5170  k += 1;
5171  }
5172  if (fFixBy == false) {
5173  a = Derby(fNPeaks, i1, i2,
5174  working_space, working_space[peak_vel],
5175  working_space[peak_vel + 1],
5176  working_space[peak_vel + 6],
5177  working_space[peak_vel + 8],
5178  working_space[peak_vel + 12],
5179  working_space[peak_vel + 13]);
5180  if (yw != 0) {
5181  working_space[2 * shift + k] += chi_opt; //der[k]
5182  b = a * a / yw;
5183  working_space[4 * shift + k] += b; //temp_xk[k]
5184  }
5185  k += 1;
5186  }
5187  }
5188  }
5189  }
5190  b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
5191  chi_er = chi_cel / b;
5192  for (i = 0, j = 0; i < fNPeaks; i++) {
5193  fVolume[i] =
5194  Volume(working_space[7 * i], working_space[peak_vel],
5195  working_space[peak_vel + 1], working_space[peak_vel + 2]);
5196  if (fVolume[i] > 0) {
5197  c = 0;
5198  if (fFixAmp[i] == false) {
5199  a = Derpa2(working_space[peak_vel],
5200  working_space[peak_vel + 1],
5201  working_space[peak_vel + 2]);
5202  b = working_space[4 * shift + j]; //temp_xk[j]
5203  if (b == 0)
5204  b = 1;
5205 
5206  else
5207  b = 1 / b;
5208  c = c + a * a * b;
5209  }
5210  if (fFixSigmaX == false) {
5211  a = Derpsigmax(working_space[shift + j],
5212  working_space[peak_vel + 1],
5213  working_space[peak_vel + 2]);
5214  b = working_space[4 * shift + peak_vel]; //temp_xk[j]
5215  if (b == 0)
5216  b = 1;
5217 
5218  else
5219  b = 1 / b;
5220  c = c + a * a * b;
5221  }
5222  if (fFixSigmaY == false) {
5223  a = Derpsigmay(working_space[shift + j],
5224  working_space[peak_vel],
5225  working_space[peak_vel + 2]);
5226  b = working_space[4 * shift + peak_vel + 1]; //temp_xk[j]
5227  if (b == 0)
5228  b = 1;
5229 
5230  else
5231  b = 1 / b;
5232  c = c + a * a * b;
5233  }
5234  if (fFixRo == false) {
5235  a = Derpro(working_space[shift + j], working_space[peak_vel],
5236  working_space[peak_vel + 1],
5237  working_space[peak_vel + 2]);
5238  b = working_space[4 * shift + peak_vel + 2]; //temp_xk[j]
5239  if (b == 0)
5240  b = 1;
5241 
5242  else
5243  b = 1 / b;
5244  c = c + a * a * b;
5245  }
5246  fVolumeErr[i] = TMath::Sqrt(TMath::Abs(chi_er * c));
5247  }
5248 
5249  else {
5250  fVolumeErr[i] = 0;
5251  }
5252  if (fFixAmp[i] == false) {
5253  fAmpCalc[i] = working_space[shift + j]; //xk[j]
5254  if (working_space[3 * shift + j] != 0)
5255  fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5256  j += 1;
5257  }
5258 
5259  else {
5260  fAmpCalc[i] = fAmpInit[i];
5261  fAmpErr[i] = 0;
5262  }
5263  if (fFixPositionX[i] == false) {
5264  fPositionCalcX[i] = working_space[shift + j]; //xk[j]
5265  if (working_space[3 * shift + j] != 0)
5266  fPositionErrX[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5267  j += 1;
5268  }
5269 
5270  else {
5271  fPositionCalcX[i] = fPositionInitX[i];
5272  fPositionErrX[i] = 0;
5273  }
5274  if (fFixPositionY[i] == false) {
5275  fPositionCalcY[i] = working_space[shift + j]; //xk[j]
5276  if (working_space[3 * shift + j] != 0)
5277  fPositionErrY[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5278  j += 1;
5279  }
5280 
5281  else {
5282  fPositionCalcY[i] = fPositionInitY[i];
5283  fPositionErrY[i] = 0;
5284  }
5285  if (fFixAmpX1[i] == false) {
5286  fAmpCalcX1[i] = working_space[shift + j]; //xk[j]
5287  if (working_space[3 * shift + j] != 0)
5288  fAmpErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5289  j += 1;
5290  }
5291 
5292  else {
5293  fAmpCalcX1[i] = fAmpInitX1[i];
5294  fAmpErrX1[i] = 0;
5295  }
5296  if (fFixAmpY1[i] == false) {
5297  fAmpCalcY1[i] = working_space[shift + j]; //xk[j]
5298  if (working_space[3 * shift + j] != 0)
5299  fAmpErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5300  j += 1;
5301  }
5302 
5303  else {
5304  fAmpCalcY1[i] = fAmpInitY1[i];
5305  fAmpErrY1[i] = 0;
5306  }
5307  if (fFixPositionX1[i] == false) {
5308  fPositionCalcX1[i] = working_space[shift + j]; //xk[j]
5309  if (working_space[3 * shift + j] != 0)
5310  fPositionErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5311  j += 1;
5312  }
5313 
5314  else {
5315  fPositionCalcX1[i] = fPositionInitX1[i];
5316  fPositionErrX1[i] = 0;
5317  }
5318  if (fFixPositionY1[i] == false) {
5319  fPositionCalcY1[i] = working_space[shift + j]; //xk[j]
5320  if (working_space[3 * shift + j] != 0)
5321  fPositionErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5322  j += 1;
5323  }
5324 
5325  else {
5326  fPositionCalcY1[i] = fPositionInitY1[i];
5327  fPositionErrY1[i] = 0;
5328  }
5329  }
5330  if (fFixSigmaX == false) {
5331  fSigmaCalcX = working_space[shift + j]; //xk[j]
5332  if (working_space[3 * shift + j] != 0) //temp[j]
5333  fSigmaErrX = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5334  j += 1;
5335  }
5336 
5337  else {
5338  fSigmaCalcX = fSigmaInitX;
5339  fSigmaErrX = 0;
5340  }
5341  if (fFixSigmaY == false) {
5342  fSigmaCalcY = working_space[shift + j]; //xk[j]
5343  if (working_space[3 * shift + j] != 0) //temp[j]
5344  fSigmaErrY = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5345  j += 1;
5346  }
5347 
5348  else {
5349  fSigmaCalcY = fSigmaInitY;
5350  fSigmaErrY = 0;
5351  }
5352  if (fFixRo == false) {
5353  fRoCalc = working_space[shift + j]; //xk[j]
5354  if (working_space[3 * shift + j] != 0) //temp[j]
5355  fRoErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5356  j += 1;
5357  }
5358 
5359  else {
5360  fRoCalc = fRoInit;
5361  fRoErr = 0;
5362  }
5363  if (fFixA0 == false) {
5364  fA0Calc = working_space[shift + j]; //xk[j]
5365  if (working_space[3 * shift + j] != 0) //temp[j]
5366  fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5367  j += 1;
5368  }
5369 
5370  else {
5371  fA0Calc = fA0Init;
5372  fA0Err = 0;
5373  }
5374  if (fFixAx == false) {
5375  fAxCalc = working_space[shift + j]; //xk[j]
5376  if (working_space[3 * shift + j] != 0) //temp[j]
5377  fAxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5378  j += 1;
5379  }
5380 
5381  else {
5382  fAxCalc = fAxInit;
5383  fAxErr = 0;
5384  }
5385  if (fFixAy == false) {
5386  fAyCalc = working_space[shift + j]; //xk[j]
5387  if (working_space[3 * shift + j] != 0) //temp[j]
5388  fAyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5389  j += 1;
5390  }
5391 
5392  else {
5393  fAyCalc = fAyInit;
5394  fAyErr = 0;
5395  }
5396  if (fFixTxy == false) {
5397  fTxyCalc = working_space[shift + j]; //xk[j]
5398  if (working_space[3 * shift + j] != 0) //temp[j]
5399  fTxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5400  j += 1;
5401  }
5402 
5403  else {
5404  fTxyCalc = fTxyInit;
5405  fTxyErr = 0;
5406  }
5407  if (fFixSxy == false) {
5408  fSxyCalc = working_space[shift + j]; //xk[j]
5409  if (working_space[3 * shift + j] != 0) //temp[j]
5410  fSxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5411  j += 1;
5412  }
5413 
5414  else {
5415  fSxyCalc = fSxyInit;
5416  fSxyErr = 0;
5417  }
5418  if (fFixTx == false) {
5419  fTxCalc = working_space[shift + j]; //xk[j]
5420  if (working_space[3 * shift + j] != 0) //temp[j]
5421  fTxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5422  j += 1;
5423  }
5424 
5425  else {
5426  fTxCalc = fTxInit;
5427  fTxErr = 0;
5428  }
5429  if (fFixTy == false) {
5430  fTyCalc = working_space[shift + j]; //xk[j]
5431  if (working_space[3 * shift + j] != 0) //temp[j]
5432  fTyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5433  j += 1;
5434  }
5435 
5436  else {
5437  fTyCalc = fTyInit;
5438  fTyErr = 0;
5439  }
5440  if (fFixSx == false) {
5441  fSxCalc = working_space[shift + j]; //xk[j]
5442  if (working_space[3 * shift + j] != 0) //temp[j]
5443  fSxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5444  j += 1;
5445  }
5446 
5447  else {
5448  fSxCalc = fSxInit;
5449  fSxErr = 0;
5450  }
5451  if (fFixSy == false) {
5452  fSyCalc = working_space[shift + j]; //xk[j]
5453  if (working_space[3 * shift + j] != 0) //temp[j]
5454  fSyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5455  j += 1;
5456  }
5457 
5458  else {
5459  fSyCalc = fSyInit;
5460  fSyErr = 0;
5461  }
5462  if (fFixBx == false) {
5463  fBxCalc = working_space[shift + j]; //xk[j]
5464  if (working_space[3 * shift + j] != 0) //temp[j]
5465  fBxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5466  j += 1;
5467  }
5468 
5469  else {
5470  fBxCalc = fBxInit;
5471  fBxErr = 0;
5472  }
5473  if (fFixBy == false) {
5474  fByCalc = working_space[shift + j]; //xk[j]
5475  if (working_space[3 * shift + j] != 0) //temp[j]
5476  fByErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
5477  j += 1;
5478  }
5479 
5480  else {
5481  fByCalc = fByInit;
5482  fByErr = 0;
5483  }
5484  b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
5485  fChi = chi_cel / b;
5486  for (i1 = fXmin; i1 <= fXmax; i1++) {
5487  for (i2 = fYmin; i2 <= fYmax; i2++) {
5488  f = Shape2(fNPeaks, i1, i2,
5489  working_space, working_space[peak_vel],
5490  working_space[peak_vel + 1],
5491  working_space[peak_vel + 2],
5492  working_space[peak_vel + 3],
5493  working_space[peak_vel + 4],
5494  working_space[peak_vel + 5],
5495  working_space[peak_vel + 6],
5496  working_space[peak_vel + 7],
5497  working_space[peak_vel + 8],
5498  working_space[peak_vel + 9],
5499  working_space[peak_vel + 10],
5500  working_space[peak_vel + 11],
5501  working_space[peak_vel + 12],
5502  working_space[peak_vel + 13]);
5503  source[i1][i2] = f;
5504 
5505  }
5506  }
5507  for (i = 0; i < size; i++) delete [] working_matrix[i];
5508  delete [] working_matrix;
5509  delete [] working_space;
5510  return;
5511 }
5512 
5513 ////////////////////////////////////////////////////////////////////////////////
5514 /// This function sets the following fitting parameters:
5515 /// - xmin, xmax, ymin, ymax - fitting region
5516 /// - numberIterations - # of desired iterations in the fit
5517 /// - alpha - convergence coefficient, it should be positive number and <=1, for details see references
5518 /// - 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
5519 /// - alphaOptim - optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal
5520 /// - power - possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting function.
5521 /// - fitTaylor - order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi fitting function.
5522 
5523 void TSpectrum2Fit::SetFitParameters(Int_t xmin,Int_t xmax,Int_t ymin,Int_t ymax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor)
5524 {
5525  if(xmin<0 || xmax <= xmin || ymin<0 || ymax <= ymin){
5526  Error("SetFitParameters", "Wrong range");
5527  return;
5528  }
5529  if (numberIterations <= 0){
5530  Error("SetFitParameters","Invalid number of iterations, must be positive");
5531  return;
5532  }
5533  if (alpha <= 0 || alpha > 1){
5534  Error ("SetFitParameters","Invalid step coefficient alpha, must be > than 0 and <=1");
5535  return;
5536  }
5537  if (statisticType != kFitOptimChiCounts
5538  && statisticType != kFitOptimChiFuncValues
5539  && statisticType != kFitOptimMaxLikelihood){
5540  Error("SetFitParameters","Wrong type of statistic");
5541  return;
5542  }
5543  if (alphaOptim != kFitAlphaHalving
5544  && alphaOptim != kFitAlphaOptimal){
5545  Error("SetFitParameters","Wrong optimization algorithm");
5546  return;
5547  }
5548  if (power != kFitPower2 && power != kFitPower4
5549  && power != kFitPower6 && power != kFitPower8
5550  && power != kFitPower10 && power != kFitPower12){
5551  Error("SetFitParameters","Wrong power");
5552  return;
5553  }
5554  if (fitTaylor != kFitTaylorOrderFirst
5555  && fitTaylor != kFitTaylorOrderSecond){
5556  Error("SetFitParameters","Wrong order of Taylor development");
5557  return;
5558  }
5559  fXmin=xmin,fXmax=xmax,fYmin=ymin,fYmax=ymax,fNumberIterations=numberIterations,fAlpha=alpha,fStatisticType=statisticType,fAlphaOptim=alphaOptim,fPower=power,fFitTaylor=fitTaylor;
5560 }
5561 
5562 ////////////////////////////////////////////////////////////////////////////////
5563 /// This function sets the following fitting parameters of peaks:
5564 /// - sigmaX - initial value of sigma x parameter
5565 /// - fixSigmaX - logical value of sigma x parameter, which allows to fix the parameter (not to fit)
5566 /// - sigmaY - initial value of sigma y parameter
5567 /// - fixSigmaY - logical value of sigma y parameter, which allows to fix the parameter (not to fit)
5568 /// - ro - initial value of ro parameter (correlation coefficient)
5569 /// - fixRo - logical value of ro parameter, which allows to fix the parameter (not to fit)
5570 /// - positionInitX - array of initial values of peaks x positions
5571 /// - fixPositionX - array of logical values which allow to fix appropriate x positions (not fit). However they are present in the estimated functional.
5572 /// - positionInitY - array of initial values of peaks y positions
5573 /// - fixPositionY - array of logical values which allow to fix appropriate y positions (not fit). However they are present in the estimated functional.
5574 /// - ampInit - array of initial values of 2D peaks amplitudes
5575 /// - fixAmp - array of logical values which allow to fix appropriate amplitudes of 2D peaks (not fit). However they are present in the estimated functional
5576 /// - ampInitX1 - array of initial values of amplitudes of 1D ridges in x direction
5577 /// - fixAmpX1 - array of logical values which allow to fix appropriate amplitudes of 1D ridges in x direction (not fit). However they are present in the estimated functional
5578 /// - ampInitY1 - array of initial values of amplitudes of 1D ridges in y direction
5579 /// - fixAmpY1 - array of logical values which allow to fix appropriate amplitudes of 1D ridges in y direction (not fit). However they are present in the estimated functional
5580 
5581 void TSpectrum2Fit::SetPeakParameters(Double_t sigmaX, Bool_t fixSigmaX, Double_t sigmaY, Bool_t fixSigmaY, Double_t ro, Bool_t fixRo, const Double_t *positionInitX, const Bool_t *fixPositionX, const Double_t *positionInitY, const Bool_t *fixPositionY, const Double_t *positionInitX1, const Bool_t *fixPositionX1, const Double_t *positionInitY1, const Bool_t *fixPositionY1, const Double_t *ampInit, const Bool_t *fixAmp, const Double_t *ampInitX1, const Bool_t *fixAmpX1, const Double_t *ampInitY1, const Bool_t *fixAmpY1)
5582 {
5583  if (sigmaX <= 0 || sigmaY <= 0){
5584  Error ("SetPeakParameters","Invalid sigma, must be > than 0");
5585  return;
5586  }
5587  if (ro < -1 || ro > 1){
5588  Error ("SetPeakParameters","Invalid ro, must be from region <-1,1>");
5589  return;
5590  }
5591  Int_t i;
5592  for(i=0; i < fNPeaks; i++){
5593  if(positionInitX[i] < fXmin || positionInitX[i] > fXmax){
5594  Error ("SetPeakParameters","Invalid peak position, must be in the range fXmin, fXmax");
5595  return;
5596  }
5597  if(positionInitY[i] < fYmin || positionInitY[i] > fYmax){
5598  Error ("SetPeakParameters","Invalid peak position, must be in the range fYmin, fYmax");
5599  return;
5600  }
5601  if(positionInitX1[i] < fXmin || positionInitX1[i] > fXmax){
5602  Error ("SetPeakParameters","Invalid ridge position, must be in the range fXmin, fXmax");
5603  return;
5604  }
5605  if(positionInitY1[i] < fYmin || positionInitY1[i] > fYmax){
5606  Error ("SetPeakParameters","Invalid ridge position, must be in the range fYmin, fYmax");
5607  return;
5608  }
5609  if(ampInit[i] < 0){
5610  Error ("SetPeakParameters","Invalid peak amplitude, must be > than 0");
5611  return;
5612  }
5613  if(ampInitX1[i] < 0){
5614  Error ("SetPeakParameters","Invalid x ridge amplitude, must be > than 0");
5615  return;
5616  }
5617  if(ampInitY1[i] < 0){
5618  Error ("SetPeakParameters","Invalid y ridge amplitude, must be > than 0");
5619  return;
5620  }
5621  }
5622  fSigmaInitX = sigmaX, fFixSigmaX = fixSigmaX, fSigmaInitY = sigmaY, fFixSigmaY = fixSigmaY, fRoInit = ro, fFixRo = fixRo;
5623  for(i=0; i < fNPeaks; i++){
5624  fPositionInitX[i] = positionInitX[i];
5625  fFixPositionX[i] = fixPositionX[i];
5626  fPositionInitY[i] = positionInitY[i];
5627  fFixPositionY[i] = fixPositionY[i];
5628  fPositionInitX1[i] = positionInitX1[i];
5629  fFixPositionX1[i] = fixPositionX1[i];
5630  fPositionInitY1[i] = positionInitY1[i];
5631  fFixPositionY1[i] = fixPositionY1[i];
5632  fAmpInit[i] = ampInit[i];
5633  fFixAmp[i] = fixAmp[i];
5634  fAmpInitX1[i] = ampInitX1[i];
5635  fFixAmpX1[i] = fixAmpX1[i];
5636  fAmpInitY1[i] = ampInitY1[i];
5637  fFixAmpY1[i] = fixAmpY1[i];
5638  }
5639 }
5640 
5641 ////////////////////////////////////////////////////////////////////////////////
5642 /// This function sets the following fitting parameters of background:
5643 /// - a0Init - initial value of a0 parameter (background is estimated as a0+ax*x+ay*y)
5644 /// - fixA0 - logical value of a0 parameter, which allows to fix the parameter (not to fit)
5645 /// - axInit - initial value of ax parameter
5646 /// - fixAx - logical value of ax parameter, which allows to fix the parameter (not to fit)
5647 /// - ayInit - initial value of ay parameter
5648 /// - fixAy - logical value of ay parameter, which allows to fix the parameter (not to fit)
5649 
5650 void TSpectrum2Fit::SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t axInit, Bool_t fixAx, Double_t ayInit, Bool_t fixAy)
5651 {
5652  fA0Init = a0Init;
5653  fFixA0 = fixA0;
5654  fAxInit = axInit;
5655  fFixAx = fixAx;
5656  fAyInit = ayInit;
5657  fFixAy = fixAy;
5658 }
5659 
5660 ////////////////////////////////////////////////////////////////////////////////
5661 /// This function sets the following fitting parameters of tails of peaks
5662 /// - tInitXY - initial value of txy parameter
5663 /// - fixTxy - logical value of txy parameter, which allows to fix the parameter (not to fit)
5664 /// - tInitX - initial value of tx parameter
5665 /// - fixTx - logical value of tx parameter, which allows to fix the parameter (not to fit)
5666 /// - tInitY - initial value of ty parameter
5667 /// - fixTy - logical value of ty parameter, which allows to fix the parameter (not to fit)
5668 /// - bInitX - initial value of bx parameter
5669 /// - fixBx - logical value of bx parameter, which allows to fix the parameter (not to fit)
5670 /// - bInitY - initial value of by parameter
5671 /// - fixBy - logical value of by parameter, which allows to fix the parameter (not to fit)
5672 /// - sInitXY - initial value of sxy parameter
5673 /// - fixSxy - logical value of sxy parameter, which allows to fix the parameter (not to fit)
5674 /// - sInitX - initial value of sx parameter
5675 /// - fixSx - logical value of sx parameter, which allows to fix the parameter (not to fit)
5676 /// - sInitY - initial value of sy parameter
5677 /// - fixSy - logical value of sy parameter, which allows to fix the parameter (not to fit)
5678 
5679 void TSpectrum2Fit::SetTailParameters(Double_t tInitXY, Bool_t fixTxy, Double_t tInitX, Bool_t fixTx, Double_t tInitY, Bool_t fixTy, Double_t bInitX, Bool_t fixBx, Double_t bInitY, Bool_t fixBy, Double_t sInitXY, Bool_t fixSxy, Double_t sInitX, Bool_t fixSx, Double_t sInitY, Bool_t fixSy)
5680 {
5681  fTxyInit = tInitXY;
5682  fFixTxy = fixTxy;
5683  fTxInit = tInitX;
5684  fFixTx = fixTx;
5685  fTyInit = tInitY;
5686  fFixTy = fixTy;
5687  fBxInit = bInitX;
5688  fFixBx = fixBx;
5689  fByInit = bInitY;
5690  fFixBy = fixBy;
5691  fSxyInit = sInitXY;
5692  fFixSxy = fixSxy;
5693  fSxInit = sInitX;
5694  fFixSx = fixSx;
5695  fSyInit = sInitY;
5696  fFixSy = fixSy;
5697 }
5698 
5699 ////////////////////////////////////////////////////////////////////////////////
5700 /// This function gets the positions of fitted 2D peaks and 1D ridges
5701 /// - positionX - gets vector of x positions of 2D peaks
5702 /// - positionY - gets vector of y positions of 2D peaks
5703 /// - positionX1 - gets vector of x positions of 1D ridges
5704 /// - positionY1 - gets vector of y positions of 1D ridges
5705 
5706 void TSpectrum2Fit::GetPositions(Double_t *positionsX, Double_t *positionsY, Double_t *positionsX1, Double_t *positionsY1)
5707 {
5708  for( Int_t i=0; i < fNPeaks; i++){
5709  positionsX[i] = fPositionCalcX[i];
5710  positionsY[i] = fPositionCalcY[i];
5711  positionsX1[i] = fPositionCalcX1[i];
5712  positionsY1[i] = fPositionCalcY1[i];
5713  }
5714 }
5715 
5716 ////////////////////////////////////////////////////////////////////////////////
5717 /// This function gets the errors of positions of fitted 2D peaks and 1D ridges
5718 /// - positionErrorsX - gets vector of errors of x positions of 2D peaks
5719 /// - positionErrorsY - gets vector of errors of y positions of 2D peaks
5720 /// - positionErrorsX1 - gets vector of errors of x positions of 1D ridges
5721 /// - positionErrorsY1 - gets vector of errors of y positions of 1D ridges
5722 
5723 void TSpectrum2Fit::GetPositionErrors(Double_t *positionErrorsX, Double_t *positionErrorsY, Double_t *positionErrorsX1, Double_t *positionErrorsY1)
5724 {
5725  for( Int_t i=0; i < fNPeaks; i++){
5726  positionErrorsX[i] = fPositionErrX[i];
5727  positionErrorsY[i] = fPositionErrY[i];
5728  positionErrorsX1[i] = fPositionErrX1[i];
5729  positionErrorsY1[i] = fPositionErrY1[i];
5730  }
5731 }
5732 
5733 ////////////////////////////////////////////////////////////////////////////////
5734 /// This function gets the amplitudes of fitted 2D peaks and 1D ridges
5735 /// - amplitudes - gets vector of amplitudes of 2D peaks
5736 /// - amplitudesX1 - gets vector of amplitudes of 1D ridges in x direction
5737 /// - amplitudesY1 - gets vector of amplitudes of 1D ridges in y direction
5738 
5739 void TSpectrum2Fit::GetAmplitudes(Double_t *amplitudes, Double_t *amplitudesX1, Double_t *amplitudesY1)
5740 {
5741  for( Int_t i=0; i < fNPeaks; i++){
5742  amplitudes[i] = fAmpCalc[i];
5743  amplitudesX1[i] = fAmpCalcX1[i];
5744  amplitudesY1[i] = fAmpCalcY1[i];
5745  }
5746 }
5747 
5748 ////////////////////////////////////////////////////////////////////////////////
5749 /// This function gets the amplitudes of fitted 2D peaks and 1D ridges
5750 /// - amplitudeErrors - gets vector of amplitudes errors of 2D peaks
5751 /// - amplitudeErrorsX1 - gets vector of amplitudes errors of 1D ridges in x direction
5752 /// - amplitudesErrorY1 - gets vector of amplitudes errors of 1D ridges in y direction
5753 
5754 void TSpectrum2Fit::GetAmplitudeErrors(Double_t *amplitudeErrors, Double_t *amplitudeErrorsX1, Double_t *amplitudeErrorsY1)
5755 {
5756  for( Int_t i=0; i < fNPeaks; i++){
5757  amplitudeErrors[i] = fAmpErr[i];
5758  amplitudeErrorsX1[i] = fAmpErrX1[i];
5759  amplitudeErrorsY1[i] = fAmpErrY1[i];
5760  }
5761 }
5762 
5763 ////////////////////////////////////////////////////////////////////////////////
5764 /// This function gets the volumes of fitted 2D peaks
5765 /// - volumes - gets vector of volumes of 2D peaks
5766 
5767 void TSpectrum2Fit::GetVolumes(Double_t *volumes)
5768 {
5769  for( Int_t i=0; i < fNPeaks; i++){
5770  volumes[i] = fVolume[i];
5771  }
5772 }
5773 
5774 ////////////////////////////////////////////////////////////////////////////////
5775 /// This function gets errors of the volumes of fitted 2D peaks
5776 /// - volumeErrors - gets vector of volumes errors of 2D peaks
5777 
5778 void TSpectrum2Fit::GetVolumeErrors(Double_t *volumeErrors)
5779 {
5780  for( Int_t i=0; i < fNPeaks; i++){
5781  volumeErrors[i] = fVolumeErr[i];
5782  }
5783 }
5784 
5785 ////////////////////////////////////////////////////////////////////////////////
5786 /// This function gets the sigma x parameter and its error
5787 /// - sigmaX - gets the fitted value of sigma x parameter
5788 /// - sigmaErrX - gets error value of sigma x parameter
5789 
5790 void TSpectrum2Fit::GetSigmaX(Double_t &sigmaX, Double_t &sigmaErrX)
5791 {
5792  sigmaX=fSigmaCalcX;
5793  sigmaErrX=fSigmaErrX;
5794 }
5795 
5796 ////////////////////////////////////////////////////////////////////////////////
5797 /// This function gets the sigma y parameter and its error
5798 /// - sigmaY - gets the fitted value of sigma y parameter
5799 /// - sigmaErrY - gets error value of sigma y parameter
5800 
5801 void TSpectrum2Fit::GetSigmaY(Double_t &sigmaY, Double_t &sigmaErrY)
5802 {
5803  sigmaY=fSigmaCalcY;
5804  sigmaErrY=fSigmaErrY;
5805 }
5806 
5807 ////////////////////////////////////////////////////////////////////////////////
5808 /// This function gets the ro parameter and its error
5809 /// - ro - gets the fitted value of ro parameter
5810 /// - roErr - gets error value of ro parameter
5811 
5812 void TSpectrum2Fit::GetRo(Double_t &ro, Double_t &roErr)
5813 {
5814  ro=fRoCalc;
5815  roErr=fRoErr;
5816 }
5817 
5818 ////////////////////////////////////////////////////////////////////////////////
5819 /// This function gets the background parameters and their errors
5820 /// - a0 - gets the fitted value of a0 parameter
5821 /// - a0Err - gets error value of a0 parameter
5822 /// - ax - gets the fitted value of ax parameter
5823 /// - axErr - gets error value of ax parameter
5824 /// - ay - gets the fitted value of ay parameter
5825 /// - ayErr - gets error value of ay parameter
5826 
5827 void TSpectrum2Fit::GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &ax, Double_t &axErr, Double_t &ay, Double_t &ayErr)
5828 {
5829  a0 = fA0Calc;
5830  a0Err = fA0Err;
5831  ax = fAxCalc;
5832  axErr = fAxErr;
5833  ay = fAyCalc;
5834  ayErr = fAyErr;
5835 }
5836 
5837 ////////////////////////////////////////////////////////////////////////////////
5838 /// This function gets the tail parameters and their errors
5839 /// - txy - gets the fitted value of txy parameter
5840 /// - txyErr - gets error value of txy parameter
5841 /// - tx - gets the fitted value of tx parameter
5842 /// - txErr - gets error value of tx parameter
5843 /// - ty - gets the fitted value of ty parameter
5844 /// - tyErr - gets error value of ty parameter
5845 /// - bx - gets the fitted value of bx parameter
5846 /// - bxErr - gets error value of bx parameter
5847 /// - by - gets the fitted value of by parameter
5848 /// - byErr - gets error value of by parameter
5849 /// - sxy - gets the fitted value of sxy parameter
5850 /// - sxyErr - gets error value of sxy parameter
5851 /// - sx - gets the fitted value of sx parameter
5852 /// - sxErr - gets error value of sx parameter
5853 /// - sy - gets the fitted value of sy parameter
5854 /// - syErr - gets error value of sy parameter
5855 
5856 void TSpectrum2Fit::GetTailParameters(Double_t &txy, Double_t &txyErr, Double_t &tx, Double_t &txErr, Double_t &ty, Double_t &tyErr, Double_t &bx, Double_t &bxErr, Double_t &by, Double_t &byErr, Double_t &sxy, Double_t &sxyErr, Double_t &sx, Double_t &sxErr, Double_t &sy, Double_t &syErr)
5857 {
5858  txy = fTxyCalc;
5859  txyErr = fTxyErr;
5860  tx = fTxCalc;
5861  txErr = fTxErr;
5862  ty = fTyCalc;
5863  tyErr = fTyErr;
5864  bx = fBxCalc;
5865  bxErr = fBxErr;
5866  by = fByCalc;
5867  byErr = fByErr;
5868  sxy = fSxyCalc;
5869  sxyErr = fSxyErr;
5870  sx = fSxCalc;
5871  sxErr = fSxErr;
5872  sy = fSyCalc;
5873  syErr = fSyErr;
5874 }
5875