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],