00001
00002
00004 #ifdef WIN32
00005 #pragma warning(disable:4786)
00006 #endif
00007
00008 #include "FittingFunctionGaussianWithBackground.h"
00009
00010
00012
00014
00015 std::string FittingFunctionGaussianWithBackground::getEquation() const
00016 {
00017 return "f(x) = Height * exp[ - (x - Position)˛ / (2 * Width˛)] + Background";
00018 }
00019
00020
00021
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
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
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
00063 PtrFunction FittingFunctionGaussianWithBackground::getFunction() const{
00064 return FittingFunctionGaussianWithBackground::Function;
00065 }
00066
00067
00068 PtrFunctionDerivate FittingFunctionGaussianWithBackground::getFunctionDerivate() const{
00069 return FittingFunctionGaussianWithBackground::FunctionDerivate;
00070 }
00071
00072
00073 PtrFunctionAndDerivate FittingFunctionGaussianWithBackground::getFunctionAndDerivate() const {
00074 return FittingFunctionGaussianWithBackground::FunctionAndDerivate;
00075 }
00076
00077
00082
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
00097
00098
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
00111 gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);
00112 }
00113
00114 return GSL_SUCCESS;
00115 }
00116
00129
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
00143
00144 double d_pos, d_width, d_height, d_background;
00145 double s,xpos2,width2,e;
00146
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;
00155 d_width = (height/width2/width * xpos2 * e) / s;
00156 d_height = e/s;
00157 d_background = 1/s;
00158
00159
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
00171
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
00180
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
00195 for (i=0; i<points; i++)
00196 {
00197 datax[i] = gsl_vector_get (x, i);
00198 datay[i] = gsl_vector_get (y, i);
00199 datas[i] = gsl_vector_get (sigma, i);
00200 }
00201
00202
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
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();
00227 f.df = getFunctionDerivate();
00228 f.fdf = getFunctionAndDerivate();
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, 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
00296
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
00313
00314 double dEstimatedWidth = _pParameters[1];
00315 _dQualityFactor = pow(chi, 2.0)/ (points - _iNbParameters);
00316
00317
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
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
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
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
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
00455
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
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
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;
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
00555
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
00566
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
00584 dSigma = lVectorSize / 6;
00585 }
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599 dCovariance = dSigma * dSigma;
00600 return dSigma;
00601 }