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