ROOT_528-00b_version: hist/spectrum/src/TSpectrumFit.cxx Source File

00001 // @(#)root/spectrum:$Id: TSpectrumFit.cxx 35405 2010-09-19 17:02:28Z brun $
00002 // Author: Miroslav Morhac   25/09/06
00003 
00004 //__________________________________________________________________________
00005 //   THIS CLASS CONTAINS ADVANCED SPECTRA FITTING FUNCTIONS.               //
00006 //                                                                         //
00007 //                                                                         //
00008 //   These functions were written by:                                      //
00009 //   Miroslav Morhac                                                       //
00010 //   Institute of Physics                                                  //
00011 //   Slovak Academy of Sciences                                            //
00012 //   Dubravska cesta 9, 842 28 BRATISLAVA                                  //
00013 //   SLOVAKIA                                                              //
00014 //                                                                         //
00015 //   email:fyzimiro@savba.sk,    fax:+421 7 54772479                       //
00016 //                                                                         //
00017 //  The original code in C has been repackaged as a C++ class by R.Brun    //
00018 //                                                                         //
00019 //  The algorithms in this class have been published in the following      //
00020 //  references:                                                            //
00021 //   [1] M. Morhac et al.: Efficient fitting algorithms applied to         //
00022 //   analysis of coincidence gamma-ray spectra. Computer Physics           //
00023 //   Communications, Vol 172/1 (2005) pp. 19-41.                           //
00024 //                                                                         //
00025 //   [2]  M. Morhac et al.: Study of fitting algorithms applied to         //
00026 //   simultaneous analysis of large number of peaks in gamma-ray spectra.  //
00027 //   Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003.              //
00028 //                                                                         //
00029 //                                                                         //
00030 //____________________________________________________________________________
00031     
00032 #include "TSpectrumFit.h"
00033 #include "TMath.h"
00034 
00035 ClassImp(TSpectrumFit)  
00036     
00037 //______________________________________________________________________________
00038 TSpectrumFit::TSpectrumFit() :TNamed("SpectrumFit", "Miroslav Morhac peak fitter") 
00039 {
00040    //default constructor
00041 
00042    fNPeaks = 0;
00043    fNumberIterations = 1;
00044    fXmin = 0;
00045    fXmax = 100; 
00046    fStatisticType = kFitOptimChiCounts;
00047    fAlphaOptim = kFitAlphaHalving;     
00048    fPower = kFitPower2;                
00049    fFitTaylor = kFitTaylorOrderFirst;  
00050    fAlpha =1; 
00051    fChi = 0;                    
00052    fPositionInit   = 0;
00053    fPositionCalc   = 0;
00054    fPositionErr   = 0;
00055    fFixPosition   = 0;
00056    fAmpInit   = 0;
00057    fAmpCalc   = 0;
00058    fAmpErr    = 0;
00059    fFixAmp    = 0;
00060    fArea      = 0;
00061    fAreaErr   = 0;
00062    fSigmaInit = 2;
00063    fSigmaCalc = 1;
00064    fSigmaErr  = 0;
00065    fTInit = 0; 
00066    fTCalc = 0;
00067    fTErr = 0; 
00068    fBInit = 1; 
00069    fBCalc = 0;
00070    fBErr = 0; 
00071    fSInit = 0; 
00072    fSCalc = 0;
00073    fSErr = 0; 
00074    fA0Init = 0; 
00075    fA0Calc = 0;
00076    fA0Err = 0;
00077    fA1Init = 0;
00078    fA1Calc = 0;
00079    fA1Err = 0;
00080    fA2Init = 0;
00081    fA2Calc = 0;
00082    fA2Err = 0;
00083    fFixSigma = false;
00084    fFixT = true;
00085    fFixB = true;
00086    fFixS = true;
00087    fFixA0 = true;
00088    fFixA1 = true;
00089    fFixA2 = true;
00090 }
00091 
00092 //______________________________________________________________________________
00093 TSpectrumFit::TSpectrumFit(Int_t numberPeaks) :TNamed("SpectrumFit", "Miroslav Morhac peak fitter") 
00094 {
00095    //numberPeaks: number of fitted peaks (must be greater than zero)
00096    //the constructor allocates arrays for all fitted parameters (peak positions, amplitudes etc) and sets the member
00097    //variables to their default values. One can change these variables by member functions (setters) of TSpectrumFit class.
00098 //Begin_Html <!--
00099 /* -->
00100 <div class=Section1>
00101 
00102 <p class=MsoNormal style='text-align:justify'>Shape function of the fitted
00103 peaks is </p>
00104 
00105 <p class=MsoNormal style='text-align:justify'>
00106 
00107 <table cellpadding=0 cellspacing=0 align=left>
00108  <tr>
00109   <td width=68 height=6></td>
00110  </tr>
00111  <tr>
00112   <td></td>
00113   <td><img width=388 height=132 src="gif/spectrumfit_constructor_image001.gif"></td>
00114  </tr>
00115 </table>
00116 
00117 <span style='font-family:Arial'>&nbsp;</span></p>
00118 
00119 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
00120 
00121 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
00122 
00123 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
00124 
00125 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
00126 
00127 <p class=MsoNormal><i>&nbsp;</i></p>
00128 
00129 <p class=MsoNormal><i>&nbsp;</i></p>
00130 
00131 <p class=MsoNormal><i>&nbsp;</i></p>
00132 
00133 <br clear=ALL>
00134 
00135 <p class=MsoNormal style='text-align:justify'>where a represents vector of
00136 fitted parameters (positions p(j), amplitudes A(j), sigma, relative amplitudes
00137 T, S and slope B).</p>
00138 
00139 <p class=MsoNormal><span style='font-size:16.0pt'>&nbsp;</span></p>
00140 
00141 </div>
00142    
00143 <!-- */
00144 // --> End_Html
00145     
00146    if (numberPeaks <= 0){
00147       Error ("TSpectrumFit","Invalid number of peaks, must be > than 0");
00148       return;
00149    }
00150    fNPeaks = numberPeaks;  
00151    fNumberIterations = 1;
00152    fXmin = 0;
00153    fXmax = 100; 
00154    fStatisticType = kFitOptimChiCounts;
00155    fAlphaOptim = kFitAlphaHalving;     
00156    fPower = kFitPower2;                
00157    fFitTaylor = kFitTaylorOrderFirst;  
00158    fAlpha =1; 
00159    fChi = 0;                    
00160    fPositionInit   = new Double_t[numberPeaks];
00161    fPositionCalc   = new Double_t[numberPeaks];   
00162    fPositionErr   = new Double_t[numberPeaks];   
00163    fFixPosition   = new Bool_t[numberPeaks];   
00164    fAmpInit   = new Double_t[numberPeaks];   
00165    fAmpCalc   = new Double_t[numberPeaks];   
00166    fAmpErr    = new Double_t[numberPeaks];   
00167    fFixAmp    = new Bool_t[numberPeaks];     
00168    fArea      = new Double_t[numberPeaks];   
00169    fAreaErr   = new Double_t[numberPeaks];      
00170    fSigmaInit = 2;
00171    fSigmaCalc = 1;
00172    fSigmaErr  = 0;
00173    fTInit = 0; 
00174    fTCalc = 0;
00175    fTErr = 0; 
00176    fBInit = 1; 
00177    fBCalc = 0;
00178    fBErr = 0; 
00179    fSInit = 0; 
00180    fSCalc = 0;
00181    fSErr = 0; 
00182    fA0Init = 0; 
00183    fA0Calc = 0;
00184    fA0Err = 0;
00185    fA1Init = 0;
00186    fA1Calc = 0;
00187    fA1Err = 0;
00188    fA2Init = 0;
00189    fA2Calc = 0;
00190    fA2Err = 0;
00191    fFixSigma = false;
00192    fFixT = true;
00193    fFixB = true;
00194    fFixS = true;
00195    fFixA0 = true;
00196    fFixA1 = true;
00197    fFixA2 = true;
00198 }
00199     
00200     
00201 
00202 //______________________________________________________________________________
00203 TSpectrumFit::~TSpectrumFit() 
00204 {
00205    //destructor
00206    delete [] fPositionInit;   
00207    delete [] fPositionCalc;   
00208    delete [] fPositionErr;
00209    delete [] fFixPosition;      
00210    delete [] fAmpInit;  
00211    delete [] fAmpCalc;  
00212    delete [] fAmpErr;   
00213    delete [] fFixAmp;
00214    delete [] fArea;  
00215    delete [] fAreaErr;         
00216 }
00217 
00218 //_____________________________________________________________________________
00219 /////////////////BEGINNING OF AUXILIARY FUNCTIONS USED BY FITTING FUNCTION Fit1//////////////////////////
00220 Double_t TSpectrumFit::Erfc(Double_t x) 
00221 {   
00222 //////////////////////////////////////////////////////////////////////////////
00223 //   AUXILIARY FUNCTION                                                      //
00224 //                                                                          //
00225 //   This function calculates error function of x.                           //
00226 //                                                                          //
00227 //////////////////////////////////////////////////////////////////////////////
00228    Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap = 0.47047;
00229    Double_t a, t, c, w;
00230    a = TMath::Abs(x);
00231    w = 1. + dap * a;
00232    t = 1. / w;
00233    w = a * a;
00234    if (w < 700)
00235       c = exp(-w);
00236    
00237    else {
00238       c = 0;
00239    }
00240    c = c * t * (da1 + t * (da2 + t * da3));
00241    if (x < 0)
00242       c = 1. - c;
00243    return (c);
00244 }
00245 
00246 //____________________________________________________________________________
00247 Double_t TSpectrumFit::Derfc(Double_t x) 
00248 {
00249 //////////////////////////////////////////////////////////////////////////////
00250 //   AUXILIARY FUNCTION                                                      //
00251 //                                                                          //
00252 //   This function calculates derivative of error function of x.             //
00253 //                                                                          //
00254 //////////////////////////////////////////////////////////////////////////////
00255    Double_t a, t, c, w;
00256    Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap = 0.47047;
00257    a = TMath::Abs(x);
00258    w = 1. + dap * a;
00259    t = 1. / w;
00260    w = a * a;
00261    if (w < 700)
00262       c = exp(-w);
00263    
00264    else {
00265       c = 0;
00266    }
00267    c = (-1.) * dap * c * t * t * (da1 + t * (2. * da2 + t * 3. * da3)) -
00268        2. * a * Erfc(a);
00269    return (c);
00270 }
00271 
00272 //____________________________________________________________________________
00273 Double_t TSpectrumFit::Deramp(Double_t i, Double_t i0, Double_t sigma, Double_t t,
00274                            Double_t s, Double_t b) 
00275 {   
00276 //////////////////////////////////////////////////////////////////////////////
00277 //   AUXILIARY FUNCTION                                                      //
00278 //                                                                          //
00279 //   This function calculates derivative of peak shape function (see manual) //
00280 //   according to amplitude of peak.                                        //
00281 //      Function parameters:                                                //
00282 //              -i-channel                                                  //
00283 //              -i0-position of peak                                        //
00284 //              -sigma-sigma of peak                                        //
00285 //              -t, s-relative amplitudes                                   //
00286 //              -b-slope                                                    //
00287 //                                                                          //
00288 //////////////////////////////////////////////////////////////////////////////
00289    Double_t p, q, r, a;
00290    p = (i - i0) / sigma;
00291    if ((p * p) < 700)
00292       q = exp(-p * p);
00293    
00294    else {
00295       q = 0;
00296    }
00297    r = 0;
00298    if (t != 0) {
00299       a = p / b;
00300       if (a > 700)
00301          a = 700;
00302       r = t * exp(a) / 2.;
00303    }
00304    if (r != 0)
00305       r = r * Erfc(p + 1. / (2. * b));
00306    q = q + r;
00307    if (s != 0)
00308       q = q + s * Erfc(p) / 2.;
00309    return (q);
00310 }
00311 
00312 //____________________________________________________________________________
00313 Double_t TSpectrumFit::Deri0(Double_t i, Double_t amp, Double_t i0, Double_t sigma,
00314                           Double_t t, Double_t s, Double_t b) 
00315 {  
00316 //////////////////////////////////////////////////////////////////////////////
00317 //   AUXILIARY FUNCTION                                                      //
00318 //                                                                          //
00319 //   This function calculates derivative of peak shape function (see manual) //
00320 //   according to peak position.                                            //
00321 //      Function parameters:                                                //
00322 //              -i-channel                                                  //
00323 //              -amp-amplitude of peak                                      //
00324 //              -i0-position of peak                                        //
00325 //              -sigma-sigma of peak                                        //
00326 //              -t, s-relative amplitudes                                   //
00327 //              -b-slope                                                    //
00328 //                                                                          //
00329 //////////////////////////////////////////////////////////////////////////////
00330    Double_t p, r1, r2, r3, r4, c, d, e;
00331    p = (i - i0) / sigma;
00332    d = 2. * sigma;
00333    if ((p * p) < 700)
00334       r1 = 2. * p * exp(-p * p) / sigma;
00335    
00336    else {
00337       r1 = 0;
00338    }
00339    r2 = 0, r3 = 0;
00340    if (t != 0) {
00341       c = p + 1. / (2. * b);
00342       e = p / b;
00343       if (e > 700)
00344          e = 700;
00345       r2 = -t * exp(e) * Erfc(c) / (d * b);
00346       r3 = -t * exp(e) * Derfc(c) / d;
00347    }
00348    r4 = 0;
00349    if (s != 0)
00350       r4 = -s * Derfc(p) / d;
00351    r1 = amp * (r1 + r2 + r3 + r4);
00352    return (r1);
00353 }
00354 
00355 //____________________________________________________________________________
00356 Double_t TSpectrumFit::Derderi0(Double_t i, Double_t amp, Double_t i0,
00357                              Double_t sigma) 
00358 {  
00359 //////////////////////////////////////////////////////////////////////////////
00360 //   AUXILIARY FUNCTION                                                      //
00361 //                                                                          //
00362 //   This function calculates second derivative of peak shape function       //
00363 //   (see manual) according to peak position.                               //
00364 //      Function parameters:                                                //
00365 //              -i-channel                                                  //
00366 //              -amp-amplitude of peak                                      //
00367 //              -i0-position of peak                                        //
00368 //              -sigma-width of peak                                        //
00369 //                                                                          //
00370 //////////////////////////////////////////////////////////////////////////////
00371    Double_t p, r1, r2, r3, r4;
00372    p = (i - i0) / sigma;
00373    if ((p * p) < 700)
00374       r1 = exp(-p * p);
00375    
00376    else {
00377       r1 = 0;
00378    }
00379    r1 = r1 * (4 * p * p - 2) / (sigma * sigma);
00380    r2 = 0, r3 = 0, r4 = 0;
00381    r1 = amp * (r1 + r2 + r3 + r4);
00382    return (r1);
00383 }
00384 
00385 //____________________________________________________________________________
00386 Double_t TSpectrumFit::Dersigma(Int_t num_of_fitted_peaks, Double_t i,
00387                              const Double_t *parameter, Double_t sigma,
00388                              Double_t t, Double_t s, Double_t b) 
00389 {  
00390 //////////////////////////////////////////////////////////////////////////////////
00391 //   AUXILIARY FUNCTION                                                          //
00392 //                                                                              //
00393 //   This function calculates derivative of peaks shape function (see manual)    //
00394 //   according to sigma of peaks.                                               //
00395 //      Function parameters:                                                    //
00396 //              -num_of_fitted_peaks-number of fitted peaks                     //
00397 //              -i-channel                                                      //
00398 //              -parameter-array of peaks parameters (amplitudes and positions) //
00399 //              -sigma-sigma of peak                                            //
00400 //              -t, s-relative amplitudes                                       //
00401 //              -b-slope                                                        //
00402 //                                                                              //
00403 //////////////////////////////////////////////////////////////////////////////////
00404    Int_t j;
00405    Double_t r, p, r1, r2, r3, r4, c, d, e;
00406    r = 0;
00407    d = 2. * sigma;
00408    for (j = 0; j < num_of_fitted_peaks; j++) {
00409       p = (i - parameter[2 * j + 1]) / sigma;
00410       r1 = 0;
00411       if (TMath::Abs(p) < 3) {
00412          if ((p * p) < 700)
00413             r1 = 2. * p * p * exp(-p * p) / sigma;
00414          
00415          else {
00416             r1 = 0;
00417          }
00418       }
00419       r2 = 0, r3 = 0;
00420       if (t != 0) {
00421          c = p + 1. / (2. * b);
00422          e = p / b;
00423          if (e > 700)
00424             e = 700;
00425          r2 = -t * p * exp(e) * Erfc(c) / (d * b);
00426          r3 = -t * p * exp(e) * Derfc(c) / d;
00427       }
00428       r4 = 0;
00429       if (s != 0)
00430          r4 = -s * p * Derfc(p) / d;
00431       r = r + parameter[2 * j] * (r1 + r2 + r3 + r4);
00432    }
00433    return (r);
00434 }
00435 
00436 //____________________________________________________________________________
00437 Double_t TSpectrumFit::Derdersigma(Int_t num_of_fitted_peaks, Double_t i,
00438                                const Double_t *parameter, Double_t sigma) 
00439 {  
00440 //////////////////////////////////////////////////////////////////////////////////
00441 //   AUXILIARY FUNCTION                                                          //
00442 //                                                                              //
00443 //   This function calculates second derivative of peaks shape function          //
00444 //   (see manual) according to sigma of peaks.                                  //
00445 //      Function parameters:                                                    //
00446 //              -num_of_fitted_peaks-number of fitted peaks                     //
00447 //              -i-channel                                                      //
00448 //              -parameter-array of peaks parameters (amplitudes and positions) //
00449 //              -sigma-sigma of peak                                            //
00450 //                                                                              //
00451 //////////////////////////////////////////////////////////////////////////////////
00452    Int_t j;
00453    Double_t r, p, r1, r2, r3, r4;
00454    r = 0;
00455    for (j = 0; j < num_of_fitted_peaks; j++) {
00456       p = (i - parameter[2 * j + 1]) / sigma;
00457       r1 = 0;
00458       if (TMath::Abs(p) < 3) {
00459          if ((p * p) < 700)
00460             r1 = exp(-p * p) * p * p * (4. * p * p - 6) / (sigma * sigma);
00461          
00462          else {
00463             r1 = 0;
00464          }
00465       }
00466       r2 = 0, r3 = 0, r4 = 0;
00467       r = r + parameter[2 * j] * (r1 + r2 + r3 + r4);
00468    }
00469    return (r);
00470 }
00471 
00472 //____________________________________________________________________________
00473 Double_t TSpectrumFit::Dert(Int_t num_of_fitted_peaks, Double_t i,
00474                         const Double_t *parameter, Double_t sigma, Double_t b) 
00475 {  
00476 //////////////////////////////////////////////////////////////////////////////////
00477 //   AUXILIARY FUNCTION                                                          //
00478 //                                                                              //
00479 //   This function calculates derivative of peaks shape function (see manual)    //
00480 //   according to relative amplitude t.                                         //
00481 //      Function parameters:                                                    //
00482 //              -num_of_fitted_peaks-number of fitted peaks                     //
00483 //              -i-channel                                                      //
00484 //              -parameter-array of peaks parameters (amplitudes and positions) //
00485 //              -sigma-sigma of peak                                            //
00486 //              -b-slope                                                        //
00487 //                                                                              //
00488 //////////////////////////////////////////////////////////////////////////////////
00489    Int_t j;
00490    Double_t r, p, r1, c, e;
00491    r = 0;
00492    for (j = 0; j < num_of_fitted_peaks; j++) {
00493       p = (i - parameter[2 * j + 1]) / sigma;
00494       c = p + 1. / (2. * b);
00495       e = p / b;
00496       if (e > 700)
00497          e = 700;
00498       r1 = exp(e) * Erfc(c);
00499       r = r + parameter[2 * j] * r1;
00500    }
00501    r = r / 2.;
00502    return (r);
00503 }
00504 
00505 //____________________________________________________________________________
00506 Double_t TSpectrumFit::Ders(Int_t num_of_fitted_peaks, Double_t i,
00507                         const Double_t *parameter, Double_t sigma) 
00508 {  
00509 //////////////////////////////////////////////////////////////////////////////////
00510 //   AUXILIARY FUNCTION                                                          //
00511 //                                                                              //
00512 //   This function calculates derivative of peaks shape function (see manual)    //
00513 //   according to relative amplitude s.                                               //
00514 //      Function parameters:                                                    //
00515 //              -num_of_fitted_peaks-number of fitted peaks                     //
00516 //              -i-channel                                                      //
00517 //              -parameter-array of peaks parameters (amplitudes and positions) //
00518 //              -sigma-sigma of peak                                            //
00519 //                                                                              //
00520 //////////////////////////////////////////////////////////////////////////////////
00521    Int_t j;
00522    Double_t r, p, r1;
00523    r = 0;
00524    for (j = 0; j < num_of_fitted_peaks; j++) {
00525       p = (i - parameter[2 * j + 1]) / sigma;
00526       r1 = Erfc(p);
00527       r = r + parameter[2 * j] * r1;
00528    }
00529    r = r / 2.;
00530    return (r);
00531 }
00532 
00533 //____________________________________________________________________________
00534 Double_t TSpectrumFit::Derb(Int_t num_of_fitted_peaks, Double_t i,
00535                         const Double_t *parameter, Double_t sigma, Double_t t,
00536                         Double_t b) 
00537 {  
00538 //////////////////////////////////////////////////////////////////////////////////
00539 //   AUXILIARY FUNCTION                                                          //
00540 //                                                                              //
00541 //   This function calculates derivative of peaks shape function (see manual)    //
00542 //   according to slope b.                                                      //
00543 //      Function parameters:                                                    //
00544 //              -num_of_fitted_peaks-number of fitted peaks                     //
00545 //              -i-channel                                                      //
00546 //              -parameter-array of peaks parameters (amplitudes and positions) //
00547 //              -sigma-sigma of peak                                            //
00548 //              -t-relative amplitude                                           //
00549 //              -b-slope                                                        //
00550 //                                                                              //
00551 //////////////////////////////////////////////////////////////////////////////////
00552    Int_t j;
00553    Double_t r, p, r1, c, e;
00554    r = 0;
00555    for (j = 0; j < num_of_fitted_peaks && t != 0; j++) {
00556       p = (i - parameter[2 * j + 1]) / sigma;
00557       c = p + 1. / (2. * b);
00558       e = p / b;
00559       r1 = p * Erfc(c);
00560       r1 = r1 + Derfc(c) / 2.;
00561       if (e > 700)
00562          e = 700;
00563       if (e < -700)
00564          r1 = 0;
00565       
00566       else
00567          r1 = r1 * exp(e);
00568       r = r + parameter[2 * j] * r1;
00569    }
00570    r = -r * t / (2. * b * b);
00571    return (r);
00572 }
00573 
00574 //____________________________________________________________________________
00575 Double_t TSpectrumFit::Dera1(Double_t i)       
00576 {
00577    //derivative of background according to a1
00578    return (i);
00579 }
00580 
00581 //____________________________________________________________________________
00582 Double_t TSpectrumFit::Dera2(Double_t i)      
00583 {
00584    //derivative of background according to a2
00585    return (i * i);
00586 }
00587 
00588 //____________________________________________________________________________
00589 Double_t TSpectrumFit::Shape(Int_t num_of_fitted_peaks, Double_t i,
00590                          const Double_t *parameter, Double_t sigma, Double_t t,
00591                          Double_t s, Double_t b, Double_t a0, Double_t a1,
00592                          Double_t a2) 
00593 {  
00594 //////////////////////////////////////////////////////////////////////////////////
00595 //   AUXILIARY FUNCTION                                                          //
00596 //                                                                              //
00597 //   This function calculates peaks shape function (see manual)                  //
00598 //      Function parameters:                                                    //
00599 //              -num_of_fitted_peaks-number of fitted peaks                     //
00600 //              -i-channel                                                      //
00601 //              -parameter-array of peaks parameters (amplitudes and positions) //
00602 //              -sigma-sigma of peak                                            //
00603 //              -t, s-relative amplitudes                                       //
00604 //              -b-slope                                                        //
00605 //              -a0, a1, a2- background coefficients                            //
00606 //                                                                              //
00607 //////////////////////////////////////////////////////////////////////////////////
00608    Int_t j;
00609    Double_t r, p, r1, r2, r3, c, e;
00610    r = 0;
00611    for (j = 0; j < num_of_fitted_peaks; j++) {
00612       if (sigma > 0.0001)
00613          p = (i - parameter[2 * j + 1]) / sigma;
00614       
00615       else {
00616          if (i == parameter[2 * j + 1])
00617             p = 0;
00618          
00619          else
00620             p = 10;
00621       }
00622       r1 = 0;
00623       if (TMath::Abs(p) < 3) {
00624          if ((p * p) < 700)
00625             r1 = exp(-p * p);
00626          
00627          else {
00628             r1 = 0;
00629          }
00630       }
00631       r2 = 0;
00632       if (t != 0) {
00633          c = p + 1. / (2. * b);
00634          e = p / b;
00635          if (e > 700)
00636             e = 700;
00637          r2 = t * exp(e) * Erfc(c) / 2.;
00638       }
00639       r3 = 0;
00640       if (s != 0)
00641          r3 = s * Erfc(p) / 2.;
00642       r = r + parameter[2 * j] * (r1 + r2 + r3);
00643    }
00644    r = r + a0 + a1 * i + a2 * i * i;
00645    return (r);
00646 }
00647 
00648 //____________________________________________________________________________
00649 Double_t TSpectrumFit::Area(Double_t a, Double_t sigma, Double_t t, Double_t b) 
00650 {  
00651 //////////////////////////////////////////////////////////////////////////////////
00652 //   AUXILIARY FUNCTION                                                          //
00653 //                                                                              //
00654 //   This function calculates area of a peak                                     //
00655 //      Function parameters:                                                    //
00656 //              -a-amplitude of the peak                                        //
00657 //              -sigma-sigma of peak                                            //
00658 //              -t-relative amplitude                                           //
00659 //              -b-slope                                                        //
00660 //                                                                              //
00661 //////////////////////////////////////////////////////////////////////////////////
00662    Double_t odm_pi = 1.7724538, r = 0;
00663    if (b != 0)
00664       r = 0.5 / b;
00665    r = (-1.) * r * r;
00666    if (TMath::Abs(r) < 700)
00667       r = a * sigma * (odm_pi + t * b * exp(r));
00668    
00669    else {
00670       r = a * sigma * odm_pi;
00671    }
00672    return (r);
00673 }
00674 
00675 //____________________________________________________________________________
00676 Double_t TSpectrumFit::Derpa(Double_t sigma, Double_t t, Double_t b) 
00677 {  
00678 //////////////////////////////////////////////////////////////////////////////////
00679 //   AUXILIARY FUNCTION                                                          //
00680 //                                                                              //
00681 //   This function calculates derivative of the area of peak                     //
00682 //   according to its amplitude.                                                //
00683 //      Function parameters:                                                    //
00684 //              -sigma-sigma of peak                                            //
00685 //              -t-relative amplitudes                                          //
00686 //              -b-slope                                                        //
00687 //                                                                              //
00688 //////////////////////////////////////////////////////////////////////////////////
00689    Double_t odm_pi = 1.7724538, r;
00690    r = 0.5 / b;
00691    r = (-1.) * r * r;
00692    if (TMath::Abs(r) < 700)
00693       r = sigma * (odm_pi + t * b * exp(r));
00694    
00695    else {
00696       r = sigma * odm_pi;
00697    }
00698    return (r);
00699 }
00700 Double_t TSpectrumFit::Derpsigma(Double_t a, Double_t t, Double_t b) 
00701 {  
00702 //////////////////////////////////////////////////////////////////////////////////
00703 //   AUXILIARY FUNCTION                                                          //
00704 //                                                                              //
00705 //   This function calculates derivative of the area of peak                     //
00706 //   according to sigma of peaks.                                               //
00707 //      Function parameters:                                                    //
00708 //              -a-amplitude of peak                                            //
00709 //              -t-relative amplitudes                                          //
00710 //              -b-slope                                                        //
00711 //                                                                              //
00712 //////////////////////////////////////////////////////////////////////////////////
00713    Double_t odm_pi = 1.7724538, r;
00714    r = 0.5 / b;
00715    r = (-1.) * r * r;
00716    if (TMath::Abs(r) < 700)
00717       r = a * (odm_pi + t * b * exp(r));
00718    
00719    else {
00720       r = a * odm_pi;
00721    }
00722    return (r);
00723 }
00724 
00725 //______________________________________________________________________________
00726 Double_t TSpectrumFit::Derpt(Double_t a, Double_t sigma, Double_t b) 
00727 {  
00728 //////////////////////////////////////////////////////////////////////////////////
00729 //   AUXILIARY FUNCTION                                                          //
00730 //                                                                              //
00731 //   This function calculates derivative of the area of peak                     //
00732 //   according to t parameter.                                                  //
00733 //      Function parameters:                                                    //
00734 //              -sigma-sigma of peak                                            //
00735 //              -t-relative amplitudes                                          //
00736 //              -b-slope                                                        //
00737 //                                                                              //
00738 //////////////////////////////////////////////////////////////////////////////////
00739    Double_t r;
00740    r = 0.5 / b;
00741    r = (-1.) * r * r;
00742    if (TMath::Abs(r) < 700)
00743       r = a * sigma * b * exp(r);
00744    
00745    else {
00746       r = 0;
00747    }
00748    return (r);
00749 }
00750 
00751 //______________________________________________________________________________
00752 Double_t TSpectrumFit::Derpb(Double_t a, Double_t sigma, Double_t t, Double_t b) 
00753 {  
00754 //////////////////////////////////////////////////////////////////////////////////
00755 //   AUXILIARY FUNCTION                                                          //
00756 //                                                                              //
00757 //   This function calculates derivative of the area of peak                     //
00758 //   according to b parameter.                                                  //
00759 //      Function parameters:                                                    //
00760 //              -sigma-sigma of peak                                            //
00761 //              -t-relative amplitudes                                          //
00762 //              -b-slope                                                        //
00763 //                                                                              //
00764 //////////////////////////////////////////////////////////////////////////////////
00765    Double_t r;
00766    r = (-1) * 0.25 / (b * b);
00767    if (TMath::Abs(r) < 700)
00768       r = a * sigma * t * exp(r) * (1 - 2 * r);
00769    
00770    else {
00771       r = 0;
00772    }
00773    return (r);
00774 }
00775 
00776 //______________________________________________________________________________
00777 Double_t TSpectrumFit::Ourpowl(Double_t a, Int_t pw)
00778 {                               
00779    //power function
00780    Double_t c;
00781    Double_t a2 = a*a;
00782    c = 1;
00783    if (pw >  0) c *= a2;
00784    if (pw >  2) c *= a2;
00785    if (pw >  4) c *= a2;
00786    if (pw >  6) c *= a2;
00787    if (pw >  8) c *= a2;
00788    if (pw > 10) c *= a2;
00789    if (pw > 12) c *= a2;
00790    return (c);
00791 }
00792 
00793 /////////////////END OF AUXILIARY FUNCTIONS USED BY FITTING FUNCTIONS FitAWMI, FitStiefel//////////////////////////
00794 /////////////////FITTING FUNCTION WITHOUT MATRIX INVERSION///////////////////////////////////////
00795 
00796 //____________________________________________________________________________
00797 void TSpectrumFit::FitAwmi(Float_t *source) 
00798 {
00799 /////////////////////////////////////////////////////////////////////////////
00800 //        ONE-DIMENSIONAL FIT FUNCTION                                 
00801 //        ALGORITHM WITHOUT MATRIX INVERSION                                  
00802 //        This function fits the source spectrum. The calling program should 
00803 //        fill in input parameters of the TSpectrumFit class            
00804 //        The fitted parameters are written into           
00805 //        TSpectrumFit class output parameters and fitted data are written into 
00806 //        source spectrum.                                                    
00807 //                                                                          
00808 //        Function parameters:                                             
00809 //        source-pointer to the vector of source spectrum                  
00810 //                                                                         
00811 /////////////////////////////////////////////////////////////////////////////
00812 //
00813 //Begin_Html <!--
00814 /* -->
00815 <div class=Section2>
00816 
00817 <p class=MsoNormal><b><span style='font-size:14.0pt'>Fitting</span></b></p>
00818 
00819 <p class=MsoNormal style='text-align:justify'><i>Goal: to estimate
00820 simultaneously peak shape parameters in spectra with large number of peaks</i></p>
00821 
00822 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
00823 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
00824 </span>peaks can be fitted separately, each peak (or multiplets) in a region or
00825 together all peaks in a spectrum. To fit separately each peak one needs to
00826 determine the fitted region. However it can happen that the regions of
00827 neighboring peaks are overlapping. Then the results of fitting are very poor.
00828 On the other hand, when fitting together all peaks found in a  spectrum, one
00829 needs to have a method that is  stable (converges) and fast enough to carry out
00830 fitting in reasonable time </p>
00831 
00832 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
00833 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
00834 </span>we have implemented the nonsymmetrical semiempirical peak shape function
00835 [1]</p>
00836 
00837 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
00838 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
00839 </span>it contains the symmetrical Gaussian as well as nonsymmetrical terms.</p>
00840 
00841 <p class=MsoNormal style='text-align:justify'>
00842 
00843 <table cellpadding=0 cellspacing=0 align=left>
00844  <tr>
00845   <td width=84 height=18></td>
00846  </tr>
00847  <tr>
00848   <td></td>
00849   <td><img width=372 height=127 src="gif/spectrumfit_awni_image001.gif"></td>
00850  </tr>
00851 </table>
00852 
00853 <span style='font-size:16.0pt'>&nbsp;</span></p>
00854 
00855 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>&nbsp;</span></p>
00856 
00857 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>&nbsp;</span></p>
00858 
00859 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>&nbsp;</span></p>
00860 
00861 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>&nbsp;</span></p>
00862 
00863 <br clear=ALL>
00864 
00865 <p class=MsoNormal style='text-indent:34.2pt'>where T and S are relative amplitudes
00866 and B is slope.</p>
00867 
00868 <p class=MsoNormal>&nbsp;</p>
00869 
00870 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent:
00871 -18.0pt'>•<span style='font:7.0pt "Times New Roman"'>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
00872 </span>algorithm without matrix inversion (AWMI) allows fitting tens, hundreds
00873 of peaks simultaneously that represent sometimes thousands of parameters [2],
00874 [5]. </p>
00875 
00876 <p class=MsoNormal><i>Function:</i></p>
00877 
00878 <p class=MsoNormal style='text-align:justify'>void <a
00879 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrumFit::FitAwmi</b></a>(<a
00880 href="http://root.cern.ch/root/html/ListOfTypes.html#float"><b>float</b></a> *fSource)
00881 </p>
00882 
00883 <p class=MsoNormal style='text-align:justify'>This function fits the source
00884 spectrum using AWMI algorithm. The calling program should fill in input fitting
00885 parameters of the TSpectrumFit class using a set of TSpectrumFit setters. The
00886 fitted parameters are written into the class and the fitted data are written
00887 into source spectrum. </p>
00888 
00889 <p class=MsoNormal>&nbsp;</p>
00890 
00891 <p class=MsoNormal><i><span style='color:red'>Parameter:</span></i></p>
00892 
00893 <p class=MsoNormal style='text-align:justify'>        <b>fSource</b>-pointer to
00894 the vector of source spectrum                  </p>
00895 
00896 <p class=MsoNormal style='text-align:justify'>        </p>
00897 
00898 <p class=MsoNormal><i><span style='color:red'>Member variables of the
00899 TSpectrumFit class:</span></i></p>
00900 
00901 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00902 Int_t     fNPeaks;                    //number of peaks present in fit, input
00903 parameter, it should be &gt; 0</span></p>
00904 
00905 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00906 Int_t     fNumberIterations;          //number of iterations in fitting
00907 procedure, input parameter, it should be &gt; 0</span></p>
00908 
00909 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00910 Int_t     fXmin;                      //first fitted channel</span></p>
00911 
00912 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00913 Int_t     fXmax;                      //last fitted channel</span></p>
00914 
00915 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00916 Int_t     fStatisticType;             //type of statistics, possible values
00917 kFitOptimChiCounts (chi square statistics with counts as weighting
00918 coefficients), kFitOptimChiFuncValues (chi square statistics with function
00919 values as weighting coefficients),kFitOptimMaxLikelihood</span></p>
00920 
00921 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00922 Int_t     fAlphaOptim;                //optimization of convergence algorithm, possible
00923 values kFitAlphaHalving, kFitAlphaOptimal</span></p>
00924 
00925 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00926 Int_t     fPower;                     //possible values kFitPower2,4,6,8,10,12,
00927 for details see references. It applies only for Awmi fitting function.</span></p>
00928 
00929 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00930 Int_t     fFitTaylor;                 //order of Taylor expansion, possible
00931 values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi
00932 fitting function.</span></p>
00933 
00934 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00935 Double_t  fAlpha;                     //convergence coefficient, input
00936 parameter, it should be positive number and &lt;=1, for details see references</span></p>
00937 
00938 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00939 Double_t  fChi;                       //here the fitting functions return
00940 resulting chi square   </span></p>
00941 
00942 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00943 Double_t *fPositionInit;              //[fNPeaks] array of initial values of
00944 peaks positions, input parameters</span></p>
00945 
00946 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00947 Double_t *fPositionCalc;              //[fNPeaks] array of calculated values of
00948 fitted positions, output parameters</span></p>
00949 
00950 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00951 Double_t *fPositionErr;               //[fNPeaks] array of position errors</span></p>
00952 
00953 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00954 Double_t *fAmpInit;                   //[fNPeaks] array of initial values of
00955 peaks amplitudes, input parameters</span></p>
00956 
00957 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00958 Double_t *fAmpCalc;                   //[fNPeaks] array of calculated values of
00959 fitted amplitudes, output parameters</span></p>
00960 
00961 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00962 Double_t *fAmpErr;                    //[fNPeaks] array of amplitude errors</span></p>
00963 
00964 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00965 Double_t *fArea;                      //[fNPeaks] array of calculated areas of
00966 peaks</span></p>
00967 
00968 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00969 Double_t *fAreaErr;                   //[fNPeaks] array of errors of peak areas</span></p>
00970 
00971 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00972 Double_t  fSigmaInit;                 //initial value of sigma parameter</span></p>
00973 
00974 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00975 Double_t  fSigmaCalc;                 //calculated value of sigma parameter</span></p>
00976 
00977 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00978 Double_t  fSigmaErr;                  //error value of sigma parameter</span></p>
00979 
00980 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00981 Double_t  fTInit;                     //initial value of t parameter (relative
00982 amplitude of tail), for details see html manual and references</span></p>
00983 
00984 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00985 Double_t  fTCalc;                     //calculated value of t parameter</span></p>
00986 
00987 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00988 Double_t  fTErr;                      //error value of t parameter</span></p>
00989 
00990 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00991 Double_t  fBInit;                     //initial value of b parameter (slope),
00992 for details see html manual and references</span></p>
00993 
00994 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00995 Double_t  fBCalc;                     //calculated value of b parameter</span></p>
00996 
00997 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
00998 Double_t  fBErr;                      //error value of b parameter</span></p>
00999 
01000 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01001 Double_t  fSInit;                     //initial value of s parameter (relative
01002 amplitude of step), for details see html manual and references</span></p>
01003 
01004 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01005 Double_t  fSCalc;                     //calculated value of s parameter</span></p>
01006 
01007 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01008 Double_t  fSErr;                      //error value of s parameter</span></p>
01009 
01010 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01011 Double_t  fA0Init;                    //initial value of background a0
01012 parameter(backgroud is estimated as a0+a1*x+a2*x*x)</span></p>
01013 
01014 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01015 Double_t  fA0Calc;                    //calculated value of background a0
01016 parameter</span></p>
01017 
01018 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01019 Double_t  fA0Err;                     //error value of background a0 parameter</span></p>
01020 
01021 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01022 Double_t  fA1Init;                    //initial value of background a1
01023 parameter(backgroud is estimated as a0+a1*x+a2*x*x)</span></p>
01024 
01025 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01026 Double_t  fA1Calc;                    //calculated value of background a1
01027 parameter</span></p>
01028 
01029 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01030 Double_t  fA1Err;                     //error value of background a1 parameter</span></p>
01031 
01032 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01033 Double_t  fA2Init;                    //initial value of background a2
01034 parameter(backgroud is estimated as a0+a1*x+a2*x*x)</span></p>
01035 
01036 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01037 Double_t  fA2Calc;                    //calculated value of background a2
01038 parameter</span></p>
01039 
01040 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01041 Double_t  fA2Err;                     //error value of background a2 parameter</span></p>
01042 
01043 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01044 Bool_t   *fFixPosition;               //[fNPeaks] array of logical values which
01045 allow to fix appropriate positions (not fit). However they are present in the
01046 estimated functional   </span></p>
01047 
01048 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01049 Bool_t   *fFixAmp;                    //[fNPeaks] array of logical values which
01050 allow to fix appropriate amplitudes (not fit). However they are present in the
01051 estimated functional      </span></p>
01052 
01053 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01054 Bool_t    fFixSigma;                  //logical value of sigma parameter, which
01055 allows to fix the parameter (not to fit).   </span></p>
01056 
01057 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01058 Bool_t    fFixT;                      //logical value of t parameter, which
01059 allows to fix the parameter (not to fit).      </span></p>
01060 
01061 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01062 Bool_t    fFixB;                      //logical value of b parameter, which
01063 allows to fix the parameter (not to fit).   </span></p>
01064 
01065 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01066 Bool_t    fFixS;                      //logical value of s parameter, which
01067 allows to fix the parameter (not to fit).      </span></p>
01068 
01069 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01070 Bool_t    fFixA0;                     //logical value of a0 parameter, which
01071 allows to fix the parameter (not to fit).</span></p>
01072 
01073 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01074 Bool_t    fFixA1;                     //logical value of a1 parameter, which
01075 allows to fix the parameter (not to fit).   </span></p>
01076 
01077 <p class=MsoNormal style='text-align:justify'><span style='font-size:10.0pt'>  
01078 Bool_t    fFixA2;                     //logical value of a2 parameter, which
01079 allows to fix the parameter (not to fit).</span></p>
01080 
01081 <p class=MsoNormal style='text-align:justify'><b><i>&nbsp;</i></b></p>
01082 
01083 <p class=MsoNormal style='text-align:justify'><b><i>References:</i></b></p>
01084 
01085 <p class=MsoNormal style='text-align:justify'>[1] Phillps G.W., Marlow K.W.,
01086 NIM 137 (1976) 525.</p>
01087 
01088 <p class=MsoNormal style='text-align:justify'>[2] I. A. Slavic: Nonlinear
01089 least-squares fitting without matrix inversion applied to complex Gaussian
01090 spectra analysis. NIM 134 (1976) 285-289.</p>
01091 
01092 <p class=MsoNormal style='text-align:justify'>[3] T. Awaya: A new method for
01093 curve fitting to the data with low statistics not using chi-square method. NIM
01094 165 (1979) 317-323.</p>
01095 
01096 <p class=MsoNormal style='text-align:justify'>[4] T. Hauschild, M. Jentschel:
01097 Comparison of maximum likelihood estimation and chi-square statistics applied
01098 to counting experiments. NIM A 457 (2001) 384-401.</p>
01099 
01100 <p class=MsoNormal style='text-align:justify'> [5]  M. Morhá&#269;,  J.
01101 Kliman,  M. Jandel,  &#317;. Krupa, V. Matoušek: Study of fitting algorithms
01102 applied to simultaneous analysis of large number of peaks in -ray spectra. <span
01103 lang=EN-GB>Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003</span></p>
01104 
01105 <p class=MsoNormal style='text-align:justify'> </p>
01106 
01107 <p class=MsoNormal style='text-align:justify'><i>Example  – script FitAwmi.c:</i></p>
01108 
01109 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
01110 border=0 width=601 height=402 src="gif/spectrumfit_awni_image002.jpg"></span></i></p>
01111 
01112 <p class=MsoNormal style='text-align:justify'><b>Fig. 1 Original spectrum
01113 (black line) and fitted spectrum using AWMI algorithm (red line) and number of
01114 iteration steps = 1000. Positions of fitted peaks are denoted by markers</b></p>
01115 
01116 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
01117 
01118 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
01119 fitting function using AWMI algorithm.</span></p>
01120 
01121 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example,
01122 do</span></p>
01123 
01124 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x FitAwmi.C</span></p>
01125 
01126 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>&nbsp;</span></p>
01127 
01128 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>void FitAwmi() {</span></p>
01129 
01130 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t a;</span></p>
01131 
01132 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t
01133 i,nfound=0,bin;</span></p>
01134 
01135 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t nbins = 256;</span></p>
01136 
01137 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmin  = 0;</span></p>
01138 
01139 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmax  =
01140 nbins;</span></p>
01141 
01142 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Float_t * source =
01143 new float[nbins];</span></p>
01144 
01145 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Float_t * dest =
01146 new float[nbins];   </span></p>
01147 
01148 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TH1F *h = new
01149 TH1F(&quot;h&quot;,&quot;Fitting using AWMI algorithm&quot;,nbins,xmin,xmax);</span></p>
01150 
01151 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TH1F *d = new
01152 TH1F(&quot;d&quot;,&quot;&quot;,nbins,xmin,xmax);      </span></p>
01153 
01154 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TFile *f = new
01155 TFile(&quot;TSpectrum.root&quot;);</span></p>
01156 
01157 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   h=(TH1F*)
01158 f-&gt;Get(&quot;fit;1&quot;);   </span></p>
01159 
01160 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for (i = 0; i &lt;
01161 nbins; i++) source[i]=h-&gt;GetBinContent(i + 1);      </span></p>
01162 
01163 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TCanvas *Fit1 =
01164 gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Fit1&quot;);</span></p>
01165 
01166 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   if (!Fit1) Fit1 =
01167 new TCanvas(&quot;Fit1&quot;,&quot;Fit1&quot;,10,10,1000,700);</span></p>
01168 
01169 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01170 h-&gt;Draw(&quot;L&quot;);</span></p>
01171 
01172 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TSpectrum *s = new
01173 TSpectrum();</span></p>
01174 
01175 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   //searching for
01176 candidate peaks positions</span></p>
01177 
01178 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   nfound =
01179 s-&gt;SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);</span></p>
01180 
01181 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Bool_t *FixPos =
01182 new Bool_t[nfound];</span></p>
01183 
01184 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Bool_t *FixAmp =
01185 new Bool_t[nfound];      </span></p>
01186 
01187 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for(i = 0; i&lt;
01188 nfound ; i++){</span></p>
01189 
01190 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>      FixPos[i] =
01191 kFALSE;</span></p>
01192 
01193 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>      FixAmp[i] =
01194 kFALSE;    </span></p>
01195 
01196 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }</span></p>
01197 
01198 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   //filling in the
01199 initial estimates of the input parameters</span></p>
01200 
01201 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Float_t *PosX =
01202 new Float_t[nfound];         </span></p>
01203 
01204 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Float_t *PosY =
01205 new Float_t[nfound];</span></p>
01206 
01207 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   PosX =
01208 s-&gt;GetPositionX();</span></p>
01209 
01210 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for (i = 0; i &lt;
01211 nfound; i++) {</span></p>
01212 
01213 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>                                                a=PosX[i];</span></p>
01214 
01215 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        bin = 1 +
01216 Int_t(a + 0.5);</span></p>
01217 
01218 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        PosY[i] =
01219 h-&gt;GetBinContent(bin);</span></p>
01220 
01221 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }   </span></p>
01222 
01223 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TSpectrumFit
01224 *pfit=new TSpectrumFit(nfound);</span></p>
01225 
01226 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01227 pfit-&gt;SetFitParameters(xmin, xmax-1, 1000, 0.1, pfit-&gt;kFitOptimChiCounts,
01228 pfit-&gt;kFitAlphaHalving, pfit-&gt;kFitPower2,
01229 pfit-&gt;kFitTaylorOrderFirst);   </span></p>
01230 
01231 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   pfit-&gt;SetPeakParameters(2,
01232 kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);   </span></p>
01233 
01234 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01235 pfit-&gt;FitAwmi(source);</span></p>
01236 
01237 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t
01238 *CalcPositions = new Double_t[nfound];      </span></p>
01239 
01240 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Double_t
01241 *CalcAmplitudes = new Double_t[nfound];         </span></p>
01242 
01243 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01244 CalcPositions=pfit-&gt;GetPositions();</span></p>
01245 
01246 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01247 CalcAmplitudes=pfit-&gt;GetAmplitudes();   </span></p>
01248 
01249 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for (i = 0; i &lt;
01250 nbins; i++) d-&gt;SetBinContent(i + 1,source[i]);</span></p>
01251 
01252 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01253 d-&gt;SetLineColor(kRed);   </span></p>
01254 
01255 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01256 d-&gt;Draw(&quot;SAME L&quot;);  </span></p>
01257 
01258 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   for (i = 0; i &lt;
01259 nfound; i++) {</span></p>
01260 
01261 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>                                                a=CalcPositions[i];</span></p>
01262 
01263 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        bin = 1 +
01264 Int_t(a + 0.5);                </span></p>
01265 
01266 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        PosX[i] =
01267 d-&gt;GetBinCenter(bin);</span></p>
01268 
01269 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>        PosY[i] =
01270 d-&gt;GetBinContent(bin);</span></p>
01271 
01272 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }</span></p>
01273 
01274 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   TPolyMarker * pm =
01275 (TPolyMarker*)h-&gt;GetListOfFunctions()-&gt;FindObject(&quot;TPolyMarker&quot;);</span></p>
01276 
01277 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   if (pm) {</span></p>
01278 
01279 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>     
01280 h-&gt;GetListOfFunctions()-&gt;Remove(pm);</span></p>
01281 
01282 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>      delete pm;</span></p>
01283 
01284 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   }</span></p>
01285 
01286 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   pm = new
01287 TPolyMarker(nfound, PosX, PosY);</span></p>
01288 
01289 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01290 h-&gt;GetListOfFunctions()-&gt;Add(pm);</span></p>
01291 
01292 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01293 pm-&gt;SetMarkerStyle(23);</span></p>
01294 
01295 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01296 pm-&gt;SetMarkerColor(kRed);</span></p>
01297 
01298 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>  
01299 pm-&gt;SetMarkerSize(1);   </span></p>
01300 
01301 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>}</span></p>
01302 
01303 </div>
01304 
01305 <!-- */
01306 // --> End_Html
01307 
01308    Int_t i, j, k, shift =
01309        2 * fNPeaks + 7, peak_vel, rozmer, iter, pw, regul_cycle,
01310        flag;
01311    Double_t a, b, c, d = 0, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
01312        0, pi, pmin = 0, chi_cel = 0, chi_er;
01313    Double_t *working_space = new Double_t[5 * (2 * fNPeaks + 7)];
01314    for (i = 0, j = 0; i < fNPeaks; i++) {
01315       working_space[2 * i] = fAmpInit[i];        //vector parameter
01316       if (fFixAmp[i] == false) {
01317          working_space[shift + j] = fAmpInit[i];        //vector xk
01318          j += 1;
01319       }
01320       working_space[2 * i + 1] = fPositionInit[i];        //vector parameter
01321       if (fFixPosition[i] == false) {
01322          working_space[shift + j] = fPositionInit[i];        //vector xk
01323          j += 1;
01324       }
01325    }
01326    peak_vel = 2 * i;
01327    working_space[2 * i] = fSigmaInit;        //vector parameter
01328    if (fFixSigma == false) {
01329       working_space[shift + j] = fSigmaInit;        //vector xk
01330       j += 1;
01331    }
01332    working_space[2 * i + 1] = fTInit;        //vector parameter
01333    if (fFixT == false) {
01334       working_space[shift + j] = fTInit;        //vector xk
01335       j += 1;
01336    }
01337    working_space[2 * i + 2] = fBInit;        //vector parameter
01338    if (fFixB == false) {
01339       working_space[shift + j] = fBInit;        //vector xk
01340       j += 1;
01341    }
01342    working_space[2 * i + 3] = fSInit;        //vector parameter
01343    if (fFixS == false) {
01344       working_space[shift + j] = fSInit;        //vector xk
01345       j += 1;
01346    }
01347    working_space[2 * i + 4] = fA0Init;        //vector parameter
01348    if (fFixA0 == false) {
01349       working_space[shift + j] = fA0Init;        //vector xk
01350       j += 1;
01351    }
01352    working_space[2 * i + 5] = fA1Init;        //vector parameter
01353    if (fFixA1 == false) {
01354       working_space[shift + j] = fA1Init;        //vector xk
01355       j += 1;
01356    }
01357    working_space[2 * i + 6] = fA2Init;        //vector parameter
01358    if (fFixA2 == false) {
01359       working_space[shift + j] = fA2Init;        //vector xk
01360       j += 1;
01361    }
01362    rozmer = j;
01363    if (rozmer == 0){
01364       delete [] working_space;
01365       Error ("FitAwmi","All parameters are fixed");   
01366       return;
01367    }
01368    if (rozmer >= fXmax - fXmin + 1){
01369       delete [] working_space;
01370       Error ("FitAwmi","Number of fitted parameters is larger than # of fitted points");   
01371       return;    
01372    }
01373    for (iter = 0; iter < fNumberIterations; iter++) {
01374       for (j = 0; j < rozmer; j++) {
01375          working_space[2 * shift + j] = 0, working_space[3 * shift + j] = 0;        //der,temp
01376       }
01377       
01378           //filling vectors
01379       alpha = fAlpha;
01380       chi_opt = 0, pw = fPower - 2;
01381       for (i = fXmin; i <= fXmax; i++) {
01382          yw = source[i];
01383          ywm = yw;
01384          f = Shape(fNPeaks, (Double_t) i, working_space,
01385                     working_space[peak_vel], working_space[peak_vel + 1],
01386                     working_space[peak_vel + 3],
01387                     working_space[peak_vel + 2],
01388                     working_space[peak_vel + 4],
01389                     working_space[peak_vel + 5],
01390                     working_space[peak_vel + 6]);
01391          if (fStatisticType == kFitOptimMaxLikelihood) {
01392             if (f > 0.00001)
01393                chi_opt += yw * TMath::Log(f) - f;
01394          }
01395          
01396          else {
01397             if (ywm != 0)
01398                chi_opt += (yw - f) * (yw - f) / ywm;
01399          }
01400          if (fStatisticType == kFitOptimChiFuncValues) {
01401             ywm = f;
01402             if (f < 0.00001)
01403                ywm = 0.00001;
01404          }
01405          
01406          else if (fStatisticType == kFitOptimMaxLikelihood) {
01407             ywm = f;
01408             if (f < 0.001)
01409                ywm = 0.001;
01410          }
01411          
01412          else {
01413             if (ywm == 0)
01414                ywm = 1;
01415          }
01416          
01417              //calculation of gradient vector
01418              for (j = 0, k = 0; j < fNPeaks; j++) {
01419             if (fFixAmp[j] == false) {
01420                a = Deramp((Double_t) i, working_space[2 * j + 1],
01421                            working_space[peak_vel],
01422                            working_space[peak_vel + 1],
01423                            working_space[peak_vel + 3],
01424                            working_space[peak_vel + 2]);
01425                if (ywm != 0) {
01426                   c = Ourpowl(a, pw);
01427                   if (fStatisticType == kFitOptimChiFuncValues) {
01428                      b = a * (yw * yw - f * f) / (ywm * ywm);
01429                      working_space[2 * shift + k] += b * c;        //der
01430                      b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01431                      working_space[3 * shift + k] += b * c;        //temp
01432                   }
01433                   
01434                   else {
01435                      b = a * (yw - f) / ywm;
01436                      working_space[2 * shift + k] += b * c;        //der
01437                      b = a * a / ywm;
01438                      working_space[3 * shift + k] += b * c;        //temp
01439                   }
01440                }
01441                k += 1;
01442             }
01443             if (fFixPosition[j] == false) {
01444                a = Deri0((Double_t) i, working_space[2 * j],
01445                           working_space[2 * j + 1],
01446                           working_space[peak_vel],
01447                           working_space[peak_vel + 1],
01448                           working_space[peak_vel + 3],
01449                           working_space[peak_vel + 2]);
01450                if (fFitTaylor == kFitTaylorOrderSecond)
01451                   d = Derderi0((Double_t) i, working_space[2 * j],
01452                                 working_space[2 * j + 1],
01453                                 working_space[peak_vel]);
01454                if (ywm != 0) {
01455                   c = Ourpowl(a, pw);
01456                   if (TMath::Abs(a) > 0.00000001
01457                        && fFitTaylor == kFitTaylorOrderSecond) {
01458                      d = d * TMath::Abs(yw - f) / (2 * a * ywm);
01459                      if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
01460                         d = 0;
01461                   }
01462                   
01463                   else
01464                      d = 0;
01465                   a = a + d;
01466                   if (fStatisticType == kFitOptimChiFuncValues) {
01467                      b = a * (yw * yw - f * f) / (ywm * ywm);
01468                      working_space[2 * shift + k] += b * c;        //der
01469                      b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01470                      working_space[3 * shift + k] += b * c;        //temp
01471                   }
01472                   
01473                   else {
01474                      b = a * (yw - f) / ywm;
01475                      working_space[2 * shift + k] += b * c;        //der
01476                      b = a * a / ywm;
01477                      working_space[3 * shift + k] += b * c;        //temp
01478                   }
01479                }
01480                k += 1;
01481             }
01482          }
01483          if (fFixSigma == false) {
01484             a = Dersigma(fNPeaks, (Double_t) i, working_space,
01485                           working_space[peak_vel],
01486                           working_space[peak_vel + 1],
01487                           working_space[peak_vel + 3],
01488                           working_space[peak_vel + 2]);
01489             if (fFitTaylor == kFitTaylorOrderSecond)
01490                d = Derdersigma(fNPeaks, (Double_t) i,
01491                                 working_space, working_space[peak_vel]);
01492             if (ywm != 0) {
01493                c = Ourpowl(a, pw);
01494                if (TMath::Abs(a) > 0.00000001
01495                     && fFitTaylor == kFitTaylorOrderSecond) {
01496                   d = d * TMath::Abs(yw - f) / (2 * a * ywm);
01497                   if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
01498                      d = 0;
01499                }
01500                
01501                else
01502                   d = 0;
01503                a = a + d;
01504                if (fStatisticType == kFitOptimChiFuncValues) {
01505                   b = a * (yw * yw - f * f) / (ywm * ywm);
01506                   working_space[2 * shift + k] += b * c;        //der
01507                   b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01508                   working_space[3 * shift + k] += b * c;        //temp
01509                }
01510                
01511                else {
01512                   b = a * (yw - f) / ywm;
01513                   working_space[2 * shift + k] += b * c;        //der
01514                   b = a * a / ywm;
01515                   working_space[3 * shift + k] += b * c;        //temp
01516                }
01517             }
01518             k += 1;
01519          }
01520          if (fFixT == false) {
01521             a = Dert(fNPeaks, (Double_t) i, working_space,
01522                       working_space[peak_vel],
01523                       working_space[peak_vel + 2]);
01524             if (ywm != 0) {
01525                c = Ourpowl(a, pw);
01526                if (fStatisticType == kFitOptimChiFuncValues) {
01527                   b = a * (yw * yw - f * f) / (ywm * ywm);
01528                   working_space[2 * shift + k] += b * c;        //der
01529                   b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01530                   working_space[3 * shift + k] += b * c;        //temp
01531                }
01532                
01533                else {
01534                   b = a * (yw - f) / ywm;
01535                   working_space[2 * shift + k] += b * c;        //der
01536                   b = a * a / ywm;
01537                   working_space[3 * shift + k] += b * c;        //temp
01538                }
01539             }
01540             k += 1;
01541          }
01542          if (fFixB == false) {
01543             a = Derb(fNPeaks, (Double_t) i, working_space,
01544                       working_space[peak_vel], working_space[peak_vel + 1],
01545                       working_space[peak_vel + 2]);
01546             if (ywm != 0) {
01547                c = Ourpowl(a, pw);
01548                if (fStatisticType == kFitOptimChiFuncValues) {
01549                   b = a * (yw * yw - f * f) / (ywm * ywm);
01550                   working_space[2 * shift + k] += b * c;        //der
01551                   b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01552                   working_space[3 * shift + k] += b * c;        //temp
01553                }
01554                
01555                else {
01556                   b = a * (yw - f) / ywm;
01557                   working_space[2 * shift + k] += b * c;        //der
01558                   b = a * a / ywm;
01559                   working_space[3 * shift + k] += b * c;        //temp
01560                }
01561             }
01562             k += 1;
01563          }
01564          if (fFixS == false) {
01565             a = Ders(fNPeaks, (Double_t) i, working_space,
01566                       working_space[peak_vel]);
01567             if (ywm != 0) {
01568                c = Ourpowl(a, pw);
01569                if (fStatisticType == kFitOptimChiFuncValues) {
01570                   b = a * (yw * yw - f * f) / (ywm * ywm);
01571                   working_space[2 * shift + k] += b * c;        //der
01572                   b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01573                   working_space[3 * shift + k] += b * c;        //temp
01574                }
01575                
01576                else {
01577                   b = a * (yw - f) / ywm;
01578                   working_space[2 * shift + k] += b * c;        //der
01579                   b = a * a / ywm;
01580                   working_space[3 * shift + k] += b * c;        //temp
01581                }
01582             }
01583             k += 1;
01584          }
01585          if (fFixA0 == false) {
01586             a = 1.;
01587             if (ywm != 0) {
01588                c = Ourpowl(a, pw);
01589                if (fStatisticType == kFitOptimChiFuncValues) {
01590                   b = a * (yw * yw - f * f) / (ywm * ywm);
01591                   working_space[2 * shift + k] += b * c;        //der
01592                   b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01593                   working_space[3 * shift + k] += b * c;        //temp
01594                }
01595                
01596                else {
01597                   b = a * (yw - f) / ywm;
01598                   working_space[2 * shift + k] += b * c;        //der
01599                   b = a * a / ywm;
01600                   working_space[3 * shift + k] += b * c;        //temp
01601                }
01602             }
01603             k += 1;
01604          }
01605          if (fFixA1 == false) {
01606             a = Dera1((Double_t) i);
01607             if (ywm != 0) {
01608                c = Ourpowl(a, pw);
01609                if (fStatisticType == kFitOptimChiFuncValues) {
01610                   b = a * (yw * yw - f * f) / (ywm * ywm);
01611                   working_space[2 * shift + k] += b * c;        //der
01612                   b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01613                   working_space[3 * shift + k] += b * c;        //temp
01614                }
01615                
01616                else {
01617                   b = a * (yw - f) / ywm;
01618                   working_space[2 * shift + k] += b * c;        //der
01619                   b = a * a / ywm;
01620                   working_space[3 * shift + k] += b * c;        //temp
01621                }
01622             }
01623             k += 1;
01624          }
01625          if (fFixA2 == false) {
01626             a = Dera2((Double_t) i);
01627             if (ywm != 0) {
01628                c = Ourpowl(a, pw);
01629                if (fStatisticType == kFitOptimChiFuncValues) {
01630                   b = a * (yw * yw - f * f) / (ywm * ywm);
01631                   working_space[2 * shift + k] += b * c;        //der
01632                   b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
01633                   working_space[3 * shift + k] += b * c;        //temp
01634                }
01635                
01636                else {
01637                   b = a * (yw - f) / ywm;
01638                   working_space[2 * shift + k] += b * c;        //der
01639                   b = a * a / ywm;
01640                   working_space[3 * shift + k] += b * c;        //temp
01641                }
01642             }
01643             k += 1;
01644          }
01645       }
01646       for (j = 0; j < rozmer; j++) {
01647          if (TMath::Abs(working_space[3 * shift + j]) > 0.000001)
01648             working_space[2 * shift + j] = working_space[2 * shift + j] / TMath::Abs(working_space[3 * shift + j]);        //der[j]=der[j]/temp[j]
01649          else
01650             working_space[2 * shift + j] = 0;        //der[j]
01651       }
01652       
01653       //calculate chi_opt
01654       chi2 = chi_opt;
01655       chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
01656       
01657       //calculate new parameters
01658       regul_cycle = 0;
01659       for (j = 0; j < rozmer; j++) {
01660          working_space[4 * shift + j] = working_space[shift + j];        //temp_xk[j]=xk[j]
01661       }
01662       
01663       do {
01664          if (fAlphaOptim == kFitAlphaOptimal) {
01665             if (fStatisticType != kFitOptimMaxLikelihood)
01666                chi_min = 10000 * chi2;
01667             
01668             else
01669                chi_min = 0.1 * chi2;
01670             flag = 0;
01671             for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
01672                for (j = 0; j < rozmer; j++) {
01673                   working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j];        //xk[j]=temp_xk[j]+pi*alpha*der[j]
01674                }
01675                for (i = 0, j = 0; i < fNPeaks; i++) {
01676                   if (fFixAmp[i] == false) {
01677                      if (working_space[shift + j] < 0)        //xk[j]
01678                         working_space[shift + j] = 0;        //xk[j]
01679                      working_space[2 * i] = working_space[shift + j];        //parameter[2*i]=xk[j]
01680                      j += 1;
01681                   }
01682                   if (fFixPosition[i] == false) {
01683                      if (working_space[shift + j] < fXmin)        //xk[j]
01684                         working_space[shift + j] = fXmin;        //xk[j]
01685                      if (working_space[shift + j] > fXmax)        //xk[j]
01686                         working_space[shift + j] = fXmax;        //xk[j]
01687                      working_space[2 * i + 1] = working_space[shift + j];        //parameter[2*i+1]=xk[j]
01688                      j += 1;
01689                   }
01690                }
01691                if (fFixSigma == false) {
01692                   if (working_space[shift + j] < 0.001) {        //xk[j]
01693                      working_space[shift + j] = 0.001;        //xk[j]
01694                   }
01695                   working_space[peak_vel] = working_space[shift + j];        //parameter[peak_vel]=xk[j]
01696                   j += 1;
01697                }
01698                if (fFixT == false) {
01699                   working_space[peak_vel + 1] = working_space[shift + j];        //parameter[peak_vel+1]=xk[j]
01700                   j += 1;
01701                }
01702                if (fFixB == false) {
01703                   if (TMath::Abs(working_space[shift + j]) < 0.001) {        //xk[j]
01704                      if (working_space[shift + j] < 0)        //xk[j]
01705                         working_space[shift + j] = -0.001;        //xk[j]
01706                      else
01707                         working_space[shift + j] = 0.001;        //xk[j]
01708                   }
01709                   working_space[peak_vel + 2] = working_space[shift + j];        //parameter[peak_vel+2]=xk[j]
01710                   j += 1;
01711                }
01712                if (fFixS == false) {
01713                   working_space[peak_vel + 3] = working_space[shift + j];        //parameter[peak_vel+3]=xk[j]
01714                   j += 1;
01715                }
01716                if (fFixA0 == false) {
01717                   working_space[peak_vel + 4] = working_space[shift + j];        //parameter[peak_vel+4]=xk[j]
01718                   j += 1;
01719                }
01720                if (fFixA1 == false) {
01721                   working_space[peak_vel + 5] = working_space[shift + j];        //parameter[peak_vel+5]=xk[j]
01722                   j += 1;
01723                }
01724                if (fFixA2 == false) {
01725                   working_space[peak_vel + 6] = working_space[shift + j];        //parameter[peak_vel+6]=xk[j]
01726                   j += 1;
01727                }
01728                chi2 = 0;
01729                for (i = fXmin; i <= fXmax; i++) {
01730                   yw = source[i];
01731                   ywm = yw;
01732                   f = Shape(fNPeaks, (Double_t) i, working_space,
01733                   working_space[peak_vel],
01734                   working_space[peak_vel + 1],
01735                   working_space[peak_vel + 3],
01736                   working_space[peak_vel + 2],
01737                   working_space[peak_vel + 4],
01738                   working_space[peak_vel + 5],
01739                   working_space[peak_vel + 6]);
01740                   if (fStatisticType == kFitOptimChiFuncValues) {
01741                      ywm = f;
01742                      if (f < 0.00001)
01743                         ywm = 0.00001;
01744                   }
01745                   if (fStatisticType == kFitOptimMaxLikelihood) {
01746                      if (f > 0.00001)
01747                         chi2 += yw * TMath::Log(f) - f;
01748                   }
01749                   
01750                   else {
01751                      if (ywm != 0)
01752                         chi2 += (yw - f) * (yw - f) / ywm;
01753                   }
01754                }
01755                if ((chi2 < chi_min
01756                     && fStatisticType != kFitOptimMaxLikelihood)
01757                     || (chi2 > chi_min
01758                     && fStatisticType == kFitOptimMaxLikelihood)) {
01759                   pmin = pi, chi_min = chi2;
01760                }
01761                
01762                else
01763                   flag = 1;
01764                if (pi == 0.1)
01765                   chi_min = chi2;
01766                chi = chi_min;
01767             }
01768             if (pmin != 0.1) {
01769                for (j = 0; j < rozmer; j++) {
01770                   working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j];        //xk[j]=temp_xk[j]+pmin*alpha*der[j]
01771                }
01772                for (i = 0, j = 0; i < fNPeaks; i++) {
01773                   if (fFixAmp[i] == false) {
01774                      if (working_space[shift + j] < 0)        //xk[j]
01775                         working_space[shift + j] = 0;        //xk[j]
01776                      working_space[2 * i] = working_space[shift + j];        //parameter[2*i]=xk[j]
01777                      j += 1;
01778                   }
01779                   if (fFixPosition[i] == false) {
01780                      if (working_space[shift + j] < fXmin)        //xk[j]
01781                         working_space[shift + j] = fXmin;        //xk[j]
01782                      if (working_space[shift + j] > fXmax)        //xk[j]
01783                         working_space[shift + j] = fXmax;        //xk[j]
01784                      working_space[2 * i + 1] = working_space[shift + j];        //parameter[2*i+1]=xk[j]
01785                      j += 1;
01786                   }
01787                }
01788                if (fFixSigma == false) {
01789                   if (working_space[shift + j] < 0.001) {        //xk[j]
01790                      working_space[shift + j] = 0.001;        //xk[j]
01791                   }
01792                   working_space[peak_vel] = working_space[shift + j];        //parameter[peak_vel]=xk[j]
01793                   j += 1;
01794                }
01795                if (fFixT == false) {
01796                   working_space[peak_vel + 1] = working_space[shift + j];        //parameter[peak_vel+1]=xk[j]
01797                   j += 1;
01798                }
01799                if (fFixB == false) {
01800                   if (TMath::Abs(working_space[shift + j]) < 0.001) {        //xk[j]
01801                      if (working_space[shift + j] < 0)        //xk[j]
01802                         working_space[shift + j] = -0.001;        //xk[j]
01803                      else
01804                         working_space[shift + j] = 0.001;        //xk[j]
01805                   }
01806                   working_space[peak_vel + 2] = working_space[shift + j];        //parameter[peak_vel+2]=xk[j]
01807                   j += 1;
01808                }
01809                if (fFixS == false) {
01810                   working_space[peak_vel + 3] = working_space[shift + j];        //parameter[peak_vel+3]=xk[j]
01811                   j += 1;
01812                }
01813                if (fFixA0 == false) {
01814                   working_space[peak_vel + 4] = working_space[shift + j];        //parameter[peak_vel+4]=xk[j]
01815                   j += 1;
01816                }
01817                if (fFixA1 == false) {
01818                   working_space[peak_vel + 5] = working_space[shift + j];        //parameter[peak_vel+5]=xk[j]
01819                   j += 1;
01820                }
01821                if (fFixA2 == false) {
01822                   working_space[peak_vel + 6] = working_space[shift + j];        //parameter[peak_vel+6]=xk[j]
01823                   j += 1;
01824                }
01825                chi = chi_min;
01826             }
01827          }
01828          
01829          else {
01830             for (j = 0; j < rozmer; j++) {
01831                working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j];        //xk[j]=temp_xk[j]+pi*alpha*der[j]
01832             }
01833             for (i = 0, j = 0; i < fNPeaks; i++) {
01834                if (fFixAmp[i] == false) {
01835                   if (working_space[shift + j] < 0)        //xk[j]
01836                      working_space[shift + j] = 0;        //xk[j]
01837                   working_space[2 * i] = working_space[shift + j];        //parameter[2*i]=xk[j]
01838                   j += 1;
01839                }
01840                if (fFixPosition[i] == false) {
01841                   if (working_space[shift + j] < fXmin)        //xk[j]
01842                      working_space[shift + j] = fXmin;        //xk[j]
01843                   if (working_space[shift + j] > fXmax)        //xk[j]
01844                      working_space[shift + j] = fXmax;        //xk[j]
01845                   working_space[2 * i + 1] = working_space[shift + j];        //parameter[2*i+1]=xk[j]
01846                   j += 1;
01847                }
01848             }
01849             if (fFixSigma == false) {
01850                if (working_space[shift + j] < 0.001) {        //xk[j]
01851                   working_space[shift + j] = 0.001;        //xk[j]
01852                }
01853                working_space[peak_vel] = working_space[shift + j];        //parameter[peak_vel]=xk[j]
01854                j += 1;
01855             }
01856             if (fFixT == false) {
01857                working_space[peak_vel + 1] = working_space[shift + j];        //parameter[peak_vel+1]=xk[j]
01858                j += 1;
01859             }
01860             if (fFixB == false) {
01861                if (TMath::Abs(working_space[shift + j]) < 0.001) {        //xk[j]
01862                   if (working_space[shift + j] < 0)        //xk[j]
01863                      working_space[shift + j] = -0.001;        //xk[j]
01864                   else
01865                      working_space[shift + j] = 0.001;        //xk[j]
01866                }
01867                working_space[peak_vel + 2] = working_space[shift + j];        //parameter[peak_vel+2]=xk[j]
01868                j += 1;
01869             }
01870             if (fFixS == false) {
01871                working_space[peak_vel + 3] = working_space[shift + j];        //parameter[peak_vel+3]=xk[j]
01872                j += 1;
01873             }
01874             if (fFixA0 == false) {
01875                working_space[peak_vel + 4] = working_space[shift + j];        //parameter[peak_vel+4]=xk[j]
01876                j += 1;
01877             }
01878             if (fFixA1 == false) {
01879                working_space[peak_vel + 5] = working_space[shift + j];        //parameter[peak_vel+5]=xk[j]
01880                j += 1;
01881             }
01882             if (fFixA2 == false) {
01883                working_space[peak_vel + 6] = working_space[shift + j];        //parameter[peak_vel+6]=xk[j]
01884                j += 1;
01885             }
01886             chi = 0;
01887             for (i = fXmin; i <= fXmax; i++) {
01888                yw = source[i];
01889                ywm = yw;
01890                f = Shape(fNPeaks, (Double_t) i, working_space,
01891                working_space[peak_vel],
01892                working_space[peak_vel + 1],
01893                working_space[peak_vel + 3],
01894                working_space[peak_vel + 2],
01895                working_space[peak_vel + 4],
01896                working_space[peak_vel + 5],
01897                working_space[peak_vel + 6]);
01898                if (fStatisticType == kFitOptimChiFuncValues) {
01899                   ywm = f;
01900                   if (f < 0.00001)
01901                      ywm = 0.00001;
01902                }
01903                if (fStatisticType == kFitOptimMaxLikelihood) {
01904                   if (f > 0.00001)
01905                      chi += yw * TMath::Log(f) - f;
01906                }
01907                
01908                else {
01909                   if (ywm != 0)
01910                      chi += (yw - f) * (yw - f) / ywm;
01911                }
01912             }
01913          }
01914          chi2 = chi;
01915          chi = TMath::Sqrt(TMath::Abs(chi));
01916          if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
01917             alpha = alpha * chi_opt / (2 * chi);
01918          
01919          else if (fAlphaOptim == kFitAlphaOptimal)
01920             alpha = alpha / 10.0;
01921          iter += 1;
01922          regul_cycle += 1;
01923       } while (((chi > chi_opt
01924                  && fStatisticType != kFitOptimMaxLikelihood)
01925                  || (chi < chi_opt
01926                  && fStatisticType == kFitOptimMaxLikelihood))
01927                 && regul_cycle < kFitNumRegulCycles);
01928       for (j = 0; j < rozmer; j++) {
01929          working_space[4 * shift + j] = 0;        //temp_xk[j]
01930          working_space[2 * shift + j] = 0;        //der[j]
01931       }
01932       for (i = fXmin, chi_cel = 0; i <= fXmax; i++) {
01933          yw = source[i];
01934          if (yw == 0)
01935             yw = 1;
01936          f = Shape(fNPeaks, (Double_t) i, working_space,
01937          working_space[peak_vel], working_space[peak_vel + 1],
01938          working_space[peak_vel + 3],
01939          working_space[peak_vel + 2],
01940          working_space[peak_vel + 4],
01941          working_space[peak_vel + 5],
01942          working_space[peak_vel + 6]);
01943          chi_opt = (yw - f) * (yw - f) / yw;
01944          chi_cel += (yw - f) * (yw - f) / yw;
01945          
01946              //calculate gradient vector
01947          for (j = 0, k = 0; j < fNPeaks; j++) {
01948             if (fFixAmp[j] == false) {
01949                a = Deramp((Double_t) i, working_space[2 * j + 1],
01950                working_space[peak_vel],
01951                working_space[peak_vel + 1],
01952                working_space[peak_vel + 3],
01953                working_space[peak_vel + 2]);
01954                if (yw != 0) {
01955                   c = Ourpowl(a, pw);
01956                   working_space[2 * shift + k] += chi_opt * c;        //der[k]
01957                   b = a * a / yw;
01958                   working_space[4 * shift + k] += b * c;        //temp_xk[k]
01959                }
01960                k += 1;
01961             }
01962             if (fFixPosition[j] == false) {
01963                a = Deri0((Double_t) i, working_space[2 * j],
01964                working_space[2 * j + 1],
01965                working_space[peak_vel],
01966                working_space[peak_vel + 1],
01967                working_space[peak_vel + 3],
01968                working_space[peak_vel + 2]);
01969                if (yw != 0) {
01970                   c = Ourpowl(a, pw);
01971                   working_space[2 * shift + k] += chi_opt * c;        //der[k]
01972                   b = a * a / yw;
01973                   working_space[4 * shift + k] += b * c;        //temp_xk[k]
01974                }
01975                k += 1;
01976             }
01977          }
01978          if (fFixSigma == false) {
01979             a = Dersigma(fNPeaks, (Double_t) i, working_space,
01980             working_space[peak_vel],
01981             working_space[peak_vel + 1],
01982             working_space[peak_vel + 3],
01983             working_space[peak_vel + 2]);
01984             if (yw != 0) {
01985                c = Ourpowl(a, pw);
01986                working_space[2 * shift + k] += chi_opt * c;        //der[k]
01987                b = a * a / yw;
01988                working_space[4 * shift + k] += b * c;        //temp_xk[k]
01989             }
01990             k += 1;
01991          }
01992          if (fFixT == false) {
01993             a = Dert(fNPeaks, (Double_t) i, working_space,
01994             working_space[peak_vel],
01995             working_space[peak_vel + 2]);
01996             if (yw != 0) {
01997                c = Ourpowl(a, pw);
01998                working_space[2 * shift + k] += chi_opt * c;        //der[k]
01999                b = a * a / yw;
02000                working_space[4 * shift + k] += b * c;        //temp_xk[k]
02001             }
02002             k += 1;
02003          }
02004          if (fFixB == false) {
02005             a = Derb(fNPeaks, (Double_t) i, working_space,
02006             working_space[peak_vel], working_space[peak_vel + 1],
02007             working_space[peak_vel + 2]);
02008             if (yw != 0) {
02009                c = Ourpowl(a, pw);
02010                working_space[2 * shift + k] += chi_opt * c;        //der[k]
02011                b = a * a / yw;
02012                working_space[4 * shift + k] += b * c;        //temp_xk[k]
02013             }
02014             k += 1;
02015          }
02016          if (fFixS == false) {
02017             a = Ders(fNPeaks, (Double_t) i, working_space,
02018             working_space[peak_vel]);
02019             if (yw != 0) {
02020                c = Ourpowl(a, pw);
02021                working_space[2 * shift + k] += chi_opt * c;        //der[k]
02022                b = a * a / yw;
02023                working_space[4 * shift + k] += b * c;        //tem_xk[k]
02024             }
02025             k += 1;
02026          }
02027          if (fFixA0 == false) {
02028             a = 1.0;
02029             if (yw != 0) {
02030                c = Ourpowl(a, pw);
02031                working_space[2 * shift + k] += chi_opt * c;        //der[k]
02032                b = a * a / yw;
02033                working_space[4 * shift + k] += b * c;        //temp_xk[k]
02034             }
02035             k += 1;
02036          }
02037          if (fFixA1 == false) {
02038             a = Dera1((Double_t) i);
02039             if (yw != 0) {
02040                c = Ourpowl(a, pw);
02041                working_space[2 * shift + k] += chi_opt * c;        //der[k]
02042                b = a * a / yw;
02043                working_space[4 * shift + k] += b * c;        //temp_xk[k]
02044             }
02045             k += 1;
02046          }
02047          if (fFixA2 == false) {
02048             a = Dera2((Double_t) i);
02049             if (yw != 0) {
02050                c = Ourpowl(a, pw);
02051                working_space[2 * shift + k] += chi_opt * c;        //der[k]
02052                b = a * a / yw;
02053                working_space[4 * shift + k] += b * c;        //temp_xk[k]
02054             }
02055             k += 1;
02056          }
02057       }
02058    }
02059    b = fXmax - fXmin + 1 - rozmer;
02060    chi_er = chi_cel / b;
02061    for (i = 0, j = 0; i < fNPeaks; i++) {
02062       fArea[i] =
02063           Area(working_space[2 * i], working_space[peak_vel],
02064                working_space[peak_vel + 1], working_space[peak_vel + 2]);
02065       if (fFixAmp[i] == false) {
02066          fAmpCalc[i] = working_space[shift + j];        //xk[j]
02067          if (working_space[3 * shift + j] != 0)
02068             fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02069          if (fArea[i] > 0) {
02070             a = Derpa(working_space[peak_vel],
02071                        working_space[peak_vel + 1],
02072                        working_space[peak_vel + 2]);
02073             b = working_space[4 * shift + j];        //temp_xk[j]
02074             if (b == 0)
02075                b = 1;
02076             
02077             else
02078                b = 1 / b;
02079             fAreaErr[i] = TMath::Sqrt(TMath::Abs(a * a * b * chi_er));
02080          }
02081          
02082          else
02083             fAreaErr[i] = 0;
02084          j += 1;
02085       }
02086       
02087       else {
02088          fAmpCalc[i] = fAmpInit[i];
02089          fAmpErr[i] = 0;
02090          fAreaErr[i] = 0;
02091       }
02092       if (fFixPosition[i] == false) {
02093          fPositionCalc[i] = working_space[shift + j];        //xk[j]
02094          if (working_space[3 * shift + j] != 0)        //temp[j]
02095             fPositionErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02096          j += 1;
02097       }
02098       
02099       else {
02100          fPositionCalc[i] = fPositionInit[i];
02101          fPositionErr[i] = 0;
02102       }
02103    }
02104    if (fFixSigma == false) {
02105       fSigmaCalc = working_space[shift + j];        //xk[j]
02106       if (working_space[3 * shift + j] != 0)        //temp[j]
02107          fSigmaErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02108       j += 1;
02109    }
02110    
02111    else {
02112       fSigmaCalc = fSigmaInit;
02113       fSigmaErr = 0;
02114    }
02115    if (fFixT == false) {
02116       fTCalc = working_space[shift + j];        //xk[j]
02117       if (working_space[3 * shift + j] != 0)        //temp[j]
02118          fTErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02119       j += 1;
02120    }
02121    
02122    else {
02123       fTCalc = fTInit;
02124       fTErr = 0;
02125    }
02126    if (fFixB == false) {
02127       fBCalc = working_space[shift + j];        //xk[j]
02128       if (working_space[3 * shift + j] != 0)        //temp[j]
02129          fBErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02130       j += 1;
02131    }
02132    
02133    else {
02134       fBCalc = fBInit;
02135       fBErr = 0;
02136    }
02137    if (fFixS == false) {
02138       fSCalc = working_space[shift + j];        //xk[j]
02139       if (working_space[3 * shift + j] != 0)        //temp[j]
02140          fSErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02141       j += 1;
02142    }
02143    
02144    else {
02145       fSCalc = fSInit;
02146       fSErr = 0;
02147    }
02148    if (fFixA0 == false) {
02149       fA0Calc = working_space[shift + j];        //xk[j]
02150       if (working_space[3 * shift + j] != 0)        //temp[j]
02151          fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02152       j += 1;
02153    }
02154    
02155    else {
02156       fA0Calc = fA0Init;
02157       fA0Err = 0;
02158    }
02159    if (fFixA1 == false) {
02160       fA1Calc = working_space[shift + j];        //xk[j]
02161       if (working_space[3 * shift + j] != 0)        //temp[j]
02162          fA1Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02163       j += 1;
02164    }
02165    
02166    else {
02167       fA1Calc = fA1Init;
02168       fA1Err = 0;
02169    }
02170    if (fFixA2 == false) {
02171       fA2Calc = working_space[shift + j];        //xk[j]
02172       if (working_space[3 * shift + j] != 0)        //temp[j]
02173          fA2Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
02174       j += 1;
02175    }
02176    
02177    else {
02178       fA2Calc = fA2Init;
02179       fA2Err = 0;
02180    }
02181    b = fXmax - fXmin + 1 - rozmer;
02182    fChi = chi_cel / b;
02183    for (i = fXmin; i <= fXmax; i++) {
02184       f = Shape(fNPeaks, (Double_t) i, working_space,
02185                  working_space[peak_vel], working_space[peak_vel + 1],
02186                  working_space[peak_vel + 3], working_space[peak_vel + 2],
02187                  working_space[peak_vel + 4], working_space[peak_vel + 5],
02188                  working_space[peak_vel + 6]);
02189       source[i] = f;
02190    }
02191    delete[]working_space;
02192    return;
02193 }
02194     
02195 /////////////////FITTING FUNCTION WITH MATRIX INVERSION///////////////////////////////////////
02196 
02197 //_______________________________________________________________________________
02198 void TSpectrumFit::StiefelInversion(Double_t **a, Int_t size)
02199 {   
02200 //////////////////////////////////////////////////////////////////////////////////
02201 //   AUXILIARY FUNCTION                                                          //
02202 //                                                                              //
02203 //   This function calculates solution of the system of linear equations.        //
02204 //   The matrix a should have a dimension size*(size+4)                         //
02205 //   The calling function should fill in the matrix, the column size should     //
02206 //   contain vector y (right side of the system of equations). The result is    //
02207 //   placed into size+1 column of the matrix.                                   //
02208 //   according to sigma of peaks.                                               //
02209 //      Function parameters:                                                    //
02210 //              -a-matrix with dimension size*(size+4)                          //                                            //
02211 //              -size-number of rows of the matrix                              //
02212 //                                                                              //
02213 //////////////////////////////////////////////////////////////////////////////////
02214    Int_t i, j, k = 0;
02215    Double_t sk = 0, b, lambdak, normk, normk_old = 0;
02216    
02217    do {
02218       normk = 0;
02219       
02220           //calculation of rk and norm
02221           for (i = 0; i < size; i++) {
02222          a[i][size + 2] = -a[i][size];        //rk=-C
02223          for (j = 0; j < size; j++) {
02224             a[i][size + 2] += a[i][j] * a[j][size + 1];        //A*xk-C
02225          }
02226          normk += a[i][size + 2] * a[i][size + 2];        //calculation normk
02227       }
02228       
02229       //calculation of sk
02230       if (k != 0) {
02231          sk = normk / normk_old;
02232       }
02233       
02234       //calculation of uk
02235       for (i = 0; i < size; i++) {
02236          a[i][size + 3] = -a[i][size + 2] + sk * a[i][size + 3];        //uk=-rk+sk*uk-1
02237       }
02238       
02239       //calculation of lambdak
02240       lambdak = 0;
02241       for (i = 0; i < size; i++) {
02242          for (j = 0, b = 0; j < size; j++) {
02243             b += a[i][j] * a[j][size + 3];        //A*uk
02244          }
02245          lambdak += b * a[i][size + 3];
02246       }
02247       if (TMath::Abs(lambdak) > 1e-50)        //computer zero
02248          lambdak = normk / lambdak;
02249       
02250       else
02251          lambdak = 0;
02252       for (i = 0; i < size; i++)
02253          a[i][size + 1] += lambdak * a[i][size + 3];        //xk+1=xk+lambdak*uk
02254       normk_old = normk;
02255       k += 1;
02256    } while (k < size && TMath::Abs(normk) > 1e-50);        //computer zero
02257    return;
02258 }
02259 
02260 //_______________________________________________________________________________
02261 void TSpectrumFit::FitStiefel(Float_t *source) 
02262 {
02263 /////////////////////////////////////////////////////////////////////////////
02264 //        ONE-DIMENSIONAL FIT FUNCTION                                    
02265 //        ALGORITHM WITH MATRIX INVERSION (STIEFEL-HESTENS METHOD)            
02266 //        This function fits the source spectrum. The calling program should 
02267 //        fill in input parameters
02268 //        The fitted parameters are written into           
02269 //        output parameters and fitted data are written into 
02270 //        source spectrum.                                                    
02271 //                                                                          
02272 //        Function parameters:                                             
02273 //        source-pointer to the vector of source spectrum                  
02274 //                                                                         
02275 /////////////////////////////////////////////////////////////////////////////
02276 //Begin_Html <!--
02277 /* -->
02278 <div class=Section3>
02279 
02280 <p class=MsoNormal><b><span style='font-size:14.0pt'>Stiefel fitting algorithm</span></b></p>
02281 
02282 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>&nbsp;</span></i></p>
02283 
02284 <p class=MsoNormal><i>Function:</i></p>
02285 
02286 <p class=MsoNormal style='text-align:justify'>void <a
02287 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrumFit::</b></a>FitStiefel(<a
02288 href="http://root.cern.ch/root/html/ListOfTypes.html#float"><b>float</b></a> *fSource)
02289 </p>
02290 
02291 <p class=MsoNormal style='text-align:justify'>&nbsp;</p>
02292 
02293 <p class=MsoNormal style='text-align:justify'>This function fits the source
02294 spectrum using Stiefel-Hestens method [1] (see Awmi function).  The calling
02295 program should fill in input fitting parameters of the TSpectrumFit class using
02296 a set of TSpectrumFit setters. The fitted parameters are written into the class
02297 and the fitted data are written into source spectrum. It converges faster than
02298 Awmi method.</p>
02299 
02300 <p class=MsoNormal>&nbsp;</p>
02301 
02302 <p class=MsoNormal><i><span style='color:red'>Parameter:</span></i></p>
02303 
02304 <p class=MsoNormal style='text-align:justify'>        <b>fSource</b>-pointer to
02305 the vector of source spectrum                  </p>
02306 
02307 <p class=MsoNormal style='text-align:justify'>        </p>
02308 
02309 <p class=MsoNormal style='text-align:justify'><i>Example – script FitStiefel.c:</i></p>
02310 
02311 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'><img
02312 border=0 width=601 height=402 src="gif/spectrumfit_stiefel_image001.jpg"></span></i></p>
02313 
02314 <p class=MsoNormal style='text-align:justify'><b>Fig. 2 Original spectrum
02315 (black line) and fitted spectrum using Stiefel-Hestens method (red line) and
02316 number of iteration steps = 100. Positions of fitted peaks are denoted by
02317 markers</b></p>
02318 
02319 <p class=MsoNormal><b><span style='color:#339966'>&nbsp;</span></b></p>
02320 
02321 <p class=MsoNormal><b><span style='color:#339966'>Script:</span></b></p>
02322 
02323 <p class=MsoNormal><span style='font-size:10.0pt'>// Example to illustrate
02324 fitting function using Stiefel-Hestens method.</span></p>
02325 
02326 <p class=MsoNormal><span style='font-size:10.0pt'>// To execute this example, do</span></p>
02327 
02328 <p class=MsoNormal><span style='font-size:10.0pt'>// root &gt; .x FitStiefel.C</span></p>
02329 
02330 <p class=MsoNormal><span style='font-size:10.0pt'>&nbsp;</span></p>
02331 
02332 <p class=MsoNormal><span style='font-size:10.0pt'>void FitStiefel() {</span></p>
02333 
02334 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t a;</span></p>
02335 
02336 <p class=MsoNormal><span style='font-size:10.0pt'>   Int_t i,nfound=0,bin;</span></p>
02337 
02338 <p class=MsoNormal><span style='font-size:10.0pt'>   Int_t nbins = 256;</span></p>
02339 
02340 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
02341 style='font-size:10.0pt'>Int_t xmin  = 0;</span></p>
02342 
02343 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   Int_t xmax  =
02344 nbins;</span></p>
02345 
02346 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
02347 style='font-size:10.0pt'>Float_t * source = new float[nbins];</span></p>
02348 
02349 <p class=MsoNormal><span style='font-size:10.0pt'>   Float_t * dest = new
02350 float[nbins];   </span></p>
02351 
02352 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *h = new TH1F(&quot;h&quot;,&quot;Fitting
02353 using AWMI algorithm&quot;,nbins,xmin,xmax);</span></p>
02354 
02355 <p class=MsoNormal><span style='font-size:10.0pt'>   TH1F *d = new
02356 TH1F(&quot;d&quot;,&quot;&quot;,nbins,xmin,xmax);      </span></p>
02357 
02358 <p class=MsoNormal><span style='font-size:10.0pt'>   TFile *f = new
02359 TFile(&quot;TSpectrum.root&quot;);</span></p>
02360 
02361 <p class=MsoNormal><span style='font-size:10.0pt'>   h=(TH1F*)
02362 f-&gt;Get(&quot;fit;1&quot;);   </span></p>
02363 
02364 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
02365 i++) source[i]=h-&gt;GetBinContent(i + 1);      </span></p>
02366 
02367 <p class=MsoNormal><span style='font-size:10.0pt'>   TCanvas *Fit1 =
02368 gROOT-&gt;GetListOfCanvases()-&gt;FindObject(&quot;Fit1&quot;);</span></p>
02369 
02370 <p class=MsoNormal><span style='font-size:10.0pt'>   if (!Fit1) Fit1 = new
02371 TCanvas(&quot;Fit1&quot;,&quot;Fit1&quot;,10,10,1000,700);</span></p>
02372 
02373 <p class=MsoNormal><span style='font-size:10.0pt'>   h-&gt;Draw(&quot;L&quot;);</span></p>
02374 
02375 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrum *s = new
02376 TSpectrum();</span></p>
02377 
02378 <p class=MsoNormal><span style='font-size:10.0pt'>   //searching for candidate
02379 peaks positions</span></p>
02380 
02381 <p class=MsoNormal><span style='font-size:10.0pt'>   nfound =
02382 s-&gt;SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);</span></p>
02383 
02384 <p class=MsoNormal><span style='font-size:10.0pt'>   Bool_t *FixPos = new
02385 Bool_t[nfound];</span></p>
02386 
02387 <p class=MsoNormal><span style='font-size:10.0pt'>   Bool_t *FixAmp = new
02388 Bool_t[nfound];      </span></p>
02389 
02390 <p class=MsoNormal><span style='font-size:10.0pt'>   for(i = 0; i&lt; nfound ;
02391 i++){</span></p>
02392 
02393 <p class=MsoNormal><span style='font-size:10.0pt'>      FixPos[i] = kFALSE;</span></p>
02394 
02395 <p class=MsoNormal><span style='font-size:10.0pt'>      FixAmp[i] = kFALSE;    </span></p>
02396 
02397 <p class=MsoNormal><span style='font-size:10.0pt'>   }</span></p>
02398 
02399 <p class=MsoNormal><span style='font-size:10.0pt'>   //filling in the initial
02400 estimates of the input parameters</span></p>
02401 
02402 <p class=MsoNormal><span style='font-size:10.0pt'>   Float_t *PosX = new
02403 Float_t[nfound];         </span></p>
02404 
02405 <p class=MsoNormal><span style='font-size:10.0pt'>   Float_t *PosY = new
02406 Float_t[nfound];</span></p>
02407 
02408 <p class=MsoNormal><span style='font-size:10.0pt'>   PosX =
02409 s-&gt;GetPositionX();</span></p>
02410 
02411 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nfound;
02412 i++) {</span></p>
02413 
02414 <p class=MsoNormal><span style='font-size:10.0pt'>                                                a=PosX[i];</span></p>
02415 
02416 <p class=MsoNormal><span style='font-size:10.0pt'>        bin = 1 + Int_t(a +
02417 0.5);</span></p>
02418 
02419 <p class=MsoNormal><span style='font-size:10.0pt'>        PosY[i] =
02420 h-&gt;GetBinContent(bin);</span></p>
02421 
02422 <p class=MsoNormal><span style='font-size:10.0pt'>   }   </span></p>
02423 
02424 <p class=MsoNormal><span style='font-size:10.0pt'>   TSpectrumFit *pfit=new
02425 TSpectrumFit(nfound);</span></p>
02426 
02427 <p class=MsoNormal><span style='font-size:10.0pt'>  
02428 pfit-&gt;SetFitParameters(xmin, xmax-1, 1000, 0.1, pfit-&gt;kFitOptimChiCounts,
02429 pfit-&gt;kFitAlphaHalving, pfit-&gt;kFitPower2,
02430 pfit-&gt;kFitTaylorOrderFirst);   </span></p>
02431 
02432 <p class=MsoNormal><span style='font-size:10.0pt'>   pfit-&gt;SetPeakParameters(2,
02433 kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);   </span></p>
02434 
02435 <p class=MsoNormal><span style='font-size:10.0pt'>  
02436 pfit-&gt;FitStiefel(source);</span></p>
02437 
02438 <p class=MsoNormal><span style='font-size:10.0pt'>   Double_t *CalcPositions =
02439 new Double_t[nfound];      </span></p>
02440 
02441 <p class=MsoNormal><span style='font-size:10.0pt'>   </span><span lang=FR
02442 style='font-size:10.0pt'>Double_t *CalcAmplitudes = new
02443 Double_t[nfound];         </span></p>
02444 
02445 <p class=MsoNormal><span lang=FR style='font-size:10.0pt'>   </span><span
02446 style='font-size:10.0pt'>CalcPositions=pfit-&gt;GetPositions();</span></p>
02447 
02448 <p class=MsoNormal><span style='font-size:10.0pt'>  
02449 CalcAmplitudes=pfit-&gt;GetAmplitudes();   </span></p>
02450 
02451 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nbins;
02452 i++) d-&gt;SetBinContent(i + 1,source[i]);</span></p>
02453 
02454 <p class=MsoNormal><span style='font-size:10.0pt'>  
02455 d-&gt;SetLineColor(kRed);   </span></p>
02456 
02457 <p class=MsoNormal><span style='font-size:10.0pt'>   d-&gt;Draw(&quot;SAME
02458 L&quot;);  </span></p>
02459 
02460 <p class=MsoNormal><span style='font-size:10.0pt'>   for (i = 0; i &lt; nfound;
02461 i++) {</span></p>
02462 
02463 <p class=MsoNormal><span style='font-size:10.0pt'>                                                a=CalcPositions[i];</span></p>
02464 
02465 <p class=MsoNormal><span style='font-size:10.0pt'>        bin = 1 + Int_t(a +
02466 0.5);                </span></p>
02467 
02468 <p class=MsoNormal><span style='font-size:10.0pt'>        PosX[i] =
02469 d-&gt;GetBinCenter(bin);</span></p>
02470 
02471 <p class=MsoNormal><span style='font-size:10.0pt'>        PosY[i] =
02472 d-&gt;GetBinContent(bin);</span></p>
02473 
02474 <p class=MsoNormal><span style='font-size:10.0pt'>   }</span></p>
02475 
02476 <p class=MsoNormal><span style='font-size:10.0pt'>   TPolyMarker * pm =
02477 (TPolyMarker*)h-&gt;GetListOfFunctions()-&gt;FindObject(&quot;TPolyMarker&quot;);</span></p>
02478 
02479 <p class=MsoNormal><span style='font-size:10.0pt'>   if (pm) {</span></p>
02480 
02481 <p class=MsoNormal><span style='font-size:10.0pt'>     
02482 h-&gt;GetListOfFunctions()-&gt;Remove(pm);</span></p>
02483 
02484 <p class=MsoNormal><span style='font-size:10.0pt'>      delete pm;</span></p>
02485 
02486 <p class=MsoNormal><span style='font-size:10.0pt'>   }</span></p>
02487 
02488 <p class=MsoNormal><span style='font-size:10.0pt'>   pm = new
02489 TPolyMarker(nfound, PosX, PosY);</span></p>
02490 
02491 <p class=MsoNormal><span style='font-size:10.0pt'>  
02492 h-&gt;GetListOfFunctions()-&gt;Add(pm);</span></p>
02493 
02494 <p class=MsoNormal><span style='font-size:10.0pt'>   pm-&gt;SetMarkerStyle(23);</span></p>
02495 
02496 <p class=MsoNormal><span style='font-size:10.0pt'>  
02497 pm-&gt;SetMarkerColor(kRed);</span></p>
02498 
02499 <p class=MsoNormal><span style='font-size:10.0pt'>   pm-&gt;SetMarkerSize(1);  
02500 </span></p>
02501 
02502 <p class=MsoNormal><span style='font-size:10.0pt'>}</span></p>
02503 
02504 </div>
02505 
02506 <!-- */
02507 // --> End_Html
02508 
02509    Int_t i, j, k, shift =
02510        2 * fNPeaks + 7, peak_vel, rozmer, iter, regul_cycle,
02511        flag;
02512    Double_t a, b, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
02513        0, pi, pmin = 0, chi_cel = 0, chi_er;
02514    Double_t *working_space = new Double_t[5 * (2 * fNPeaks + 7)];
02515    for (i = 0, j = 0; i < fNPeaks; i++) {
02516       working_space[2 * i] = fAmpInit[i];        //vector parameter
02517       if (fFixAmp[i] == false) {
02518          working_space[shift + j] = fAmpInit[i];        //vector xk
02519          j += 1;
02520       }
02521       working_space[2 * i + 1] = fPositionInit[i];        //vector parameter
02522       if (fFixPosition[i] == false) {
02523          working_space[shift + j] = fPositionInit[i];        //vector xk
02524          j += 1;
02525       }
02526    }
02527    peak_vel = 2 * i;
02528    working_space[2 * i] = fSigmaInit;        //vector parameter
02529    if (fFixSigma == false) {
02530       working_space[shift + j] = fSigmaInit;        //vector xk
02531       j += 1;
02532    }
02533    working_space[2 * i + 1] = fTInit;        //vector parameter
02534    if (fFixT == false) {
02535       working_space[shift + j] = fTInit;        //vector xk
02536       j += 1;
02537    }
02538    working_space[2 * i + 2] = fBInit;        //vector parameter
02539    if (fFixB == false) {
02540       working_space[shift + j] = fBInit;        //vector xk
02541       j += 1;
02542    }
02543    working_space[2 * i + 3] = fSInit;        //vector parameter
02544    if (fFixS == false) {
02545       working_space[shift + j] = fSInit;        //vector xk
02546       j += 1;
02547    }
02548    working_space[2 * i + 4] = fA0Init;        //vector parameter
02549    if (fFixA0 == false) {
02550       working_space[shift + j] = fA0Init;        //vector xk
02551       j += 1;
02552    }
02553    working_space[2 * i + 5] = fA1Init;        //vector parameter
02554    if (fFixA1 == false) {
02555       working_space[shift + j] = fA1Init;        //vector xk
02556       j += 1;
02557    }
02558    working_space[2 * i + 6] = fA2Init;        //vector parameter
02559    if (fFixA2 == false) {
02560       working_space[shift + j] = fA2Init;        //vector xk
02561       j += 1;
02562    }
02563    rozmer = j;
02564    if (rozmer == 0){
02565       Error ("FitAwmi","All parameters are fixed");
02566       delete [] working_space;
02567       return;    
02568    }
02569    if (rozmer >= fXmax - fXmin + 1){
02570       Error ("FitAwmi","Number of fitted parameters is larger than # of fitted points");   
02571       delete [] working_space;
02572       return;    
02573    }
02574    Double_t **working_matrix = new Double_t *[rozmer];
02575    for (i = 0; i < rozmer; i++)
02576       working_matrix[i] = new Double_t[rozmer + 4];
02577    for (iter = 0; iter < fNumberIterations; iter++) {
02578       for (j = 0; j < rozmer; j++) {
02579          working_space[3 * shift + j] = 0;        //temp
02580          for (k = 0; k <= rozmer; k++) {
02581             working_matrix[j][k] = 0;
02582          }
02583       }
02584       
02585       //filling working matrix
02586       alpha = fAlpha;
02587       chi_opt = 0;
02588       for (i = fXmin; i <= fXmax; i++) {
02589          
02590              //calculation of gradient vector
02591              for (j = 0, k = 0; j < fNPeaks; j++) {
02592             if (fFixAmp[j] == false) {
02593                working_space[2 * shift + k] =
02594                    Deramp((Double_t) i, working_space[2 * j + 1],
02595                working_space[peak_vel],
02596                working_space[peak_vel + 1],
02597                working_space[peak_vel + 3],
02598                working_space[peak_vel + 2]);
02599                k += 1;
02600             }
02601             if (fFixPosition[j] == false) {
02602                working_space[2 * shift + k] =
02603                    Deri0((Double_t) i, working_space[2 * j],
02604                working_space[2 * j + 1], working_space[peak_vel],
02605                working_space[peak_vel + 1],
02606                working_space[peak_vel + 3],
02607                working_space[peak_vel + 2]);
02608                k += 1;
02609             }
02610          } if (fFixSigma == false) {
02611             working_space[2 * shift + k] =
02612                 Dersigma(fNPeaks, (Double_t) i, working_space,
02613             working_space[peak_vel],
02614             working_space[peak_vel + 1],
02615             working_space[peak_vel + 3],
02616             working_space[peak_vel + 2]);
02617             k += 1;
02618          }
02619          if (fFixT == false) {
02620             working_space[2 * shift + k] =
02621                 Dert(fNPeaks, (Double_t) i, working_space,
02622             working_space[peak_vel], working_space[peak_vel + 2]);
02623             k += 1;
02624          }
02625          if (fFixB == false) {
02626             working_space[2 * shift + k] =
02627                 Derb(fNPeaks, (Double_t) i, working_space,
02628             working_space[peak_vel], working_space[peak_vel + 1],
02629             working_space[peak_vel + 2]);
02630             k += 1;
02631          }
02632          if (fFixS == false) {
02633             working_space[2 * shift + k] =
02634                 Ders(fNPeaks, (Double_t) i, working_space,
02635             working_space[peak_vel]);
02636             k += 1;
02637          }
02638          if (fFixA0 == false) {
02639             working_space[2 * shift + k] = 1.;
02640             k += 1;
02641          }
02642          if (fFixA1 == false) {
02643             working_space[2 * shift + k] = Dera1((Double_t) i);
02644             k += 1;
02645          }
02646          if (fFixA2 == false) {
02647             working_space[2 * shift + k] = Dera2((Double_t) i);
02648             k += 1;
02649          }
02650          yw = source[i];
02651          ywm = yw;
02652          f = Shape(fNPeaks, (Double_t) i, working_space,
02653          working_space[peak_vel], working_space[peak_vel + 1],
02654          working_space[peak_vel + 3],
02655          working_space[peak_vel + 2],
02656          working_space[peak_vel + 4],
02657          working_space[peak_vel + 5],
02658          working_space[peak_vel + 6]);
02659          if (fStatisticType == kFitOptimMaxLikelihood) {
02660             if (f > 0.00001)
02661                chi_opt += yw * TMath::Log(f) - f;
02662          }
02663          
02664          else {
02665             if (ywm != 0)
02666                chi_opt += (yw - f) * (yw - f) / ywm;
02667          }
02668          if (fStatisticType == kFitOptimChiFuncValues) {
02669             ywm = f;
02670             if (f < 0.00001)
02671                ywm = 0.00001;
02672          }
02673          
02674          else if (fStatisticType == kFitOptimMaxLikelihood) {
02675             ywm = f;
02676             if (f < 0.00001)
02677                ywm = 0.00001;
02678          }
02679          
02680          else {
02681             if (ywm == 0)
02682                ywm = 1;
02683          }
02684          for (j = 0; j < rozmer; j++) {
02685             for (k = 0; k < rozmer; k++) {
02686                b = working_space[2 * shift +
02687                                   j] * working_space[2 * shift + k] / ywm;
02688                if (fStatisticType == kFitOptimChiFuncValues)
02689                   b = b * (4 * yw - 2 * f) / ywm;
02690                working_matrix[j][k] += b;
02691                if (j == k)
02692                   working_space[3 * shift + j] += b;
02693             }
02694          }
02695          if (fStatisticType == kFitOptimChiFuncValues)
02696             b = (f * f - yw * yw) / (ywm * ywm);
02697          
02698          else
02699             b = (f - yw) / ywm;
02700          for (j = 0; j < rozmer; j++) {
02701             working_matrix[j][rozmer] -= b * working_space[2 * shift + j];
02702          }
02703       }
02704       for (i = 0; i < rozmer; i++) {
02705          working_matrix[i][rozmer + 1] = 0;        //xk
02706       }
02707       StiefelInversion(working_matrix, rozmer);
02708       for (i = 0; i < rozmer; i++) {
02709          working_space[2 * shift + i] = working_matrix[i][rozmer + 1];        //der
02710       }
02711       
02712       //calculate chi_opt
02713       chi2 = chi_opt;
02714       chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
02715       
02716       //calculate new parameters
02717       regul_cycle = 0;
02718       for (j = 0; j < rozmer; j++) {
02719          working_space[4 * shift + j] = working_space[shift + j];        //temp_xk[j]=xk[j]
02720       }
02721       
02722       do {
02723          if (fAlphaOptim == kFitAlphaOptimal) {
02724             if (fStatisticType != kFitOptimMaxLikelihood)
02725                chi_min = 10000 * chi2;
02726             
02727             else
02728                chi_min = 0.1 * chi2;
02729             flag = 0;
02730             for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
02731                for (j = 0; j < rozmer; j++) {
02732                   working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j];        //xk[j]=temp_xk[j]+pi*alpha*der[j]
02733                }
02734                for (i = 0, j = 0; i < fNPeaks; i++) {
02735                   if (fFixAmp[i] == false) {
02736                      if (working_space[shift + j] < 0)        //xk[j]
02737                         working_space[shift + j] = 0;        //xk[j]
02738                      working_space[2 * i] = working_space[shift + j];        //parameter[2*i]=xk[j]
02739                      j += 1;
02740                   }
02741                   if (fFixPosition[i] == false) {
02742                      if (working_space[shift + j] < fXmin)        //xk[j]
02743                         working_space[shift + j] = fXmin;        //xk[j]
02744                      if (working_space[shift + j] > fXmax)        //xk[j]
02745                         working_space[shift + j] = fXmax;        //xk[j]
02746                      working_space[2 * i + 1] = working_space[shift + j];        //parameter[2*i+1]=xk[j]
02747                      j += 1;
02748                   }
02749                }
02750                if (fFixSigma == false) {
02751                   if (working_space[shift + j] < 0.001) {        //xk[j]
02752                      working_space[shift + j] = 0.001;        //xk[j]
02753                   }
02754                   working_space[peak_vel] = working_space[shift + j];        //parameter[peak_vel]=xk[j]
02755                   j += 1;
02756                }
02757                if (fFixT == false) {
02758                   working_space[peak_vel + 1] = working_space[shift + j];        //parameter[peak_vel+1]=xk[j]
02759                   j += 1;
02760                }
02761                if (fFixB == false) {
02762                   if (TMath::Abs(working_space[shift + j]) < 0.001) {        //xk[j]
02763                      if (working_space[shift + j] < 0)        //xk[j]
02764                         working_space[shift + j] = -0.001;        //xk[j]
02765                      else
02766                         working_space[shift + j] = 0.001;        //xk[j]
02767                   }
02768                   working_space[peak_vel + 2] = working_space[shift + j];        //parameter[peak_vel+2]=xk[j]
02769                   j += 1;
02770                }
02771                if (fFixS == false) {
02772                   working_space[peak_vel + 3] = working_space[shift + j];        //parameter[peak_vel+3]=xk[j]
02773                   j += 1;
02774                }
02775                if (fFixA0 == false) {
02776                   working_space[peak_vel + 4] = working_space[shift + j];        //parameter[peak_vel+4]=xk[j]
02777                   j += 1;
02778                }
02779                if (fFixA1 == false) {
02780                   working_space[peak_vel + 5] = working_space[shift + j];        //parameter[peak_vel+5]=xk[j]
02781                   j += 1;
02782                }
02783                if (fFixA2 == false) {
02784                   working_space[peak_vel + 6] = working_space[shift + j];        //parameter[peak_vel+6]=xk[j]
02785                   j += 1;
02786                }
02787                chi2 = 0;
02788                for (i = fXmin; i <= fXmax; i++) {
02789                   yw = source[i];
02790                   ywm = yw;
02791                   f = Shape(fNPeaks, (Double_t) i, working_space,
02792                   working_space[peak_vel],
02793                   working_space[peak_vel + 1],
02794                   working_space[peak_vel + 3],
02795                   working_space[peak_vel + 2],
02796                   working_space[peak_vel + 4],
02797                   working_space[peak_vel + 5],
02798                   working_space[peak_vel + 6]);
02799                   if (fStatisticType == kFitOptimChiFuncValues) {
02800                      ywm = f;
02801                      if (f < 0.00001)
02802                         ywm = 0.00001;
02803                   }
02804                   if (fStatisticType == kFitOptimMaxLikelihood) {
02805                      if (f > 0.00001)
02806                         chi2 += yw * TMath::Log(f) - f;
02807                   }
02808                   
02809                   else {
02810                      if (ywm != 0)
02811                         chi2 += (yw - f) * (yw - f) / ywm;
02812                   }
02813                }
02814                if ((chi2 < chi_min
02815                     && fStatisticType != kFitOptimMaxLikelihood)
02816                     || (chi2 > chi_min
02817                     && fStatisticType == kFitOptimMaxLikelihood)) {
02818                   pmin = pi, chi_min = chi2;
02819                }
02820                
02821                else
02822                   flag = 1;
02823                if (pi == 0.1)
02824                   chi_min = chi2;
02825                chi = chi_min;
02826             }
02827             if (pmin != 0.1) {
02828                for (j = 0; j < rozmer; j++) {
02829                   working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j];        //xk[j]=temp_xk[j]+pmin*alpha*der[j]
02830                }
02831                for (i = 0, j = 0; i < fNPeaks; i++) {
02832                   if (fFixAmp[i] == false) {
02833                      if (working_space[shift + j] < 0)        //xk[j]
02834                         working_space[shift + j] = 0;        //xk[j]
02835                      working_space[2 * i] = working_space[shift + j];        //parameter[2*i]=xk[j]
02836                      j += 1;
02837                   }
02838                   if (fFixPosition[i] == false) {
02839                      if (working_space[shift + j] < fXmin)        //xk[j]
02840                         working_space[shift + j] = fXmin;        //xk[j]
02841                      if (working_space[shift + j] > fXmax)        //xk[j]
02842                         working_space[shift + j] = fXmax;        //xk[j]
02843                      working_space[2 * i + 1] = working_space[shift + j];        //parameter[2*i+1]=xk[j]
02844                      j += 1;
02845                   }
02846                }
02847                if (fFixSigma == false) {
02848                   if (working_space[shift + j] < 0.001) {        //xk[j]
02849                      working_space[shift + j] = 0.001;        //xk[j]
02850                   }
02851                   working_space[peak_vel] = working_space[shift + j];        //parameter[peak_vel]=xk[j]
02852                   j += 1;
02853                }
02854                if (fFixT == false) {
02855                   working_space[peak_vel + 1] = working_space[shift + j];        //parameter[peak_vel+1]=xk[j]
02856                   j += 1;
02857                }
02858                if (fFixB == false) {
02859                   if (TMath::Abs(working_space[shift + j]) < 0.001) {        //xk[j]
02860                      if (working_space[shift + j] < 0)        //xk[j]
02861                         working_space[shift + j] = -0.001;        //xk[j]
02862                      else
02863                         working_space[shift + j] = 0.001;        //xk[j]
02864                   }
02865                   working_space[peak_vel + 2] = working_space[shift + j];        //parameter[peak_vel+2]=xk[j]
02866                   j += 1;
02867                }
02868                if (fFixS == false) {
02869                   working_space[peak_vel + 3] = working_space[shift + j];        //parameter[peak_vel+3]=xk[j]
02870                   j += 1;
02871                }
02872                if (fFixA0 == false) {
02873                   working_space[peak_vel + 4] = working_space[shift + j];        //parameter[peak_vel+4]=xk[j]
02874                   j += 1;
02875                }
02876                if (fFixA1 == false) {
02877                   working_space[peak_vel + 5] = working_space[shift + j];        //parameter[peak_vel+5]=xk[j]
02878                   j += 1;
02879                }
02880                if (fFixA2 == false) {
02881                   working_space[peak_vel + 6] = working_space[shift + j];        //parameter[peak_vel+6]=xk[j]
02882                   j += 1;
02883                }
02884                chi = chi_min;
02885             }
02886          }
02887          
02888          else {
02889             for (j = 0; j < rozmer; j++) {
02890                working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j];        //xk[j]=temp_xk[j]+alpha*der[j]
02891             }
02892             for (i = 0, j = 0; i < fNPeaks; i++) {
02893                if (fFixAmp[i] == false) {
02894                   if (working_space[shift + j] < 0)        //xk[j]
02895                      working_space[shift + j] = 0;        //xk[j]
02896                   working_space[2 * i] = working_space[shift + j];        //parameter[2*i]=xk[j]
02897                   j += 1;
02898                }
02899                if (fFixPosition[i] == false) {
02900                   if (working_space[shift + j] < fXmin)        //xk[j]
02901                      working_space[shift + j] = fXmin;        //xk[j]
02902                   if (working_space[shift + j] > fXmax)        //xk[j]
02903                      working_space[shift + j] = fXmax;        //xk[j]
02904                   working_space[2 * i + 1] = working_space[shift + j];        //parameter[2*i+1]=xk[j]
02905                   j += 1;
02906                }
02907             }
02908             if (fFixSigma == false) {
02909                if (working_space[shift + j] < 0.001) {        //xk[j]
02910                   working_space[shift + j] = 0.001;        //xk[j]
02911                }
02912                working_space[peak_vel] = working_space[shift + j];        //parameter[peak_vel]=xk[j]
02913                j += 1;
02914             }
02915             if (fFixT == false) {
02916                working_space[peak_vel + 1] = working_space[shift + j];        //parameter[peak_vel+1]=xk[j]
02917                j += 1;
02918             }
02919             if (fFixB == false) {
02920                if (TMath::Abs(working_space[shift + j]) < 0.001) {        //xk[j]
02921                   if (working_space[shift + j] < 0)        //xk[j]
02922                      working_space[shift + j] = -0.001;        //xk[j]
02923                   else
02924                      working_space[shift + j] = 0.001;        //xk[j]
02925                }
02926                working_space[peak_vel + 2] = working_space[shift + j];        //parameter[peak_vel+2]=xk[j]
02927                j += 1;
02928             }
02929             if (fFixS == false) {
02930                working_space[peak_vel + 3] = working_space[shift + j];        //parameter[peak_vel+3]=xk[j]
02931                j += 1;
02932             }
02933             if (fFixA0 == false) {
02934                working_space[peak_vel + 4] = working_space[shift + j];        //parameter[peak_vel+4]=xk[j]
02935                j += 1;
02936             }
02937             if (fFixA1 == false) {
02938                working_space[peak_vel + 5] = working_space[shift + j];        //parameter[peak_vel+5]=xk[j]
02939                j += 1;
02940             }
02941             if (fFixA2 == false) {
02942                working_space[peak_vel + 6] = working_space[shift + j];        //parameter[peak_vel+6]=xk[j]
02943                j += 1;
02944             }
02945             chi = 0;
02946             for (i = fXmin; i <= fXmax; i++) {
02947                yw = source[i];
02948                ywm = yw;
02949                f = Shape(fNPeaks, (Double_t) i, working_space,
02950                working_space[peak_vel],
02951                working_space[peak_vel + 1],
02952                working_space[peak_vel + 3],
02953                working_space[peak_vel + 2],
02954                working_space[peak_vel + 4],
02955                working_space[peak_vel + 5],
02956                working_space[peak_vel + 6]);
02957                if (fStatisticType == kFitOptimChiFuncValues) {
02958                   ywm = f;
02959                   if (f < 0.00001)
02960                      ywm = 0.00001;
02961                }
02962                if (fStatisticType == kFitOptimMaxLikelihood) {
02963                   if (f > 0.00001)
02964                      chi += yw * TMath::Log(f) - f;
02965                }
02966                
02967                else {
02968                   if (ywm != 0)
02969                      chi += (yw - f) * (yw - f) / ywm;
02970                }
02971             }
02972          }
02973          chi2 = chi;
02974          chi = TMath::Sqrt(TMath::Abs(chi));
02975          if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
02976             alpha = alpha * chi_opt / (2 * chi);
02977          
02978          else if (fAlphaOptim == kFitAlphaOptimal)
02979             alpha = alpha / 10.0;
02980          iter += 1;
02981          regul_cycle += 1;
02982       } while (((chi > chi_opt
02983                  && fStatisticType != kFitOptimMaxLikelihood)
02984                  || (chi < chi_opt
02985                  && fStatisticType == kFitOptimMaxLikelihood))
02986                 && regul_cycle < kFitNumRegulCycles);
02987       for (j = 0; j < rozmer; j++) {
02988          working_space[4 * shift + j] = 0;        //temp_xk[j]
02989          working_space[2 * shift + j] = 0;        //der[j]
02990       }
02991       for (i = fXmin, chi_cel = 0; i <= fXmax; i++) {
02992          yw = source[i];
02993          if (yw == 0)
02994             yw = 1;
02995          f = Shape(fNPeaks, (Double_t) i, working_space,
02996          working_space[peak_vel], working_space[peak_vel + 1],
02997          working_space[peak_vel + 3],
02998          working_space[peak_vel + 2],
02999          working_space[peak_vel + 4],
03000          working_space[peak_vel + 5],
03001          working_space[peak_vel + 6]);
03002          chi_opt = (yw - f) * (yw - f) / yw;
03003          chi_cel += (yw - f) * (yw - f) / yw;
03004          
03005              //calculate gradient vector
03006              for (j = 0, k = 0; j < fNPeaks; j++) {
03007             if (fFixAmp[j] == false) {
03008                a = Deramp((Double_t) i, working_space[2 * j + 1],
03009                working_space[peak_vel],
03010                working_space[peak_vel + 1],
03011                working_space[peak_vel + 3],
03012                working_space[peak_vel + 2]);
03013                if (yw != 0) {
03014                   working_space[2 * shift + k] += chi_opt;        //der[k]
03015                   b = a * a / yw;
03016                   working_space[4 * shift + k] += b;        //temp_xk[k]
03017                }
03018                k += 1;
03019             }
03020             if (fFixPosition[j] == false) {
03021                a = Deri0((Double_t) i, working_space[2 * j],
03022                working_space[2 * j + 1],
03023                working_space[peak_vel],
03024                working_space[peak_vel + 1],
03025                working_space[peak_vel + 3],
03026                working_space[peak_vel + 2]);
03027                if (yw != 0) {
03028                   working_space[2 * shift + k] += chi_opt;        //der[k]
03029                   b = a * a / yw;
03030                   working_space[4 * shift + k] += b;        //temp_xk[k]
03031                }
03032                k += 1;
03033             }
03034          }
03035          if (fFixSigma == false) {
03036             a = Dersigma(fNPeaks, (Double_t) i, working_space,
03037             working_space[peak_vel],
03038             working_space[peak_vel + 1],
03039             working_space[peak_vel + 3],
03040            working_space[peak_vel + 2]);
03041             if (yw != 0) {
03042                working_space[2 * shift + k] += chi_opt;        //der[k]
03043                b = a * a / yw;
03044                working_space[4 * shift + k] += b;        //temp_xk[k]
03045             }
03046             k += 1;
03047          }
03048          if (fFixT == false) {
03049             a = Dert(fNPeaks, (Double_t) i, working_space,
03050             working_space[peak_vel],
03051             working_space[peak_vel + 2]);
03052             if (yw != 0) {
03053                working_space[2 * shift + k] += chi_opt;        //der[k]
03054                b = a * a / yw;
03055                working_space[4 * shift + k] += b;        //temp_xk[k]
03056             }
03057             k += 1;
03058          }
03059          if (fFixB == false) {
03060             a = Derb(fNPeaks, (Double_t) i, working_space,
03061             working_space[peak_vel], working_space[peak_vel + 1],
03062             working_space[peak_vel + 2]);
03063             if (yw != 0) {
03064                working_space[2 * shift + k] += chi_opt;        //der[k]
03065                b = a * a / yw;
03066                working_space[4 * shift + k] += b;        //temp_xk[k]
03067             }
03068             k += 1;
03069          }
03070          if (fFixS == false) {
03071             a = Ders(fNPeaks, (Double_t) i, working_space,
03072             working_space[peak_vel]);
03073             if (yw != 0) {
03074                working_space[2 * shift + k] += chi_opt;        //der[k]
03075                b = a * a / yw;
03076                working_space[4 * shift + k] += b;        //tem_xk[k]
03077             }
03078             k += 1;
03079          }
03080          if (fFixA0 == false) {
03081             a = 1.0;
03082             if (yw != 0) {
03083                working_space[2 * shift + k] += chi_opt;        //der[k]
03084                b = a * a / yw;
03085                working_space[4 * shift + k] += b;        //temp_xk[k]
03086             }
03087             k += 1;
03088          }
03089          if (fFixA1 == false) {
03090             a = Dera1((Double_t) i);
03091             if (yw != 0) {
03092                working_space[2 * shift + k] += chi_opt;        //der[k]
03093                b = a * a / yw;
03094                working_space[4 * shift + k] += b;        //temp_xk[k]
03095             }
03096             k += 1;
03097          }
03098          if (fFixA2 == false) {
03099             a = Dera2((Double_t) i);
03100             if (yw != 0) {
03101                working_space[2 * shift + k] += chi_opt;        //der[k]
03102                b = a * a / yw;
03103                working_space[4 * shift + k] += b;        //temp_xk[k]
03104             }
03105             k += 1;
03106          }
03107       }
03108    }
03109    b = fXmax - fXmin + 1 - rozmer;
03110    chi_er = chi_cel / b;
03111    for (i = 0, j = 0; i < fNPeaks; i++) {
03112       fArea[i] =
03113           Area(working_space[2 * i], working_space[peak_vel],
03114                working_space[peak_vel + 1], working_space[peak_vel + 2]);
03115       if (fFixAmp[i] == false) {
03116          fAmpCalc[i] = working_space[shift + j];        //xk[j]
03117          if (working_space[3 * shift + j] != 0)
03118             fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
03119          if (fArea[i] > 0) {
03120             a = Derpa(working_space[peak_vel],
03121             working_space[peak_vel + 1],
03122             working_space[peak_vel + 2]);
03123             b = working_space[4 * shift + j];        //temp_xk[j]
03124             if (b == 0)
03125                b = 1;
03126             
03127             else
03128                b = 1 / b;
03129             fAreaErr[i] = TMath::Sqrt(TMath::Abs(a * a * b * chi_er));
03130          }
03131          
03132          else
03133             fAreaErr[i] = 0;
03134          j += 1;
03135       }
03136       
03137       else {
03138          fAmpCalc[i] = fAmpInit[i];
03139          fAmpErr[i] = 0;
03140          fAreaErr[i] = 0;
03141       }
03142       if (fFixPosition[i] == false) {
03143          fPositionCalc[i] = working_space[shift + j];        //xk[j]
03144          if (working_space[3 * shift + j] != 0)        //temp[j]
03145             fPositionErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //Der[j]/temp[j]
03146          j += 1;
03147       }
03148       
03149       else {
03150          fPositionCalc[i] = fPositionInit[i];
03151          fPositionErr[i] = 0;
03152       }
03153    }
03154    if (fFixSigma == false) {
03155       fSigmaCalc = working_space[shift + j];        //xk[j]
03156       if (working_space[3 * shift + j] != 0)        //temp[j]
03157          fSigmaErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
03158       j += 1;
03159    }
03160    
03161    else {
03162       fSigmaCalc = fSigmaInit;
03163       fSigmaErr = 0;
03164    }
03165    if (fFixT == false) {
03166       fTCalc = working_space[shift + j];        //xk[j]
03167       if (working_space[3 * shift + j] != 0)        //temp[j]
03168          fTErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
03169       j += 1;
03170    }
03171    
03172    else {
03173       fTCalc = fTInit;
03174       fTErr = 0;
03175    }
03176    if (fFixB == false) {
03177       fBCalc = working_space[shift + j];        //xk[j]
03178       if (working_space[3 * shift + j] != 0)        //temp[j]
03179          fBErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
03180       j += 1;
03181    }
03182    
03183    else {
03184       fBCalc = fBInit;
03185       fBErr = 0;
03186    }
03187    if (fFixS == false) {
03188       fSCalc = working_space[shift + j];        //xk[j]
03189       if (working_space[3 * shift + j] != 0)        //temp[j]
03190          fSErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
03191       j += 1;
03192    }
03193    
03194    else {
03195       fSCalc = fSInit;
03196       fSErr = 0;
03197    }
03198    if (fFixA0 == false) {
03199       fA0Calc = working_space[shift + j];        //xk[j]
03200       if (working_space[3 * shift + j] != 0)        //temp[j]
03201          fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
03202       j += 1;
03203    }
03204    
03205    else {
03206       fA0Calc = fA0Init;
03207       fA0Err = 0;
03208    }
03209    if (fFixA1 == false) {
03210       fA1Calc = working_space[shift + j];        //xk[j]
03211       if (working_space[3 * shift + j] != 0)        //temp[j]
03212          fA1Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
03213       j += 1;
03214    }
03215    
03216    else {
03217       fA1Calc = fA1Init;
03218       fA1Err = 0;
03219    }
03220    if (fFixA2 == false) {
03221       fA2Calc = working_space[shift + j];        //xk[j]
03222       if (working_space[3 * shift + j] != 0)        //temp[j]
03223          fA2Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j]));        //der[j]/temp[j]
03224       j += 1;
03225    }
03226    
03227    else {
03228       fA2Calc = fA2Init;
03229       fA2Err = 0;
03230    }
03231    b = fXmax - fXmin + 1 - rozmer;
03232    fChi = chi_cel / b;
03233    for (i = fXmin; i <= fXmax; i++) {
03234       f = Shape(fNPeaks, (Double_t) i, working_space,
03235       working_space[peak_vel], working_space[peak_vel + 1],
03236       working_space[peak_vel + 3], working_space[peak_vel + 2],
03237       working_space[peak_vel + 4], working_space[peak_vel + 5],
03238       working_space[peak_vel + 6]);
03239       source[i] = f;
03240    } 
03241    for (i = 0; i < rozmer; i++)
03242       delete [] working_matrix[i];
03243    delete [] working_matrix;
03244    delete [] working_space;
03245    return;
03246 }
03247 
03248 //____________________________________________________________________________
03249 void TSpectrumFit::SetFitParameters(Int_t xmin,Int_t xmax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor)
03250 {
03251 //////////////////////////////////////////////////////////////////////////////
03252 //   SETTER FUNCTION                                                      
03253 //                                                     
03254 //   This function sets the following fitting parameters:
03255 //         -xmin, xmax - fitting region
03256 //         -numberIterations - # of desired iterations in the fit
03257 //         -alpha - convergence coefficient, it should be positive number and <=1, for details see references
03258 //         -statisticType - type of statistics, possible values kFitOptimChiCounts (chi square statistics with counts as weighting coefficients), kFitOptimChiFuncValues (chi square statistics with function values as weighting coefficients),kFitOptimMaxLikelihood 
03259 //         -alphaOptim - optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal
03260 //         -power - possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting function.
03261 //         -fitTaylor - order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi fitting function.
03262 //////////////////////////////////////////////////////////////////////////////   
03263    if(xmin<0 || xmax <= xmin){ 
03264       Error("SetFitParameters", "Wrong range");
03265       return;
03266    }      
03267    if (numberIterations <= 0){
03268       Error("SetFitParameters","Invalid number of iterations, must be positive");
03269       return;
03270    }      
03271    if (alpha <= 0 || alpha > 1){
03272       Error ("SetFitParameters","Invalid step coefficient alpha, must be > than 0 and <=1");
03273       return;
03274    }      
03275    if (statisticType != kFitOptimChiCounts
03276         && statisticType != kFitOptimChiFuncValues
03277         && statisticType != kFitOptimMaxLikelihood){
03278       Error("SetFitParameters","Wrong type of statistic");
03279       return;
03280    }      
03281    if (alphaOptim != kFitAlphaHalving
03282         && alphaOptim != kFitAlphaOptimal){
03283       Error("SetFitParameters","Wrong optimization algorithm");
03284       return;
03285    }      
03286    if (power != kFitPower2 && power != kFitPower4
03287         && power != kFitPower6 && power != kFitPower8
03288         && power != kFitPower10 && power != kFitPower12){
03289       Error("SetFitParameters","Wrong power");
03290       return;
03291    }      
03292    if (fitTaylor != kFitTaylorOrderFirst
03293         && fitTaylor != kFitTaylorOrderSecond){
03294       Error("SetFitParameters","Wrong order of Taylor development");
03295       return;
03296    }      
03297    fXmin=xmin,fXmax=xmax,fNumberIterations=numberIterations,fAlpha=alpha,fStatisticType=statisticType,fAlphaOptim=alphaOptim,fPower=power,fFitTaylor=fitTaylor;
03298 }
03299 
03300 //_______________________________________________________________________________
03301 void TSpectrumFit::SetPeakParameters(Double_t sigma, Bool_t fixSigma, const Float_t *positionInit, const Bool_t *fixPosition, const Float_t *ampInit, const Bool_t *fixAmp)
03302 {
03303 //////////////////////////////////////////////////////////////////////////////
03304 //   SETTER FUNCTION                                                      
03305 //                                                     
03306 //   This function sets the following fitting parameters of peaks:
03307 //         -sigma - initial value of sigma parameter
03308 //         -fixSigma - logical value of sigma parameter, which allows to fix the parameter (not to fit)   
03309 //         -positionInit - aray of initial values of peaks positions
03310 //         -fixPosition - array of logical values which allow to fix appropriate positions (not fit). However they are present in the estimated functional.
03311 //         -ampInit - aray of initial values of peaks amplitudes
03312 //         -fixAmp - aray of logical values which allow to fix appropriate amplitudes (not fit). However they are present in the estimated functional      
03313 //////////////////////////////////////////////////////////////////////////////   
03314    
03315    Int_t i;
03316    if (sigma <= 0){
03317       Error ("SetPeakParameters","Invalid sigma, must be > than 0");
03318       return;
03319    }      
03320    for(i=0; i < fNPeaks; i++){
03321       if((Int_t) positionInit[i] < fXmin || (Int_t) positionInit[i] > fXmax){
03322          Error ("SetPeakParameters","Invalid peak position, must be in the range fXmin, fXmax");
03323          return;
03324       }         
03325       if(ampInit[i] < 0){
03326          Error ("SetPeakParameters","Invalid peak amplitude, must be > than 0");        
03327          return;
03328       }         
03329    }
03330    fSigmaInit = sigma,fFixSigma = fixSigma;
03331    for(i=0; i < fNPeaks; i++){
03332       fPositionInit[i] = (Double_t) positionInit[i];
03333       fFixPosition[i] = fixPosition[i];
03334       fAmpInit[i] = (Double_t) ampInit[i];
03335       fFixAmp[i] = fixAmp[i]; 
03336    }
03337 }
03338 
03339 //_______________________________________________________________________________
03340 void TSpectrumFit::SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t a1Init, Bool_t fixA1, Double_t a2Init, Bool_t fixA2)
03341 {
03342 //////////////////////////////////////////////////////////////////////////////
03343 //   SETTER FUNCTION                                                      
03344 //                                                     
03345 //   This function sets the following fitting parameters of background:
03346 //         -a0Init - initial value of a0 parameter (backgroud is estimated as a0+a1*x+a2*x*x)
03347 //         -fixA0 - logical value of a0 parameter, which allows to fix the parameter (not to fit)  
03348 //         -a1Init - initial value of a1 parameter
03349 //         -fixA1 - logical value of a1 parameter, which allows to fix the parameter (not to fit)   
03350 //         -a2Init - initial value of a2 parameter
03351 //         -fixA2 - logical value of a2 parameter, which allows to fix the parameter (not to fit)     
03352 //////////////////////////////////////////////////////////////////////////////   
03353    
03354    fA0Init = a0Init;
03355    fFixA0 = fixA0;
03356    fA1Init = a1Init;
03357    fFixA1 = fixA1;
03358    fA2Init = a2Init;
03359    fFixA2 = fixA2;    
03360 }
03361 
03362 //_______________________________________________________________________________
03363 void TSpectrumFit::SetTailParameters(Double_t tInit, Bool_t fixT, Double_t bInit, Bool_t fixB, Double_t sInit, Bool_t fixS)
03364 {
03365 //////////////////////////////////////////////////////////////////////////////
03366 //   SETTER FUNCTION                                                      
03367 //                                                     
03368 //   This function sets the following fitting parameters of tails of peaks
03369 //         -tInit - initial value of t parameter
03370 //         -fixT - logical value of t parameter, which allows to fix the parameter (not to fit)  
03371 //         -bInit - initial value of b parameter
03372 //         -fixB - logical value of b parameter, which allows to fix the parameter (not to fit)  
03373 //         -sInit - initial value of s parameter
03374 //         -fixS - logical value of s parameter, which allows to fix the parameter (not to fit)  
03375 //////////////////////////////////////////////////////////////////////////////   
03376    
03377    fTInit = tInit;
03378    fFixT = fixT;
03379    fBInit = bInit;
03380    fFixB = fixB;
03381    fSInit = sInit;
03382    fFixS = fixS;     
03383 }
03384 
03385 //_______________________________________________________________________________
03386 void TSpectrumFit::GetSigma(Double_t &sigma, Double_t &sigmaErr)
03387 {
03388 //////////////////////////////////////////////////////////////////////////////
03389 //   GETTER FUNCTION                                                      
03390 //                                                     
03391 //   This function gets the sigma parameter and its error
03392 //         -sigma - gets the fitted value of sigma parameter
03393 //         -sigmaErr - gets error value of sigma parameter
03394 //////////////////////////////////////////////////////////////////////////////   
03395    sigma=fSigmaCalc;
03396    sigmaErr=fSigmaErr;
03397 }
03398 
03399 //_______________________________________________________________________________
03400 void TSpectrumFit::GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &a1, Double_t &a1Err, Double_t &a2, Double_t &a2Err)
03401 {
03402 //////////////////////////////////////////////////////////////////////////////
03403 //   GETTER FUNCTION                                                      
03404 //                                                     
03405 //   This function gets the background parameters and their errors
03406 //         -a0 - gets the fitted value of a0 parameter
03407 //         -a0Err - gets error value of a0 parameter
03408 //         -a1 - gets the fitted value of a1 parameter
03409 //         -a1Err - gets error value of a1 parameter
03410 //         -a2 - gets the fitted value of a2 parameter
03411 //         -a2Err - gets error value of a2 parameter
03412 //////////////////////////////////////////////////////////////////////////////      
03413    a0 = fA0Calc;
03414    a0Err = fA0Err;
03415    a1 = fA1Calc;
03416    a1Err = fA1Err;
03417    a2 = fA2Calc;
03418    a2Err = fA2Err;    
03419 }
03420 
03421 //_______________________________________________________________________________
03422 void TSpectrumFit::GetTailParameters(Double_t &t, Double_t &tErr, Double_t &b, Double_t &bErr, Double_t &s, Double_t &sErr)
03423 {
03424 //////////////////////////////////////////////////////////////////////////////
03425 //   GETTER FUNCTION                                                      
03426 //                                                     
03427 //   This function gets the tail parameters and their errors
03428 //         -t - gets the fitted value of t parameter
03429 //         -tErr - gets error value of t parameter
03430 //         -b - gets the fitted value of b parameter
03431 //         -bErr - gets error value of b parameter
03432 //         -s - gets the fitted value of s parameter
03433 //         -sErr - gets error value of s parameter
03434 //////////////////////////////////////////////////////////////////////////////         
03435    t = fTCalc;
03436    tErr = fTErr;
03437    b = fBCalc;
03438    bErr = fBErr;  
03439    s = fSCalc;
03440    sErr = fSErr;  
03441 }
03442