FittingFunctionGaussianWithBackground.cpp

Go to the documentation of this file.
00001 // FittingFunctionGaussianWithBackground.cpp: implementation of the FittingFunctionGaussianWithBackground class.
00002 //
00004 #ifdef WIN32 
00005 #pragma warning(disable:4786)
00006 #endif
00007 
00008 #include "FittingFunctionGaussianWithBackground.h"
00009 
00010 
00012 // Construction/Destruction
00014 
00015 std::string     FittingFunctionGaussianWithBackground::getEquation() const
00016 {
00017         return "f(x) = Height * exp[ - (x - Position)˛ / (2 * Width˛)] + Background";
00018 }
00019 
00020 
00021 //##ModelId=43849A020392
00022 FittingFunctionGaussianWithBackground::FittingFunctionGaussianWithBackground()
00023 {
00024         _iNbParameters                          =       4;
00025         _lNbMaxIteration                        =       NB_ITERATION_MAX;
00026         _dSearchLimit                           =       SEARCH_LIMIT;
00027         _pParameters                            =       new double[_iNbParameters];
00028         _pParametersErrors                      =       new double[_iNbParameters];
00029         _mFunctionMap["Position"]       =       0;
00030         _mFunctionMap["Width"]          =       1;
00031         _mFunctionMap["Height"]         =       2;
00032         _mFunctionMap["Background"]     =       3;
00033         _vInitialGuess = gsl_vector_alloc(4);
00034 }
00035 
00036 //##ModelId=43849A02039F
00037 FittingFunctionGaussianWithBackground::FittingFunctionGaussianWithBackground(long lNbIteration,double dLimitSearch) 
00038 {
00039         _iNbParameters                          =       4;
00040         _lNbMaxIteration                        =       lNbIteration;
00041         _dSearchLimit                           =       dLimitSearch;
00042         _pParameters                            =       new double[_iNbParameters];
00043         _pParametersErrors                      =       new double[_iNbParameters];
00044         _mFunctionMap["Position"]       =       0;
00045         _mFunctionMap["Width"]          =       1;
00046         _mFunctionMap["Height"]         =       2;
00047         _mFunctionMap["Background"]     =       3;
00048         _vInitialGuess = gsl_vector_alloc(4);
00049 }
00050 
00051 //##ModelId=43849A0203A2
00052 FittingFunctionGaussianWithBackground::~FittingFunctionGaussianWithBackground()
00053 {
00054         printf ("Destruction FittingFunctionGaussianWithBackground !!!\n");
00055         delete[] _pParameters;
00056         delete[] _pParametersErrors;
00057         _pParameters = NULL;
00058         _pParametersErrors = NULL;
00059         gsl_vector_free(_vInitialGuess);
00060 }
00061 
00062 //##ModelId=43849A0203A4
00063 PtrFunction FittingFunctionGaussianWithBackground::getFunction() const{
00064         return FittingFunctionGaussianWithBackground::Function;
00065 }
00066 
00067 //##ModelId=43849A0203A6
00068 PtrFunctionDerivate FittingFunctionGaussianWithBackground::getFunctionDerivate() const{
00069         return FittingFunctionGaussianWithBackground::FunctionDerivate;
00070 }
00071 
00072 //##ModelId=43849A0203A8
00073 PtrFunctionAndDerivate FittingFunctionGaussianWithBackground::getFunctionAndDerivate() const {
00074         return FittingFunctionGaussianWithBackground::FunctionAndDerivate;
00075 }
00076 
00077 
00082 //##ModelId=43849A0203AF
00083 int FittingFunctionGaussianWithBackground::Function(const gsl_vector *peak_params, void *calc_data, gsl_vector *func)
00084 {
00085         size_t i;
00086         size_t iNbPoints        = ((struct FittingData *)calc_data)->n;
00087         double *x                       = ((struct FittingData *)calc_data)->x;
00088         double *y                       = ((struct FittingData *)calc_data)->y;
00089         double *sigma           = ((struct FittingData *)calc_data)->sigma;
00090 
00091         double pos                      = gsl_vector_get (peak_params, 0);
00092         double width            = gsl_vector_get (peak_params, 1);
00093         double height           = gsl_vector_get (peak_params, 2);
00094         double background       = gsl_vector_get (peak_params, 3);
00095 
00096 //#ifdef DEBUG_FIT      
00097 //  printf("pos = %g, width = %g, height = %g, background = %g\n", pos, width, height,background);
00098 //#endif        
00099 
00100 double dXPos,dXPos2,dWidth2,dYi;
00101 
00102         for (i=0; i < iNbPoints; i++)
00103         {       
00104                 dXPos   = x[i] - pos;
00105                 dXPos2  = dXPos * dXPos;
00106                 dWidth2 = width * width;
00107         
00108                 dYi = background + height * exp(-0.5 * dXPos2 / dWidth2);
00109         
00110                 //gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);
00111                 gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);                              
00112     }
00113 
00114   return GSL_SUCCESS;
00115 }
00116 
00129 //##ModelId=43849A0203B4
00130 int FittingFunctionGaussianWithBackground::FunctionDerivate(const gsl_vector *peak_params, void *calc_data, gsl_matrix *Jacobian)
00131 {
00132 
00133         size_t i;
00134         
00135         size_t  n                       = ((struct FittingData *)calc_data)->n;
00136         double *x                       = ((struct FittingData *)calc_data)->x;
00137         double *sigma           = ((struct FittingData *)calc_data)->sigma;
00138 
00139         double pos                      = gsl_vector_get (peak_params, 0);
00140         double width            = gsl_vector_get (peak_params, 1);
00141         double height           = gsl_vector_get (peak_params, 2);
00142         //double background     = gsl_vector_get (peak_params, 3);
00143 
00144         double d_pos, d_width, d_height, d_background;
00145         double s,xpos2,width2,e;
00146         // Fill the Jacobian Matrix
00147         for (i=0; i< n; i++)
00148     {
00149                 s                               = sigma[i];
00150                 xpos2                   = (x[i]-pos)*(x[i]-pos);
00151                 width2                  = width * width;
00152                 e                               = exp( -0.5 * xpos2 / width2);          
00153                  
00154                 d_pos                   = (height/width2 * (x[i] - pos) * e) / s;       // derivate dYi/dpos    
00155                 d_width                 = (height/width2/width * xpos2 * e) / s;        // derivate dYi/dwidth 
00156                 d_height                = e/s;                                                                          // derivate dYi/dheight                  
00157                 d_background    = 1/s;                                                                          // derivate dYi/dbackground
00158 
00159                 // fill the jacobian matrix J(points, parameters) 
00160         gsl_matrix_set (Jacobian, i, 0, d_pos); 
00161                 gsl_matrix_set (Jacobian, i, 1, d_width);
00162                 gsl_matrix_set (Jacobian, i, 2, d_height);
00163                 gsl_matrix_set (Jacobian, i, 3, d_background);
00164     }
00165          
00166         return GSL_SUCCESS;
00167 }
00168 
00169 
00170 // Interface function for the gsl_multifit_fdfsolver
00171 //##ModelId=43849A0203B9
00172 int FittingFunctionGaussianWithBackground::FunctionAndDerivate(const gsl_vector *peak_params, void *calc_data,gsl_vector *func, gsl_matrix *Jacobian) 
00173 {       
00174         Function  (peak_params, calc_data, func);
00175         FunctionDerivate (peak_params, calc_data, Jacobian);
00176         return GSL_SUCCESS;
00177 }
00178 
00179 //The size of the 3 gsl_vector must be the same ! Not checked here !!!
00180 //##ModelId=43849A0203C4
00181 int FittingFunctionGaussianWithBackground::doFit(long points, gsl_vector *x, gsl_vector *y,gsl_vector *sigma,gsl_vector *guess,  bool bEstimateGuess,int iSearchStoppingMethod,bool bUseScaled)
00182 {
00183 
00184 initializeInitialsParameters(x,y,guess,bEstimateGuess);
00185 
00186         long            i;
00187         long            iter;
00188         long            status;
00189         
00190         double* datax = new double[points];
00191         double* datay = new double[points];
00192         double* datas = new double[points];
00193                 
00194         // fill the data arrays 
00195         for (i=0; i<points; i++)
00196         {               
00197                 datax[i] = gsl_vector_get (x, i);               // x = index            
00198                 datay[i] = gsl_vector_get (y, i);               // y = data to be fit           
00199                 datas[i] = gsl_vector_get (sigma, i);   // s = sigma (weight) of data values 
00200         }
00201         
00202         //gsl data structure    
00203         struct FittingData gsl_data;
00204         gsl_data.n      = points;
00205         gsl_data.x      = datax;
00206         gsl_data.y              = datay;
00207         gsl_data.sigma  = datas;
00208 
00209         cout << "\n---------------------------------------------" << endl;
00210         cout << "Initials Parameters For The Fit" << endl;
00211         // set-up the solver 
00212         const gsl_multifit_fdfsolver_type *T;
00213         if (bUseScaled)
00214         {
00215                 T       = gsl_multifit_fdfsolver_lmsder;
00216                 cout << "Jacobian Data Scaled" << endl;
00217         }
00218         else
00219         {
00220                 T       = gsl_multifit_fdfsolver_lmder;
00221                 cout << "Jacobian Data Not Scaled" << endl;
00222         }
00223         gsl_multifit_fdfsolver *s                               = gsl_multifit_fdfsolver_alloc (T, points, _iNbParameters);
00224         gsl_multifit_function_fdf f;
00225 
00226         f.f             = getFunction();              // gauss_f;
00227         f.df            = getFunctionDerivate();      // gauss_df;
00228         f.fdf           = getFunctionAndDerivate();   // gauss_fdf;
00229         f.n             = points;
00230         f.p             = _iNbParameters;
00231         f.params        = &gsl_data;
00232   
00233         gsl_multifit_fdfsolver_set ( s,  &f,  _vInitialGuess );
00234 
00235         switch(iSearchStoppingMethod)
00236         {
00237                 case 1 :
00238                         cout << "Test Delta Search Stopping Function Used" << endl;
00239                         break;
00240 
00241                 case 2 :
00242                         cout << "Test Gradient Search Stopping Function Used" << endl;
00243                         break;
00244         }
00245 
00246         cout << "\n----------------------------------------------" << endl;
00247         cout << "Initials Values For The Fit " << endl;
00248         cout << "Position   = " << gsl_vector_get (_vInitialGuess, 0) << endl;
00249         cout << "Width      = " << gsl_vector_get (_vInitialGuess ,1) << endl;
00250         cout << "Height     = " << gsl_vector_get (_vInitialGuess, 2) << endl;
00251         cout << "Background = " << gsl_vector_get (_vInitialGuess, 3) << endl;
00252 
00253 
00254                         cout << "\n---------------------------------------------" << endl;
00255                         cout << "Iterations" << endl;
00256         iter = 0;
00257         #ifdef DEBUG_FIT
00258           printState (iter, s);
00259         #endif
00260 
00261         do
00262     {
00263                 iter++;
00264                 status = gsl_multifit_fdfsolver_iterate (s);
00265 
00266                 #ifdef DEBUG_FIT
00267                           printf ("status = %s\n", gsl_strerror (status));
00268                           printState (iter, s);
00269                 #endif
00270                 
00271                 if (status)
00272                 {
00273                 break;
00274                 }
00275 
00276 
00277                 //
00278                 switch(iSearchStoppingMethod)
00279                 {
00280                 case 1 :
00281                         status = gsl_multifit_test_delta (s->dx, s->x,  /*_dSearchLimit*/0.0, _dSearchLimit);
00282                         break;
00283 
00284                 case 2 :
00285                         gsl_vector* gradient = gsl_vector_alloc(_iNbParameters);
00286                         gsl_multifit_gradient(s->J,s->f,gradient);
00287                         status = gsl_multifit_test_gradient (gradient, _dSearchLimit);
00288                         gsl_vector_free(gradient);
00289                         break;
00290                 }
00291 
00292     }
00293         while (status == GSL_CONTINUE && iter < _lNbMaxIteration);
00294         _lNbIterations = iter;
00295         // Check whether the fitting has converged or failed  
00296         //if ( status == GSL_SUCCESS )
00297         //      {
00298         gsl_matrix *covar = gsl_matrix_alloc (_iNbParameters, _iNbParameters);
00299         gsl_multifit_covar (s->J, 0.0, covar);
00300         double chi                      = gsl_blas_dnrm2(s->f);
00301 
00302         _pParameters[0]                 =       FIT(0);
00303         _pParameters[1]                 =       FIT(1);
00304         _pParameters[2]                 =       FIT(2);
00305         _pParameters[3]                 =       FIT(3);
00306 
00307         _pParametersErrors[0]   =       ERR(0);
00308         _pParametersErrors[1]   =       ERR(1);
00309         _pParametersErrors[2]   =       ERR(2);
00310         _pParametersErrors[3]   =       ERR(3);
00311 
00312         //double dHeight = _pParameters[2];
00313 
00314         double dEstimatedWidth  =       _pParameters[1];
00315         _dQualityFactor                 =       pow(chi, 2.0)/ (points - _iNbParameters);
00316 //      _dSigma                                 =       1.0/ (sqrt(2*M_PI)*dHeight);
00317 //      _dFHWM                                  =       2.355*_dSigma;
00318         _dHWHM                                  =       sqrt(2*log(2.0))*dEstimatedWidth;               
00319         _dFHWM                                  =       2*_dHWHM;               
00320         cout << "2*sqrt(2*log(2.0))=" << 2*sqrt(2*log(2.0)) <<endl;
00321 
00322         #ifdef DEBUG_FIT        
00323                         cout << "\n---------------------------------------------" << endl;
00324                         cout << "Results For Gaussian With Background" << endl;
00325                         printResults();
00326         #endif          
00327 
00328         gsl_matrix_free(covar);         
00329         //      }
00330 
00331         #ifdef DEBUG_FIT  
00332                 printf ("status = %s\n", gsl_strerror (status));
00333         #endif
00334         
00335         // Free data
00336         gsl_multifit_fdfsolver_free (s);
00337         delete[] datax;
00338         delete[] datay;
00339         delete[] datas;
00340         
00341         datax = NULL;
00342         datax = NULL;
00343         datax = NULL;
00344 
00345         return status;
00346 }
00347 
00348 
00349 
00350 void    FittingFunctionGaussianWithBackground::estimateInitialGuess(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY)
00351 {
00352         // fill the initials parameters values estimation vector for the initialisation
00353         // [0] = position 
00354         // [1] = width  
00355         // [2] = height
00356         // [3] = background
00357         // Thanks to the max index of all the points, we can get the
00358         // position_value       = x_value at index_max
00359         // width                        = standard deviation
00360         // height_value         = y_value - min_value (=background)     at index_max
00361         // background_value = minimal value
00362         //min_index = gsl_vector_min_index (vExperimentalDataY);
00363         //max_index = gsl_vector_max_index (vExperimentalDataY);
00364         double dInitialPosition = estimateInitialPosition(vExperimentalDataX,vExperimentalDataY);
00365         double dInitialSigma    = estimateInitialSigma(vExperimentalDataX,vExperimentalDataY);
00366         double dInitialHeight   = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00367         double dInitialBackground = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00368 
00369         gsl_vector_set(_vInitialGuess, 0, dInitialPosition);
00370         gsl_vector_set(_vInitialGuess, 1, dInitialSigma);
00371         gsl_vector_set(_vInitialGuess, 2, dInitialHeight);
00372         gsl_vector_set(_vInitialGuess, 3, dInitialBackground);  
00373 }
00374 
00375 
00376 void    FittingFunctionGaussianWithBackground::initializeInitialsParameters(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY,gsl_vector* vInitialGuess,bool bEstimate)
00377 {
00378         if (bEstimate) 
00379         {
00380                 estimateInitialGuess(vExperimentalDataX,vExperimentalDataY);
00381         }
00382         else
00383         {
00384                 gsl_vector_set(_vInitialGuess,0,gsl_vector_get(vInitialGuess,0));
00385                 gsl_vector_set(_vInitialGuess,1,gsl_vector_get(vInitialGuess,1));
00386                 gsl_vector_set(_vInitialGuess,2,gsl_vector_get(vInitialGuess,2));
00387                 gsl_vector_set(_vInitialGuess,3,gsl_vector_get(vInitialGuess,3));
00388         }
00389 }
00390 
00391 
00392 
00393 
00394 
00395 //##ModelId=43849A0203D0
00396 int FittingFunctionGaussianWithBackground::generateFunctionFit(double dValMinX,double dResolutionX,long lNbPoints,gsl_vector* funcX,gsl_vector* funcY) {
00397 
00398         int i;
00399         double dPos                     =       _pParameters[0];
00400         double dWidth           =       _pParameters[1];
00401         double dHeight          =       _pParameters[2];
00402         double dBackground      =       _pParameters[3];
00403         double dXi                      =       0;
00404         double dXPos            =       0;
00405         double dXPos2           =       0;
00406         double dWidth2          =       1;
00407         double dYi                      =       0;
00408 
00409         for (i=0;i<lNbPoints;i++) 
00410         {
00411                 dXi             = dValMinX+i*dResolutionX;
00412                 dXPos   = dXi - dPos;
00413                 dXPos2  = dXPos * dXPos;
00414                 dWidth2 = dWidth * dWidth;
00415                 dYi             = dBackground + dHeight * exp(-0.5 * dXPos2 / dWidth2);
00416         
00417                 gsl_vector_set(funcX,i,dXi);
00418                 gsl_vector_set(funcY,i,dYi);
00419         }
00420         return 1;
00421 }
00422 
00423 
00424 int FittingFunctionGaussianWithBackground::generateFunctionFit(long lNbPoints, gsl_vector *x, gsl_vector* funcX,gsl_vector* funcY)
00425 {
00426 
00427         int i;
00428         double dPos                     =       _pParameters[0];
00429         double dWidth           =       _pParameters[1];
00430         double dHeight          =       _pParameters[2];
00431         double dBackground      =       _pParameters[3];
00432         double dXi                      =       0;
00433         double dXPos            =       0;
00434         double dXPos2           =       0;
00435         double dWidth2          =       1;
00436         double dYi                      =       0;
00437 
00438         for (i=0;i<lNbPoints;i++) 
00439         {
00440                 dXi             = gsl_vector_get(x,i);
00441                 dXPos   = dXi - dPos;
00442                 dXPos2  = dXPos * dXPos;
00443                 dWidth2 = dWidth * dWidth;
00444                 dYi             = dBackground + dHeight * exp(-0.5 * dXPos2 / dWidth2);
00445         
00446                 gsl_vector_set(funcX,i,dXi);
00447                 gsl_vector_set(funcY,i,dYi);
00448         }
00449         return 1;
00450 }
00451 
00452 
00453 
00454 // print the state of the fitting during iterations 
00455 //##ModelId=43849A0203D7
00456 void FittingFunctionGaussianWithBackground::printState (size_t iter, gsl_multifit_fdfsolver * s)
00457 {
00458   printf ("iter: %3u x : %g %g %g %g|f(x)|=%g\n",
00459           iter,
00460           gsl_vector_get (s->x, 0), 
00461           gsl_vector_get (s->x, 1),
00462           gsl_vector_get (s->x, 2),
00463                   gsl_vector_get (s->x, 3),
00464           gsl_blas_dnrm2 (s->f));
00465 }
00466 
00467 //##ModelId=43849A0203E0
00468 void FittingFunctionGaussianWithBackground::printResults()
00469 {
00470                 printf ("pos            = %g \t+/- %g\n",       getParameter("Position")        ,_pParametersErrors[0]);
00471                 printf ("width          = %g \t+/- %g\n",       getParameter("Width")           ,_pParametersErrors[1]);
00472                 printf ("height         = %g \t+/- %g\n",       getParameter("Height")          ,_pParametersErrors[2]);
00473                 printf ("background     = %g \t+/- %g\n",       getParameter("Background")      ,_pParametersErrors[3]);
00474         printf ("chisq/dof      = %g\n",                                _dQualityFactor);
00475                 printf ("Sigma          = %g\n",                                _dSigma);
00476                 printf ("HWHM           = %g\n",                                _dHWHM);
00477                 printf ("FHWM           = %g\n",                                _dFHWM);
00478 }
00479 
00480 
00481 double  FittingFunctionGaussianWithBackground::computeValue(double dX,double dPosition, double dWidth, double dHeight, double dBackground)
00482 {
00483         double dXPos,dXPos2,dWidth2,dYi;
00484 
00485         dXPos   = dX - dPosition;
00486         dXPos2  = dXPos * dXPos;
00487         dWidth2 = dWidth * dWidth;
00488         
00489         dYi = dBackground + dHeight * exp(-0.5 * dXPos2 / dWidth2);
00490         return dYi;
00491 }
00492 
00494 double  FittingFunctionGaussianWithBackground::computeDerivateValue(double dX,double dPosition, double dWidth, double dHeight)
00495 {
00496         return 0.0;
00497 }
00498 
00499 double FittingFunctionGaussianWithBackground::estimateInitialPosition(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00500 {
00501         size_t max_index        = gsl_vector_max_index (vExperimentalDataY);
00502         double dPosition        = gsl_vector_get(vExperimentalDataX,max_index);
00503 
00504         return dPosition;
00505 }
00506 
00507 
00508 
00509 
00510 double FittingFunctionGaussianWithBackground::estimateInitialHeight(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00511 {
00512         size_t max_index        = gsl_vector_max_index (vExperimentalDataY);
00513 
00514         double background       = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00515         double maxValue         = gsl_vector_get(vExperimentalDataY, max_index);
00516         double height           = maxValue - background;
00517         return height;
00518 }
00519 
00520 
00521 double FittingFunctionGaussianWithBackground::estimateInitialBackground(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00522 {
00523         double dBackground = 0.0;
00524         size_t min_index = gsl_vector_min_index (vExperimentalDataY);
00525         //TODO do a mean with more points to better estimate background
00526         dBackground =  gsl_vector_get(vExperimentalDataY, min_index);
00527 
00528         return dBackground;
00529 }
00530 
00531 
00532 double FittingFunctionGaussianWithBackground::estimateInitialSigma(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00533 {
00534         double dSigma;
00535         double dCovariance;
00536         double dHeight                          = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00537         double dBackground                      = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00538         
00539         size_t iIndexMax                        = gsl_vector_max_index(vExperimentalDataY);
00540         double dYMaxValue                       = gsl_vector_get(vExperimentalDataY, iIndexMax);
00541   
00542         double  dYCurrentValue          = dYMaxValue;
00543 
00544         size_t iCurrentIndex            = iIndexMax;
00545         size_t lVectorSize                      = vExperimentalDataX->size;//lSizeX;
00546         double lHalfVectorSize          = 0.5 * lVectorSize;
00547         double dHalfMaximum                     = 0.5 * dHeight + dBackground;
00548   
00549         size_t iIndexLow                        = iCurrentIndex;
00550         size_t iIndexHigh                       = iCurrentIndex;
00551 
00552     if (iIndexMax < lHalfVectorSize)
00553     {
00554 //              cout << "\n\t--> Look in increasing direction" << endl;
00555                 //- look in increasing direction
00556                 while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue > dHalfMaximum)
00557                 {
00558                         dYCurrentValue  = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00559                         iIndexHigh              = iCurrentIndex;
00560                         iIndexLow               = iIndexMax - iCurrentIndex;
00561                 }
00562     }
00563     else
00564     {
00565 //              cout << "\n\t--> Look in decreasing direction" << endl;
00566                 //- look in decreasing direction
00567                 while( iCurrentIndex > 0 && dYCurrentValue > dHalfMaximum)
00568                 {
00569                         dYCurrentValue =  gsl_vector_get(vExperimentalDataY, --iCurrentIndex);
00570                         iIndexLow               = iCurrentIndex;
00571                         iIndexHigh              = iIndexMax + iCurrentIndex;
00572                 }
00573     }
00574 
00575     if (iCurrentIndex != iIndexMax)
00576     {
00577                 double dPosition        = gsl_vector_get(vExperimentalDataX, iIndexMax);
00578                 double dX2                      = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00579                 dSigma = ::fabs(dX2 - dPosition);
00580     }
00581     else
00582     {
00583                 //- try an arbitrary value that should be not so far from reality
00584                 dSigma = lVectorSize / 6;
00585     }
00586 
00587 /*      cout << "----------------------------" << endl;
00588         cout << "Peak : (x,y) --> ("    << gsl_vector_get(vExperimentalDataX, iIndexMax) << "," 
00589                                                                         << gsl_vector_get(vExperimentalDataY, iIndexMax) << ")" << endl;
00590                 
00591         
00592         cout << "Low : index=" << iIndexLow << endl;
00593         cout << "\txlow=" << gsl_vector_get(vExperimentalDataX, iIndexLow) << " | ylow=" << gsl_vector_get(vExperimentalDataY, iIndexLow) << endl;
00594         cout << "\txhigh=" << gsl_vector_get(vExperimentalDataX, iIndexHigh) << " | yhigh=" << gsl_vector_get(vExperimentalDataY, iIndexHigh) << endl;
00595         
00596         cout << "\txhigh - xlow= " << gsl_vector_get(vExperimentalDataX, iIndexHigh) - gsl_vector_get(vExperimentalDataX, iIndexLow) << endl;
00597         cout << "\tsigma=" << dSigma << endl;
00598 */
00599     dCovariance = dSigma * dSigma;
00600         return dSigma;
00601 }

Generated on Tue Apr 14 09:34:27 2009 for Data Fitting Library by  doxygen 1.4.5