00001
00002
00004 #ifdef WIN32
00005 #pragma warning(disable:4786)
00006 #endif
00007
00008 #include "FittingFunctionGaussian.h"
00010
00012
00013
00014 std::string FittingFunctionGaussian::getEquation() const
00015 {
00016 return "f(x) = Height * exp[ - (x - Position)˛ / (2 * Width˛)]";
00017 }
00018
00019
00020
00021
00022 FittingFunctionGaussian::FittingFunctionGaussian()
00023 {
00024 _iNbParameters = 3;
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 _vInitialGuess = gsl_vector_alloc(3);
00033 }
00034
00035
00036 FittingFunctionGaussian::FittingFunctionGaussian(long lNbIteration,double dLimitSearch)
00037 {
00038 _iNbParameters = 3;
00039 _lNbMaxIteration = lNbIteration;
00040 _dSearchLimit = dLimitSearch;
00041 _pParameters = new double[_iNbParameters];
00042 _pParametersErrors = new double[_iNbParameters];
00043 _mFunctionMap["Position"] = 0;
00044 _mFunctionMap["Width"] = 1;
00045 _mFunctionMap["Height"] = 2;
00046 _vInitialGuess = gsl_vector_alloc(3);
00047 }
00048
00049
00050 FittingFunctionGaussian::~FittingFunctionGaussian()
00051 {
00052 printf ("Destruction FittingFunctionGaussian !!!\n");
00053 delete[] _pParameters;
00054 delete[] _pParametersErrors;
00055 _pParameters = NULL;
00056 _pParametersErrors = NULL;
00057 gsl_vector_free(_vInitialGuess);
00058 }
00059
00060
00061 PtrFunction FittingFunctionGaussian::getFunction() const{
00062 return FittingFunctionGaussian::Function;
00063 }
00064
00065
00066 PtrFunctionDerivate FittingFunctionGaussian::getFunctionDerivate() const{
00067 return FittingFunctionGaussian::FunctionDerivate;
00068 }
00069
00070
00071 PtrFunctionAndDerivate FittingFunctionGaussian::getFunctionAndDerivate() const {
00072 return FittingFunctionGaussian::FunctionAndDerivate;
00073 }
00074
00075
00080
00081 int FittingFunctionGaussian::Function(const gsl_vector *peak_params, void *calc_data, gsl_vector *func)
00082 {
00083 size_t i;
00084 size_t iNbPoints = ((struct FittingData *)calc_data)->n;
00085 double *x = ((struct FittingData *)calc_data)->x;
00086 double *y = ((struct FittingData *)calc_data)->y;
00087 double *sigma = ((struct FittingData *)calc_data)->sigma;
00088
00089 double pos = gsl_vector_get (peak_params, 0);
00090 double width = gsl_vector_get (peak_params, 1);
00091 double height = gsl_vector_get (peak_params, 2);
00092
00093
00094
00095
00096
00097 double dXPos,dXPos2,dWidth2,dYi;
00098
00099 for (i=0; i < iNbPoints; i++)
00100 {
00101 dXPos = x[i] - pos;
00102 dXPos2 = dXPos * dXPos;
00103 dWidth2 = width * width;
00104
00105 dYi = height * exp(-0.5 * dXPos2 / dWidth2);
00106
00107 gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);
00108 }
00109 return GSL_SUCCESS;
00110 }
00111
00122
00123 int FittingFunctionGaussian::FunctionDerivate(const gsl_vector *peak_params, void *calc_data, gsl_matrix *Jacobian)
00124 {
00125 size_t i;
00126
00127 size_t n = ((struct FittingData *)calc_data)->n;
00128 double *x = ((struct FittingData *)calc_data)->x;
00129 double *sigma = ((struct FittingData *)calc_data)->sigma;
00130
00131 double pos = gsl_vector_get (peak_params, 0);
00132 double width = gsl_vector_get (peak_params, 1);
00133 double height = gsl_vector_get (peak_params, 2);
00134
00135 double d_pos, d_width, d_height;
00136
00137 double s, xpos2, width2,e;
00138
00139 for (i=0; i< n; i++)
00140 {
00141 s = sigma[i];
00142 xpos2 = (x[i]-pos)*(x[i]-pos);
00143 width2 = width * width;
00144 e = exp( -0.5 * xpos2 / width2);
00145
00146 d_pos = (height/width2 * (x[i] - pos) * e) / s;
00147 d_width = (height/width2/width * xpos2 * e) / s;
00148 d_height = e/s;
00149
00150
00151 gsl_matrix_set (Jacobian, i, 0, d_pos);
00152 gsl_matrix_set (Jacobian, i, 1, d_width);
00153 gsl_matrix_set (Jacobian, i, 2, d_height);
00154 }
00155 return GSL_SUCCESS;
00156 }
00157
00158
00159
00160 int FittingFunctionGaussian::FunctionAndDerivate(const gsl_vector *peak_params, void *calc_data,gsl_vector *func, gsl_matrix *Jacobian)
00161 {
00162 Function (peak_params, calc_data, func);
00163 FunctionDerivate(peak_params, calc_data, Jacobian);
00164 return GSL_SUCCESS;
00165 }
00166
00167
00168
00169 int FittingFunctionGaussian::doFit(long points, gsl_vector *x, gsl_vector *y,gsl_vector *sigma,gsl_vector* guess, bool bEstimateGuess, int iSearchStoppingMethod,bool bUseScaled)
00170 {
00171
00172 initializeInitialsParameters(x,y,guess,bEstimateGuess);
00173
00174 long i;
00175 long iter;
00176 long status;
00177
00178 double* datax = new double[points];
00179 double* datay = new double[points];
00180 double* datas = new double[points];
00181
00182
00183 for (i=0; i<points; i++)
00184 {
00185 datax[i] = gsl_vector_get (x, i);
00186 datay[i] = gsl_vector_get (y, i);
00187 datas[i] = gsl_vector_get (sigma, i);
00188 }
00189
00190
00191 struct FittingData gsl_data;
00192 gsl_data.n = points;
00193 gsl_data.x = datax;
00194 gsl_data.y = datay;
00195 gsl_data.sigma = datas;
00196
00197 cout << "\n---------------------------------------------" << endl;
00198 cout << "Initials Parameters For The Fit" << endl;
00199
00200 const gsl_multifit_fdfsolver_type *T;
00201 if (bUseScaled)
00202 {
00203 T = gsl_multifit_fdfsolver_lmsder;
00204 cout << "Jacobian Data Scaled" << endl;
00205 }
00206 else
00207 {
00208 T = gsl_multifit_fdfsolver_lmder;
00209 cout << "Jacobian Data Not Scaled" << endl;
00210 }
00211 gsl_multifit_fdfsolver *s = gsl_multifit_fdfsolver_alloc (T, points, _iNbParameters);
00212 gsl_multifit_function_fdf f;
00213
00214 f.f = getFunction();
00215 f.df = getFunctionDerivate();
00216 f.fdf = getFunctionAndDerivate();
00217 f.n = points;
00218 f.p = _iNbParameters;
00219 f.params = &gsl_data;
00220
00221 switch(iSearchStoppingMethod)
00222 {
00223 case 1 :
00224 cout << "Test Delta Search Stopping Function Used" << endl;
00225 break;
00226
00227 case 2 :
00228 cout << "Test Gradient Search Stopping Function Used" << endl;
00229 break;
00230 }
00231
00232 cout << "\n----------------------------------------------" << endl;
00233 cout << "Initials Values For The Fit " << endl;
00234 cout << "Position = " << gsl_vector_get (_vInitialGuess, 0) << endl;
00235 cout << "Width = " << gsl_vector_get (_vInitialGuess, 1) << endl;
00236 cout << "Height = " << gsl_vector_get (_vInitialGuess, 2) << endl;
00237
00238
00239 gsl_multifit_fdfsolver_set ( s, &f, _vInitialGuess );
00240
00241 cout << "\n---------------------------------------------" << endl;
00242 cout << "Iterations" << endl;
00243 iter = 0;
00244 #ifdef DEBUG_FIT
00245 printState (iter, s);
00246 #endif
00247
00248 do
00249 {
00250 iter++;
00251 status = gsl_multifit_fdfsolver_iterate (s);
00252
00253 #ifdef DEBUG_FIT
00254 printf ("status = %s\n", gsl_strerror (status));
00255 printState (iter, s);
00256 #endif
00257
00258 if (status)
00259 {
00260 break;
00261 }
00262
00263
00264 switch(iSearchStoppingMethod)
00265 {
00266 case 1 :
00267 status = gsl_multifit_test_delta (s->dx, s->x, 0.0, _dSearchLimit);
00268 break;
00269
00270 case 2 :
00271 gsl_vector* gradient = gsl_vector_alloc(_iNbParameters);
00272 gsl_multifit_gradient(s->J,s->f,gradient);
00273 status = gsl_multifit_test_gradient (gradient, _dSearchLimit);
00274 gsl_vector_free(gradient);
00275 break;
00276 }
00277 }
00278 while (status == GSL_CONTINUE && iter < _lNbMaxIteration);
00279
00280 _lNbIterations = iter;
00281
00282
00283
00284
00285 gsl_matrix *covar = gsl_matrix_alloc (_iNbParameters, _iNbParameters);
00286 gsl_multifit_covar (s->J, 0.0, covar);
00287 double chi = gsl_blas_dnrm2(s->f);
00288
00289 _pParameters[0] = FIT(0);
00290 _pParameters[1] = FIT(1);
00291 _pParameters[2] = FIT(2);
00292
00293 _pParametersErrors[0] = ERR(0);
00294 _pParametersErrors[1] = ERR(1);
00295 _pParametersErrors[2] = ERR(2);
00296
00297
00298
00299
00300
00301
00302 double dEstimatedWidth = _pParameters[1];
00303 _dQualityFactor = pow(chi, 2.0)/ (points - _iNbParameters);
00304
00305
00306 _dHWHM = dEstimatedWidth * sqrt(2*log(2.0));
00307 _dFHWM = 2*_dHWHM;
00308 cout << "2*sqrt(2*log(2.0))=" << 2*sqrt(2*log(2.0)) <<endl;
00309
00310 #ifdef DEBUG_FIT
00311 cout << "\n---------------------------------------------" << endl;
00312 cout << "Results For Gaussian" << endl;
00313 printResults();
00314 #endif
00315
00316 gsl_matrix_free(covar);
00317
00318
00319 #ifdef DEBUG_FIT
00320 printf ("status = %s\n", gsl_strerror (status));
00321 #endif
00322
00323
00324 gsl_multifit_fdfsolver_free (s);
00325
00326 delete[] datax;
00327 delete[] datay;
00328 delete[] datas;
00329
00330 datax = NULL;
00331 datax = NULL;
00332 datax = NULL;
00333
00334 return status;
00335 }
00336
00337
00338
00339 int FittingFunctionGaussian::generateFunctionFit(double dValMinX,double dResolutionX,long lNbPoints,gsl_vector* funcX,gsl_vector* funcY)
00340 {
00341
00342 int i;
00343 double dPos = _pParameters[0];
00344 double dWidth = _pParameters[1];
00345 double dHeight = _pParameters[2];
00346 double dXi = 0;
00347 double dXPos = 0;
00348 double dXPos2 = 0;
00349 double dWidth2 = 1;
00350 double dYi = 0;
00351
00352 for (i=0;i<lNbPoints;i++)
00353 {
00354 dXi = dValMinX+i*dResolutionX;
00355 dXPos = dXi - dPos;
00356 dXPos2 = dXPos * dXPos;
00357 dWidth2 = dWidth * dWidth;
00358 dYi = dHeight * exp(-0.5 * dXPos2 / dWidth2);
00359
00360 gsl_vector_set(funcX,i,dXi);
00361 gsl_vector_set(funcY,i,dYi);
00362 }
00363 return 1;
00364 }
00365
00366 int FittingFunctionGaussian::generateFunctionFit(long lNbPoints, gsl_vector *x, gsl_vector* funcX,gsl_vector* funcY)
00367 {
00368 int i;
00369 double dPos = _pParameters[0];
00370 double dWidth = _pParameters[1];
00371 double dHeight = _pParameters[2];
00372 double dXi = 0;
00373 double dXPos = 0;
00374 double dXPos2 = 0;
00375 double dWidth2 = 1;
00376 double dYi = 0;
00377
00378 for (i=0;i<lNbPoints;i++)
00379 {
00380 dXi = gsl_vector_get(x,i);
00381
00382 dXPos = dXi - dPos;
00383 dXPos2 = dXPos * dXPos;
00384 dWidth2 = dWidth * dWidth;
00385 dYi = dHeight * exp(-0.5 * dXPos2 / dWidth2);
00386
00387 gsl_vector_set(funcX,i,dXi);
00388 gsl_vector_set(funcY,i,dYi);
00389 }
00390 return 1;
00391 }
00392
00393
00394
00395 void FittingFunctionGaussian::printState (size_t iter, gsl_multifit_fdfsolver * s)
00396 {
00397 printf ("iter: %3u x : %g %g %g |f(x)|=%g\n",
00398 iter,
00399 gsl_vector_get (s->x, 0),
00400 gsl_vector_get (s->x, 1),
00401 gsl_vector_get (s->x, 2),
00402 gsl_blas_dnrm2 (s->f));
00403 }
00404
00405
00406 void FittingFunctionGaussian::printResults()
00407 {
00408 printf ("pos = %g \t+/- %g\n", getParameter("Position"), _pParametersErrors[0]);
00409 printf ("width = %g \t+/- %g\n", getParameter("Width") , _pParametersErrors[1]);
00410 printf ("height = %g \t+/- %g\n", getParameter("Height") , _pParametersErrors[2]);
00411 printf ("chisq/dof = %g\n", _dQualityFactor);
00412 printf ("Sigma = %g\n", _dSigma);
00413 printf ("HWHM = %g\n", _dHWHM);
00414 printf ("FHWM = %g\n", _dFHWM);
00415 }
00416
00417 double FittingFunctionGaussian::computeValue(double dX,double dPosition, double dWidth, double dHeight, double dBackground)
00418 {
00419 double dXPos,dXPos2,dWidth2,dYi;
00420
00421 dXPos = dX - dPosition;
00422 dXPos2 = dXPos * dXPos;
00423 dWidth2 = dWidth * dWidth;
00424
00425 dYi = dHeight * exp(-0.5 * dXPos2 / dWidth2);
00426 return dYi;
00427 }
00428
00430 double FittingFunctionGaussian::computeDerivateValue(double dX,double dPosition, double dWidth, double dHeight)
00431 {
00432 return 0.0;
00433 }
00434
00435
00436
00437 void FittingFunctionGaussian::estimateInitialGuess(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY)
00438 {
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449 double dInitialPosition = estimateInitialPosition(vExperimentalDataX,vExperimentalDataY);
00450 double dInitialSigma = estimateInitialSigma(vExperimentalDataX,vExperimentalDataY);
00451 double dInitialHeight = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00452
00453 gsl_vector_set(_vInitialGuess, 0, dInitialPosition);
00454 gsl_vector_set(_vInitialGuess, 1, dInitialSigma);
00455 gsl_vector_set(_vInitialGuess, 2, dInitialHeight);
00456 }
00457
00458
00459
00460 void FittingFunctionGaussian::initializeInitialsParameters(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY,gsl_vector* vInitialGuess,bool bEstimate)
00461 {
00462
00463 if (bEstimate)
00464 {
00465 estimateInitialGuess(vExperimentalDataX,vExperimentalDataY);
00466 }
00467 else
00468 {
00469 gsl_vector_set(_vInitialGuess,0,gsl_vector_get(vInitialGuess,0));
00470 gsl_vector_set(_vInitialGuess,1,gsl_vector_get(vInitialGuess,1));
00471 gsl_vector_set(_vInitialGuess,2,gsl_vector_get(vInitialGuess,2));
00472 }
00473 }
00474
00475
00476 double FittingFunctionGaussian::estimateInitialPosition(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00477 {
00478 size_t max_index = gsl_vector_max_index (vExperimentalDataY);
00479 double dPosition = gsl_vector_get(vExperimentalDataX,max_index);
00480
00481 return dPosition;
00482 }
00483
00484
00485
00486
00487 double FittingFunctionGaussian::estimateInitialHeight(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00488 {
00489 size_t max_index = gsl_vector_max_index (vExperimentalDataY);
00490
00491 double background = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00492 double maxValue = gsl_vector_get(vExperimentalDataY, max_index);
00493 double height = maxValue - background;
00494 return height;
00495 }
00496
00497
00498 double FittingFunctionGaussian::estimateInitialBackground(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00499 {
00500 double dBackground = 0.0;
00501 return dBackground;
00502 }
00503
00504
00505 double FittingFunctionGaussian::estimateInitialSigma(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00506 {
00507 double dSigma;
00508 double dCovariance;
00509 double dHeight = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00510 double dBackground = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00511
00512 size_t iIndexMax = gsl_vector_max_index(vExperimentalDataY);
00513 double dYMaxValue = gsl_vector_get(vExperimentalDataY, iIndexMax);
00514
00515 double dYCurrentValue = dYMaxValue;
00516
00517 size_t iCurrentIndex = iIndexMax;
00518 size_t lVectorSize = vExperimentalDataX->size;
00519 double lHalfVectorSize = 0.5 * lVectorSize;
00520 double dHalfMaximum = 0.5 * dHeight + dBackground;
00521
00522 size_t iIndexLow = iCurrentIndex;
00523 size_t iIndexHigh = iCurrentIndex;
00524
00525 if (iIndexMax < lHalfVectorSize)
00526 {
00527
00528
00529 while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue > dHalfMaximum)
00530 {
00531 dYCurrentValue = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00532 iIndexHigh = iCurrentIndex;
00533 iIndexLow = iIndexMax - (iIndexHigh - iIndexMax);
00534 }
00535 }
00536 else
00537 {
00538
00539
00540 while( iCurrentIndex > 0 && dYCurrentValue > dHalfMaximum)
00541 {
00542 dYCurrentValue = gsl_vector_get(vExperimentalDataY, --iCurrentIndex);
00543 iIndexLow = iCurrentIndex;
00544 iIndexHigh = iIndexMax + (iIndexMax - iIndexLow);
00545 }
00546 }
00547
00548 if (iCurrentIndex != iIndexMax)
00549 {
00550 double dPosition = gsl_vector_get(vExperimentalDataX, iIndexMax);
00551 double dX2 = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00552 dSigma = ::fabs(dX2 - dPosition);
00553 }
00554 else
00555 {
00556
00557 dSigma = lVectorSize / 6;
00558 }
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571 dCovariance = dSigma * dSigma;
00572 return dSigma;
00573 }
00574