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