00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "TSpectrum2Fit.h"
00033 #include "TMath.h"
00034
00035 ClassImp(TSpectrum2Fit)
00036
00037
00038 TSpectrum2Fit::TSpectrum2Fit() :TNamed("Spectrum2Fit", "Miroslav Morhac peak fitter")
00039 {
00040
00041 fNPeaks = 0;
00042 fNumberIterations = 1;
00043 fXmin = 0;
00044 fXmax = 100;
00045 fYmin = 0;
00046 fYmax = 100;
00047 fStatisticType = kFitOptimChiCounts;
00048 fAlphaOptim = kFitAlphaHalving;
00049 fPower = kFitPower2;
00050 fFitTaylor = kFitTaylorOrderFirst;
00051 fAlpha = 1;
00052 fChi = 0;
00053 fPositionInitX = 0;
00054 fPositionCalcX = 0;
00055 fPositionErrX = 0;
00056 fPositionInitY = 0;
00057 fPositionCalcY = 0;
00058 fPositionErrY = 0;
00059 fPositionInitX1 = 0;
00060 fPositionCalcX1 = 0;
00061 fPositionErrX1 = 0;
00062 fPositionInitY1 = 0;
00063 fPositionCalcY1 = 0;
00064 fPositionErrY1 = 0;
00065 fAmpInit = 0;
00066 fAmpCalc = 0;
00067 fAmpErr = 0;
00068 fAmpInitX1 = 0;
00069 fAmpCalcX1 = 0;
00070 fAmpErrX1 = 0;
00071 fAmpInitY1 = 0;
00072 fAmpCalcY1 = 0;
00073 fAmpErrY1 = 0;
00074 fVolume = 0;
00075 fVolumeErr = 0;
00076 fSigmaInitX = 2;
00077 fSigmaCalcX = 0;
00078 fSigmaErrX = 0;
00079 fSigmaInitY = 2;
00080 fSigmaCalcY = 0;
00081 fSigmaErrY = 0;
00082 fRoInit = 0;
00083 fRoCalc = 0;
00084 fRoErr = 0;
00085 fTxyInit = 0;
00086 fTxyCalc = 0;
00087 fTxyErr = 0;
00088 fTxInit = 0;
00089 fTxCalc = 0;
00090 fTxErr = 0;
00091 fTyInit = 0;
00092 fTyCalc = 0;
00093 fTyErr = 0;
00094 fBxInit = 1;
00095 fBxCalc = 0;
00096 fBxErr = 0;
00097 fByInit = 1;
00098 fByCalc = 0;
00099 fByErr = 0;
00100 fSxyInit = 0;
00101 fSxyCalc = 0;
00102 fSxyErr = 0;
00103 fSxInit = 0;
00104 fSxCalc = 0;
00105 fSxErr = 0;
00106 fSyInit = 0;
00107 fSyCalc = 0;
00108 fSyErr = 0;
00109 fA0Init = 0;
00110 fA0Calc = 0;
00111 fA0Err = 0;
00112 fAxInit = 0;
00113 fAxCalc = 0;
00114 fAxErr = 0;
00115 fAyInit = 0;
00116 fAyCalc = 0;
00117 fAyErr = 0;
00118 fFixPositionX = 0;
00119 fFixPositionY = 0;
00120 fFixPositionX1 = 0;
00121 fFixPositionY1 = 0;
00122 fFixAmp = 0;
00123 fFixAmpX1 = 0;
00124 fFixAmpY1 = 0;
00125 fFixSigmaX = false;
00126 fFixSigmaY = false;
00127 fFixRo = true;
00128 fFixTxy = true;
00129 fFixTx = true;
00130 fFixTy = true;
00131 fFixBx = true;
00132 fFixBy = true;
00133 fFixSxy = true;
00134 fFixSx = true;
00135 fFixSy = true;
00136 fFixA0 = true;
00137 fFixAx = true;
00138 fFixAy = true;
00139
00140 }
00141
00142 TSpectrum2Fit::TSpectrum2Fit(Int_t numberPeaks) :TNamed("Spectrum2Fit", "Miroslav Morhac peak fitter")
00143 {
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 if (numberPeaks <= 0){
00163 Error ("TSpectrum2Fit","Invalid number of peaks, must be > than 0");
00164 return;
00165 }
00166 fNPeaks = numberPeaks;
00167 fNumberIterations = 1;
00168 fXmin = 0;
00169 fXmax = 100;
00170 fYmin = 0;
00171 fYmax = 100;
00172 fStatisticType = kFitOptimChiCounts;
00173 fAlphaOptim = kFitAlphaHalving;
00174 fPower = kFitPower2;
00175 fFitTaylor = kFitTaylorOrderFirst;
00176 fAlpha = 1;
00177 fChi = 0;
00178 fPositionInitX = new Double_t[numberPeaks];
00179 fPositionCalcX = new Double_t[numberPeaks];
00180 fPositionErrX = new Double_t[numberPeaks];
00181 fPositionInitY = new Double_t[numberPeaks];
00182 fPositionCalcY = new Double_t[numberPeaks];
00183 fPositionErrY = new Double_t[numberPeaks];
00184 fPositionInitX1 = new Double_t[numberPeaks];
00185 fPositionCalcX1 = new Double_t[numberPeaks];
00186 fPositionErrX1 = new Double_t[numberPeaks];
00187 fPositionInitY1 = new Double_t[numberPeaks];
00188 fPositionCalcY1 = new Double_t[numberPeaks];
00189 fPositionErrY1 = new Double_t[numberPeaks];
00190 fAmpInit = new Double_t[numberPeaks];
00191 fAmpCalc = new Double_t[numberPeaks];
00192 fAmpErr = new Double_t[numberPeaks];
00193 fAmpInitX1 = new Double_t[numberPeaks];
00194 fAmpCalcX1 = new Double_t[numberPeaks];
00195 fAmpErrX1 = new Double_t[numberPeaks];
00196 fAmpInitY1 = new Double_t[numberPeaks];
00197 fAmpCalcY1 = new Double_t[numberPeaks];
00198 fAmpErrY1 = new Double_t[numberPeaks];
00199 fVolume = new Double_t[numberPeaks];
00200 fVolumeErr = new Double_t[numberPeaks];
00201 fSigmaInitX = 2;
00202 fSigmaCalcX = 0;
00203 fSigmaErrX = 0;
00204 fSigmaInitY = 2;
00205 fSigmaCalcY = 0;
00206 fSigmaErrY = 0;
00207 fRoInit = 0;
00208 fRoCalc = 0;
00209 fRoErr = 0;
00210 fTxyInit = 0;
00211 fTxyCalc = 0;
00212 fTxyErr = 0;
00213 fTxInit = 0;
00214 fTxCalc = 0;
00215 fTxErr = 0;
00216 fTyInit = 0;
00217 fTyCalc = 0;
00218 fTyErr = 0;
00219 fBxInit = 1;
00220 fBxCalc = 0;
00221 fBxErr = 0;
00222 fByInit = 1;
00223 fByCalc = 0;
00224 fByErr = 0;
00225 fSxyInit = 0;
00226 fSxyCalc = 0;
00227 fSxyErr = 0;
00228 fSxInit = 0;
00229 fSxCalc = 0;
00230 fSxErr = 0;
00231 fSyInit = 0;
00232 fSyCalc = 0;
00233 fSyErr = 0;
00234 fA0Init = 0;
00235 fA0Calc = 0;
00236 fA0Err = 0;
00237 fAxInit = 0;
00238 fAxCalc = 0;
00239 fAxErr = 0;
00240 fAyInit = 0;
00241 fAyCalc = 0;
00242 fAyErr = 0;
00243 fFixPositionX = new Bool_t[numberPeaks];
00244 fFixPositionY = new Bool_t[numberPeaks];
00245 fFixPositionX1 = new Bool_t[numberPeaks];
00246 fFixPositionY1 = new Bool_t[numberPeaks];
00247 fFixAmp = new Bool_t[numberPeaks];
00248 fFixAmpX1 = new Bool_t[numberPeaks];
00249 fFixAmpY1 = new Bool_t[numberPeaks];
00250 fFixSigmaX = false;
00251 fFixSigmaY = false;
00252 fFixRo = true;
00253 fFixTxy = true;
00254 fFixTx = true;
00255 fFixTy = true;
00256 fFixBx = true;
00257 fFixBy = true;
00258 fFixSxy = true;
00259 fFixSx = true;
00260 fFixSy = true;
00261 fFixA0 = true;
00262 fFixAx = true;
00263 fFixAy = true;
00264 }
00265
00266
00267 TSpectrum2Fit::~TSpectrum2Fit()
00268 {
00269
00270 delete [] fPositionInitX;
00271 delete [] fPositionCalcX;
00272 delete [] fPositionErrX;
00273 delete [] fFixPositionX;
00274 delete [] fPositionInitY;
00275 delete [] fPositionCalcY;
00276 delete [] fPositionErrY;
00277 delete [] fFixPositionY;
00278 delete [] fPositionInitX1;
00279 delete [] fPositionCalcX1;
00280 delete [] fPositionErrX1;
00281 delete [] fFixPositionX1;
00282 delete [] fPositionInitY1;
00283 delete [] fPositionCalcY1;
00284 delete [] fPositionErrY1;
00285 delete [] fFixPositionY1;
00286 delete [] fAmpInit;
00287 delete [] fAmpCalc;
00288 delete [] fAmpErr;
00289 delete [] fFixAmp;
00290 delete [] fAmpInitX1;
00291 delete [] fAmpCalcX1;
00292 delete [] fAmpErrX1;
00293 delete [] fFixAmpX1;
00294 delete [] fAmpInitY1;
00295 delete [] fAmpCalcY1;
00296 delete [] fAmpErrY1;
00297 delete [] fFixAmpY1;
00298 delete [] fVolume;
00299 delete [] fVolumeErr;
00300 }
00301
00302
00303
00304
00305 Double_t TSpectrum2Fit::Erfc(Double_t x)
00306 {
00307
00308
00309
00310
00311
00312
00313 Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap =
00314 0.47047;
00315 Double_t a, t, c, w;
00316 a = TMath::Abs(x);
00317 w = 1. + dap * a;
00318 t = 1. / w;
00319 w = a * a;
00320 if (w < 700)
00321 c = exp(-w);
00322
00323 else {
00324 c = 0;
00325 }
00326 c = c * t * (da1 + t * (da2 + t * da3));
00327 if (x < 0)
00328 c = 1. - c;
00329 return (c);
00330 }
00331
00332
00333 Double_t TSpectrum2Fit::Derfc(Double_t x)
00334 {
00335
00336
00337
00338
00339
00340
00341 Double_t a, t, c, w;
00342 Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap =
00343 0.47047;
00344 a = TMath::Abs(x);
00345 w = 1. + dap * a;
00346 t = 1. / w;
00347 w = a * a;
00348 if (w < 700)
00349 c = exp(-w);
00350
00351 else {
00352 c = 0;
00353 }
00354 c = (-1.) * dap * c * t * t * (da1 + t * (2. * da2 + t * 3. * da3)) -
00355 2. * a * Erfc(a);
00356 return (c);
00357 }
00358
00359
00360 Double_t TSpectrum2Fit::Ourpowl(Double_t a, Int_t pw)
00361 {
00362
00363 Double_t c;
00364 Double_t a2 = a*a;
00365 c = 1;
00366 if (pw > 0) c *= a2;
00367 if (pw > 2) c *= a2;
00368 if (pw > 4) c *= a2;
00369 if (pw > 6) c *= a2;
00370 if (pw > 8) c *= a2;
00371 if (pw > 10) c *= a2;
00372 if (pw > 12) c *= a2;
00373 return (c);
00374 }
00375 void TSpectrum2Fit::StiefelInversion(Double_t **a, Int_t size)
00376 {
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391 Int_t i, j, k = 0;
00392 Double_t sk = 0, b, lambdak, normk, normk_old = 0;
00393
00394 do {
00395 normk = 0;
00396
00397
00398 for (i = 0; i < size; i++) {
00399 a[i][size + 2] = -a[i][size];
00400 for (j = 0; j < size; j++) {
00401 a[i][size + 2] += a[i][j] * a[j][size + 1];
00402 }
00403 normk += a[i][size + 2] * a[i][size + 2];
00404 }
00405
00406
00407 if (k != 0) {
00408 sk = normk / normk_old;
00409 }
00410
00411
00412 for (i = 0; i < size; i++) {
00413 a[i][size + 3] = -a[i][size + 2] + sk * a[i][size + 3];
00414 }
00415
00416
00417 lambdak = 0;
00418 for (i = 0; i < size; i++) {
00419 for (j = 0, b = 0; j < size; j++) {
00420 b += a[i][j] * a[j][size + 3];
00421 }
00422 lambdak += b * a[i][size + 3];
00423 }
00424 if (TMath::Abs(lambdak) > 1e-50)
00425 lambdak = normk / lambdak;
00426
00427 else
00428 lambdak = 0;
00429 for (i = 0; i < size; i++)
00430 a[i][size + 1] += lambdak * a[i][size + 3];
00431 normk_old = normk;
00432 k += 1;
00433 } while (k < size && TMath::Abs(normk) > 1e-50);
00434 return;
00435 }
00436
00437
00438 Double_t TSpectrum2Fit::Shape2(Int_t numOfFittedPeaks, Double_t x, Double_t y,
00439 const Double_t *parameter, Double_t sigmax,
00440 Double_t sigmay, Double_t ro, Double_t a0, Double_t ax,
00441 Double_t ay, Double_t txy, Double_t sxy, Double_t tx,
00442 Double_t ty, Double_t sx, Double_t sy, Double_t bx,
00443 Double_t by)
00444 {
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462 Int_t j;
00463 Double_t r, p, r1, e, ex, ey, vx, s2, px, py, rx, ry, erx, ery;
00464 vx = 0;
00465 s2 = TMath::Sqrt(2.0);
00466 for (j = 0; j < numOfFittedPeaks; j++) {
00467 p = (x - parameter[7 * j + 1]) / sigmax;
00468 r = (y - parameter[7 * j + 2]) / sigmay;
00469 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
00470 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
00471 if (e < 700)
00472 r1 = exp(-e);
00473
00474 else {
00475 r1 = 0;
00476 }
00477 if (txy != 0) {
00478 px = 0, py = 0;
00479 erx = Erfc(p / s2 + 1 / (2 * bx)), ery =
00480 Erfc(r / s2 + 1 / (2 * by));
00481 ex = p / (s2 * bx), ey = r / (s2 * by);
00482 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
00483 px = exp(ex) * erx, py = exp(ey) * ery;
00484 }
00485 r1 += 0.5 * txy * px * py;
00486 }
00487 if (sxy != 0) {
00488 rx = Erfc(p / s2), ry = Erfc(r / s2);
00489 r1 += 0.5 * sxy * rx * ry;
00490 }
00491 vx = vx + parameter[7 * j] * r1;
00492 }
00493 p = (x - parameter[7 * j + 5]) / sigmax;
00494 if (TMath::Abs(p) < 3) {
00495 e = p * p / 2;
00496 if (e < 700)
00497 r1 = exp(-e);
00498
00499 else {
00500 r1 = 0;
00501 }
00502 if (tx != 0) {
00503 px = 0;
00504 erx = Erfc(p / s2 + 1 / (2 * bx));
00505 ex = p / (s2 * bx);
00506 if (TMath::Abs(ex) < 9) {
00507 px = exp(ex) * erx;
00508 }
00509 r1 += 0.5 * tx * px;
00510 }
00511 if (sx != 0) {
00512 rx = Erfc(p / s2);
00513 r1 += 0.5 * sx * rx;
00514 }
00515 vx = vx + parameter[7 * j + 3] * r1;
00516 }
00517 r = (y - parameter[7 * j + 6]) / sigmay;
00518 if (TMath::Abs(r) < 3) {
00519 e = r * r / 2;
00520 if (e < 700)
00521 r1 = exp(-e);
00522
00523 else {
00524 r1 = 0;
00525 }
00526 if (ty != 0) {
00527 py = 0;
00528 ery = Erfc(r / s2 + 1 / (2 * by));
00529 ey = r / (s2 * by);
00530 if (TMath::Abs(ey) < 9) {
00531 py = exp(ey) * ery;
00532 }
00533 r1 += 0.5 * ty * py;
00534 }
00535 if (sy != 0) {
00536 ry = Erfc(r / s2);
00537 r1 += 0.5 * sy * ry;
00538 }
00539 vx = vx + parameter[7 * j + 4] * r1;
00540 }
00541 }
00542 vx = vx + a0 + ax * x + ay * y;
00543 return (vx);
00544 }
00545
00546
00547 Double_t TSpectrum2Fit::Deramp2(Double_t x, Double_t y, Double_t x0, Double_t y0,
00548 Double_t sigmax, Double_t sigmay, Double_t ro,
00549 Double_t txy, Double_t sxy, Double_t bx, Double_t by)
00550 {
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568 Double_t p, r, r1 = 0, e, ex, ey, px, py, rx, ry, erx, ery, s2;
00569 p = (x - x0) / sigmax;
00570 r = (y - y0) / sigmay;
00571 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
00572 s2 = TMath::Sqrt(2.0);
00573 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
00574 if (e < 700)
00575 r1 = exp(-e);
00576
00577 else {
00578 r1 = 0;
00579 }
00580 if (txy != 0) {
00581 px = 0, py = 0;
00582 erx = Erfc(p / s2 + 1 / (2 * bx)), ery =
00583 Erfc(r / s2 + 1 / (2 * by));
00584 ex = p / (s2 * bx), ey = r / (s2 * by);
00585 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
00586 px = exp(ex) * erx, py = exp(ey) * ery;
00587 }
00588 r1 += 0.5 * txy * px * py;
00589 }
00590 if (sxy != 0) {
00591 rx = Erfc(p / s2), ry = Erfc(r / s2);
00592 r1 += 0.5 * sxy * rx * ry;
00593 }
00594 }
00595 return (r1);
00596 }
00597
00598
00599 Double_t TSpectrum2Fit::Derampx(Double_t x, Double_t x0, Double_t sigmax, Double_t tx,
00600 Double_t sx, Double_t bx)
00601 {
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617 Double_t p, r1 = 0, px, erx, rx, ex, s2;
00618 p = (x - x0) / sigmax;
00619 if (TMath::Abs(p) < 3) {
00620 s2 = TMath::Sqrt(2.0);
00621 p = p * p / 2;
00622 if (p < 700)
00623 r1 = exp(-p);
00624
00625 else {
00626 r1 = 0;
00627 }
00628 if (tx != 0) {
00629 px = 0;
00630 erx = Erfc(p / s2 + 1 / (2 * bx));
00631 ex = p / (s2 * bx);
00632 if (TMath::Abs(ex) < 9) {
00633 px = exp(ex) * erx;
00634 }
00635 r1 += 0.5 * tx * px;
00636 }
00637 if (sx != 0) {
00638 rx = Erfc(p / s2);
00639 r1 += 0.5 * sx * rx;
00640 }
00641 }
00642 return (r1);
00643 }
00644
00645
00646 Double_t TSpectrum2Fit::Deri02(Double_t x, Double_t y, Double_t a, Double_t x0,
00647 Double_t y0, Double_t sigmax, Double_t sigmay,
00648 Double_t ro, Double_t txy, Double_t sxy, Double_t bx,
00649 Double_t by)
00650 {
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669 Double_t p, r, r1 = 0, e, ex, ey, px, py, rx, ry, erx, ery, s2;
00670 p = (x - x0) / sigmax;
00671 r = (y - y0) / sigmay;
00672 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
00673 s2 = TMath::Sqrt(2.0);
00674 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
00675 if (e < 700)
00676 r1 = exp(-e);
00677
00678 else {
00679 r1 = 0;
00680 }
00681 e = -(ro * r - p) / sigmax;
00682 e = e / (1 - ro * ro);
00683 r1 = r1 * e;
00684 if (txy != 0) {
00685 px = 0, py = 0;
00686 erx =
00687 (-Erfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * sigmax) -
00688 Derfc(p / s2 + 1 / (2 * bx)) / (s2 * sigmax)), ery =
00689 Erfc(r / s2 + 1 / (2 * by));
00690 ex = p / (s2 * bx), ey = r / (s2 * by);
00691 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
00692 px = exp(ex) * erx, py = exp(ey) * ery;
00693 }
00694 r1 += 0.5 * txy * px * py;
00695 }
00696 if (sxy != 0) {
00697 rx = -Derfc(p / s2) / (s2 * sigmax), ry = Erfc(r / s2);
00698 r1 += 0.5 * sxy * rx * ry;
00699 }
00700 r1 = a * r1;
00701 }
00702 return (r1);
00703 }
00704
00705
00706 Double_t TSpectrum2Fit::Derderi02(Double_t x, Double_t y, Double_t a, Double_t x0,
00707 Double_t y0, Double_t sigmax, Double_t sigmay,
00708 Double_t ro)
00709 {
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726 Double_t p, r, r1 = 0, e;
00727 p = (x - x0) / sigmax;
00728 r = (y - y0) / sigmay;
00729 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
00730 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
00731 if (e < 700)
00732 r1 = exp(-e);
00733
00734 else {
00735 r1 = 0;
00736 }
00737 e = -(ro * r - p) / sigmax;
00738 e = e / (1 - ro * ro);
00739 r1 = r1 * (e * e - 1 / ((1 - ro * ro) * sigmax * sigmax));
00740 r1 = a * r1;
00741 }
00742 return (r1);
00743 }
00744 Double_t TSpectrum2Fit::Derj02(Double_t x, Double_t y, Double_t a, Double_t x0,
00745 Double_t y0, Double_t sigmax, Double_t sigmay,
00746 Double_t ro, Double_t txy, Double_t sxy, Double_t bx,
00747 Double_t by)
00748 {
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767 Double_t p, r, r1 = 0, e, ex, ey, px, py, rx, ry, erx, ery, s2;
00768 p = (x - x0) / sigmax;
00769 r = (y - y0) / sigmay;
00770 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
00771 s2 = TMath::Sqrt(2.0);
00772 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
00773 if (e < 700)
00774 r1 = exp(-e);
00775
00776 else {
00777 r1 = 0;
00778 }
00779 e = -(ro * p - r) / sigmay;
00780 e = e / (1 - ro * ro);
00781 r1 = r1 * e;
00782 if (txy != 0) {
00783 px = 0, py = 0;
00784 ery =
00785 (-Erfc(r / s2 + 1 / (2 * by)) / (s2 * by * sigmay) -
00786 Derfc(r / s2 + 1 / (2 * by)) / (s2 * sigmay)), erx =
00787 Erfc(p / s2 + 1 / (2 * bx));
00788 ex = p / (s2 * bx), ey = r / (s2 * by);
00789 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
00790 px = exp(ex) * erx, py = exp(ey) * ery;
00791 }
00792 r1 += 0.5 * txy * px * py;
00793 }
00794 if (sxy != 0) {
00795 ry = -Derfc(r / s2) / (s2 * sigmay), rx = Erfc(p / s2);
00796 r1 += 0.5 * sxy * rx * ry;
00797 }
00798 r1 = a * r1;
00799 }
00800 return (r1);
00801 }
00802
00803
00804 Double_t TSpectrum2Fit::Derderj02(Double_t x, Double_t y, Double_t a, Double_t x0,
00805 Double_t y0, Double_t sigmax, Double_t sigmay,
00806 Double_t ro)
00807 {
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824 Double_t p, r, r1 = 0, e;
00825 p = (x - x0) / sigmax;
00826 r = (y - y0) / sigmay;
00827 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
00828 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
00829 if (e < 700)
00830 r1 = exp(-e);
00831
00832 else {
00833 r1 = 0;
00834 }
00835 e = -(ro * p - r) / sigmay;
00836 e = e / (1 - ro * ro);
00837 r1 = r1 * (e * e - 1 / ((1 - ro * ro) * sigmay * sigmay));
00838 r1 = a * r1;
00839 }
00840 return (r1);
00841 }
00842
00843
00844 Double_t TSpectrum2Fit::Deri01(Double_t x, Double_t ax, Double_t x0, Double_t sigmax,
00845 Double_t tx, Double_t sx, Double_t bx)
00846 {
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862 Double_t p, e, r1 = 0, px, rx, erx, ex, s2;
00863 p = (x - x0) / sigmax;
00864 if (TMath::Abs(p) < 3) {
00865 s2 = TMath::Sqrt(2.0);
00866 e = p * p / 2;
00867 if (e < 700)
00868 r1 = exp(-e);
00869
00870 else {
00871 r1 = 0;
00872 }
00873 r1 = r1 * p / sigmax;
00874 if (tx != 0) {
00875 px = 0;
00876 erx =
00877 (-Erfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * sigmax) -
00878 Derfc(p / s2 + 1 / (2 * bx)) / (s2 * sigmax));
00879 ex = p / (s2 * bx);
00880 if (TMath::Abs(ex) < 9)
00881 px = exp(ex) * erx;
00882 r1 += 0.5 * tx * px;
00883 }
00884 if (sx != 0) {
00885 rx = -Derfc(p / s2) / (s2 * sigmax);
00886 r1 += 0.5 * sx * rx;
00887 }
00888 r1 = ax * r1;
00889 }
00890 return (r1);
00891 }
00892
00893
00894 Double_t TSpectrum2Fit::Derderi01(Double_t x, Double_t ax, Double_t x0,
00895 Double_t sigmax)
00896 {
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909 Double_t p, e, r1 = 0;
00910 p = (x - x0) / sigmax;
00911 if (TMath::Abs(p) < 3) {
00912 e = p * p / 2;
00913 if (e < 700)
00914 r1 = exp(-e);
00915
00916 else {
00917 r1 = 0;
00918 }
00919 r1 = r1 * (p * p / (sigmax * sigmax) - 1 / (sigmax * sigmax));
00920 r1 = ax * r1;
00921 }
00922 return (r1);
00923 }
00924
00925
00926 Double_t TSpectrum2Fit::Dersigmax(Int_t numOfFittedPeaks, Double_t x, Double_t y,
00927 const Double_t *parameter, Double_t sigmax,
00928 Double_t sigmay, Double_t ro, Double_t txy,
00929 Double_t sxy, Double_t tx, Double_t sx, Double_t bx,
00930 Double_t by)
00931 {
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948 Double_t p, r, r1 =
00949 0, e, a, b, x0, y0, s2, px, py, rx, ry, erx, ery, ex, ey;
00950 Int_t j;
00951 s2 = TMath::Sqrt(2.0);
00952 for (j = 0; j < numOfFittedPeaks; j++) {
00953 a = parameter[7 * j];
00954 x0 = parameter[7 * j + 1];
00955 y0 = parameter[7 * j + 2];
00956 p = (x - x0) / sigmax;
00957 r = (y - y0) / sigmay;
00958 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
00959 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
00960 if (e < 700)
00961 e = exp(-e);
00962
00963 else {
00964 e = 0;
00965 }
00966 b = -(ro * p * r - p * p) / sigmax;
00967 e = e * b / (1 - ro * ro);
00968 if (txy != 0) {
00969 px = 0, py = 0;
00970 erx =
00971 -Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * sigmax) -
00972 Derfc(p / s2 + 1 / (2 * bx)) * p / (s2 * sigmax), ery =
00973 Erfc(r / s2 + 1 / (2 * by));
00974 ex = p / (s2 * bx), ey = r / (s2 * by);
00975 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
00976 px = exp(ex) * erx, py = exp(ey) * ery;
00977 }
00978 e += 0.5 * txy * px * py;
00979 }
00980 if (sxy != 0) {
00981 rx = -Derfc(p / s2) * p / (s2 * sigmax), ry = Erfc(r / s2);
00982 e += 0.5 * sxy * rx * ry;
00983 }
00984 r1 = r1 + a * e;
00985 }
00986 if (TMath::Abs(p) < 3) {
00987 x0 = parameter[7 * j + 5];
00988 p = (x - x0) / sigmax;
00989 b = p * p / 2;
00990 if (b < 700)
00991 e = exp(-b);
00992
00993 else {
00994 e = 0;
00995 }
00996 e = 2 * b * e / sigmax;
00997 if (tx != 0) {
00998 px = 0;
00999 erx =
01000 (-Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * sigmax) -
01001 Derfc(p / s2 + 1 / (2 * bx)) * p / (s2 * sigmax));
01002 ex = p / (s2 * bx);
01003 if (TMath::Abs(ex) < 9)
01004 px = exp(ex) * erx;
01005 e += 0.5 * tx * px;
01006 }
01007 if (sx != 0) {
01008 rx = -Derfc(p / s2) * p / (s2 * sigmax);
01009 e += 0.5 * sx * rx;
01010 }
01011 r1 += parameter[7 * j + 3] * e;
01012 }
01013 }
01014 return (r1);
01015 }
01016
01017
01018 Double_t TSpectrum2Fit::Derdersigmax(Int_t numOfFittedPeaks, Double_t x,
01019 Double_t y, const Double_t *parameter,
01020 Double_t sigmax, Double_t sigmay,
01021 Double_t ro)
01022 {
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037 Double_t p, r, r1 = 0, e, a, b, x0, y0;
01038 Int_t j;
01039 for (j = 0; j < numOfFittedPeaks; j++) {
01040 a = parameter[7 * j];
01041 x0 = parameter[7 * j + 1];
01042 y0 = parameter[7 * j + 2];
01043 p = (x - x0) / sigmax;
01044 r = (y - y0) / sigmay;
01045 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
01046 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
01047 if (e < 700)
01048 e = exp(-e);
01049
01050 else {
01051 e = 0;
01052 }
01053 b = -(ro * p * r - p * p) / sigmax;
01054 e = e * (b * b / (1 - ro * ro) -
01055 (3 * p * p - 2 * ro * p * r) / (sigmax * sigmax)) / (1 -
01056 ro
01057 *
01058 ro);
01059 r1 = r1 + a * e;
01060 }
01061 if (TMath::Abs(p) < 3) {
01062 x0 = parameter[7 * j + 5];
01063 p = (x - x0) / sigmax;
01064 b = p * p / 2;
01065 if (b < 700)
01066 e = exp(-b);
01067
01068 else {
01069 e = 0;
01070 }
01071 e = e * (4 * b * b - 6 * b) / (sigmax * sigmax);
01072 r1 += parameter[7 * j + 3] * e;
01073 }
01074 }
01075 return (r1);
01076 }
01077
01078
01079 Double_t TSpectrum2Fit::Dersigmay(Int_t numOfFittedPeaks, Double_t x, Double_t y,
01080 const Double_t *parameter, Double_t sigmax,
01081 Double_t sigmay, Double_t ro, Double_t txy,
01082 Double_t sxy, Double_t ty, Double_t sy, Double_t bx,
01083 Double_t by)
01084 {
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101 Double_t p, r, r1 =
01102 0, e, a, b, x0, y0, s2, px, py, rx, ry, erx, ery, ex, ey;
01103 Int_t j;
01104 s2 = TMath::Sqrt(2.0);
01105 for (j = 0; j < numOfFittedPeaks; j++) {
01106 a = parameter[7 * j];
01107 x0 = parameter[7 * j + 1];
01108 y0 = parameter[7 * j + 2];
01109 p = (x - x0) / sigmax;
01110 r = (y - y0) / sigmay;
01111 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
01112 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
01113 if (e < 700)
01114 e = exp(-e);
01115
01116 else {
01117 e = 0;
01118 }
01119 b = -(ro * p * r - r * r) / sigmay;
01120 e = e * b / (1 - ro * ro);
01121 if (txy != 0) {
01122 px = 0, py = 0;
01123 ery =
01124 -Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * sigmay) -
01125 Derfc(r / s2 + 1 / (2 * by)) * r / (s2 * sigmay), erx =
01126 Erfc(p / s2 + 1 / (2 * bx));
01127 ex = p / (s2 * bx), ey = r / (s2 * by);
01128 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
01129 px = exp(ex) * erx, py = exp(ey) * ery;
01130 }
01131 e += 0.5 * txy * px * py;
01132 }
01133 if (sxy != 0) {
01134 ry = -Derfc(r / s2) * r / (s2 * sigmay), rx = Erfc(p / s2);
01135 e += 0.5 * sxy * rx * ry;
01136 }
01137 r1 = r1 + a * e;
01138 }
01139 if (TMath::Abs(r) < 3) {
01140 y0 = parameter[7 * j + 6];
01141 r = (y - y0) / sigmay;
01142 b = r * r / 2;
01143 if (b < 700)
01144 e = exp(-b);
01145
01146 else {
01147 e = 0;
01148 }
01149 e = 2 * b * e / sigmay;
01150 if (ty != 0) {
01151 py = 0;
01152 ery =
01153 (-Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * sigmay) -
01154 Derfc(r / s2 + 1 / (2 * by)) * r / (s2 * sigmay));
01155 ey = r / (s2 * by);
01156 if (TMath::Abs(ey) < 9)
01157 py = exp(ey) * ery;
01158 e += 0.5 * ty * py;
01159 }
01160 if (sy != 0) {
01161 ry = -Derfc(r / s2) * r / (s2 * sigmay);
01162 e += 0.5 * sy * ry;
01163 }
01164 r1 += parameter[7 * j + 4] * e;
01165 }
01166 }
01167 return (r1);
01168 }
01169
01170
01171 Double_t TSpectrum2Fit::Derdersigmay(Int_t numOfFittedPeaks, Double_t x,
01172 Double_t y, const Double_t *parameter,
01173 Double_t sigmax, Double_t sigmay,
01174 Double_t ro)
01175 {
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190 Double_t p, r, r1 = 0, e, a, b, x0, y0;
01191 Int_t j;
01192 for (j = 0; j < numOfFittedPeaks; j++) {
01193 a = parameter[7 * j];
01194 x0 = parameter[7 * j + 1];
01195 y0 = parameter[7 * j + 2];
01196 p = (x - x0) / sigmax;
01197 r = (y - y0) / sigmay;
01198 if (TMath::Abs(p) < 3 && TMath::Abs(r) < 3) {
01199 e = (p * p - 2 * ro * p * r + r * r) / (2 * (1 - ro * ro));
01200 if (e < 700)
01201 e = exp(-e);
01202
01203 else {
01204 e = 0;
01205 }
01206 b = -(ro * p * r - r * r) / sigmay;
01207 e = e * (b * b / (1 - ro * ro) -
01208 (3 * r * r - 2 * ro * r * p) / (sigmay * sigmay)) / (1 -
01209 ro
01210 *
01211 ro);
01212 r1 = r1 + a * e;
01213 }
01214 if (TMath::Abs(r) < 3) {
01215 y0 = parameter[7 * j + 6];
01216 r = (y - y0) / sigmay;
01217 b = r * r / 2;
01218 if (b < 700)
01219 e = exp(-b);
01220
01221 else {
01222 e = 0;
01223 }
01224 e = e * (4 * b * b - 6 * b) / (sigmay * sigmay);
01225 r1 += parameter[7 * j + 4] * e;
01226 }
01227 }
01228 return (r1);
01229 }
01230
01231
01232 Double_t TSpectrum2Fit::Derro(Int_t numOfFittedPeaks, Double_t x, Double_t y,
01233 const Double_t *parameter, Double_t sx, Double_t sy,
01234 Double_t r)
01235 {
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250 Double_t px, qx, rx, vx, x0, y0, a, ex, tx;
01251 Int_t j;
01252 vx = 0;
01253 for (j = 0; j < numOfFittedPeaks; j++) {
01254 a = parameter[7 * j];
01255 x0 = parameter[7 * j + 1];
01256 y0 = parameter[7 * j + 2];
01257 px = (x - x0) / sx;
01258 qx = (y - y0) / sy;
01259 if (TMath::Abs(px) < 3 && TMath::Abs(qx) < 3) {
01260 rx = (px * px - 2 * r * px * qx + qx * qx);
01261 ex = rx / (2 * (1 - r * r));
01262 if ((ex) < 700)
01263 ex = exp(-ex);
01264
01265 else {
01266 ex = 0;
01267 }
01268 tx = px * qx / (1 - r * r);
01269 tx = tx - r * rx / ((1 - r * r) * (1 - r * r));
01270 vx = vx + a * ex * tx;
01271 }
01272 }
01273 return (vx);
01274 }
01275
01276
01277 Double_t TSpectrum2Fit::Dertxy(Int_t numOfFittedPeaks, Double_t x, Double_t y,
01278 const Double_t *parameter, Double_t sigmax,
01279 Double_t sigmay, Double_t bx, Double_t by)
01280 {
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295 Double_t p, r, r1 = 0, ex, ey, px, py, erx, ery, s2, x0, y0, a;
01296 Int_t j;
01297 s2 = TMath::Sqrt(2.0);
01298 for (j = 0; j < numOfFittedPeaks; j++) {
01299 a = parameter[7 * j];
01300 x0 = parameter[7 * j + 1];
01301 y0 = parameter[7 * j + 2];
01302 p = (x - x0) / sigmax;
01303 r = (y - y0) / sigmay;
01304 px = 0, py = 0;
01305 erx = Erfc(p / s2 + 1 / (2 * bx)), ery =
01306 Erfc(r / s2 + 1 / (2 * by));
01307 ex = p / (s2 * bx), ey = r / (s2 * by);
01308 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
01309 px = exp(ex) * erx, py = exp(ey) * ery;
01310 }
01311 r1 += 0.5 * a * px * py;
01312 }
01313 return (r1);
01314 }
01315
01316
01317 Double_t TSpectrum2Fit::Dersxy(Int_t numOfFittedPeaks, Double_t x, Double_t y,
01318 const Double_t *parameter, Double_t sigmax,
01319 Double_t sigmay)
01320 {
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334 Double_t p, r, r1 = 0, rx, ry, x0, y0, a, s2;
01335 Int_t j;
01336 s2 = TMath::Sqrt(2.0);
01337 for (j = 0; j < numOfFittedPeaks; j++) {
01338 a = parameter[7 * j];
01339 x0 = parameter[7 * j + 1];
01340 y0 = parameter[7 * j + 2];
01341 p = (x - x0) / sigmax;
01342 r = (y - y0) / sigmay;
01343 rx = Erfc(p / s2), ry = Erfc(r / s2);
01344 r1 += 0.5 * a * rx * ry;
01345 }
01346 return (r1);
01347 }
01348
01349
01350 Double_t TSpectrum2Fit::Dertx(Int_t numOfFittedPeaks, Double_t x,
01351 const Double_t *parameter, Double_t sigmax,
01352 Double_t bx)
01353 {
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367 Double_t p, r1 = 0, ex, px, erx, s2, ax, x0;
01368 Int_t j;
01369 s2 = TMath::Sqrt(2.0);
01370 for (j = 0; j < numOfFittedPeaks; j++) {
01371 ax = parameter[7 * j + 3];
01372 x0 = parameter[7 * j + 5];
01373 p = (x - x0) / sigmax;
01374 px = 0;
01375 erx = Erfc(p / s2 + 1 / (2 * bx));
01376 ex = p / (s2 * bx);
01377 if (TMath::Abs(ex) < 9) {
01378 px = exp(ex) * erx;
01379 }
01380 r1 += 0.5 * ax * px;
01381 }
01382 return (r1);
01383 }
01384
01385
01386 Double_t TSpectrum2Fit::Derty(Int_t numOfFittedPeaks, Double_t x,
01387 const Double_t *parameter, Double_t sigmax,
01388 Double_t bx)
01389 {
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403 Double_t p, r1 = 0, ex, px, erx, s2, ax, x0;
01404 Int_t j;
01405 s2 = TMath::Sqrt(2.0);
01406 for (j = 0; j < numOfFittedPeaks; j++) {
01407 ax = parameter[7 * j + 4];
01408 x0 = parameter[7 * j + 6];
01409 p = (x - x0) / sigmax;
01410 px = 0;
01411 erx = Erfc(p / s2 + 1 / (2 * bx));
01412 ex = p / (s2 * bx);
01413 if (TMath::Abs(ex) < 9) {
01414 px = exp(ex) * erx;
01415 }
01416 r1 += 0.5 * ax * px;
01417 }
01418 return (r1);
01419 }
01420
01421
01422 Double_t TSpectrum2Fit::Dersx(Int_t numOfFittedPeaks, Double_t x,
01423 const Double_t *parameter, Double_t sigmax)
01424 {
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437 Double_t p, r1 = 0, rx, ax, x0, s2;
01438 Int_t j;
01439 s2 = TMath::Sqrt(2.0);
01440 for (j = 0; j < numOfFittedPeaks; j++) {
01441 ax = parameter[7 * j + 3];
01442 x0 = parameter[7 * j + 5];
01443 p = (x - x0) / sigmax;
01444 s2 = TMath::Sqrt(2.0);
01445 rx = Erfc(p / s2);
01446 r1 += 0.5 * ax * rx;
01447 }
01448 return (r1);
01449 }
01450
01451
01452 Double_t TSpectrum2Fit::Dersy(Int_t numOfFittedPeaks, Double_t x,
01453 const Double_t *parameter, Double_t sigmax)
01454 {
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467 Double_t p, r1 = 0, rx, ax, x0, s2;
01468 Int_t j;
01469 s2 = TMath::Sqrt(2.0);
01470 for (j = 0; j < numOfFittedPeaks; j++) {
01471 ax = parameter[7 * j + 4];
01472 x0 = parameter[7 * j + 6];
01473 p = (x - x0) / sigmax;
01474 s2 = TMath::Sqrt(2.0);
01475 rx = Erfc(p / s2);
01476 r1 += 0.5 * ax * rx;
01477 }
01478 return (r1);
01479 }
01480
01481
01482 Double_t TSpectrum2Fit::Derbx(Int_t numOfFittedPeaks, Double_t x, Double_t y,
01483 const Double_t *parameter, Double_t sigmax,
01484 Double_t sigmay, Double_t txy, Double_t tx, Double_t bx,
01485 Double_t by)
01486 {
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502 Double_t p, r, r1 = 0, a, x0, y0, s2, px, py, erx, ery, ex, ey;
01503 Int_t j;
01504 s2 = TMath::Sqrt(2.0);
01505 for (j = 0; j < numOfFittedPeaks; j++) {
01506 a = parameter[7 * j];
01507 x0 = parameter[7 * j + 1];
01508 y0 = parameter[7 * j + 2];
01509 p = (x - x0) / sigmax;
01510 r = (y - y0) / sigmay;
01511 if (txy != 0) {
01512 px = 0, py = 0;
01513 erx =
01514 -Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * bx) -
01515 Derfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * bx), ery =
01516 Erfc(r / s2 + 1 / (2 * by));
01517 ex = p / (s2 * bx), ey = r / (s2 * by);
01518 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
01519 px = exp(ex) * erx, py = exp(ey) * ery;
01520 }
01521 r1 += 0.5 * a * txy * px * py;
01522 }
01523 a = parameter[7 * j + 3];
01524 x0 = parameter[7 * j + 5];
01525 p = (x - x0) / sigmax;
01526 if (tx != 0) {
01527 px = 0;
01528 erx =
01529 (-Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * bx) -
01530 Derfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * bx));
01531 ex = p / (s2 * bx);
01532 if (TMath::Abs(ex) < 9)
01533 px = exp(ex) * erx;
01534 r1 += 0.5 * a * tx * px;
01535 }
01536 }
01537 return (r1);
01538 }
01539
01540
01541 Double_t TSpectrum2Fit::Derby(Int_t numOfFittedPeaks, Double_t x, Double_t y,
01542 const Double_t *parameter, Double_t sigmax,
01543 Double_t sigmay, Double_t txy, Double_t ty, Double_t bx,
01544 Double_t by)
01545 {
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561 Double_t p, r, r1 = 0, a, x0, y0, s2, px, py, erx, ery, ex, ey;
01562 Int_t j;
01563 s2 = TMath::Sqrt(2.0);
01564 for (j = 0; j < numOfFittedPeaks; j++) {
01565 a = parameter[7 * j];
01566 x0 = parameter[7 * j + 1];
01567 y0 = parameter[7 * j + 2];
01568 p = (x - x0) / sigmax;
01569 r = (y - y0) / sigmay;
01570 if (txy != 0) {
01571 px = 0, py = 0;
01572 ery =
01573 -Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * by) -
01574 Derfc(r / s2 + 1 / (2 * by)) / (s2 * by * by), erx =
01575 Erfc(p / s2 + 1 / (2 * bx));
01576 ex = p / (s2 * bx), ey = r / (s2 * by);
01577 if (TMath::Abs(ex) < 9 && TMath::Abs(ey) < 9) {
01578 px = exp(ex) * erx, py = exp(ey) * ery;
01579 }
01580 r1 += 0.5 * a * txy * px * py;
01581 }
01582 a = parameter[7 * j + 4];
01583 y0 = parameter[7 * j + 6];
01584 r = (y - y0) / sigmay;
01585 if (ty != 0) {
01586 py = 0;
01587 ery =
01588 (-Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * by) -
01589 Derfc(r / s2 + 1 / (2 * by)) / (s2 * by * by));
01590 ey = r / (s2 * by);
01591 if (TMath::Abs(ey) < 9)
01592 py = exp(ey) * ery;
01593 r1 += 0.5 * a * ty * py;
01594 }
01595 }
01596 return (r1);
01597 }
01598
01599
01600 Double_t TSpectrum2Fit::Volume(Double_t a, Double_t sx, Double_t sy, Double_t ro)
01601 {
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01612 Double_t pi = 3.1415926535, r;
01613 r = 1 - ro * ro;
01614 if (r > 0)
01615 r = TMath::Sqrt(r);
01616
01617 else {
01618 return (0);
01619 }
01620 r = 2 * a * pi * sx * sy * r;
01621 return (r);
01622 }
01623
01624
01625 Double_t TSpectrum2Fit::Derpa2(Double_t sx, Double_t sy, Double_t ro)
01626 {
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637 Double_t pi = 3.1415926535, r;
01638 r = 1 - ro * ro;
01639 if (r > 0)
01640 r = TMath::Sqrt(r);
01641
01642 else {
01643 return (0);
01644 }
01645 r = 2 * pi * sx * sy * r;
01646 return (r);
01647 }
01648
01649
01650 Double_t TSpectrum2Fit::Derpsigmax(Double_t a, Double_t sy, Double_t ro)
01651 {
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663 Double_t pi = 3.1415926535, r;
01664 r = 1 - ro * ro;
01665 if (r > 0)
01666 r = TMath::Sqrt(r);
01667
01668 else {
01669 return (0);
01670 }
01671 r = a * 2 * pi * sy * r;
01672 return (r);
01673 }
01674
01675
01676 Double_t TSpectrum2Fit::Derpsigmay(Double_t a, Double_t sx, Double_t ro)
01677 {
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689 Double_t pi = 3.1415926535, r;
01690 r = 1 - ro * ro;
01691 if (r > 0)
01692 r = TMath::Sqrt(r);
01693
01694 else {
01695 return (0);
01696 }
01697 r = a * 2 * pi * sx * r;
01698 return (r);
01699 }
01700
01701
01702 Double_t TSpectrum2Fit::Derpro(Double_t a, Double_t sx, Double_t sy, Double_t ro)
01703 {
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715 Double_t pi = 3.1415926535, r;
01716 r = 1 - ro * ro;
01717 if (r > 0)
01718 r = TMath::Sqrt(r);
01719
01720 else {
01721 return (0);
01722 }
01723 r = -a * 2 * pi * sx * sy * ro / r;
01724 return (r);
01725 }
01726
01727
01728
01729
01730
01731 void TSpectrum2Fit::FitAwmi(Float_t **source)
01732 {
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098
02099
02100
02101
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221
02222
02223
02224
02225
02226
02227
02228
02229
02230
02231
02232
02233
02234
02235
02236
02237
02238
02239
02240
02241
02242
02243
02244
02245
02246
02247
02248
02249
02250
02251
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368
02369
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386
02387
02388
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460
02461
02462
02463
02464
02465
02466
02467
02468
02469
02470
02471
02472
02473
02474
02475
02476
02477
02478
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488
02489
02490
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562
02563
02564
02565
02566
02567
02568
02569
02570
02571
02572
02573
02574
02575
02576
02577
02578
02579
02580
02581
02582
02583
02584
02585
02586
02587
02588
02589
02590
02591
02592
02593
02594
02595
02596
02597
02598
02599
02600
02601
02602
02603
02604
02605
02606
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618
02619
02620
02621
02622
02623
02624
02625
02626
02627
02628
02629
02630
02631
02632
02633
02634
02635
02636
02637
02638
02639
02640 Int_t i, i1, i2, j, k, shift =
02641 7 * fNPeaks + 14, peak_vel, size, iter, pw,
02642 regul_cycle, flag;
02643 Double_t a, b, c, d = 0, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
02644 0, pi, pmin = 0, chi_cel = 0, chi_er;
02645 Double_t *working_space = new Double_t[5 * (7 * fNPeaks + 14)];
02646 for (i = 0, j = 0; i < fNPeaks; i++) {
02647 working_space[7 * i] = fAmpInit[i];
02648 if (fFixAmp[i] == false) {
02649 working_space[shift + j] = fAmpInit[i];
02650 j += 1;
02651 }
02652 working_space[7 * i + 1] = fPositionInitX[i];
02653 if (fFixPositionX[i] == false) {
02654 working_space[shift + j] = fPositionInitX[i];
02655 j += 1;
02656 }
02657 working_space[7 * i + 2] = fPositionInitY[i];
02658 if (fFixPositionY[i] == false) {
02659 working_space[shift + j] = fPositionInitY[i];
02660 j += 1;
02661 }
02662 working_space[7 * i + 3] = fAmpInitX1[i];
02663 if (fFixAmpX1[i] == false) {
02664 working_space[shift + j] = fAmpInitX1[i];
02665 j += 1;
02666 }
02667 working_space[7 * i + 4] = fAmpInitY1[i];
02668 if (fFixAmpY1[i] == false) {
02669 working_space[shift + j] = fAmpInitY1[i];
02670 j += 1;
02671 }
02672 working_space[7 * i + 5] = fPositionInitX1[i];
02673 if (fFixPositionX1[i] == false) {
02674 working_space[shift + j] = fPositionInitX1[i];
02675 j += 1;
02676 }
02677 working_space[7 * i + 6] = fPositionInitY1[i];
02678 if (fFixPositionY1[i] == false) {
02679 working_space[shift + j] = fPositionInitY1[i];
02680 j += 1;
02681 }
02682 }
02683 peak_vel = 7 * i;
02684 working_space[7 * i] = fSigmaInitX;
02685 if (fFixSigmaX == false) {
02686 working_space[shift + j] = fSigmaInitX;
02687 j += 1;
02688 }
02689 working_space[7 * i + 1] = fSigmaInitY;
02690 if (fFixSigmaY == false) {
02691 working_space[shift + j] = fSigmaInitY;
02692 j += 1;
02693 }
02694 working_space[7 * i + 2] = fRoInit;
02695 if (fFixRo == false) {
02696 working_space[shift + j] = fRoInit;
02697 j += 1;
02698 }
02699 working_space[7 * i + 3] = fA0Init;
02700 if (fFixA0 == false) {
02701 working_space[shift + j] = fA0Init;
02702 j += 1;
02703 }
02704 working_space[7 * i + 4] = fAxInit;
02705 if (fFixAx == false) {
02706 working_space[shift + j] = fAxInit;
02707 j += 1;
02708 }
02709 working_space[7 * i + 5] = fAyInit;
02710 if (fFixAy == false) {
02711 working_space[shift + j] = fAyInit;
02712 j += 1;
02713 }
02714 working_space[7 * i + 6] = fTxyInit;
02715 if (fFixTxy == false) {
02716 working_space[shift + j] = fTxyInit;
02717 j += 1;
02718 }
02719 working_space[7 * i + 7] = fSxyInit;
02720 if (fFixSxy == false) {
02721 working_space[shift + j] = fSxyInit;
02722 j += 1;
02723 }
02724 working_space[7 * i + 8] = fTxInit;
02725 if (fFixTx == false) {
02726 working_space[shift + j] = fTxInit;
02727 j += 1;
02728 }
02729 working_space[7 * i + 9] = fTyInit;
02730 if (fFixTy == false) {
02731 working_space[shift + j] = fTyInit;
02732 j += 1;
02733 }
02734 working_space[7 * i + 10] = fSxyInit;
02735 if (fFixSx == false) {
02736 working_space[shift + j] = fSxInit;
02737 j += 1;
02738 }
02739 working_space[7 * i + 11] = fSyInit;
02740 if (fFixSy == false) {
02741 working_space[shift + j] = fSyInit;
02742 j += 1;
02743 }
02744 working_space[7 * i + 12] = fBxInit;
02745 if (fFixBx == false) {
02746 working_space[shift + j] = fBxInit;
02747 j += 1;
02748 }
02749 working_space[7 * i + 13] = fByInit;
02750 if (fFixBy == false) {
02751 working_space[shift + j] = fByInit;
02752 j += 1;
02753 }
02754 size = j;
02755 for (iter = 0; iter < fNumberIterations; iter++) {
02756 for (j = 0; j < size; j++) {
02757 working_space[2 * shift + j] = 0, working_space[3 * shift + j] = 0;
02758 }
02759
02760
02761 alpha = fAlpha;
02762 chi_opt = 0, pw = fPower - 2;
02763 for (i1 = fXmin; i1 <= fXmax; i1++) {
02764 for (i2 = fYmin; i2 <= fYmax; i2++) {
02765 yw = source[i1][i2];
02766 ywm = yw;
02767 f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
02768 working_space, working_space[peak_vel],
02769 working_space[peak_vel + 1],
02770 working_space[peak_vel + 2],
02771 working_space[peak_vel + 3],
02772 working_space[peak_vel + 4],
02773 working_space[peak_vel + 5],
02774 working_space[peak_vel + 6],
02775 working_space[peak_vel + 7],
02776 working_space[peak_vel + 8],
02777 working_space[peak_vel + 9],
02778 working_space[peak_vel + 10],
02779 working_space[peak_vel + 11],
02780 working_space[peak_vel + 12],
02781 working_space[peak_vel + 13]);
02782 if (fStatisticType == kFitOptimMaxLikelihood) {
02783 if (f > 0.00001)
02784 chi_opt += yw * TMath::Log(f) - f;
02785 }
02786
02787 else {
02788 if (ywm != 0)
02789 chi_opt += (yw - f) * (yw - f) / ywm;
02790 }
02791 if (fStatisticType == kFitOptimChiFuncValues) {
02792 ywm = f;
02793 if (f < 0.00001)
02794 ywm = 0.00001;
02795 }
02796
02797 else if (fStatisticType == kFitOptimMaxLikelihood) {
02798 ywm = f;
02799 if (f < 0.00001)
02800 ywm = 0.00001;
02801 }
02802
02803 else {
02804 if (ywm == 0)
02805 ywm = 1;
02806 }
02807
02808
02809 for (j = 0, k = 0; j < fNPeaks; j++) {
02810 if (fFixAmp[j] == false) {
02811 a = Deramp2((Double_t) i1, (Double_t) i2,
02812 working_space[7 * j + 1],
02813 working_space[7 * j + 2],
02814 working_space[peak_vel],
02815 working_space[peak_vel + 1],
02816 working_space[peak_vel + 2],
02817 working_space[peak_vel + 6],
02818 working_space[peak_vel + 7],
02819 working_space[peak_vel + 12],
02820 working_space[peak_vel + 13]);
02821 if (ywm != 0) {
02822 c = Ourpowl(a, pw);
02823 if (fStatisticType == kFitOptimChiFuncValues) {
02824 b = a * (yw * yw - f * f) / (ywm * ywm);
02825 working_space[2 * shift + k] += b * c;
02826 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
02827 working_space[3 * shift + k] += b * c;
02828 }
02829
02830 else {
02831 b = a * (yw - f) / ywm;
02832 working_space[2 * shift + k] += b * c;
02833 b = a * a / ywm;
02834 working_space[3 * shift + k] += b * c;
02835 }
02836 }
02837 k += 1;
02838 }
02839 if (fFixPositionX[j] == false) {
02840 a = Deri02((Double_t) i1, (Double_t) i2,
02841 working_space[7 * j],
02842 working_space[7 * j + 1],
02843 working_space[7 * j + 2],
02844 working_space[peak_vel],
02845 working_space[peak_vel + 1],
02846 working_space[peak_vel + 2],
02847 working_space[peak_vel + 6],
02848 working_space[peak_vel + 7],
02849 working_space[peak_vel + 12],
02850 working_space[peak_vel + 13]);
02851 if (fFitTaylor == kFitTaylorOrderSecond)
02852 d = Derderi02((Double_t) i1, (Double_t) i2,
02853 working_space[7 * j],
02854 working_space[7 * j + 1],
02855 working_space[7 * j + 2],
02856 working_space[peak_vel],
02857 working_space[peak_vel + 1],
02858 working_space[peak_vel + 2]);
02859 if (ywm != 0) {
02860 c = Ourpowl(a, pw);
02861 if (TMath::Abs(a) > 0.00000001
02862 && fFitTaylor == kFitTaylorOrderSecond) {
02863 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
02864 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
02865 && a <= 0))
02866 d = 0;
02867 }
02868
02869 else
02870 d = 0;
02871 a = a + d;
02872 if (fStatisticType == kFitOptimChiFuncValues) {
02873 b = a * (yw * yw - f * f) / (ywm * ywm);
02874 working_space[2 * shift + k] += b * c;
02875 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
02876 working_space[3 * shift + k] += b * c;
02877 }
02878
02879 else {
02880 b = a * (yw - f) / ywm;
02881 working_space[2 * shift + k] += b * c;
02882 b = a * a / ywm;
02883 working_space[3 * shift + k] += b * c;
02884 }
02885 }
02886 k += 1;
02887 }
02888 if (fFixPositionY[j] == false) {
02889 a = Derj02((Double_t) i1, (Double_t) i2,
02890 working_space[7 * j],
02891 working_space[7 * j + 1],
02892 working_space[7 * j + 2],
02893 working_space[peak_vel],
02894 working_space[peak_vel + 1],
02895 working_space[peak_vel + 2],
02896 working_space[peak_vel + 6],
02897 working_space[peak_vel + 7],
02898 working_space[peak_vel + 12],
02899 working_space[peak_vel + 13]);
02900 if (fFitTaylor == kFitTaylorOrderSecond)
02901 d = Derderj02((Double_t) i1, (Double_t) i2,
02902 working_space[7 * j],
02903 working_space[7 * j + 1],
02904 working_space[7 * j + 2],
02905 working_space[peak_vel],
02906 working_space[peak_vel + 1],
02907 working_space[peak_vel + 2]);
02908 if (ywm != 0) {
02909 c = Ourpowl(a, pw);
02910 if (TMath::Abs(a) > 0.00000001
02911 && fFitTaylor == kFitTaylorOrderSecond) {
02912 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
02913 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
02914 && a <= 0))
02915 d = 0;
02916 }
02917
02918 else
02919 d = 0;
02920 a = a + d;
02921 if (fStatisticType == kFitOptimChiFuncValues) {
02922 b = a * (yw * yw - f * f) / (ywm * ywm);
02923 working_space[2 * shift + k] += b * c;
02924 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
02925 working_space[3 * shift + k] += b * c;
02926 }
02927
02928 else {
02929 b = a * (yw - f) / ywm;
02930 working_space[2 * shift + k] += b * c;
02931 b = a * a / ywm;
02932 working_space[3 * shift + k] += b * c;
02933 }
02934 }
02935 k += 1;
02936 }
02937 if (fFixAmpX1[j] == false) {
02938 a = Derampx((Double_t) i1, working_space[7 * j + 5],
02939 working_space[peak_vel],
02940 working_space[peak_vel + 8],
02941 working_space[peak_vel + 10],
02942 working_space[peak_vel + 12]);
02943 if (ywm != 0) {
02944 c = Ourpowl(a, pw);
02945 if (fStatisticType == kFitOptimChiFuncValues) {
02946 b = a * (yw * yw - f * f) / (ywm * ywm);
02947 working_space[2 * shift + k] += b * c;
02948 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
02949 working_space[3 * shift + k] += b * c;
02950 }
02951
02952 else {
02953 b = a * (yw - f) / ywm;
02954 working_space[2 * shift + k] += b * c;
02955 b = a * a / ywm;
02956 working_space[3 * shift + k] += b * c;
02957 }
02958 }
02959 k += 1;
02960 }
02961 if (fFixAmpY1[j] == false) {
02962 a = Derampx((Double_t) i2, working_space[7 * j + 6],
02963 working_space[peak_vel + 1],
02964 working_space[peak_vel + 9],
02965 working_space[peak_vel + 11],
02966 working_space[peak_vel + 13]);
02967 if (ywm != 0) {
02968 c = Ourpowl(a, pw);
02969 if (fStatisticType == kFitOptimChiFuncValues) {
02970 b = a * (yw * yw - f * f) / (ywm * ywm);
02971 working_space[2 * shift + k] += b * c;
02972 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
02973 working_space[3 * shift + k] += b * c;
02974 }
02975
02976 else {
02977 b = a * (yw - f) / ywm;
02978 working_space[2 * shift + k] += b * c;
02979 b = a * a / ywm;
02980 working_space[3 * shift + k] += b * c;
02981 }
02982 }
02983 k += 1;
02984 }
02985 if (fFixPositionX1[j] == false) {
02986 a = Deri01((Double_t) i1, working_space[7 * j + 3],
02987 working_space[7 * j + 5],
02988 working_space[peak_vel],
02989 working_space[peak_vel + 8],
02990 working_space[peak_vel + 10],
02991 working_space[peak_vel + 12]);
02992 if (fFitTaylor == kFitTaylorOrderSecond)
02993 d = Derderi01((Double_t) i1, working_space[7 * j + 3],
02994 working_space[7 * j + 5],
02995 working_space[peak_vel]);
02996 if (ywm != 0) {
02997 c = Ourpowl(a, pw);
02998 if (TMath::Abs(a) > 0.00000001
02999 && fFitTaylor == kFitTaylorOrderSecond) {
03000 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
03001 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
03002 && a <= 0))
03003 d = 0;
03004 }
03005
03006 else
03007 d = 0;
03008 a = a + d;
03009 if (fStatisticType == kFitOptimChiFuncValues) {
03010 b = a * (yw * yw - f * f) / (ywm * ywm);
03011 working_space[2 * shift + k] += b * c;
03012 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03013 working_space[3 * shift + k] += b * c;
03014 }
03015
03016 else {
03017 b = a * (yw - f) / ywm;
03018 working_space[2 * shift + k] += b * c;
03019 b = a * a / ywm;
03020 working_space[3 * shift + k] += b * c;
03021 }
03022 }
03023 k += 1;
03024 }
03025 if (fFixPositionY1[j] == false) {
03026 a = Deri01((Double_t) i2, working_space[7 * j + 4],
03027 working_space[7 * j + 6],
03028 working_space[peak_vel + 1],
03029 working_space[peak_vel + 9],
03030 working_space[peak_vel + 11],
03031 working_space[peak_vel + 13]);
03032 if (fFitTaylor == kFitTaylorOrderSecond)
03033 d = Derderi01((Double_t) i2, working_space[7 * j + 4],
03034 working_space[7 * j + 6],
03035 working_space[peak_vel + 1]);
03036 if (ywm != 0) {
03037 c = Ourpowl(a, pw);
03038 if (TMath::Abs(a) > 0.00000001
03039 && fFitTaylor == kFitTaylorOrderSecond) {
03040 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
03041 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
03042 && a <= 0))
03043 d = 0;
03044 }
03045
03046 else
03047 d = 0;
03048 a = a + d;
03049 if (fStatisticType == kFitOptimChiFuncValues) {
03050 b = a * (yw * yw - f * f) / (ywm * ywm);
03051 working_space[2 * shift + k] += b * c;
03052 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03053 working_space[3 * shift + k] += b * c;
03054 }
03055
03056 else {
03057 b = a * (yw - f) / ywm;
03058 working_space[2 * shift + k] += b * c;
03059 b = a * a / ywm;
03060 working_space[3 * shift + k] += b * c;
03061 }
03062 }
03063 k += 1;
03064 }
03065 }
03066 if (fFixSigmaX == false) {
03067 a = Dersigmax(fNPeaks, (Double_t) i1, (Double_t) i2,
03068 working_space, working_space[peak_vel],
03069 working_space[peak_vel + 1],
03070 working_space[peak_vel + 2],
03071 working_space[peak_vel + 6],
03072 working_space[peak_vel + 7],
03073 working_space[peak_vel + 8],
03074 working_space[peak_vel + 10],
03075 working_space[peak_vel + 12],
03076 working_space[peak_vel + 13]);
03077 if (fFitTaylor == kFitTaylorOrderSecond)
03078 d = Derdersigmax(fNPeaks, (Double_t) i1,
03079 (Double_t) i2, working_space,
03080 working_space[peak_vel],
03081 working_space[peak_vel + 1],
03082 working_space[peak_vel + 2]);
03083 if (ywm != 0) {
03084 c = Ourpowl(a, pw);
03085 if (TMath::Abs(a) > 0.00000001
03086 && fFitTaylor == kFitTaylorOrderSecond) {
03087 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
03088 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
03089 d = 0;
03090 }
03091
03092 else
03093 d = 0;
03094 a = a + d;
03095 if (fStatisticType == kFitOptimChiFuncValues) {
03096 b = a * (yw * yw - f * f) / (ywm * ywm);
03097 working_space[2 * shift + k] += b * c;
03098 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03099 working_space[3 * shift + k] += b * c;
03100 }
03101
03102 else {
03103 b = a * (yw - f) / ywm;
03104 working_space[2 * shift + k] += b * c;
03105 b = a * a / ywm;
03106 working_space[3 * shift + k] += b * c;
03107 }
03108 }
03109 k += 1;
03110 }
03111 if (fFixSigmaY == false) {
03112 a = Dersigmay(fNPeaks, (Double_t) i1, (Double_t) i2,
03113 working_space, working_space[peak_vel],
03114 working_space[peak_vel + 1],
03115 working_space[peak_vel + 2],
03116 working_space[peak_vel + 6],
03117 working_space[peak_vel + 7],
03118 working_space[peak_vel + 9],
03119 working_space[peak_vel + 11],
03120 working_space[peak_vel + 12],
03121 working_space[peak_vel + 13]);
03122 if (fFitTaylor == kFitTaylorOrderSecond)
03123 d = Derdersigmay(fNPeaks, (Double_t) i1,
03124 (Double_t) i2, working_space,
03125 working_space[peak_vel],
03126 working_space[peak_vel + 1],
03127 working_space[peak_vel + 2]);
03128 if (ywm != 0) {
03129 c = Ourpowl(a, pw);
03130 if (TMath::Abs(a) > 0.00000001
03131 && fFitTaylor == kFitTaylorOrderSecond) {
03132 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
03133 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
03134 d = 0;
03135 }
03136
03137 else
03138 d = 0;
03139 a = a + d;
03140 if (fStatisticType == kFitOptimChiFuncValues) {
03141 b = a * (yw * yw - f * f) / (ywm * ywm);
03142 working_space[2 * shift + k] += b * c;
03143 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03144 working_space[3 * shift + k] += b * c;
03145 }
03146
03147 else {
03148 b = a * (yw - f) / ywm;
03149 working_space[2 * shift + k] += b * c;
03150 b = a * a / ywm;
03151 working_space[3 * shift + k] += b * c;
03152 }
03153 }
03154 k += 1;
03155 }
03156 if (fFixRo == false) {
03157 a = Derro(fNPeaks, (Double_t) i1, (Double_t) i2,
03158 working_space, working_space[peak_vel],
03159 working_space[peak_vel + 1],
03160 working_space[peak_vel + 2]);
03161 if (ywm != 0) {
03162 c = Ourpowl(a, pw);
03163 if (TMath::Abs(a) > 0.00000001
03164 && fFitTaylor == kFitTaylorOrderSecond) {
03165 d = d * TMath::Abs(yw - f) / (2 * a * ywm);
03166 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
03167 d = 0;
03168 }
03169
03170 else
03171 d = 0;
03172 a = a + d;
03173 if (fStatisticType == kFitOptimChiFuncValues) {
03174 b = a * (yw * yw - f * f) / (ywm * ywm);
03175 working_space[2 * shift + k] += b * c;
03176 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03177 working_space[3 * shift + k] += b * c;
03178 }
03179
03180 else {
03181 b = a * (yw - f) / ywm;
03182 working_space[2 * shift + k] += b * c;
03183 b = a * a / ywm;
03184 working_space[3 * shift + k] += b * c;
03185 }
03186 }
03187 k += 1;
03188 }
03189 if (fFixA0 == false) {
03190 a = 1.;
03191 if (ywm != 0) {
03192 c = Ourpowl(a, pw);
03193 if (fStatisticType == kFitOptimChiFuncValues) {
03194 b = a * (yw * yw - f * f) / (ywm * ywm);
03195 working_space[2 * shift + k] += b * c;
03196 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03197 working_space[3 * shift + k] += b * c;
03198 }
03199
03200 else {
03201 b = a * (yw - f) / ywm;
03202 working_space[2 * shift + k] += b * c;
03203 b = a * a / ywm;
03204 working_space[3 * shift + k] += b * c;
03205 }
03206 }
03207 k += 1;
03208 }
03209 if (fFixAx == false) {
03210 a = i1;
03211 if (ywm != 0) {
03212 c = Ourpowl(a, pw);
03213 if (fStatisticType == kFitOptimChiFuncValues) {
03214 b = a * (yw * yw - f * f) / (ywm * ywm);
03215 working_space[2 * shift + k] += b * c;
03216 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03217 working_space[3 * shift + k] += b * c;
03218 }
03219
03220 else {
03221 b = a * (yw - f) / ywm;
03222 working_space[2 * shift + k] += b * c;
03223 b = a * a / ywm;
03224 working_space[3 * shift + k] += b * c;
03225 }
03226 }
03227 k += 1;
03228 }
03229 if (fFixAy == false) {
03230 a = i2;
03231 if (ywm != 0) {
03232 c = Ourpowl(a, pw);
03233 if (fStatisticType == kFitOptimChiFuncValues) {
03234 b = a * (yw * yw - f * f) / (ywm * ywm);
03235 working_space[2 * shift + k] += b * c;
03236 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03237 working_space[3 * shift + k] += b * c;
03238 }
03239
03240 else {
03241 b = a * (yw - f) / ywm;
03242 working_space[2 * shift + k] += b * c;
03243 b = a * a / ywm;
03244 working_space[3 * shift + k] += b * c;
03245 }
03246 }
03247 k += 1;
03248 }
03249 if (fFixTxy == false) {
03250 a = Dertxy(fNPeaks, (Double_t) i1, (Double_t) i2,
03251 working_space, working_space[peak_vel],
03252 working_space[peak_vel + 1],
03253 working_space[peak_vel + 12],
03254 working_space[peak_vel + 13]);
03255 if (ywm != 0) {
03256 c = Ourpowl(a, pw);
03257 if (fStatisticType == kFitOptimChiFuncValues) {
03258 b = a * (yw * yw - f * f) / (ywm * ywm);
03259 working_space[2 * shift + k] += b * c;
03260 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03261 working_space[3 * shift + k] += b * c;
03262 }
03263
03264 else {
03265 b = a * (yw - f) / ywm;
03266 working_space[2 * shift + k] += b * c;
03267 b = a * a / ywm;
03268 working_space[3 * shift + k] += b * c;
03269 }
03270 }
03271 k += 1;
03272 }
03273 if (fFixSxy == false) {
03274 a = Dersxy(fNPeaks, (Double_t) i1, (Double_t) i2,
03275 working_space, working_space[peak_vel],
03276 working_space[peak_vel + 1]);
03277 if (ywm != 0) {
03278 c = Ourpowl(a, pw);
03279 if (fStatisticType == kFitOptimChiFuncValues) {
03280 b = a * (yw * yw - f * f) / (ywm * ywm);
03281 working_space[2 * shift + k] += b * c;
03282 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03283 working_space[3 * shift + k] += b * c;
03284 }
03285
03286 else {
03287 b = a * (yw - f) / ywm;
03288 working_space[2 * shift + k] += b * c;
03289 b = a * a / ywm;
03290 working_space[3 * shift + k] += b * c;
03291 }
03292 }
03293 k += 1;
03294 }
03295 if (fFixTx == false) {
03296 a = Dertx(fNPeaks, (Double_t) i1, working_space,
03297 working_space[peak_vel],
03298 working_space[peak_vel + 12]);
03299 if (ywm != 0) {
03300 c = Ourpowl(a, pw);
03301 if (fStatisticType == kFitOptimChiFuncValues) {
03302 b = a * (yw * yw - f * f) / (ywm * ywm);
03303 working_space[2 * shift + k] += b * c;
03304 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03305 working_space[3 * shift + k] += b * c;
03306 }
03307
03308 else {
03309 b = a * (yw - f) / ywm;
03310 working_space[2 * shift + k] += b * c;
03311 b = a * a / ywm;
03312 working_space[3 * shift + k] += b * c;
03313 }
03314 }
03315 k += 1;
03316 }
03317 if (fFixTy == false) {
03318 a = Derty(fNPeaks, (Double_t) i2, working_space,
03319 working_space[peak_vel + 1],
03320 working_space[peak_vel + 13]);
03321 if (ywm != 0) {
03322 c = Ourpowl(a, pw);
03323 if (fStatisticType == kFitOptimChiFuncValues) {
03324 b = a * (yw * yw - f * f) / (ywm * ywm);
03325 working_space[2 * shift + k] += b * c;
03326 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03327 working_space[3 * shift + k] += b * c;
03328 }
03329
03330 else {
03331 b = a * (yw - f) / ywm;
03332 working_space[2 * shift + k] += b * c;
03333 b = a * a / ywm;
03334 working_space[3 * shift + k] += b * c;
03335 }
03336 }
03337 k += 1;
03338 }
03339 if (fFixSx == false) {
03340 a = Dersx(fNPeaks, (Double_t) i1, working_space,
03341 working_space[peak_vel]);
03342 if (ywm != 0) {
03343 c = Ourpowl(a, pw);
03344 if (fStatisticType == kFitOptimChiFuncValues) {
03345 b = a * (yw * yw - f * f) / (ywm * ywm);
03346 working_space[2 * shift + k] += b * c;
03347 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03348 working_space[3 * shift + k] += b * c;
03349 }
03350
03351 else {
03352 b = a * (yw - f) / ywm;
03353 working_space[2 * shift + k] += b * c;
03354 b = a * a / ywm;
03355 working_space[3 * shift + k] += b * c;
03356 }
03357 }
03358 k += 1;
03359 }
03360 if (fFixSy == false) {
03361 a = Dersy(fNPeaks, (Double_t) i2, working_space,
03362 working_space[peak_vel + 1]);
03363 if (ywm != 0) {
03364 c = Ourpowl(a, pw);
03365 if (fStatisticType == kFitOptimChiFuncValues) {
03366 b = a * (yw * yw - f * f) / (ywm * ywm);
03367 working_space[2 * shift + k] += b * c;
03368 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03369 working_space[3 * shift + k] += b * c;
03370 }
03371
03372 else {
03373 b = a * (yw - f) / ywm;
03374 working_space[2 * shift + k] += b * c;
03375 b = a * a / ywm;
03376 working_space[3 * shift + k] += b * c;
03377 }
03378 }
03379 k += 1;
03380 }
03381 if (fFixBx == false) {
03382 a = Derbx(fNPeaks, (Double_t) i1, (Double_t) i2,
03383 working_space, working_space[peak_vel],
03384 working_space[peak_vel + 1],
03385 working_space[peak_vel + 6],
03386 working_space[peak_vel + 8],
03387 working_space[peak_vel + 12],
03388 working_space[peak_vel + 13]);
03389 if (ywm != 0) {
03390 c = Ourpowl(a, pw);
03391 if (fStatisticType == kFitOptimChiFuncValues) {
03392 b = a * (yw * yw - f * f) / (ywm * ywm);
03393 working_space[2 * shift + k] += b * c;
03394 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03395 working_space[3 * shift + k] += b * c;
03396 }
03397
03398 else {
03399 b = a * (yw - f) / ywm;
03400 working_space[2 * shift + k] += b * c;
03401 b = a * a / ywm;
03402 working_space[3 * shift + k] += b * c;
03403 }
03404 }
03405 k += 1;
03406 }
03407 if (fFixBy == false) {
03408 a = Derby(fNPeaks, (Double_t) i1, (Double_t) i2,
03409 working_space, working_space[peak_vel],
03410 working_space[peak_vel + 1],
03411 working_space[peak_vel + 6],
03412 working_space[peak_vel + 8],
03413 working_space[peak_vel + 12],
03414 working_space[peak_vel + 13]);
03415 if (ywm != 0) {
03416 c = Ourpowl(a, pw);
03417 if (fStatisticType == kFitOptimChiFuncValues) {
03418 b = a * (yw * yw - f * f) / (ywm * ywm);
03419 working_space[2 * shift + k] += b * c;
03420 b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
03421 working_space[3 * shift + k] += b * c;
03422 }
03423
03424 else {
03425 b = a * (yw - f) / ywm;
03426 working_space[2 * shift + k] += b * c;
03427 b = a * a / ywm;
03428 working_space[3 * shift + k] += b * c;
03429 }
03430 }
03431 k += 1;
03432 }
03433 }
03434 }
03435 for (j = 0; j < size; j++) {
03436 if (TMath::Abs(working_space[3 * shift + j]) > 0.000001)
03437 working_space[2 * shift + j] = working_space[2 * shift + j] / TMath::Abs(working_space[3 * shift + j]);
03438 else
03439 working_space[2 * shift + j] = 0;
03440 }
03441
03442
03443 chi2 = chi_opt;
03444 chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
03445
03446
03447 regul_cycle = 0;
03448 for (j = 0; j < size; j++) {
03449 working_space[4 * shift + j] = working_space[shift + j];
03450 }
03451
03452 do {
03453 if (fAlphaOptim == kFitAlphaOptimal) {
03454 if (fStatisticType != kFitOptimMaxLikelihood)
03455 chi_min = 10000 * chi2;
03456
03457 else
03458 chi_min = 0.1 * chi2;
03459 flag = 0;
03460 for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
03461 for (j = 0; j < size; j++) {
03462 working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j];
03463 }
03464 for (i = 0, j = 0; i < fNPeaks; i++) {
03465 if (fFixAmp[i] == false) {
03466 if (working_space[shift + j] < 0)
03467 working_space[shift + j] = 0;
03468 working_space[7 * i] = working_space[shift + j];
03469 j += 1;
03470 }
03471 if (fFixPositionX[i] == false) {
03472 if (working_space[shift + j] < fXmin)
03473 working_space[shift + j] = fXmin;
03474 if (working_space[shift + j] > fXmax)
03475 working_space[shift + j] = fXmax;
03476 working_space[7 * i + 1] = working_space[shift + j];
03477 j += 1;
03478 }
03479 if (fFixPositionY[i] == false) {
03480 if (working_space[shift + j] < fYmin)
03481 working_space[shift + j] = fYmin;
03482 if (working_space[shift + j] > fYmax)
03483 working_space[shift + j] = fYmax;
03484 working_space[7 * i + 2] = working_space[shift + j];
03485 j += 1;
03486 }
03487 if (fFixAmpX1[i] == false) {
03488 if (working_space[shift + j] < 0)
03489 working_space[shift + j] = 0;
03490 working_space[7 * i + 3] = working_space[shift + j];
03491 j += 1;
03492 }
03493 if (fFixAmpY1[i] == false) {
03494 if (working_space[shift + j] < 0)
03495 working_space[shift + j] = 0;
03496 working_space[7 * i + 4] = working_space[shift + j];
03497 j += 1;
03498 }
03499 if (fFixPositionX1[i] == false) {
03500 if (working_space[shift + j] < fXmin)
03501 working_space[shift + j] = fXmin;
03502 if (working_space[shift + j] > fXmax)
03503 working_space[shift + j] = fXmax;
03504 working_space[7 * i + 5] = working_space[shift + j];
03505 j += 1;
03506 }
03507 if (fFixPositionY1[i] == false) {
03508 if (working_space[shift + j] < fYmin)
03509 working_space[shift + j] = fYmin;
03510 if (working_space[shift + j] > fYmax)
03511 working_space[shift + j] = fYmax;
03512 working_space[7 * i + 6] = working_space[shift + j];
03513 j += 1;
03514 }
03515 }
03516 if (fFixSigmaX == false) {
03517 if (working_space[shift + j] < 0.001) {
03518 working_space[shift + j] = 0.001;
03519 }
03520 working_space[peak_vel] = working_space[shift + j];
03521 j += 1;
03522 }
03523 if (fFixSigmaY == false) {
03524 if (working_space[shift + j] < 0.001) {
03525 working_space[shift + j] = 0.001;
03526 }
03527 working_space[peak_vel + 1] = working_space[shift + j];
03528 j += 1;
03529 }
03530 if (fFixRo == false) {
03531 if (working_space[shift + j] < -1) {
03532 working_space[shift + j] = -1;
03533 }
03534 if (working_space[shift + j] > 1) {
03535 working_space[shift + j] = 1;
03536 }
03537 working_space[peak_vel + 2] = working_space[shift + j];
03538 j += 1;
03539 }
03540 if (fFixA0 == false) {
03541 working_space[peak_vel + 3] = working_space[shift + j];
03542 j += 1;
03543 }
03544 if (fFixAx == false) {
03545 working_space[peak_vel + 4] = working_space[shift + j];
03546 j += 1;
03547 }
03548 if (fFixAy == false) {
03549 working_space[peak_vel + 5] = working_space[shift + j];
03550 j += 1;
03551 }
03552 if (fFixTxy == false) {
03553 working_space[peak_vel + 6] = working_space[shift + j];
03554 j += 1;
03555 }
03556 if (fFixSxy == false) {
03557 working_space[peak_vel + 7] = working_space[shift + j];
03558 j += 1;
03559 }
03560 if (fFixTx == false) {
03561 working_space[peak_vel + 8] = working_space[shift + j];
03562 j += 1;
03563 }
03564 if (fFixTy == false) {
03565 working_space[peak_vel + 9] = working_space[shift + j];
03566 j += 1;
03567 }
03568 if (fFixSx == false) {
03569 working_space[peak_vel + 10] = working_space[shift + j];
03570 j += 1;
03571 }
03572 if (fFixSy == false) {
03573 working_space[peak_vel + 11] = working_space[shift + j];
03574 j += 1;
03575 }
03576 if (fFixBx == false) {
03577 if (TMath::Abs(working_space[shift + j]) < 0.001) {
03578 if (working_space[shift + j] < 0)
03579 working_space[shift + j] = -0.001;
03580 else
03581 working_space[shift + j] = 0.001;
03582 }
03583 working_space[peak_vel + 12] = working_space[shift + j];
03584 j += 1;
03585 }
03586 if (fFixBy == false) {
03587 if (TMath::Abs(working_space[shift + j]) < 0.001) {
03588 if (working_space[shift + j] < 0)
03589 working_space[shift + j] = -0.001;
03590 else
03591 working_space[shift + j] = 0.001;
03592 }
03593 working_space[peak_vel + 13] = working_space[shift + j];
03594 j += 1;
03595 }
03596 chi2 = 0;
03597 for (i1 = fXmin; i1 <= fXmax; i1++) {
03598 for (i2 = fYmin; i2 <= fYmax; i2++) {
03599 yw = source[i1][i2];
03600 ywm = yw;
03601 f = Shape2(fNPeaks, (Double_t) i1,
03602 (Double_t) i2, working_space,
03603 working_space[peak_vel],
03604 working_space[peak_vel + 1],
03605 working_space[peak_vel + 2],
03606 working_space[peak_vel + 3],
03607 working_space[peak_vel + 4],
03608 working_space[peak_vel + 5],
03609 working_space[peak_vel + 6],
03610 working_space[peak_vel + 7],
03611 working_space[peak_vel + 8],
03612 working_space[peak_vel + 9],
03613 working_space[peak_vel + 10],
03614 working_space[peak_vel + 11],
03615 working_space[peak_vel + 12],
03616 working_space[peak_vel + 13]);
03617 if (fStatisticType == kFitOptimChiFuncValues) {
03618 ywm = f;
03619 if (f < 0.00001)
03620 ywm = 0.00001;
03621 }
03622 if (fStatisticType == kFitOptimMaxLikelihood) {
03623 if (f > 0.00001)
03624 chi2 += yw * TMath::Log(f) - f;
03625 }
03626
03627 else {
03628 if (ywm != 0)
03629 chi2 += (yw - f) * (yw - f) / ywm;
03630 }
03631 }
03632 }
03633 if ((chi2 < chi_min
03634 && fStatisticType != kFitOptimMaxLikelihood)
03635 || (chi2 > chi_min
03636 && fStatisticType == kFitOptimMaxLikelihood)) {
03637 pmin = pi, chi_min = chi2;
03638 }
03639
03640 else
03641 flag = 1;
03642 if (pi == 0.1)
03643 chi_min = chi2;
03644 chi = chi_min;
03645 }
03646 if (pmin != 0.1) {
03647 for (j = 0; j < size; j++) {
03648 working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j];
03649 }
03650 for (i = 0, j = 0; i < fNPeaks; i++) {
03651 if (fFixAmp[i] == false) {
03652 if (working_space[shift + j] < 0)
03653 working_space[shift + j] = 0;
03654 working_space[7 * i] = working_space[shift + j];
03655 j += 1;
03656 }
03657 if (fFixPositionX[i] == false) {
03658 if (working_space[shift + j] < fXmin)
03659 working_space[shift + j] = fXmin;
03660 if (working_space[shift + j] > fXmax)
03661 working_space[shift + j] = fXmax;
03662 working_space[7 * i + 1] = working_space[shift + j];
03663 j += 1;
03664 }
03665 if (fFixPositionY[i] == false) {
03666 if (working_space[shift + j] < fYmin)
03667 working_space[shift + j] = fYmin;
03668 if (working_space[shift + j] > fYmax)
03669 working_space[shift + j] = fYmax;
03670 working_space[7 * i + 2] = working_space[shift + j];
03671 j += 1;
03672 }
03673 if (fFixAmpX1[i] == false) {
03674 if (working_space[shift + j] < 0)
03675 working_space[shift + j] = 0;
03676 working_space[7 * i + 3] = working_space[shift + j];
03677 j += 1;
03678 }
03679 if (fFixAmpY1[i] == false) {
03680 if (working_space[shift + j] < 0)
03681 working_space[shift + j] = 0;
03682 working_space[7 * i + 4] = working_space[shift + j];
03683 j += 1;
03684 }
03685 if (fFixPositionX1[i] == false) {
03686 if (working_space[shift + j] < fXmin)
03687 working_space[shift + j] = fXmin;
03688 if (working_space[shift + j] > fXmax)
03689 working_space[shift + j] = fXmax;
03690 working_space[7 * i + 5] = working_space[shift + j];
03691 j += 1;
03692 }
03693 if (fFixPositionY1[i] == false) {
03694 if (working_space[shift + j] < fYmin)
03695 working_space[shift + j] = fYmin;
03696 if (working_space[shift + j] > fYmax)
03697 working_space[shift + j] = fYmax;
03698 working_space[7 * i + 6] = working_space[shift + j];
03699 j += 1;
03700 }
03701 }
03702 if (fFixSigmaX == false) {
03703 if (working_space[shift + j] < 0.001) {
03704 working_space[shift + j] = 0.001;
03705 }
03706 working_space[peak_vel] = working_space[shift + j];
03707 j += 1;
03708 }
03709 if (fFixSigmaY == false) {
03710 if (working_space[shift + j] < 0.001) {
03711 working_space[shift + j] = 0.001;
03712 }
03713 working_space[peak_vel + 1] = working_space[shift + j];
03714 j += 1;
03715 }
03716 if (fFixRo == false) {
03717 if (working_space[shift + j] < -1) {
03718 working_space[shift + j] = -1;
03719 }
03720 if (working_space[shift + j] > 1) {
03721 working_space[shift + j] = 1;
03722 }
03723 working_space[peak_vel + 2] = working_space[shift + j];
03724 j += 1;
03725 }
03726 if (fFixA0 == false) {
03727 working_space[peak_vel + 3] = working_space[shift + j];
03728 j += 1;
03729 }
03730 if (fFixAx == false) {
03731 working_space[peak_vel + 4] = working_space[shift + j];
03732 j += 1;
03733 }
03734 if (fFixAy == false) {
03735 working_space[peak_vel + 5] = working_space[shift + j];
03736 j += 1;
03737 }
03738 if (fFixTxy == false) {
03739 working_space[peak_vel + 6] = working_space[shift + j];
03740 j += 1;
03741 }
03742 if (fFixSxy == false) {
03743 working_space[peak_vel + 7] = working_space[shift + j];
03744 j += 1;
03745 }
03746 if (fFixTx == false) {
03747 working_space[peak_vel + 8] = working_space[shift + j];
03748 j += 1;
03749 }
03750 if (fFixTy == false) {
03751 working_space[peak_vel + 9] = working_space[shift + j];
03752 j += 1;
03753 }
03754 if (fFixSx == false) {
03755 working_space[peak_vel + 10] = working_space[shift + j];
03756 j += 1;
03757 }
03758 if (fFixSy == false) {
03759 working_space[peak_vel + 11] = working_space[shift + j];
03760 j += 1;
03761 }
03762 if (fFixBx == false) {
03763 if (TMath::Abs(working_space[shift + j]) < 0.001) {
03764 if (working_space[shift + j] < 0)
03765 working_space[shift + j] = -0.001;
03766 else
03767 working_space[shift + j] = 0.001;
03768 }
03769 working_space[peak_vel + 12] = working_space[shift + j];
03770 j += 1;
03771 }
03772 if (fFixBy == false) {
03773 if (TMath::Abs(working_space[shift + j]) < 0.001) {
03774 if (working_space[shift + j] < 0)
03775 working_space[shift + j] = -0.001;
03776 else
03777 working_space[shift + j] = 0.001;
03778 }
03779 working_space[peak_vel + 13] = working_space[shift + j];
03780 j += 1;
03781 }
03782 chi = chi_min;
03783 }
03784 }
03785
03786 else {
03787 for (j = 0; j < size; j++) {
03788 working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j];
03789 }
03790 for (i = 0, j = 0; i < fNPeaks; i++) {
03791 if (fFixAmp[i] == false) {
03792 if (working_space[shift + j] < 0)
03793 working_space[shift + j] = 0;
03794 working_space[7 * i] = working_space[shift + j];
03795 j += 1;
03796 }
03797 if (fFixPositionX[i] == false) {
03798 if (working_space[shift + j] < fXmin)
03799 working_space[shift + j] = fXmin;
03800 if (working_space[shift + j] > fXmax)
03801 working_space[shift + j] = fXmax;
03802 working_space[7 * i + 1] = working_space[shift + j];
03803 j += 1;
03804 }
03805 if (fFixPositionY[i] == false) {
03806 if (working_space[shift + j] < fYmin)
03807 working_space[shift + j] = fYmin;
03808 if (working_space[shift + j] > fYmax)
03809 working_space[shift + j] = fYmax;
03810 working_space[7 * i + 2] = working_space[shift + j];
03811 j += 1;
03812 }
03813 if (fFixAmpX1[i] == false) {
03814 if (working_space[shift + j] < 0)
03815 working_space[shift + j] = 0;
03816 working_space[7 * i + 3] = working_space[shift + j];
03817 j += 1;
03818 }
03819 if (fFixAmpY1[i] == false) {
03820 if (working_space[shift + j] < 0)
03821 working_space[shift + j] = 0;
03822 working_space[7 * i + 4] = working_space[shift + j];
03823 j += 1;
03824 }
03825 if (fFixPositionX1[i] == false) {
03826 if (working_space[shift + j] < fXmin)
03827 working_space[shift + j] = fXmin;
03828 if (working_space[shift + j] > fXmax)
03829 working_space[shift + j] = fXmax;
03830 working_space[7 * i + 5] = working_space[shift + j];
03831 j += 1;
03832 }
03833 if (fFixPositionY1[i] == false) {
03834 if (working_space[shift + j] < fYmin)
03835 working_space[shift + j] = fYmin;
03836 if (working_space[shift + j] > fYmax)
03837 working_space[shift + j] = fYmax;
03838 working_space[7 * i + 6] = working_space[shift + j];
03839 j += 1;
03840 }
03841 }
03842 if (fFixSigmaX == false) {
03843 if (working_space[shift + j] < 0.001) {
03844 working_space[shift + j] = 0.001;
03845 }
03846 working_space[peak_vel] = working_space[shift + j];
03847 j += 1;
03848 }
03849 if (fFixSigmaY == false) {
03850 if (working_space[shift + j] < 0.001) {
03851 working_space[shift + j] = 0.001;
03852 }
03853 working_space[peak_vel + 1] = working_space[shift + j];
03854 j += 1;
03855 }
03856 if (fFixRo == false) {
03857 if (working_space[shift + j] < -1) {
03858 working_space[shift + j] = -1;
03859 }
03860 if (working_space[shift + j] > 1) {
03861 working_space[shift + j] = 1;
03862 }
03863 working_space[peak_vel + 2] = working_space[shift + j];
03864 j += 1;
03865 }
03866 if (fFixA0 == false) {
03867 working_space[peak_vel + 3] = working_space[shift + j];
03868 j += 1;
03869 }
03870 if (fFixAx == false) {
03871 working_space[peak_vel + 4] = working_space[shift + j];
03872 j += 1;
03873 }
03874 if (fFixAy == false) {
03875 working_space[peak_vel + 5] = working_space[shift + j];
03876 j += 1;
03877 }
03878 if (fFixTxy == false) {
03879 working_space[peak_vel + 6] = working_space[shift + j];
03880 j += 1;
03881 }
03882 if (fFixSxy == false) {
03883 working_space[peak_vel + 7] = working_space[shift + j];
03884 j += 1;
03885 }
03886 if (fFixTx == false) {
03887 working_space[peak_vel + 8] = working_space[shift + j];
03888 j += 1;
03889 }
03890 if (fFixTy == false) {
03891 working_space[peak_vel + 9] = working_space[shift + j];
03892 j += 1;
03893 }
03894 if (fFixSx == false) {
03895 working_space[peak_vel + 10] = working_space[shift + j];
03896 j += 1;
03897 }
03898 if (fFixSy == false) {
03899 working_space[peak_vel + 11] = working_space[shift + j];
03900 j += 1;
03901 }
03902 if (fFixBx == false) {
03903 if (TMath::Abs(working_space[shift + j]) < 0.001) {
03904 if (working_space[shift + j] < 0)
03905 working_space[shift + j] = -0.001;
03906 else
03907 working_space[shift + j] = 0.001;
03908 }
03909 working_space[peak_vel + 12] = working_space[shift + j];
03910 j += 1;
03911 }
03912 if (fFixBy == false) {
03913 if (TMath::Abs(working_space[shift + j]) < 0.001) {
03914 if (working_space[shift + j] < 0)
03915 working_space[shift + j] = -0.001;
03916 else
03917 working_space[shift + j] = 0.001;
03918 }
03919 working_space[peak_vel + 13] = working_space[shift + j];
03920 j += 1;
03921 }
03922 chi = 0;
03923 for (i1 = fXmin; i1 <= fXmax; i1++) {
03924 for (i2 = fYmin; i2 <= fYmax; i2++) {
03925 yw = source[i1][i2];
03926 ywm = yw;
03927 f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
03928 working_space, working_space[peak_vel],
03929 working_space[peak_vel + 1],
03930 working_space[peak_vel + 2],
03931 working_space[peak_vel + 3],
03932 working_space[peak_vel + 4],
03933 working_space[peak_vel + 5],
03934 working_space[peak_vel + 6],
03935 working_space[peak_vel + 7],
03936 working_space[peak_vel + 8],
03937 working_space[peak_vel + 9],
03938 working_space[peak_vel + 10],
03939 working_space[peak_vel + 11],
03940 working_space[peak_vel + 12],
03941 working_space[peak_vel + 13]);
03942 if (fStatisticType == kFitOptimChiFuncValues) {
03943 ywm = f;
03944 if (f < 0.00001)
03945 ywm = 0.00001;
03946 }
03947 if (fStatisticType == kFitOptimMaxLikelihood) {
03948 if (f > 0.00001)
03949 chi += yw * TMath::Log(f) - f;
03950 }
03951
03952 else {
03953 if (ywm != 0)
03954 chi += (yw - f) * (yw - f) / ywm;
03955 }
03956 }
03957 }
03958 }
03959 chi2 = chi;
03960 chi = TMath::Sqrt(TMath::Abs(chi));
03961 if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
03962 alpha = alpha * chi_opt / (2 * chi);
03963
03964 else if (fAlphaOptim == kFitAlphaOptimal)
03965 alpha = alpha / 10.0;
03966 iter += 1;
03967 regul_cycle += 1;
03968 } while (((chi > chi_opt
03969 && fStatisticType != kFitOptimMaxLikelihood)
03970 || (chi < chi_opt
03971 && fStatisticType == kFitOptimMaxLikelihood))
03972 && regul_cycle < kFitNumRegulCycles);
03973 for (j = 0; j < size; j++) {
03974 working_space[4 * shift + j] = 0;
03975 working_space[2 * shift + j] = 0;
03976 }
03977 for (i1 = fXmin, chi_cel = 0; i1 <= fXmax; i1++) {
03978 for (i2 = fYmin; i2 <= fYmax; i2++) {
03979 yw = source[i1][i2];
03980 if (yw == 0)
03981 yw = 1;
03982 f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
03983 working_space, working_space[peak_vel],
03984 working_space[peak_vel + 1],
03985 working_space[peak_vel + 2],
03986 working_space[peak_vel + 3],
03987 working_space[peak_vel + 4],
03988 working_space[peak_vel + 5],
03989 working_space[peak_vel + 6],
03990 working_space[peak_vel + 7],
03991 working_space[peak_vel + 8],
03992 working_space[peak_vel + 9],
03993 working_space[peak_vel + 10],
03994 working_space[peak_vel + 11],
03995 working_space[peak_vel + 12],
03996 working_space[peak_vel + 13]);
03997 chi_opt = (yw - f) * (yw - f) / yw;
03998 chi_cel += (yw - f) * (yw - f) / yw;
03999
04000
04001 for (j = 0, k = 0; j < fNPeaks; j++) {
04002 if (fFixAmp[j] == false) {
04003 a = Deramp2((Double_t) i1, (Double_t) i2,
04004 working_space[7 * j + 1],
04005 working_space[7 * j + 2],
04006 working_space[peak_vel],
04007 working_space[peak_vel + 1],
04008 working_space[peak_vel + 2],
04009 working_space[peak_vel + 6],
04010 working_space[peak_vel + 7],
04011 working_space[peak_vel + 12],
04012 working_space[peak_vel + 13]);
04013 if (yw != 0) {
04014 c = Ourpowl(a, pw);
04015 working_space[2 * shift + k] += chi_opt * c;
04016 b = a * a / yw;
04017 working_space[4 * shift + k] += b * c;
04018 }
04019 k += 1;
04020 }
04021 if (fFixPositionX[j] == false) {
04022 a = Deri02((Double_t) i1, (Double_t) i2,
04023 working_space[7 * j],
04024 working_space[7 * j + 1],
04025 working_space[7 * j + 2],
04026 working_space[peak_vel],
04027 working_space[peak_vel + 1],
04028 working_space[peak_vel + 2],
04029 working_space[peak_vel + 6],
04030 working_space[peak_vel + 7],
04031 working_space[peak_vel + 12],
04032 working_space[peak_vel + 13]);
04033 if (yw != 0) {
04034 c = Ourpowl(a, pw);
04035 working_space[2 * shift + k] += chi_opt * c;
04036 b = a * a / yw;
04037 working_space[4 * shift + k] += b * c;
04038 }
04039 k += 1;
04040 }
04041 if (fFixPositionY[j] == false) {
04042 a = Derj02((Double_t) i1, (Double_t) i2,
04043 working_space[7 * j],
04044 working_space[7 * j + 1],
04045 working_space[7 * j + 2],
04046 working_space[peak_vel],
04047 working_space[peak_vel + 1],
04048 working_space[peak_vel + 2],
04049 working_space[peak_vel + 6],
04050 working_space[peak_vel + 7],
04051 working_space[peak_vel + 12],
04052 working_space[peak_vel + 13]);
04053 if (yw != 0) {
04054 c = Ourpowl(a, pw);
04055 working_space[2 * shift + k] += chi_opt * c;
04056 b = a * a / yw;
04057 working_space[4 * shift + k] += b * c;
04058 }
04059 k += 1;
04060 }
04061 if (fFixAmpX1[j] == false) {
04062 a = Derampx((Double_t) i1, working_space[7 * j + 5],
04063 working_space[peak_vel],
04064 working_space[peak_vel + 8],
04065 working_space[peak_vel + 10],
04066 working_space[peak_vel + 12]);
04067 if (yw != 0) {
04068 c = Ourpowl(a, pw);
04069 working_space[2 * shift + k] += chi_opt * c;
04070 b = a * a / yw;
04071 working_space[4 * shift + k] += b * c;
04072 }
04073 k += 1;
04074 }
04075 if (fFixAmpY1[j] == false) {
04076 a = Derampx((Double_t) i2, working_space[7 * j + 6],
04077 working_space[peak_vel + 1],
04078 working_space[peak_vel + 9],
04079 working_space[peak_vel + 11],
04080 working_space[peak_vel + 13]);
04081 if (yw != 0) {
04082 c = Ourpowl(a, pw);
04083 working_space[2 * shift + k] += chi_opt * c;
04084 b = a * a / yw;
04085 working_space[4 * shift + k] += b * c;
04086 }
04087 k += 1;
04088 }
04089 if (fFixPositionX1[j] == false) {
04090 a = Deri01((Double_t) i1, working_space[7 * j + 3],
04091 working_space[7 * j + 5],
04092 working_space[peak_vel],
04093 working_space[peak_vel + 8],
04094 working_space[peak_vel + 10],
04095 working_space[peak_vel + 12]);
04096 if (yw != 0) {
04097 c = Ourpowl(a, pw);
04098 working_space[2 * shift + k] += chi_opt * c;
04099 b = a * a / yw;
04100 working_space[4 * shift + k] += b * c;
04101 }
04102 k += 1;
04103 }
04104 if (fFixPositionY1[j] == false) {
04105 a = Deri01((Double_t) i2, working_space[7 * j + 4],
04106 working_space[7 * j + 6],
04107 working_space[peak_vel + 1],
04108 working_space[peak_vel + 9],
04109 working_space[peak_vel + 11],
04110 working_space[peak_vel + 13]);
04111 if (yw != 0) {
04112 c = Ourpowl(a, pw);
04113 working_space[2 * shift + k] += chi_opt * c;
04114 b = a * a / yw;
04115 working_space[4 * shift + k] += b * c;
04116 }
04117 k += 1;
04118 }
04119 }
04120 if (fFixSigmaX == false) {
04121 a = Dersigmax(fNPeaks, (Double_t) i1, (Double_t) i2,
04122 working_space, working_space[peak_vel],
04123 working_space[peak_vel + 1],
04124 working_space[peak_vel + 2],
04125 working_space[peak_vel + 6],
04126 working_space[peak_vel + 7],
04127 working_space[peak_vel + 8],
04128 working_space[peak_vel + 10],
04129 working_space[peak_vel + 12],
04130 working_space[peak_vel + 13]);
04131 if (yw != 0) {
04132 c = Ourpowl(a, pw);
04133 working_space[2 * shift + k] += chi_opt * c;
04134 b = a * a / yw;
04135 working_space[4 * shift + k] += b * c;
04136 }
04137 k += 1;
04138 }
04139 if (fFixSigmaY == false) {
04140 a = Dersigmay(fNPeaks, (Double_t) i1, (Double_t) i2,
04141 working_space, working_space[peak_vel],
04142 working_space[peak_vel + 1],
04143 working_space[peak_vel + 2],
04144 working_space[peak_vel + 6],
04145 working_space[peak_vel + 7],
04146 working_space[peak_vel + 9],
04147 working_space[peak_vel + 11],
04148 working_space[peak_vel + 12],
04149 working_space[peak_vel + 13]);
04150 if (yw != 0) {
04151 c = Ourpowl(a, pw);
04152 working_space[2 * shift + k] += chi_opt * c;
04153 b = a * a / yw;
04154 working_space[4 * shift + k] += b * c;
04155 }
04156 k += 1;
04157 }
04158 if (fFixRo == false) {
04159 a = Derro(fNPeaks, (Double_t) i1, (Double_t) i2,
04160 working_space, working_space[peak_vel],
04161 working_space[peak_vel + 1],
04162 working_space[peak_vel + 2]);
04163 if (yw != 0) {
04164 c = Ourpowl(a, pw);
04165 working_space[2 * shift + k] += chi_opt * c;
04166 b = a * a / yw;
04167 working_space[4 * shift + k] += b * c;
04168 }
04169 k += 1;
04170 }
04171 if (fFixA0 == false) {
04172 a = 1.;
04173 if (yw != 0) {
04174 c = Ourpowl(a, pw);
04175 working_space[2 * shift + k] += chi_opt * c;
04176 b = a * a / yw;
04177 working_space[4 * shift + k] += b * c;
04178 }
04179 k += 1;
04180 }
04181 if (fFixAx == false) {
04182 a = i1;
04183 if (yw != 0) {
04184 c = Ourpowl(a, pw);
04185 working_space[2 * shift + k] += chi_opt * c;
04186 b = a * a / yw;
04187 working_space[4 * shift + k] += b * c;
04188 }
04189 k += 1;
04190 }
04191 if (fFixAy == false) {
04192 a = i2;
04193 if (yw != 0) {
04194 c = Ourpowl(a, pw);
04195 working_space[2 * shift + k] += chi_opt * c;
04196 b = a * a / yw;
04197 working_space[4 * shift + k] += b * c;
04198 }
04199 k += 1;
04200 }
04201 if (fFixTxy == false) {
04202 a = Dertxy(fNPeaks, (Double_t) i1, (Double_t) i2,
04203 working_space, working_space[peak_vel],
04204 working_space[peak_vel + 1],
04205 working_space[peak_vel + 12],
04206 working_space[peak_vel + 13]);
04207 if (yw != 0) {
04208 c = Ourpowl(a, pw);
04209 working_space[2 * shift + k] += chi_opt * c;
04210 b = a * a / yw;
04211 working_space[4 * shift + k] += b * c;
04212 }
04213 k += 1;
04214 }
04215 if (fFixSxy == false) {
04216 a = Dersxy(fNPeaks, (Double_t) i1, (Double_t) i2,
04217 working_space, working_space[peak_vel],
04218 working_space[peak_vel + 1]);
04219 if (yw != 0) {
04220 c = Ourpowl(a, pw);
04221 working_space[2 * shift + k] += chi_opt * c;
04222 b = a * a / yw;
04223 working_space[4 * shift + k] += b * c;
04224 }
04225 k += 1;
04226 }
04227 if (fFixTx == false) {
04228 a = Dertx(fNPeaks, (Double_t) i1, working_space,
04229 working_space[peak_vel],
04230 working_space[peak_vel + 12]);
04231 if (yw != 0) {
04232 c = Ourpowl(a, pw);
04233 working_space[2 * shift + k] += chi_opt * c;
04234 b = a * a / yw;
04235 working_space[4 * shift + k] += b * c;
04236 }
04237 k += 1;
04238 }
04239 if (fFixTy == false) {
04240 a = Derty(fNPeaks, (Double_t) i2, working_space,
04241 working_space[peak_vel + 1],
04242 working_space[peak_vel + 13]);
04243 if (yw != 0) {
04244 c = Ourpowl(a, pw);
04245 working_space[2 * shift + k] += chi_opt * c;
04246 b = a * a / yw;
04247 working_space[4 * shift + k] += b * c;
04248 }
04249 k += 1;
04250 }
04251 if (fFixSx == false) {
04252 a = Dersx(fNPeaks, (Double_t) i1, working_space,
04253 working_space[peak_vel]);
04254 if (yw != 0) {
04255 c = Ourpowl(a, pw);
04256 working_space[2 * shift + k] += chi_opt * c;
04257 b = a * a / yw;
04258 working_space[4 * shift + k] += b * c;
04259 }
04260 k += 1;
04261 }
04262 if (fFixSy == false) {
04263 a = Dersy(fNPeaks, (Double_t) i2, working_space,
04264 working_space[peak_vel + 1]);
04265 if (yw != 0) {
04266 c = Ourpowl(a, pw);
04267 working_space[2 * shift + k] += chi_opt * c;
04268 b = a * a / yw;
04269 working_space[4 * shift + k] += b * c;
04270 }
04271 k += 1;
04272 }
04273 if (fFixBx == false) {
04274 a = Derbx(fNPeaks, (Double_t) i1, (Double_t) i2,
04275 working_space, working_space[peak_vel],
04276 working_space[peak_vel + 1],
04277 working_space[peak_vel + 6],
04278 working_space[peak_vel + 8],
04279 working_space[peak_vel + 12],
04280 working_space[peak_vel + 13]);
04281 if (yw != 0) {
04282 c = Ourpowl(a, pw);
04283 working_space[2 * shift + k] += chi_opt * c;
04284 b = a * a / yw;
04285 working_space[4 * shift + k] += b * c;
04286 }
04287 k += 1;
04288 }
04289 if (fFixBy == false) {
04290 a = Derby(fNPeaks, (Double_t) i1, (Double_t) i2,
04291 working_space, working_space[peak_vel],
04292 working_space[peak_vel + 1],
04293 working_space[peak_vel + 6],
04294 working_space[peak_vel + 8],
04295 working_space[peak_vel + 12],
04296 working_space[peak_vel + 13]);
04297 if (yw != 0) {
04298 c = Ourpowl(a, pw);
04299 working_space[2 * shift + k] += chi_opt * c;
04300 b = a * a / yw;
04301 working_space[4 * shift + k] += b * c;
04302 }
04303 k += 1;
04304 }
04305 }
04306 }
04307 }
04308 b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
04309 chi_er = chi_cel / b;
04310 for (i = 0, j = 0; i < fNPeaks; i++) {
04311 fVolume[i] =
04312 Volume(working_space[7 * i], working_space[peak_vel],
04313 working_space[peak_vel + 1], working_space[peak_vel + 2]);
04314 if (fVolume[i] > 0) {
04315 c = 0;
04316 if (fFixAmp[i] == false) {
04317 a = Derpa2(working_space[peak_vel],
04318 working_space[peak_vel + 1],
04319 working_space[peak_vel + 2]);
04320 b = working_space[4 * shift + j];
04321 if (b == 0)
04322 b = 1;
04323
04324 else
04325 b = 1 / b;
04326 c = c + a * a * b;
04327 }
04328 if (fFixSigmaX == false) {
04329 a = Derpsigmax(working_space[shift + j],
04330 working_space[peak_vel + 1],
04331 working_space[peak_vel + 2]);
04332 b = working_space[4 * shift + peak_vel];
04333 if (b == 0)
04334 b = 1;
04335
04336 else
04337 b = 1 / b;
04338 c = c + a * a * b;
04339 }
04340 if (fFixSigmaY == false) {
04341 a = Derpsigmay(working_space[shift + j],
04342 working_space[peak_vel],
04343 working_space[peak_vel + 2]);
04344 b = working_space[4 * shift + peak_vel + 1];
04345 if (b == 0)
04346 b = 1;
04347
04348 else
04349 b = 1 / b;
04350 c = c + a * a * b;
04351 }
04352 if (fFixRo == false) {
04353 a = Derpro(working_space[shift + j], working_space[peak_vel],
04354 working_space[peak_vel + 1],
04355 working_space[peak_vel + 2]);
04356 b = working_space[4 * shift + peak_vel + 2];
04357 if (b == 0)
04358 b = 1;
04359
04360 else
04361 b = 1 / b;
04362 c = c + a * a * b;
04363 }
04364 fVolumeErr[i] = TMath::Sqrt(TMath::Abs(chi_er * c));
04365 }
04366
04367 else {
04368 fVolumeErr[i] = 0;
04369 }
04370 if (fFixAmp[i] == false) {
04371 fAmpCalc[i] = working_space[shift + j];
04372 if (working_space[3 * shift + j] != 0)
04373 fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04374 j += 1;
04375 }
04376
04377 else {
04378 fAmpCalc[i] = fAmpInit[i];
04379 fAmpErr[i] = 0;
04380 }
04381 if (fFixPositionX[i] == false) {
04382 fPositionCalcX[i] = working_space[shift + j];
04383 if (working_space[3 * shift + j] != 0)
04384 fPositionErrX[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04385 j += 1;
04386 }
04387
04388 else {
04389 fPositionCalcX[i] = fPositionInitX[i];
04390 fPositionErrX[i] = 0;
04391 }
04392 if (fFixPositionY[i] == false) {
04393 fPositionCalcY[i] = working_space[shift + j];
04394 if (working_space[3 * shift + j] != 0)
04395 fPositionErrY[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04396 j += 1;
04397 }
04398
04399 else {
04400 fPositionCalcY[i] = fPositionInitY[i];
04401 fPositionErrY[i] = 0;
04402 }
04403 if (fFixAmpX1[i] == false) {
04404 fAmpCalcX1[i] = working_space[shift + j];
04405 if (working_space[3 * shift + j] != 0)
04406 fAmpErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04407 j += 1;
04408 }
04409
04410 else {
04411 fAmpCalcX1[i] = fAmpInitX1[i];
04412 fAmpErrX1[i] = 0;
04413 }
04414 if (fFixAmpY1[i] == false) {
04415 fAmpCalcY1[i] = working_space[shift + j];
04416 if (working_space[3 * shift + j] != 0)
04417 fAmpErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04418 j += 1;
04419 }
04420
04421 else {
04422 fAmpCalcY1[i] = fAmpInitY1[i];
04423 fAmpErrY1[i] = 0;
04424 }
04425 if (fFixPositionX1[i] == false) {
04426 fPositionCalcX1[i] = working_space[shift + j];
04427 if (working_space[3 * shift + j] != 0)
04428 fPositionErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04429 j += 1;
04430 }
04431
04432 else {
04433 fPositionCalcX1[i] = fPositionInitX1[i];
04434 fPositionErrX1[i] = 0;
04435 }
04436 if (fFixPositionY1[i] == false) {
04437 fPositionCalcY1[i] = working_space[shift + j];
04438 if (working_space[3 * shift + j] != 0)
04439 fPositionErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04440 j += 1;
04441 }
04442
04443 else {
04444 fPositionCalcY1[i] = fPositionInitY1[i];
04445 fPositionErrY1[i] = 0;
04446 }
04447 }
04448 if (fFixSigmaX == false) {
04449 fSigmaCalcX = working_space[shift + j];
04450 if (working_space[3 * shift + j] != 0)
04451 fSigmaErrX = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04452 j += 1;
04453 }
04454
04455 else {
04456 fSigmaCalcX = fSigmaInitX;
04457 fSigmaErrX = 0;
04458 }
04459 if (fFixSigmaY == false) {
04460 fSigmaCalcY = working_space[shift + j];
04461 if (working_space[3 * shift + j] != 0)
04462 fSigmaErrY = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04463 j += 1;
04464 }
04465
04466 else {
04467 fSigmaCalcY = fSigmaInitY;
04468 fSigmaErrY = 0;
04469 }
04470 if (fFixRo == false) {
04471 fRoCalc = working_space[shift + j];
04472 if (working_space[3 * shift + j] != 0)
04473 fRoErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04474 j += 1;
04475 }
04476
04477 else {
04478 fRoCalc = fRoInit;
04479 fRoErr = 0;
04480 }
04481 if (fFixA0 == false) {
04482 fA0Calc = working_space[shift + j];
04483 if (working_space[3 * shift + j] != 0)
04484 fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04485 j += 1;
04486 }
04487
04488 else {
04489 fA0Calc = fA0Init;
04490 fA0Err = 0;
04491 }
04492 if (fFixAx == false) {
04493 fAxCalc = working_space[shift + j];
04494 if (working_space[3 * shift + j] != 0)
04495 fAxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04496 j += 1;
04497 }
04498
04499 else {
04500 fAxCalc = fAxInit;
04501 fAxErr = 0;
04502 }
04503 if (fFixAy == false) {
04504 fAyCalc = working_space[shift + j];
04505 if (working_space[3 * shift + j] != 0)
04506 fAyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04507 j += 1;
04508 }
04509
04510 else {
04511 fAyCalc = fAyInit;
04512 fAyErr = 0;
04513 }
04514 if (fFixTxy == false) {
04515 fTxyCalc = working_space[shift + j];
04516 if (working_space[3 * shift + j] != 0)
04517 fTxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04518 j += 1;
04519 }
04520
04521 else {
04522 fTxyCalc = fTxyInit;
04523 fTxyErr = 0;
04524 }
04525 if (fFixSxy == false) {
04526 fSxyCalc = working_space[shift + j];
04527 if (working_space[3 * shift + j] != 0)
04528 fSxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04529 j += 1;
04530 }
04531
04532 else {
04533 fSxyCalc = fSxyInit;
04534 fSxyErr = 0;
04535 }
04536 if (fFixTx == false) {
04537 fTxCalc = working_space[shift + j];
04538 if (working_space[3 * shift + j] != 0)
04539 fTxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04540 j += 1;
04541 }
04542
04543 else {
04544 fTxCalc = fTxInit;
04545 fTxErr = 0;
04546 }
04547 if (fFixTy == false) {
04548 fTyCalc = working_space[shift + j];
04549 if (working_space[3 * shift + j] != 0)
04550 fTyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04551 j += 1;
04552 }
04553
04554 else {
04555 fTyCalc = fTyInit;
04556 fTyErr = 0;
04557 }
04558 if (fFixSx == false) {
04559 fSxCalc = working_space[shift + j];
04560 if (working_space[3 * shift + j] != 0)
04561 fSxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04562 j += 1;
04563 }
04564
04565 else {
04566 fSxCalc = fSxInit;
04567 fSxErr = 0;
04568 }
04569 if (fFixSy == false) {
04570 fSyCalc = working_space[shift + j];
04571 if (working_space[3 * shift + j] != 0)
04572 fSyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04573 j += 1;
04574 }
04575
04576 else {
04577 fSyCalc = fSyInit;
04578 fSyErr = 0;
04579 }
04580 if (fFixBx == false) {
04581 fBxCalc = working_space[shift + j];
04582 if (working_space[3 * shift + j] != 0)
04583 fBxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04584 j += 1;
04585 }
04586
04587 else {
04588 fBxCalc = fBxInit;
04589 fBxErr = 0;
04590 }
04591 if (fFixBy == false) {
04592 fByCalc = working_space[shift + j];
04593 if (working_space[3 * shift + j] != 0)
04594 fByErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
04595 j += 1;
04596 }
04597
04598 else {
04599 fByCalc = fByInit;
04600 fByErr = 0;
04601 }
04602 b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
04603 fChi = chi_cel / b;
04604 for (i1 = fXmin; i1 <= fXmax; i1++) {
04605 for (i2 = fYmin; i2 <= fYmax; i2++) {
04606 f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
04607 working_space, working_space[peak_vel],
04608 working_space[peak_vel + 1],
04609 working_space[peak_vel + 2],
04610 working_space[peak_vel + 3],
04611 working_space[peak_vel + 4],
04612 working_space[peak_vel + 5],
04613 working_space[peak_vel + 6],
04614 working_space[peak_vel + 7],
04615 working_space[peak_vel + 8],
04616 working_space[peak_vel + 9],
04617 working_space[peak_vel + 10],
04618 working_space[peak_vel + 11],
04619 working_space[peak_vel + 12],
04620 working_space[peak_vel + 13]);
04621 source[i1][i2] = f;
04622 }
04623 }
04624 delete [] working_space;
04625 return;
04626 }
04627
04628
04629
04630
04631
04632 void TSpectrum2Fit::FitStiefel(Float_t **source)
04633 {
04634
04635
04636
04637
04638
04639
04640
04641
04642
04643
04644
04645
04646
04647
04648
04649
04650
04651
04652
04653
04654
04655
04656
04657
04658
04659
04660
04661
04662
04663
04664
04665
04666
04667
04668
04669
04670
04671
04672
04673
04674
04675
04676
04677
04678
04679
04680
04681
04682
04683
04684
04685
04686
04687
04688
04689
04690
04691
04692
04693
04694
04695
04696
04697
04698
04699
04700
04701
04702
04703
04704
04705
04706
04707
04708
04709
04710
04711
04712
04713
04714
04715
04716
04717
04718
04719
04720
04721
04722
04723
04724
04725
04726
04727
04728
04729
04730
04731
04732
04733
04734
04735
04736
04737
04738
04739
04740
04741
04742
04743
04744
04745
04746
04747
04748
04749
04750
04751
04752
04753
04754
04755
04756
04757
04758
04759
04760
04761
04762
04763
04764
04765
04766
04767
04768
04769
04770
04771
04772
04773
04774
04775
04776
04777
04778
04779
04780
04781
04782
04783
04784
04785
04786
04787
04788
04789
04790
04791
04792
04793
04794
04795
04796
04797
04798
04799
04800
04801
04802
04803
04804
04805
04806
04807
04808
04809
04810
04811
04812
04813
04814
04815
04816
04817
04818
04819
04820
04821
04822
04823
04824
04825
04826
04827
04828
04829
04830
04831
04832
04833
04834
04835
04836
04837
04838
04839
04840
04841
04842
04843
04844
04845
04846
04847
04848
04849
04850
04851
04852
04853
04854
04855
04856
04857
04858
04859
04860
04861
04862
04863
04864
04865
04866
04867
04868
04869
04870
04871
04872
04873
04874
04875
04876
04877
04878
04879
04880
04881
04882 Int_t i, i1, i2, j, k, shift =
04883 7 * fNPeaks + 14, peak_vel, size, iter, regul_cycle,
04884 flag;
04885 Double_t a, b, c, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi = 0
04886 , pi, pmin = 0, chi_cel = 0, chi_er;
04887 Double_t *working_space = new Double_t[5 * (7 * fNPeaks + 14)];
04888 for (i = 0, j = 0; i < fNPeaks; i++) {
04889 working_space[7 * i] = fAmpInit[i];
04890 if (fFixAmp[i] == false) {
04891 working_space[shift + j] = fAmpInit[i];
04892 j += 1;
04893 }
04894 working_space[7 * i + 1] = fPositionInitX[i];
04895 if (fFixPositionX[i] == false) {
04896 working_space[shift + j] = fPositionInitX[i];
04897 j += 1;
04898 }
04899 working_space[7 * i + 2] = fPositionInitY[i];
04900 if (fFixPositionY[i] == false) {
04901 working_space[shift + j] = fPositionInitY[i];
04902 j += 1;
04903 }
04904 working_space[7 * i + 3] = fAmpInitX1[i];
04905 if (fFixAmpX1[i] == false) {
04906 working_space[shift + j] = fAmpInitX1[i];
04907 j += 1;
04908 }
04909 working_space[7 * i + 4] = fAmpInitY1[i];
04910 if (fFixAmpY1[i] == false) {
04911 working_space[shift + j] = fAmpInitY1[i];
04912 j += 1;
04913 }
04914 working_space[7 * i + 5] = fPositionInitX1[i];
04915 if (fFixPositionX1[i] == false) {
04916 working_space[shift + j] = fPositionInitX1[i];
04917 j += 1;
04918 }
04919 working_space[7 * i + 6] = fPositionInitY1[i];
04920 if (fFixPositionY1[i] == false) {
04921 working_space[shift + j] = fPositionInitY1[i];
04922 j += 1;
04923 }
04924 }
04925 peak_vel = 7 * i;
04926 working_space[7 * i] = fSigmaInitX;
04927 if (fFixSigmaX == false) {
04928 working_space[shift + j] = fSigmaInitX;
04929 j += 1;
04930 }
04931 working_space[7 * i + 1] = fSigmaInitY;
04932 if (fFixSigmaY == false) {
04933 working_space[shift + j] = fSigmaInitY;
04934 j += 1;
04935 }
04936 working_space[7 * i + 2] = fRoInit;
04937 if (fFixRo == false) {
04938 working_space[shift + j] = fRoInit;
04939 j += 1;
04940 }
04941 working_space[7 * i + 3] = fA0Init;
04942 if (fFixA0 == false) {
04943 working_space[shift + j] = fA0Init;
04944 j += 1;
04945 }
04946 working_space[7 * i + 4] = fAxInit;
04947 if (fFixAx == false) {
04948 working_space[shift + j] = fAxInit;
04949 j += 1;
04950 }
04951 working_space[7 * i + 5] = fAyInit;
04952 if (fFixAy == false) {
04953 working_space[shift + j] = fAyInit;
04954 j += 1;
04955 }
04956 working_space[7 * i + 6] = fTxyInit;
04957 if (fFixTxy == false) {
04958 working_space[shift + j] = fTxyInit;
04959 j += 1;
04960 }
04961 working_space[7 * i + 7] = fSxyInit;
04962 if (fFixSxy == false) {
04963 working_space[shift + j] = fSxyInit;
04964 j += 1;
04965 }
04966 working_space[7 * i + 8] = fTxInit;
04967 if (fFixTx == false) {
04968 working_space[shift + j] = fTxInit;
04969 j += 1;
04970 }
04971 working_space[7 * i + 9] = fTyInit;
04972 if (fFixTy == false) {
04973 working_space[shift + j] = fTyInit;
04974 j += 1;
04975 }
04976 working_space[7 * i + 10] = fSxyInit;
04977 if (fFixSx == false) {
04978 working_space[shift + j] = fSxInit;
04979 j += 1;
04980 }
04981 working_space[7 * i + 11] = fSyInit;
04982 if (fFixSy == false) {
04983 working_space[shift + j] = fSyInit;
04984 j += 1;
04985 }
04986 working_space[7 * i + 12] = fBxInit;
04987 if (fFixBx == false) {
04988 working_space[shift + j] = fBxInit;
04989 j += 1;
04990 }
04991 working_space[7 * i + 13] = fByInit;
04992 if (fFixBy == false) {
04993 working_space[shift + j] = fByInit;
04994 j += 1;
04995 }
04996 size = j;
04997 Double_t **working_matrix = new Double_t *[size];
04998 for (i = 0; i < size; i++)
04999 working_matrix[i] = new Double_t[size + 4];
05000 for (iter = 0; iter < fNumberIterations; iter++) {
05001 for (j = 0; j < size; j++) {
05002 working_space[3 * shift + j] = 0;
05003 for (k = 0; k <= size; k++) {
05004 working_matrix[j][k] = 0;
05005 }
05006 }
05007
05008
05009 alpha = fAlpha;
05010 chi_opt = 0;
05011 for (i1 = fXmin; i1 <= fXmax; i1++) {
05012 for (i2 = fYmin; i2 <= fYmax; i2++) {
05013
05014 for (j = 0, k = 0; j < fNPeaks; j++) {
05015 if (fFixAmp[j] == false) {
05016 working_space[2 * shift + k] =
05017 Deramp2((Double_t) i1, (Double_t) i2,
05018 working_space[7 * j + 1],
05019 working_space[7 * j + 2],
05020 working_space[peak_vel],
05021 working_space[peak_vel + 1],
05022 working_space[peak_vel + 2],
05023 working_space[peak_vel + 6],
05024 working_space[peak_vel + 7],
05025 working_space[peak_vel + 12],
05026 working_space[peak_vel + 13]);
05027 k += 1;
05028 }
05029 if (fFixPositionX[j] == false) {
05030 working_space[2 * shift + k] =
05031 Deri02((Double_t) i1, (Double_t) i2,
05032 working_space[7 * j],
05033 working_space[7 * j + 1],
05034 working_space[7 * j + 2],
05035 working_space[peak_vel],
05036 working_space[peak_vel + 1],
05037 working_space[peak_vel + 2],
05038 working_space[peak_vel + 6],
05039 working_space[peak_vel + 7],
05040 working_space[peak_vel + 12],
05041 working_space[peak_vel + 13]);
05042 k += 1;
05043 }
05044 if (fFixPositionY[j] == false) {
05045 working_space[2 * shift + k] =
05046 Derj02((Double_t) i1, (Double_t) i2,
05047 working_space[7 * j],
05048 working_space[7 * j + 1],
05049 working_space[7 * j + 2],
05050 working_space[peak_vel],
05051 working_space[peak_vel + 1],
05052 working_space[peak_vel + 2],
05053 working_space[peak_vel + 6],
05054 working_space[peak_vel + 7],
05055 working_space[peak_vel + 12],
05056 working_space[peak_vel + 13]);
05057 k += 1;
05058 }
05059 if (fFixAmpX1[j] == false) {
05060 working_space[2 * shift + k] =
05061 Derampx((Double_t) i1, working_space[7 * j + 5],
05062 working_space[peak_vel],
05063 working_space[peak_vel + 8],
05064 working_space[peak_vel + 10],
05065 working_space[peak_vel + 12]);
05066 k += 1;
05067 }
05068 if (fFixAmpY1[j] == false) {
05069 working_space[2 * shift + k] =
05070 Derampx((Double_t) i2, working_space[7 * j + 6],
05071 working_space[peak_vel + 1],
05072 working_space[peak_vel + 9],
05073 working_space[peak_vel + 11],
05074 working_space[peak_vel + 13]);
05075 k += 1;
05076 }
05077 if (fFixPositionX1[j] == false) {
05078 working_space[2 * shift + k] =
05079 Deri01((Double_t) i1, working_space[7 * j + 3],
05080 working_space[7 * j + 5],
05081 working_space[peak_vel],
05082 working_space[peak_vel + 8],
05083 working_space[peak_vel + 10],
05084 working_space[peak_vel + 12]);
05085 k += 1;
05086 }
05087 if (fFixPositionY1[j] == false) {
05088 working_space[2 * shift + k] =
05089 Deri01((Double_t) i2, working_space[7 * j + 4],
05090 working_space[7 * j + 6],
05091 working_space[peak_vel + 1],
05092 working_space[peak_vel + 9],
05093 working_space[peak_vel + 11],
05094 working_space[peak_vel + 13]);
05095 k += 1;
05096 }
05097 } if (fFixSigmaX == false) {
05098 working_space[2 * shift + k] =
05099 Dersigmax(fNPeaks, (Double_t) i1, (Double_t) i2,
05100 working_space, working_space[peak_vel],
05101 working_space[peak_vel + 1],
05102 working_space[peak_vel + 2],
05103 working_space[peak_vel + 6],
05104 working_space[peak_vel + 7],
05105 working_space[peak_vel + 8],
05106 working_space[peak_vel + 10],
05107 working_space[peak_vel + 12],
05108 working_space[peak_vel + 13]);
05109 k += 1;
05110 }
05111 if (fFixSigmaY == false) {
05112 working_space[2 * shift + k] =
05113 Dersigmay(fNPeaks, (Double_t) i1, (Double_t) i2,
05114 working_space, working_space[peak_vel],
05115 working_space[peak_vel + 1],
05116 working_space[peak_vel + 2],
05117 working_space[peak_vel + 6],
05118 working_space[peak_vel + 7],
05119 working_space[peak_vel + 9],
05120 working_space[peak_vel + 11],
05121 working_space[peak_vel + 12],
05122 working_space[peak_vel + 13]);
05123 k += 1;
05124 }
05125 if (fFixRo == false) {
05126 working_space[2 * shift + k] =
05127 Derro(fNPeaks, (Double_t) i1, (Double_t) i2,
05128 working_space, working_space[peak_vel],
05129 working_space[peak_vel + 1],
05130 working_space[peak_vel + 2]);
05131 k += 1;
05132 }
05133 if (fFixA0 == false) {
05134 working_space[2 * shift + k] = 1.;
05135 k += 1;
05136 }
05137 if (fFixAx == false) {
05138 working_space[2 * shift + k] = i1;
05139 k += 1;
05140 }
05141 if (fFixAy == false) {
05142 working_space[2 * shift + k] = i2;
05143 k += 1;
05144 }
05145 if (fFixTxy == false) {
05146 working_space[2 * shift + k] =
05147 Dertxy(fNPeaks, (Double_t) i1, (Double_t) i2,
05148 working_space, working_space[peak_vel],
05149 working_space[peak_vel + 1],
05150 working_space[peak_vel + 12],
05151 working_space[peak_vel + 13]);
05152 k += 1;
05153 }
05154 if (fFixSxy == false) {
05155 working_space[2 * shift + k] =
05156 Dersxy(fNPeaks, (Double_t) i1, (Double_t) i2,
05157 working_space, working_space[peak_vel],
05158 working_space[peak_vel + 1]);
05159 k += 1;
05160 }
05161 if (fFixTx == false) {
05162 working_space[2 * shift + k] =
05163 Dertx(fNPeaks, (Double_t) i1, working_space,
05164 working_space[peak_vel],
05165 working_space[peak_vel + 12]);
05166 k += 1;
05167 }
05168 if (fFixTy == false) {
05169 working_space[2 * shift + k] =
05170 Derty(fNPeaks, (Double_t) i2, working_space,
05171 working_space[peak_vel + 1],
05172 working_space[peak_vel + 13]);
05173 k += 1;
05174 }
05175 if (fFixSx == false) {
05176 working_space[2 * shift + k] =
05177 Dersx(fNPeaks, (Double_t) i1, working_space,
05178 working_space[peak_vel]);
05179 k += 1;
05180 }
05181 if (fFixSy == false) {
05182 working_space[2 * shift + k] =
05183 Dersy(fNPeaks, (Double_t) i2, working_space,
05184 working_space[peak_vel + 1]);
05185 k += 1;
05186 }
05187 if (fFixBx == false) {
05188 working_space[2 * shift + k] =
05189 Derbx(fNPeaks, (Double_t) i1, (Double_t) i2,
05190 working_space, working_space[peak_vel],
05191 working_space[peak_vel + 1],
05192 working_space[peak_vel + 6],
05193 working_space[peak_vel + 8],
05194 working_space[peak_vel + 12],
05195 working_space[peak_vel + 13]);
05196 k += 1;
05197 }
05198 if (fFixBy == false) {
05199 working_space[2 * shift + k] =
05200 Derby(fNPeaks, (Double_t) i1, (Double_t) i2,
05201 working_space, working_space[peak_vel],
05202 working_space[peak_vel + 1],
05203 working_space[peak_vel + 6],
05204 working_space[peak_vel + 8],
05205 working_space[peak_vel + 12],
05206 working_space[peak_vel + 13]);
05207 k += 1;
05208 }
05209 yw = source[i1][i2];
05210 ywm = yw;
05211 f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
05212 working_space, working_space[peak_vel],
05213 working_space[peak_vel + 1],
05214 working_space[peak_vel + 2],
05215 working_space[peak_vel + 3],
05216 working_space[peak_vel + 4],
05217 working_space[peak_vel + 5],
05218 working_space[peak_vel + 6],
05219 working_space[peak_vel + 7],
05220 working_space[peak_vel + 8],
05221 working_space[peak_vel + 9],
05222 working_space[peak_vel + 10],
05223 working_space[peak_vel + 11],
05224 working_space[peak_vel + 12],
05225 working_space[peak_vel + 13]);
05226 if (fStatisticType == kFitOptimMaxLikelihood) {
05227 if (f > 0.00001)
05228 chi_opt += yw * TMath::Log(f) - f;
05229 }
05230
05231 else {
05232 if (ywm != 0)
05233 chi_opt += (yw - f) * (yw - f) / ywm;
05234 }
05235 if (fStatisticType == kFitOptimChiFuncValues) {
05236 ywm = f;
05237 if (f < 0.00001)
05238 ywm = 0.00001;
05239 }
05240
05241 else if (fStatisticType == kFitOptimMaxLikelihood) {
05242 ywm = f;
05243 if (f < 0.00001)
05244 ywm = 0.00001;
05245 }
05246
05247 else {
05248 if (ywm == 0)
05249 ywm = 1;
05250 }
05251 for (j = 0; j < size; j++) {
05252 for (k = 0; k < size; k++) {
05253 b = working_space[2 * shift +
05254 j] * working_space[2 * shift +
05255 k] / ywm;
05256 if (fStatisticType == kFitOptimChiFuncValues)
05257 b = b * (4 * yw - 2 * f) / ywm;
05258 working_matrix[j][k] += b;
05259 if (j == k)
05260 working_space[3 * shift + j] += b;
05261 }
05262 }
05263 if (fStatisticType == kFitOptimChiFuncValues)
05264 b = (f * f - yw * yw) / (ywm * ywm);
05265
05266 else
05267 b = (f - yw) / ywm;
05268 for (j = 0; j < size; j++) {
05269 working_matrix[j][size] -=
05270 b * working_space[2 * shift + j];
05271 }
05272 }
05273 }
05274 for (i = 0; i < size; i++) {
05275 working_matrix[i][size + 1] = 0;
05276 }
05277 StiefelInversion(working_matrix, size);
05278 for (i = 0; i < size; i++) {
05279 working_space[2 * shift + i] = working_matrix[i][size + 1];
05280 }
05281
05282
05283 chi2 = chi_opt;
05284 chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
05285
05286
05287 regul_cycle = 0;
05288 for (j = 0; j < size; j++) {
05289 working_space[4 * shift + j] = working_space[shift + j];
05290 }
05291
05292 do {
05293 if (fAlphaOptim == kFitAlphaOptimal) {
05294 if (fStatisticType != kFitOptimMaxLikelihood)
05295 chi_min = 10000 * chi2;
05296
05297 else
05298 chi_min = 0.1 * chi2;
05299 flag = 0;
05300 for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
05301 for (j = 0; j < size; j++) {
05302 working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j];
05303 }
05304 for (i = 0, j = 0; i < fNPeaks; i++) {
05305 if (fFixAmp[i] == false) {
05306 if (working_space[shift + j] < 0)
05307 working_space[shift + j] = 0;
05308 working_space[7 * i] = working_space[shift + j];
05309 j += 1;
05310 }
05311 if (fFixPositionX[i] == false) {
05312 if (working_space[shift + j] < fXmin)
05313 working_space[shift + j] = fXmin;
05314 if (working_space[shift + j] > fXmax)
05315 working_space[shift + j] = fXmax;
05316 working_space[7 * i + 1] = working_space[shift + j];
05317 j += 1;
05318 }
05319 if (fFixPositionY[i] == false) {
05320 if (working_space[shift + j] < fYmin)
05321 working_space[shift + j] = fYmin;
05322 if (working_space[shift + j] > fYmax)
05323 working_space[shift + j] = fYmax;
05324 working_space[7 * i + 2] = working_space[shift + j];
05325 j += 1;
05326 }
05327 if (fFixAmpX1[i] == false) {
05328 if (working_space[shift + j] < 0)
05329 working_space[shift + j] = 0;
05330 working_space[7 * i + 3] = working_space[shift + j];
05331 j += 1;
05332 }
05333 if (fFixAmpY1[i] == false) {
05334 if (working_space[shift + j] < 0)
05335 working_space[shift + j] = 0;
05336 working_space[7 * i + 4] = working_space[shift + j];
05337 j += 1;
05338 }
05339 if (fFixPositionX1[i] == false) {
05340 if (working_space[shift + j] < fXmin)
05341 working_space[shift + j] = fXmin;
05342 if (working_space[shift + j] > fXmax)
05343 working_space[shift + j] = fXmax;
05344 working_space[7 * i + 5] = working_space[shift + j];
05345 j += 1;
05346 }
05347 if (fFixPositionY1[i] == false) {
05348 if (working_space[shift + j] < fYmin)
05349 working_space[shift + j] = fYmin;
05350 if (working_space[shift + j] > fYmax)
05351 working_space[shift + j] = fYmax;
05352 working_space[7 * i + 6] = working_space[shift + j];
05353 j += 1;
05354 }
05355 }
05356 if (fFixSigmaX == false) {
05357 if (working_space[shift + j] < 0.001) {
05358 working_space[shift + j] = 0.001;
05359 }
05360 working_space[peak_vel] = working_space[shift + j];
05361 j += 1;
05362 }
05363 if (fFixSigmaY == false) {
05364 if (working_space[shift + j] < 0.001) {
05365 working_space[shift + j] = 0.001;
05366 }
05367 working_space[peak_vel + 1] = working_space[shift + j];
05368 j += 1;
05369 }
05370 if (fFixRo == false) {
05371 if (working_space[shift + j] < -1) {
05372 working_space[shift + j] = -1;
05373 }
05374 if (working_space[shift + j] > 1) {
05375 working_space[shift + j] = 1;
05376 }
05377 working_space[peak_vel + 2] = working_space[shift + j];
05378 j += 1;
05379 }
05380 if (fFixA0 == false) {
05381 working_space[peak_vel + 3] = working_space[shift + j];
05382 j += 1;
05383 }
05384 if (fFixAx == false) {
05385 working_space[peak_vel + 4] = working_space[shift + j];
05386 j += 1;
05387 }
05388 if (fFixAy == false) {
05389 working_space[peak_vel + 5] = working_space[shift + j];
05390 j += 1;
05391 }
05392 if (fFixTxy == false) {
05393 working_space[peak_vel + 6] = working_space[shift + j];
05394 j += 1;
05395 }
05396 if (fFixSxy == false) {
05397 working_space[peak_vel + 7] = working_space[shift + j];
05398 j += 1;
05399 }
05400 if (fFixTx == false) {
05401 working_space[peak_vel + 8] = working_space[shift + j];
05402 j += 1;
05403 }
05404 if (fFixTy == false) {
05405 working_space[peak_vel + 9] = working_space[shift + j];
05406 j += 1;
05407 }
05408 if (fFixSx == false) {
05409 working_space[peak_vel + 10] = working_space[shift + j];
05410 j += 1;
05411 }
05412 if (fFixSy == false) {
05413 working_space[peak_vel + 11] = working_space[shift + j];
05414 j += 1;
05415 }
05416 if (fFixBx == false) {
05417 if (TMath::Abs(working_space[shift + j]) < 0.001) {
05418 if (working_space[shift + j] < 0)
05419 working_space[shift + j] = -0.001;
05420 else
05421 working_space[shift + j] = 0.001;
05422 }
05423 working_space[peak_vel + 12] = working_space[shift + j];
05424 j += 1;
05425 }
05426 if (fFixBy == false) {
05427 if (TMath::Abs(working_space[shift + j]) < 0.001) {
05428 if (working_space[shift + j] < 0)
05429 working_space[shift + j] = -0.001;
05430 else
05431 working_space[shift + j] = 0.001;
05432 }
05433 working_space[peak_vel + 13] = working_space[shift + j];
05434 j += 1;
05435 }
05436 chi2 = 0;
05437 for (i1 = fXmin; i1 <= fXmax; i1++) {
05438 for (i2 = fYmin; i2 <= fYmax; i2++) {
05439 yw = source[i1][i2];
05440 ywm = yw;
05441 f = Shape2(fNPeaks, (Double_t) i1,
05442 (Double_t) i2, working_space,
05443 working_space[peak_vel],
05444 working_space[peak_vel + 1],
05445 working_space[peak_vel + 2],
05446 working_space[peak_vel + 3],
05447 working_space[peak_vel + 4],
05448 working_space[peak_vel + 5],
05449 working_space[peak_vel + 6],
05450 working_space[peak_vel + 7],
05451 working_space[peak_vel + 8],
05452 working_space[peak_vel + 9],
05453 working_space[peak_vel + 10],
05454 working_space[peak_vel + 11],
05455 working_space[peak_vel + 12],
05456 working_space[peak_vel + 13]);
05457 if (fStatisticType == kFitOptimChiFuncValues) {
05458 ywm = f;
05459 if (f < 0.00001)
05460 ywm = 0.00001;
05461 }
05462 if (fStatisticType == kFitOptimMaxLikelihood) {
05463 if (f > 0.00001)
05464 chi2 += yw * TMath::Log(f) - f;
05465 }
05466
05467 else {
05468 if (ywm != 0)
05469 chi2 += (yw - f) * (yw - f) / ywm;
05470 }
05471 }
05472 }
05473 if ((chi2 < chi_min
05474 && fStatisticType != kFitOptimMaxLikelihood)
05475 || (chi2 > chi_min
05476 && fStatisticType == kFitOptimMaxLikelihood)) {
05477 pmin = pi, chi_min = chi2;
05478 }
05479
05480 else
05481 flag = 1;
05482 if (pi == 0.1)
05483 chi_min = chi2;
05484 chi = chi_min;
05485 }
05486 if (pmin != 0.1) {
05487 for (j = 0; j < size; j++) {
05488 working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j];
05489 }
05490 for (i = 0, j = 0; i < fNPeaks; i++) {
05491 if (fFixAmp[i] == false) {
05492 if (working_space[shift + j] < 0)
05493 working_space[shift + j] = 0;
05494 working_space[7 * i] = working_space[shift + j];
05495 j += 1;
05496 }
05497 if (fFixPositionX[i] == false) {
05498 if (working_space[shift + j] < fXmin)
05499 working_space[shift + j] = fXmin;
05500 if (working_space[shift + j] > fXmax)
05501 working_space[shift + j] = fXmax;
05502 working_space[7 * i + 1] = working_space[shift + j];
05503 j += 1;
05504 }
05505 if (fFixPositionY[i] == false) {
05506 if (working_space[shift + j] < fYmin)
05507 working_space[shift + j] = fYmin;
05508 if (working_space[shift + j] > fYmax)
05509 working_space[shift + j] = fYmax;
05510 working_space[7 * i + 2] = working_space[shift + j];
05511 j += 1;
05512 }
05513 if (fFixAmpX1[i] == false) {
05514 if (working_space[shift + j] < 0)
05515 working_space[shift + j] = 0;
05516 working_space[7 * i + 3] = working_space[shift + j];
05517 j += 1;
05518 }
05519 if (fFixAmpY1[i] == false) {
05520 if (working_space[shift + j] < 0)
05521 working_space[shift + j] = 0;
05522 working_space[7 * i + 4] = working_space[shift + j];
05523 j += 1;
05524 }
05525 if (fFixPositionX1[i] == false) {
05526 if (working_space[shift + j] < fXmin)
05527 working_space[shift + j] = fXmin;
05528 if (working_space[shift + j] > fXmax)
05529 working_space[shift + j] = fXmax;
05530 working_space[7 * i + 5] = working_space[shift + j];
05531 j += 1;
05532 }
05533 if (fFixPositionY1[i] == false) {
05534 if (working_space[shift + j] < fYmin)
05535 working_space[shift + j] = fYmin;
05536 if (working_space[shift + j] > fYmax)
05537 working_space[shift + j] = fYmax;
05538 working_space[7 * i + 6] = working_space[shift + j];
05539 j += 1;
05540 }
05541 }
05542 if (fFixSigmaX == false) {
05543 if (working_space[shift + j] < 0.001) {
05544 working_space[shift + j] = 0.001;
05545 }
05546 working_space[peak_vel] = working_space[shift + j];
05547 j += 1;
05548 }
05549 if (fFixSigmaY == false) {
05550 if (working_space[shift + j] < 0.001) {
05551 working_space[shift + j] = 0.001;
05552 }
05553 working_space[peak_vel + 1] = working_space[shift + j];
05554 j += 1;
05555 }
05556 if (fFixRo == false) {
05557 if (working_space[shift + j] < -1) {
05558 working_space[shift + j] = -1;
05559 }
05560 if (working_space[shift + j] > 1) {
05561 working_space[shift + j] = 1;
05562 }
05563 working_space[peak_vel + 2] = working_space[shift + j];
05564 j += 1;
05565 }
05566 if (fFixA0 == false) {
05567 working_space[peak_vel + 3] = working_space[shift + j];
05568 j += 1;
05569 }
05570 if (fFixAx == false) {
05571 working_space[peak_vel + 4] = working_space[shift + j];
05572 j += 1;
05573 }
05574 if (fFixAy == false) {
05575 working_space[peak_vel + 5] = working_space[shift + j];
05576 j += 1;
05577 }
05578 if (fFixTxy == false) {
05579 working_space[peak_vel + 6] = working_space[shift + j];
05580 j += 1;
05581 }
05582 if (fFixSxy == false) {
05583 working_space[peak_vel + 7] = working_space[shift + j];
05584 j += 1;
05585 }
05586 if (fFixTx == false) {
05587 working_space[peak_vel + 8] = working_space[shift + j];
05588 j += 1;
05589 }
05590 if (fFixTy == false) {
05591 working_space[peak_vel + 9] = working_space[shift + j];
05592 j += 1;
05593 }
05594 if (fFixSx == false) {
05595 working_space[peak_vel + 10] = working_space[shift + j];
05596 j += 1;
05597 }
05598 if (fFixSy == false) {
05599 working_space[peak_vel + 11] = working_space[shift + j];
05600 j += 1;
05601 }
05602 if (fFixBx == false) {
05603 if (TMath::Abs(working_space[shift + j]) < 0.001) {
05604 if (working_space[shift + j] < 0)
05605 working_space[shift + j] = -0.001;
05606 else
05607 working_space[shift + j] = 0.001;
05608 }
05609 working_space[peak_vel + 12] = working_space[shift + j];
05610 j += 1;
05611 }
05612 if (fFixBy == false) {
05613 if (TMath::Abs(working_space[shift + j]) < 0.001) {
05614 if (working_space[shift + j] < 0)
05615 working_space[shift + j] = -0.001;
05616 else
05617 working_space[shift + j] = 0.001;
05618 }
05619 working_space[peak_vel + 13] = working_space[shift + j];
05620 j += 1;
05621 }
05622 chi = chi_min;
05623 }
05624 }
05625
05626 else {
05627 for (j = 0; j < size; j++) {
05628 working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j];
05629 }
05630 for (i = 0, j = 0; i < fNPeaks; i++) {
05631 if (fFixAmp[i] == false) {
05632 if (working_space[shift + j] < 0)
05633 working_space[shift + j] = 0;
05634 working_space[7 * i] = working_space[shift + j];
05635 j += 1;
05636 }
05637 if (fFixPositionX[i] == false) {
05638 if (working_space[shift + j] < fXmin)
05639 working_space[shift + j] = fXmin;
05640 if (working_space[shift + j] > fXmax)
05641 working_space[shift + j] = fXmax;
05642 working_space[7 * i + 1] = working_space[shift + j];
05643 j += 1;
05644 }
05645 if (fFixPositionY[i] == false) {
05646 if (working_space[shift + j] < fYmin)
05647 working_space[shift + j] = fYmin;
05648 if (working_space[shift + j] > fYmax)
05649 working_space[shift + j] = fYmax;
05650 working_space[7 * i + 2] = working_space[shift + j];
05651 j += 1;
05652 }
05653 if (fFixAmpX1[i] == false) {
05654 if (working_space[shift + j] < 0)
05655 working_space[shift + j] = 0;
05656 working_space[7 * i + 3] = working_space[shift + j];
05657 j += 1;
05658 }
05659 if (fFixAmpY1[i] == false) {
05660 if (working_space[shift + j] < 0)
05661 working_space[shift + j] = 0;
05662 working_space[7 * i + 4] = working_space[shift + j];
05663 j += 1;
05664 }
05665 if (fFixPositionX1[i] == false) {
05666 if (working_space[shift + j] < fXmin)
05667 working_space[shift + j] = fXmin;
05668 if (working_space[shift + j] > fXmax)
05669 working_space[shift + j] = fXmax;
05670 working_space[7 * i + 5] = working_space[shift + j];
05671 j += 1;
05672 }
05673 if (fFixPositionY1[i] == false) {
05674 if (working_space[shift + j] < fYmin)
05675 working_space[shift + j] = fYmin;
05676 if (working_space[shift + j] > fYmax)
05677 working_space[shift + j] = fYmax;
05678 working_space[7 * i + 6] = working_space[shift + j];
05679 j += 1;
05680 }
05681 }
05682 if (fFixSigmaX == false) {
05683 if (working_space[shift + j] < 0.001) {
05684 working_space[shift + j] = 0.001;
05685 }
05686 working_space[peak_vel] = working_space[shift + j];
05687 j += 1;
05688 }
05689 if (fFixSigmaY == false) {
05690 if (working_space[shift + j] < 0.001) {
05691 working_space[shift + j] = 0.001;
05692 }
05693 working_space[peak_vel + 1] = working_space[shift + j];
05694 j += 1;
05695 }
05696 if (fFixRo == false) {
05697 if (working_space[shift + j] < -1) {
05698 working_space[shift + j] = -1;
05699 }
05700 if (working_space[shift + j] > 1) {
05701 working_space[shift + j] = 1;
05702 }
05703 working_space[peak_vel + 2] = working_space[shift + j];
05704 j += 1;
05705 }
05706 if (fFixA0 == false) {
05707 working_space[peak_vel + 3] = working_space[shift + j];
05708 j += 1;
05709 }
05710 if (fFixAx == false) {
05711 working_space[peak_vel + 4] = working_space[shift + j];
05712 j += 1;
05713 }
05714 if (fFixAy == false) {
05715 working_space[peak_vel + 5] = working_space[shift + j];
05716 j += 1;
05717 }
05718 if (fFixTxy == false) {
05719 working_space[peak_vel + 6] = working_space[shift + j];
05720 j += 1;
05721 }
05722 if (fFixSxy == false) {
05723 working_space[peak_vel + 7] = working_space[shift + j];
05724 j += 1;
05725 }
05726 if (fFixTx == false) {
05727 working_space[peak_vel + 8] = working_space[shift + j];
05728 j += 1;
05729 }
05730 if (fFixTy == false) {
05731 working_space[peak_vel + 9] = working_space[shift + j];
05732 j += 1;
05733 }
05734 if (fFixSx == false) {
05735 working_space[peak_vel + 10] = working_space[shift + j];
05736 j += 1;
05737 }
05738 if (fFixSy == false) {
05739 working_space[peak_vel + 11] = working_space[shift + j];
05740 j += 1;
05741 }
05742 if (fFixBx == false) {
05743 if (TMath::Abs(working_space[shift + j]) < 0.001) {
05744 if (working_space[shift + j] < 0)
05745 working_space[shift + j] = -0.001;
05746 else
05747 working_space[shift + j] = 0.001;
05748 }
05749 working_space[peak_vel + 12] = working_space[shift + j];
05750 j += 1;
05751 }
05752 if (fFixBy == false) {
05753 if (TMath::Abs(working_space[shift + j]) < 0.001) {
05754 if (working_space[shift + j] < 0)
05755 working_space[shift + j] = -0.001;
05756 else
05757 working_space[shift + j] = 0.001;
05758 }
05759 working_space[peak_vel + 13] = working_space[shift + j];
05760 j += 1;
05761 }
05762 chi = 0;
05763 for (i1 = fXmin; i1 <= fXmax; i1++) {
05764 for (i2 = fYmin; i2 <= fYmax; i2++) {
05765 yw = source[i1][i2];
05766 ywm = yw;
05767 f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
05768 working_space, working_space[peak_vel],
05769 working_space[peak_vel + 1],
05770 working_space[peak_vel + 2],
05771 working_space[peak_vel + 3],
05772 working_space[peak_vel + 4],
05773 working_space[peak_vel + 5],
05774 working_space[peak_vel + 6],
05775 working_space[peak_vel + 7],
05776 working_space[peak_vel + 8],
05777 working_space[peak_vel + 9],
05778 working_space[peak_vel + 10],
05779 working_space[peak_vel + 11],
05780 working_space[peak_vel + 12],
05781 working_space[peak_vel + 13]);
05782 if (fStatisticType == kFitOptimChiFuncValues) {
05783 ywm = f;
05784 if (f < 0.00001)
05785 ywm = 0.00001;
05786 }
05787 if (fStatisticType == kFitOptimMaxLikelihood) {
05788 if (f > 0.00001)
05789 chi += yw * TMath::Log(f) - f;
05790 }
05791
05792 else {
05793 if (ywm != 0)
05794 chi += (yw - f) * (yw - f) / ywm;
05795 }
05796 }
05797 }
05798 }
05799 chi2 = chi;
05800 chi = TMath::Sqrt(TMath::Abs(chi));
05801 if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
05802 alpha = alpha * chi_opt / (2 * chi);
05803
05804 else if (fAlphaOptim == kFitAlphaOptimal)
05805 alpha = alpha / 10.0;
05806 iter += 1;
05807 regul_cycle += 1;
05808 } while (((chi > chi_opt
05809 && fStatisticType != kFitOptimMaxLikelihood)
05810 || (chi < chi_opt
05811 && fStatisticType == kFitOptimMaxLikelihood))
05812 && regul_cycle < kFitNumRegulCycles);
05813 for (j = 0; j < size; j++) {
05814 working_space[4 * shift + j] = 0;
05815 working_space[2 * shift + j] = 0;
05816 }
05817 for (i1 = fXmin, chi_cel = 0; i1 <= fXmax; i1++) {
05818 for (i2 = fYmin; i2 <= fYmax; i2++) {
05819 yw = source[i1][i2];
05820 if (yw == 0)
05821 yw = 1;
05822 f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
05823 working_space, working_space[peak_vel],
05824 working_space[peak_vel + 1],
05825 working_space[peak_vel + 2],
05826 working_space[peak_vel + 3],
05827 working_space[peak_vel + 4],
05828 working_space[peak_vel + 5],
05829 working_space[peak_vel + 6],
05830 working_space[peak_vel + 7],
05831 working_space[peak_vel + 8],
05832 working_space[peak_vel + 9],
05833 working_space[peak_vel + 10],
05834 working_space[peak_vel + 11],
05835 working_space[peak_vel + 12],
05836 working_space[peak_vel + 13]);
05837 chi_opt = (yw - f) * (yw - f) / yw;
05838 chi_cel += (yw - f) * (yw - f) / yw;
05839
05840
05841 for (j = 0, k = 0; j < fNPeaks; j++) {
05842 if (fFixAmp[j] == false) {
05843 a = Deramp2((Double_t) i1, (Double_t) i2,
05844 working_space[7 * j + 1],
05845 working_space[7 * j + 2],
05846 working_space[peak_vel],
05847 working_space[peak_vel + 1],
05848 working_space[peak_vel + 2],
05849 working_space[peak_vel + 6],
05850 working_space[peak_vel + 7],
05851 working_space[peak_vel + 12],
05852 working_space[peak_vel + 13]);
05853 if (yw != 0) {
05854 working_space[2 * shift + k] += chi_opt;
05855 b = a * a / yw;
05856 working_space[4 * shift + k] += b;
05857 }
05858 k += 1;
05859 }
05860 if (fFixPositionX[j] == false) {
05861 a = Deri02((Double_t) i1, (Double_t) i2,
05862 working_space[7 * j],
05863 working_space[7 * j + 1],
05864 working_space[7 * j + 2],
05865 working_space[peak_vel],
05866 working_space[peak_vel + 1],
05867 working_space[peak_vel + 2],
05868 working_space[peak_vel + 6],
05869 working_space[peak_vel + 7],
05870 working_space[peak_vel + 12],
05871 working_space[peak_vel + 13]);
05872 if (yw != 0) {
05873 working_space[2 * shift + k] += chi_opt;
05874 b = a * a / yw;
05875 working_space[4 * shift + k] += b;
05876 }
05877 k += 1;
05878 }
05879 if (fFixPositionY[j] == false) {
05880 a = Derj02((Double_t) i1, (Double_t) i2,
05881 working_space[7 * j],
05882 working_space[7 * j + 1],
05883 working_space[7 * j + 2],
05884 working_space[peak_vel],
05885 working_space[peak_vel + 1],
05886 working_space[peak_vel + 2],
05887 working_space[peak_vel + 6],
05888 working_space[peak_vel + 7],
05889 working_space[peak_vel + 12],
05890 working_space[peak_vel + 13]);
05891 if (yw != 0) {
05892 working_space[2 * shift + k] += chi_opt;
05893 b = a * a / yw;
05894 working_space[4 * shift + k] += b;
05895 }
05896 k += 1;
05897 }
05898 if (fFixAmpX1[j] == false) {
05899 a = Derampx((Double_t) i1, working_space[7 * j + 5],
05900 working_space[peak_vel],
05901 working_space[peak_vel + 8],
05902 working_space[peak_vel + 10],
05903 working_space[peak_vel + 12]);
05904 if (yw != 0) {
05905 working_space[2 * shift + k] += chi_opt;
05906 b = a * a / yw;
05907 working_space[4 * shift + k] += b;
05908 }
05909 k += 1;
05910 }
05911 if (fFixAmpY1[j] == false) {
05912 a = Derampx((Double_t) i2, working_space[7 * j + 6],
05913 working_space[peak_vel + 1],
05914 working_space[peak_vel + 9],
05915 working_space[peak_vel + 11],
05916 working_space[peak_vel + 13]);
05917 if (yw != 0) {
05918 working_space[2 * shift + k] += chi_opt;
05919 b = a * a / yw;
05920 working_space[4 * shift + k] += b;
05921 }
05922 k += 1;
05923 }
05924 if (fFixPositionX1[j] == false) {
05925 a = Deri01((Double_t) i1, working_space[7 * j + 3],
05926 working_space[7 * j + 5],
05927 working_space[peak_vel],
05928 working_space[peak_vel + 8],
05929 working_space[peak_vel + 10],
05930 working_space[peak_vel + 12]);
05931 if (yw != 0) {
05932 working_space[2 * shift + k] += chi_opt;
05933 b = a * a / yw;
05934 working_space[4 * shift + k] += b;
05935 }
05936 k += 1;
05937 }
05938 if (fFixPositionY1[j] == false) {
05939 a = Deri01((Double_t) i2, working_space[7 * j + 4],
05940 working_space[7 * j + 6],
05941 working_space[peak_vel + 1],
05942 working_space[peak_vel + 9],
05943 working_space[peak_vel + 11],
05944 working_space[peak_vel + 13]);
05945 if (yw != 0) {
05946 working_space[2 * shift + k] += chi_opt;
05947 b = a * a / yw;
05948 working_space[4 * shift + k] += b;
05949 }
05950 k += 1;
05951 }
05952 }
05953 if (fFixSigmaX == false) {
05954 a = Dersigmax(fNPeaks, (Double_t) i1, (Double_t) i2,
05955 working_space, working_space[peak_vel],
05956 working_space[peak_vel + 1],
05957 working_space[peak_vel + 2],
05958 working_space[peak_vel + 6],
05959 working_space[peak_vel + 7],
05960 working_space[peak_vel + 8],
05961 working_space[peak_vel + 10],
05962 working_space[peak_vel + 12],
05963 working_space[peak_vel + 13]);
05964 if (yw != 0) {
05965 working_space[2 * shift + k] += chi_opt;
05966 b = a * a / yw;
05967 working_space[4 * shift + k] += b;
05968 }
05969 k += 1;
05970 }
05971 if (fFixSigmaY == false) {
05972 a = Dersigmay(fNPeaks, (Double_t) i1, (Double_t) i2,
05973 working_space, working_space[peak_vel],
05974 working_space[peak_vel + 1],
05975 working_space[peak_vel + 2],
05976 working_space[peak_vel + 6],
05977 working_space[peak_vel + 7],
05978 working_space[peak_vel + 9],
05979 working_space[peak_vel + 11],
05980 working_space[peak_vel + 12],
05981 working_space[peak_vel + 13]);
05982 if (yw != 0) {
05983 working_space[2 * shift + k] += chi_opt;
05984 b = a * a / yw;
05985 working_space[4 * shift + k] += b;
05986 }
05987 k += 1;
05988 }
05989 if (fFixRo == false) {
05990 a = Derro(fNPeaks, (Double_t) i1, (Double_t) i2,
05991 working_space, working_space[peak_vel],
05992 working_space[peak_vel + 1],
05993 working_space[peak_vel + 2]);
05994 if (yw != 0) {
05995 working_space[2 * shift + k] += chi_opt;
05996 b = a * a / yw;
05997 working_space[4 * shift + k] += b;
05998 }
05999 k += 1;
06000 }
06001 if (fFixA0 == false) {
06002 a = 1.;
06003 if (yw != 0) {
06004 working_space[2 * shift + k] += chi_opt;
06005 b = a * a / yw;
06006 working_space[4 * shift + k] += b;
06007 }
06008 k += 1;
06009 }
06010 if (fFixAx == false) {
06011 a = i1;
06012 if (yw != 0) {
06013 working_space[2 * shift + k] += chi_opt;
06014 b = a * a / yw;
06015 working_space[4 * shift + k] += b;
06016 }
06017 k += 1;
06018 }
06019 if (fFixAy == false) {
06020 a = i2;
06021 if (yw != 0) {
06022 working_space[2 * shift + k] += chi_opt;
06023 b = a * a / yw;
06024 working_space[4 * shift + k] += b;
06025 }
06026 k += 1;
06027 }
06028 if (fFixTxy == false) {
06029 a = Dertxy(fNPeaks, (Double_t) i1, (Double_t) i2,
06030 working_space, working_space[peak_vel],
06031 working_space[peak_vel + 1],
06032 working_space[peak_vel + 12],
06033 working_space[peak_vel + 13]);
06034 if (yw != 0) {
06035 working_space[2 * shift + k] += chi_opt;
06036 b = a * a / yw;
06037 working_space[4 * shift + k] += b;
06038 }
06039 k += 1;
06040 }
06041 if (fFixSxy == false) {
06042 a = Dersxy(fNPeaks, (Double_t) i1, (Double_t) i2,
06043 working_space, working_space[peak_vel],
06044 working_space[peak_vel + 1]);
06045 if (yw != 0) {
06046 working_space[2 * shift + k] += chi_opt;
06047 b = a * a / yw;
06048 working_space[4 * shift + k] += b;
06049 }
06050 k += 1;
06051 }
06052 if (fFixTx == false) {
06053 a = Dertx(fNPeaks, (Double_t) i1, working_space,
06054 working_space[peak_vel],
06055 working_space[peak_vel + 12]);
06056 if (yw != 0) {
06057 working_space[2 * shift + k] += chi_opt;
06058 b = a * a / yw;
06059 working_space[4 * shift + k] += b;
06060 }
06061 k += 1;
06062 }
06063 if (fFixTy == false) {
06064 a = Derty(fNPeaks, (Double_t) i2, working_space,
06065 working_space[peak_vel + 1],
06066 working_space[peak_vel + 13]);
06067 if (yw != 0) {
06068 working_space[2 * shift + k] += chi_opt;
06069 b = a * a / yw;
06070 working_space[4 * shift + k] += b;
06071 }
06072 k += 1;
06073 }
06074 if (fFixSx == false) {
06075 a = Dersx(fNPeaks, (Double_t) i1, working_space,
06076 working_space[peak_vel]);
06077 if (yw != 0) {
06078 working_space[2 * shift + k] += chi_opt;
06079 b = a * a / yw;
06080 working_space[4 * shift + k] += b;
06081 }
06082 k += 1;
06083 }
06084 if (fFixSy == false) {
06085 a = Dersy(fNPeaks, (Double_t) i2, working_space,
06086 working_space[peak_vel + 1]);
06087 if (yw != 0) {
06088 working_space[2 * shift + k] += chi_opt;
06089 b = a * a / yw;
06090 working_space[4 * shift + k] += b;
06091 }
06092 k += 1;
06093 }
06094 if (fFixBx == false) {
06095 a = Derbx(fNPeaks, (Double_t) i1, (Double_t) i2,
06096 working_space, working_space[peak_vel],
06097 working_space[peak_vel + 1],
06098 working_space[peak_vel + 6],
06099 working_space[peak_vel + 8],
06100 working_space[peak_vel + 12],
06101 working_space[peak_vel + 13]);
06102 if (yw != 0) {
06103 working_space[2 * shift + k] += chi_opt;
06104 b = a * a / yw;
06105 working_space[4 * shift + k] += b;
06106 }
06107 k += 1;
06108 }
06109 if (fFixBy == false) {
06110 a = Derby(fNPeaks, (Double_t) i1, (Double_t) i2,
06111 working_space, working_space[peak_vel],
06112 working_space[peak_vel + 1],
06113 working_space[peak_vel + 6],
06114 working_space[peak_vel + 8],
06115 working_space[peak_vel + 12],
06116 working_space[peak_vel + 13]);
06117 if (yw != 0) {
06118 working_space[2 * shift + k] += chi_opt;
06119 b = a * a / yw;
06120 working_space[4 * shift + k] += b;
06121 }
06122 k += 1;
06123 }
06124 }
06125 }
06126 }
06127 b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
06128 chi_er = chi_cel / b;
06129 for (i = 0, j = 0; i < fNPeaks; i++) {
06130 fVolume[i] =
06131 Volume(working_space[7 * i], working_space[peak_vel],
06132 working_space[peak_vel + 1], working_space[peak_vel + 2]);
06133 if (fVolume[i] > 0) {
06134 c = 0;
06135 if (fFixAmp[i] == false) {
06136 a = Derpa2(working_space[peak_vel],
06137 working_space[peak_vel + 1],
06138 working_space[peak_vel + 2]);
06139 b = working_space[4 * shift + j];
06140 if (b == 0)
06141 b = 1;
06142
06143 else
06144 b = 1 / b;
06145 c = c + a * a * b;
06146 }
06147 if (fFixSigmaX == false) {
06148 a = Derpsigmax(working_space[shift + j],
06149 working_space[peak_vel + 1],
06150 working_space[peak_vel + 2]);
06151 b = working_space[4 * shift + peak_vel];
06152 if (b == 0)
06153 b = 1;
06154
06155 else
06156 b = 1 / b;
06157 c = c + a * a * b;
06158 }
06159 if (fFixSigmaY == false) {
06160 a = Derpsigmay(working_space[shift + j],
06161 working_space[peak_vel],
06162 working_space[peak_vel + 2]);
06163 b = working_space[4 * shift + peak_vel + 1];
06164 if (b == 0)
06165 b = 1;
06166
06167 else
06168 b = 1 / b;
06169 c = c + a * a * b;
06170 }
06171 if (fFixRo == false) {
06172 a = Derpro(working_space[shift + j], working_space[peak_vel],
06173 working_space[peak_vel + 1],
06174 working_space[peak_vel + 2]);
06175 b = working_space[4 * shift + peak_vel + 2];
06176 if (b == 0)
06177 b = 1;
06178
06179 else
06180 b = 1 / b;
06181 c = c + a * a * b;
06182 }
06183 fVolumeErr[i] = TMath::Sqrt(TMath::Abs(chi_er * c));
06184 }
06185
06186 else {
06187 fVolumeErr[i] = 0;
06188 }
06189 if (fFixAmp[i] == false) {
06190 fAmpCalc[i] = working_space[shift + j];
06191 if (working_space[3 * shift + j] != 0)
06192 fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06193 j += 1;
06194 }
06195
06196 else {
06197 fAmpCalc[i] = fAmpInit[i];
06198 fAmpErr[i] = 0;
06199 }
06200 if (fFixPositionX[i] == false) {
06201 fPositionCalcX[i] = working_space[shift + j];
06202 if (working_space[3 * shift + j] != 0)
06203 fPositionErrX[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06204 j += 1;
06205 }
06206
06207 else {
06208 fPositionCalcX[i] = fPositionInitX[i];
06209 fPositionErrX[i] = 0;
06210 }
06211 if (fFixPositionY[i] == false) {
06212 fPositionCalcY[i] = working_space[shift + j];
06213 if (working_space[3 * shift + j] != 0)
06214 fPositionErrY[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06215 j += 1;
06216 }
06217
06218 else {
06219 fPositionCalcY[i] = fPositionInitY[i];
06220 fPositionErrY[i] = 0;
06221 }
06222 if (fFixAmpX1[i] == false) {
06223 fAmpCalcX1[i] = working_space[shift + j];
06224 if (working_space[3 * shift + j] != 0)
06225 fAmpErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06226 j += 1;
06227 }
06228
06229 else {
06230 fAmpCalcX1[i] = fAmpInitX1[i];
06231 fAmpErrX1[i] = 0;
06232 }
06233 if (fFixAmpY1[i] == false) {
06234 fAmpCalcY1[i] = working_space[shift + j];
06235 if (working_space[3 * shift + j] != 0)
06236 fAmpErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06237 j += 1;
06238 }
06239
06240 else {
06241 fAmpCalcY1[i] = fAmpInitY1[i];
06242 fAmpErrY1[i] = 0;
06243 }
06244 if (fFixPositionX1[i] == false) {
06245 fPositionCalcX1[i] = working_space[shift + j];
06246 if (working_space[3 * shift + j] != 0)
06247 fPositionErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06248 j += 1;
06249 }
06250
06251 else {
06252 fPositionCalcX1[i] = fPositionInitX1[i];
06253 fPositionErrX1[i] = 0;
06254 }
06255 if (fFixPositionY1[i] == false) {
06256 fPositionCalcY1[i] = working_space[shift + j];
06257 if (working_space[3 * shift + j] != 0)
06258 fPositionErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06259 j += 1;
06260 }
06261
06262 else {
06263 fPositionCalcY1[i] = fPositionInitY1[i];
06264 fPositionErrY1[i] = 0;
06265 }
06266 }
06267 if (fFixSigmaX == false) {
06268 fSigmaCalcX = working_space[shift + j];
06269 if (working_space[3 * shift + j] != 0)
06270 fSigmaErrX = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06271 j += 1;
06272 }
06273
06274 else {
06275 fSigmaCalcX = fSigmaInitX;
06276 fSigmaErrX = 0;
06277 }
06278 if (fFixSigmaY == false) {
06279 fSigmaCalcY = working_space[shift + j];
06280 if (working_space[3 * shift + j] != 0)
06281 fSigmaErrY = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06282 j += 1;
06283 }
06284
06285 else {
06286 fSigmaCalcY = fSigmaInitY;
06287 fSigmaErrY = 0;
06288 }
06289 if (fFixRo == false) {
06290 fRoCalc = working_space[shift + j];
06291 if (working_space[3 * shift + j] != 0)
06292 fRoErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06293 j += 1;
06294 }
06295
06296 else {
06297 fRoCalc = fRoInit;
06298 fRoErr = 0;
06299 }
06300 if (fFixA0 == false) {
06301 fA0Calc = working_space[shift + j];
06302 if (working_space[3 * shift + j] != 0)
06303 fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06304 j += 1;
06305 }
06306
06307 else {
06308 fA0Calc = fA0Init;
06309 fA0Err = 0;
06310 }
06311 if (fFixAx == false) {
06312 fAxCalc = working_space[shift + j];
06313 if (working_space[3 * shift + j] != 0)
06314 fAxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06315 j += 1;
06316 }
06317
06318 else {
06319 fAxCalc = fAxInit;
06320 fAxErr = 0;
06321 }
06322 if (fFixAy == false) {
06323 fAyCalc = working_space[shift + j];
06324 if (working_space[3 * shift + j] != 0)
06325 fAyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06326 j += 1;
06327 }
06328
06329 else {
06330 fAyCalc = fAyInit;
06331 fAyErr = 0;
06332 }
06333 if (fFixTxy == false) {
06334 fTxyCalc = working_space[shift + j];
06335 if (working_space[3 * shift + j] != 0)
06336 fTxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06337 j += 1;
06338 }
06339
06340 else {
06341 fTxyCalc = fTxyInit;
06342 fTxyErr = 0;
06343 }
06344 if (fFixSxy == false) {
06345 fSxyCalc = working_space[shift + j];
06346 if (working_space[3 * shift + j] != 0)
06347 fSxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06348 j += 1;
06349 }
06350
06351 else {
06352 fSxyCalc = fSxyInit;
06353 fSxyErr = 0;
06354 }
06355 if (fFixTx == false) {
06356 fTxCalc = working_space[shift + j];
06357 if (working_space[3 * shift + j] != 0)
06358 fTxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06359 j += 1;
06360 }
06361
06362 else {
06363 fTxCalc = fTxInit;
06364 fTxErr = 0;
06365 }
06366 if (fFixTy == false) {
06367 fTyCalc = working_space[shift + j];
06368 if (working_space[3 * shift + j] != 0)
06369 fTyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06370 j += 1;
06371 }
06372
06373 else {
06374 fTyCalc = fTyInit;
06375 fTyErr = 0;
06376 }
06377 if (fFixSx == false) {
06378 fSxCalc = working_space[shift + j];
06379 if (working_space[3 * shift + j] != 0)
06380 fSxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06381 j += 1;
06382 }
06383
06384 else {
06385 fSxCalc = fSxInit;
06386 fSxErr = 0;
06387 }
06388 if (fFixSy == false) {
06389 fSyCalc = working_space[shift + j];
06390 if (working_space[3 * shift + j] != 0)
06391 fSyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06392 j += 1;
06393 }
06394
06395 else {
06396 fSyCalc = fSyInit;
06397 fSyErr = 0;
06398 }
06399 if (fFixBx == false) {
06400 fBxCalc = working_space[shift + j];
06401 if (working_space[3 * shift + j] != 0)
06402 fBxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06403 j += 1;
06404 }
06405
06406 else {
06407 fBxCalc = fBxInit;
06408 fBxErr = 0;
06409 }
06410 if (fFixBy == false) {
06411 fByCalc = working_space[shift + j];
06412 if (working_space[3 * shift + j] != 0)
06413 fByErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));
06414 j += 1;
06415 }
06416
06417 else {
06418 fByCalc = fByInit;
06419 fByErr = 0;
06420 }
06421 b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
06422 fChi = chi_cel / b;
06423 for (i1 = fXmin; i1 <= fXmax; i1++) {
06424 for (i2 = fYmin; i2 <= fYmax; i2++) {
06425 f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
06426 working_space, working_space[peak_vel],
06427 working_space[peak_vel + 1],
06428 working_space[peak_vel + 2],
06429 working_space[peak_vel + 3],
06430 working_space[peak_vel + 4],
06431 working_space[peak_vel + 5],
06432 working_space[peak_vel + 6],
06433 working_space[peak_vel + 7],
06434 working_space[peak_vel + 8],
06435 working_space[peak_vel + 9],
06436 working_space[peak_vel + 10],
06437 working_space[peak_vel + 11],
06438 working_space[peak_vel + 12],
06439 working_space[peak_vel + 13]);
06440 source[i1][i2] = f;
06441
06442 }
06443 }
06444 for (i = 0; i < size; i++) delete [] working_matrix[i];
06445 delete [] working_matrix;
06446 delete [] working_space;
06447 return;
06448 }
06449
06450 void TSpectrum2Fit::SetFitParameters(Int_t xmin,Int_t xmax,Int_t ymin,Int_t ymax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor)
06451 {
06452
06453
06454
06455
06456
06457
06458
06459
06460
06461
06462
06463
06464
06465 if(xmin<0 || xmax <= xmin || ymin<0 || ymax <= ymin){
06466 Error("SetFitParameters", "Wrong range");
06467 return;
06468 }
06469 if (numberIterations <= 0){
06470 Error("SetFitParameters","Invalid number of iterations, must be positive");
06471 return;
06472 }
06473 if (alpha <= 0 || alpha > 1){
06474 Error ("SetFitParameters","Invalid step coefficient alpha, must be > than 0 and <=1");
06475 return;
06476 }
06477 if (statisticType != kFitOptimChiCounts
06478 && statisticType != kFitOptimChiFuncValues
06479 && statisticType != kFitOptimMaxLikelihood){
06480 Error("SetFitParameters","Wrong type of statistic");
06481 return;
06482 }
06483 if (alphaOptim != kFitAlphaHalving
06484 && alphaOptim != kFitAlphaOptimal){
06485 Error("SetFitParameters","Wrong optimization algorithm");
06486 return;
06487 }
06488 if (power != kFitPower2 && power != kFitPower4
06489 && power != kFitPower6 && power != kFitPower8
06490 && power != kFitPower10 && power != kFitPower12){
06491 Error("SetFitParameters","Wrong power");
06492 return;
06493 }
06494 if (fitTaylor != kFitTaylorOrderFirst
06495 && fitTaylor != kFitTaylorOrderSecond){
06496 Error("SetFitParameters","Wrong order of Taylor development");
06497 return;
06498 }
06499 fXmin=xmin,fXmax=xmax,fYmin=ymin,fYmax=ymax,fNumberIterations=numberIterations,fAlpha=alpha,fStatisticType=statisticType,fAlphaOptim=alphaOptim,fPower=power,fFitTaylor=fitTaylor;
06500 }
06501
06502
06503 void TSpectrum2Fit::SetPeakParameters(Double_t sigmaX, Bool_t fixSigmaX, Double_t sigmaY, Bool_t fixSigmaY, Double_t ro, Bool_t fixRo, const Float_t *positionInitX, const Bool_t *fixPositionX, const Float_t *positionInitY, const Bool_t *fixPositionY, const Float_t *positionInitX1, const Bool_t *fixPositionX1, const Float_t *positionInitY1, const Bool_t *fixPositionY1, const Float_t *ampInit, const Bool_t *fixAmp, const Float_t *ampInitX1, const Bool_t *fixAmpX1, const Float_t *ampInitY1, const Bool_t *fixAmpY1)
06504 {
06505
06506
06507
06508
06509
06510
06511
06512
06513
06514
06515
06516
06517
06518
06519
06520
06521
06522
06523
06524
06525
06526
06527 if (sigmaX <= 0 || sigmaY <= 0){
06528 Error ("SetPeakParameters","Invalid sigma, must be > than 0");
06529 return;
06530 }
06531 if (ro < -1 || ro > 1){
06532 Error ("SetPeakParameters","Invalid ro, must be from region <-1,1>");
06533 return;
06534 }
06535 Int_t i;
06536 for(i=0; i < fNPeaks; i++){
06537 if(positionInitX[i] < fXmin || positionInitX[i] > fXmax){
06538 Error ("SetPeakParameters","Invalid peak position, must be in the range fXmin, fXmax");
06539 return;
06540 }
06541 if(positionInitY[i] < fYmin || positionInitY[i] > fYmax){
06542 Error ("SetPeakParameters","Invalid peak position, must be in the range fYmin, fYmax");
06543 return;
06544 }
06545 if(positionInitX1[i] < fXmin || positionInitX1[i] > fXmax){
06546 Error ("SetPeakParameters","Invalid ridge position, must be in the range fXmin, fXmax");
06547 return;
06548 }
06549 if(positionInitY1[i] < fYmin || positionInitY1[i] > fYmax){
06550 Error ("SetPeakParameters","Invalid ridge position, must be in the range fYmin, fYmax");
06551 return;
06552 }
06553 if(ampInit[i] < 0){
06554 Error ("SetPeakParameters","Invalid peak amplitude, must be > than 0");
06555 return;
06556 }
06557 if(ampInitX1[i] < 0){
06558 Error ("SetPeakParameters","Invalid x ridge amplitude, must be > than 0");
06559 return;
06560 }
06561 if(ampInitY1[i] < 0){
06562 Error ("SetPeakParameters","Invalid y ridge amplitude, must be > than 0");
06563 return;
06564 }
06565 }
06566 fSigmaInitX = sigmaX, fFixSigmaX = fixSigmaX, fSigmaInitY = sigmaY, fFixSigmaY = fixSigmaY, fRoInit = ro, fFixRo = fixRo;
06567 for(i=0; i < fNPeaks; i++){
06568 fPositionInitX[i] = (Double_t) positionInitX[i];
06569 fFixPositionX[i] = fixPositionX[i];
06570 fPositionInitY[i] = (Double_t) positionInitY[i];
06571 fFixPositionY[i] = fixPositionY[i];
06572 fPositionInitX1[i] = (Double_t) positionInitX1[i];
06573 fFixPositionX1[i] = fixPositionX1[i];
06574 fPositionInitY1[i] = (Double_t) positionInitY1[i];
06575 fFixPositionY1[i] = fixPositionY1[i];
06576 fAmpInit[i] = (Double_t) ampInit[i];
06577 fFixAmp[i] = fixAmp[i];
06578 fAmpInitX1[i] = (Double_t) ampInitX1[i];
06579 fFixAmpX1[i] = fixAmpX1[i];
06580 fAmpInitY1[i] = (Double_t) ampInitY1[i];
06581 fFixAmpY1[i] = fixAmpY1[i];
06582 }
06583 }
06584
06585
06586 void TSpectrum2Fit::SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t axInit, Bool_t fixAx, Double_t ayInit, Bool_t fixAy)
06587 {
06588
06589
06590
06591
06592
06593
06594
06595
06596
06597
06598
06599
06600 fA0Init = a0Init;
06601 fFixA0 = fixA0;
06602 fAxInit = axInit;
06603 fFixAx = fixAx;
06604 fAyInit = ayInit;
06605 fFixAy = fixAy;
06606 }
06607
06608
06609 void TSpectrum2Fit::SetTailParameters(Double_t tInitXY, Bool_t fixTxy, Double_t tInitX, Bool_t fixTx, Double_t tInitY, Bool_t fixTy, Double_t bInitX, Bool_t fixBx, Double_t bInitY, Bool_t fixBy, Double_t sInitXY, Bool_t fixSxy, Double_t sInitX, Bool_t fixSx, Double_t sInitY, Bool_t fixSy)
06610 {
06611
06612
06613
06614
06615
06616
06617
06618
06619
06620
06621
06622
06623
06624
06625
06626
06627
06628
06629
06630
06631
06632 fTxyInit = tInitXY;
06633 fFixTxy = fixTxy;
06634 fTxInit = tInitX;
06635 fFixTx = fixTx;
06636 fTyInit = tInitY;
06637 fFixTy = fixTy;
06638 fBxInit = bInitX;
06639 fFixBx = fixBx;
06640 fByInit = bInitY;
06641 fFixBy = fixBy;
06642 fSxyInit = sInitXY;
06643 fFixSxy = fixSxy;
06644 fSxInit = sInitX;
06645 fFixSx = fixSx;
06646 fSyInit = sInitY;
06647 fFixSy = fixSy;
06648 }
06649
06650
06651 void TSpectrum2Fit::GetPositions(Float_t *positionsX, Float_t *positionsY, Float_t *positionsX1, Float_t *positionsY1)
06652 {
06653
06654
06655
06656
06657
06658
06659
06660
06661
06662 for( Int_t i=0; i < fNPeaks; i++){
06663 positionsX[i] = (Float_t) fPositionCalcX[i];
06664 positionsY[i] = (Float_t) fPositionCalcY[i];
06665 positionsX1[i] = (Float_t) fPositionCalcX1[i];
06666 positionsY1[i] = (Float_t) fPositionCalcY1[i];
06667 }
06668 }
06669
06670
06671 void TSpectrum2Fit::GetPositionErrors(Float_t *positionErrorsX, Float_t *positionErrorsY, Float_t *positionErrorsX1, Float_t *positionErrorsY1)
06672 {
06673
06674
06675
06676
06677
06678
06679
06680
06681
06682
06683 for( Int_t i=0; i < fNPeaks; i++){
06684 positionErrorsX[i] = (Float_t) fPositionErrX[i];
06685 positionErrorsY[i] = (Float_t) fPositionErrY[i];
06686 positionErrorsX1[i] = (Float_t) fPositionErrX1[i];
06687 positionErrorsY1[i] = (Float_t) fPositionErrY1[i];
06688 }
06689 }
06690
06691
06692 void TSpectrum2Fit::GetAmplitudes(Float_t *amplitudes, Float_t *amplitudesX1, Float_t *amplitudesY1)
06693 {
06694
06695
06696
06697
06698
06699
06700
06701
06702
06703 for( Int_t i=0; i < fNPeaks; i++){
06704 amplitudes[i] = (Float_t) fAmpCalc[i];
06705 amplitudesX1[i] = (Float_t) fAmpCalcX1[i];
06706 amplitudesY1[i] = (Float_t) fAmpCalcY1[i];
06707 }
06708 }
06709
06710
06711 void TSpectrum2Fit::GetAmplitudeErrors(Float_t *amplitudeErrors, Float_t *amplitudeErrorsX1, Float_t *amplitudeErrorsY1)
06712 {
06713
06714
06715
06716
06717
06718
06719
06720
06721
06722 for( Int_t i=0; i < fNPeaks; i++){
06723 amplitudeErrors[i] = (Float_t) fAmpErr[i];
06724 amplitudeErrorsX1[i] = (Float_t) fAmpErrX1[i];
06725 amplitudeErrorsY1[i] = (Float_t) fAmpErrY1[i];
06726 }
06727 }
06728
06729
06730 void TSpectrum2Fit::GetVolumes(Float_t *volumes)
06731 {
06732
06733
06734
06735
06736
06737
06738 for( Int_t i=0; i < fNPeaks; i++){
06739 volumes[i] = (Float_t) fVolume[i];
06740 }
06741 }
06742
06743
06744 void TSpectrum2Fit::GetVolumeErrors(Float_t *volumeErrors)
06745 {
06746
06747
06748
06749
06750
06751
06752 for( Int_t i=0; i < fNPeaks; i++){
06753 volumeErrors[i] = (Float_t) fVolumeErr[i];
06754 }
06755 }
06756
06757
06758 void TSpectrum2Fit::GetSigmaX(Double_t &sigmaX, Double_t &sigmaErrX)
06759 {
06760
06761
06762
06763
06764
06765
06766
06767 sigmaX=fSigmaCalcX;
06768 sigmaErrX=fSigmaErrX;
06769 }
06770
06771
06772 void TSpectrum2Fit::GetSigmaY(Double_t &sigmaY, Double_t &sigmaErrY)
06773 {
06774
06775
06776
06777
06778
06779
06780
06781 sigmaY=fSigmaCalcY;
06782 sigmaErrY=fSigmaErrY;
06783 }
06784
06785
06786 void TSpectrum2Fit::GetRo(Double_t &ro, Double_t &roErr)
06787 {
06788
06789
06790
06791
06792
06793
06794
06795 ro=fRoCalc;
06796 roErr=fRoErr;
06797 }
06798
06799
06800 void TSpectrum2Fit::GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &ax, Double_t &axErr, Double_t &ay, Double_t &ayErr)
06801 {
06802
06803
06804
06805
06806
06807
06808
06809
06810
06811
06812
06813
06814 a0 = fA0Calc;
06815 a0Err = fA0Err;
06816 ax = fAxCalc;
06817 axErr = fAxErr;
06818 ay = fAyCalc;
06819 ayErr = fAyErr;
06820 }
06821
06822
06823 void TSpectrum2Fit::GetTailParameters(Double_t &txy, Double_t &txyErr, Double_t &tx, Double_t &txErr, Double_t &ty, Double_t &tyErr, Double_t &bx, Double_t &bxErr, Double_t &by, Double_t &byErr, Double_t &sxy, Double_t &sxyErr, Double_t &sx, Double_t &sxErr, Double_t &sy, Double_t &syErr)
06824 {
06825
06826
06827
06828
06829
06830
06831
06832
06833
06834
06835
06836
06837
06838
06839
06840
06841
06842
06843
06844
06845
06846 txy = fTxyCalc;
06847 txyErr = fTxyErr;
06848 tx = fTxCalc;
06849 txErr = fTxErr;
06850 ty = fTyCalc;
06851 tyErr = fTyErr;
06852 bx = fBxCalc;
06853 bxErr = fBxErr;
06854 by = fByCalc;
06855 byErr = fByErr;
06856 sxy = fSxyCalc;
06857 sxyErr = fSxyErr;
06858 sx = fSxCalc;
06859 sxErr = fSxErr;
06860 sy = fSyCalc;
06861 syErr = fSyErr;
06862 }
06863