Logo ROOT  

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