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