00001
00002
00004 #ifdef WIN32
00005 #pragma warning(disable:4786)
00006 #endif
00007
00008 #include "FittingFunctionSigmoidWithBackground.h"
00009
00010
00012
00014
00015 std::string FittingFunctionSigmoidWithBackground::getEquation() const
00016 {
00017 return "f(x) = Height / [1 + exp[ - (x - Position) / Width] + Background";
00018 }
00019
00020
00021 FittingFunctionSigmoidWithBackground::FittingFunctionSigmoidWithBackground() : _dXLow(0),_dXHigh(0)
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
00037 FittingFunctionSigmoidWithBackground::FittingFunctionSigmoidWithBackground(long lNbIteration,double dLimitSearch) : _dXLow(0),_dXHigh(0)
00038 {
00039 _lNbMaxIteration = lNbIteration;
00040 _dSearchLimit = dLimitSearch;
00041 _iNbParameters = 4;
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 FittingFunctionSigmoidWithBackground::~FittingFunctionSigmoidWithBackground()
00053 {
00054 printf ("Destruction FittingFunctionSigmoidWithBackground !!!\n");
00055 delete[] _pParameters;
00056 delete[] _pParametersErrors;
00057 _pParameters = NULL;
00058 _pParametersErrors = NULL;
00059 gsl_vector_free(_vInitialGuess);
00060 }
00061
00062
00063 PtrFunction FittingFunctionSigmoidWithBackground::getFunction() const{
00064 return FittingFunctionSigmoidWithBackground::Function;
00065 }
00066
00067
00068 PtrFunctionDerivate FittingFunctionSigmoidWithBackground::getFunctionDerivate() const{
00069 return FittingFunctionSigmoidWithBackground::FunctionDerivate;
00070 }
00071
00072
00073 PtrFunctionAndDerivate FittingFunctionSigmoidWithBackground::getFunctionAndDerivate() const {
00074 return FittingFunctionSigmoidWithBackground::FunctionAndDerivate;
00075 }
00076
00077
00082
00083 int FittingFunctionSigmoidWithBackground::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 if (width == 0) width = 1E-20;
00097
00098 #ifdef DEBUG_FIT
00099
00100 #endif
00101 double dXPos,dYi;
00102 for (i=0; i < iNbPoints; i++)
00103 {
00104 dXPos = x[i] - pos;
00105 dYi = background + height / (1 + exp(-dXPos/width));
00106
00107 gsl_vector_set (func, i,( dYi - y[i]) / sigma[i]);
00108 }
00109
00110 return GSL_SUCCESS;
00111 }
00112
00113
00126
00127 int FittingFunctionSigmoidWithBackground::FunctionDerivate(const gsl_vector *peak_params, void *calc_data, gsl_matrix *Jacobian)
00128 {
00129
00130 size_t i;
00131
00132 size_t n = ((struct FittingData *)calc_data)->n;
00133 double *x = ((struct FittingData *)calc_data)->x;
00134 double *sigma = ((struct FittingData *)calc_data)->sigma;
00135
00136 double pos = gsl_vector_get (peak_params, 0);
00137 double width = gsl_vector_get (peak_params, 1);
00138 double height = gsl_vector_get (peak_params, 2);
00139
00140
00141 double d_pos, d_width, d_height, d_background;
00142 double s, alpha, e, e2,expo,xpos;
00143
00144
00145
00146
00147
00148
00149 for (i=0; i< n; i++)
00150 {
00151 s = sigma[i];
00152 xpos = x[i] - pos;
00153 alpha = -xpos/width;
00154 expo = exp(alpha);
00155 e = 1 + expo;
00156 e2 = e * e;
00157
00158
00159 d_pos = -((height*expo)/(e2*width)) / s;
00160
00161 d_width = -(((height*xpos)*expo)/(e2*width*width)) / s;
00162
00163 d_height = (1/e)/s;
00164
00165 d_background = 1/s;
00166
00167
00168 gsl_matrix_set (Jacobian, i, 0, d_pos);
00169 gsl_matrix_set (Jacobian, i, 1, d_width);
00170 gsl_matrix_set (Jacobian, i, 2, d_height);
00171 gsl_matrix_set (Jacobian, i, 3, d_background);
00172 }
00173 return GSL_SUCCESS;
00174 }
00175
00176
00177
00178 int FittingFunctionSigmoidWithBackground::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 FittingFunctionSigmoidWithBackground::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 cout << "\n---------------------------------------------" << endl;
00215 cout << "Initials Parameters For The Fit" << endl;
00216
00217 const gsl_multifit_fdfsolver_type *T;
00218 if (bUseScaled)
00219 {
00220 T = gsl_multifit_fdfsolver_lmsder;
00221 cout << "Jacobian Data Scaled" << endl;
00222 }
00223 else
00224 {
00225 T = gsl_multifit_fdfsolver_lmder;
00226 cout << "Jacobian Data Not Scaled" << endl;
00227 }
00228 gsl_multifit_fdfsolver *s = gsl_multifit_fdfsolver_alloc (T, points, _iNbParameters);
00229 gsl_multifit_function_fdf f;
00230
00231 f.f = getFunction();
00232 f.df = getFunctionDerivate();
00233 f.fdf = getFunctionAndDerivate();
00234 f.n = points;
00235 f.p = _iNbParameters;
00236 f.params = &gsl_data;
00237
00238 gsl_multifit_fdfsolver_set ( s, &f, _vInitialGuess );
00239
00240 switch(iSearchStoppingMethod)
00241 {
00242 case 1 :
00243 cout << "Test Delta Search Stopping Function Used" << endl;
00244 break;
00245
00246 case 2 :
00247 cout << "Test Gradient Search Stopping Function Used" << endl;
00248 break;
00249 }
00250
00251 cout << "\n------------------------------------------" << endl;
00252 cout << "Initials Values For The Fit " << endl;
00253 cout << "Position = " << gsl_vector_get (_vInitialGuess, 0) << endl;
00254 cout << "Width = " << gsl_vector_get (_vInitialGuess, 1) << endl;
00255 cout << "Height = " << gsl_vector_get (_vInitialGuess, 2) << endl;
00256 cout << "Background= " << gsl_vector_get (_vInitialGuess, 3) << endl;
00257
00258
00259 iter = 0;
00260
00261 cout << "\n---------------------------------------------" << endl;
00262 cout << "Iterations" << endl;
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 double dEstimatedWidth = _pParameters[1];
00317 _dQualityFactor = pow(chi, 2.0)/ (points - _iNbParameters);
00318
00319
00320
00321
00322
00323 _dFHWM = _pParameters[1];
00324 _dHWHM = 0.5*_dFHWM;
00325
00326 #ifdef DEBUG_FIT
00327 cout << "\n---------------------------------------------" << endl;
00328 printf ("Results For Sigmoid With Background\n");
00329 printResults();
00330 #endif
00331
00332 gsl_matrix_free(covar);
00333
00334
00335 #ifdef DEBUG_FIT
00336 printf ("status = %s\n", gsl_strerror (status));
00337 #endif
00338
00339
00340 gsl_multifit_fdfsolver_free (s);
00341
00342 estimateSigmoidIntersections();
00343
00344
00345 delete[] datax;
00346 delete[] datay;
00347 delete[] datas;
00348
00349 datax = NULL;
00350 datax = NULL;
00351 datax = NULL;
00352
00353 return status;
00354 }
00355
00356
00357
00358 int FittingFunctionSigmoidWithBackground::generateFunctionFit(double dValMinX,double dResolutionX,long lNbPoints,gsl_vector* funcX,gsl_vector* funcY){
00359
00360 int i;
00361
00362 double dPos = _pParameters[0];
00363 double dWidth = _pParameters[1];
00364 double dHeight = _pParameters[2];
00365 double dBackground = _pParameters[3];
00366
00367 double dXi,dXPos,dExp,dYi;
00368
00369 for (i=0; i < lNbPoints; i++)
00370 {
00371 dXi = dValMinX + i*dResolutionX;
00372 dXPos = dXi - dPos;
00373 dExp = exp(-dXPos/dWidth);
00374
00375 dYi = dBackground + dHeight / (1.0 + dExp);
00376
00377 gsl_vector_set (funcX, i, dXi);
00378 gsl_vector_set (funcY, i, dYi);
00379 }
00380 return 1;
00381 }
00382
00383
00384
00385 int FittingFunctionSigmoidWithBackground::generateFunctionFit(long lNbPoints, gsl_vector *x, gsl_vector* funcX,gsl_vector* funcY)
00386 {
00387 int i;
00388
00389 double dPos = _pParameters[0];
00390 double dWidth = _pParameters[1];
00391 double dHeight = _pParameters[2];
00392 double dBackground = _pParameters[3];
00393
00394 double dXi,dXPos,dExp,dYi;
00395
00396 for (i=0; i < lNbPoints; i++)
00397 {
00398 dXi = gsl_vector_get(x,i);
00399 dXPos = dXi - dPos;
00400 dExp = exp(-dXPos/dWidth);
00401
00402 dYi = dBackground + dHeight / (1.0 + dExp);
00403
00404 gsl_vector_set (funcX, i, dXi);
00405 gsl_vector_set (funcY, i, dYi);
00406 }
00407 return 1;
00408 }
00409
00410
00411
00412
00413
00414
00415 void FittingFunctionSigmoidWithBackground::printState (size_t iter, gsl_multifit_fdfsolver * s)
00416 {
00417 printf ("iter: %3u x : %g %g %g %g |f(x)|=%g\n",
00418 iter,
00419 gsl_vector_get (s->x, 0),
00420 gsl_vector_get (s->x, 1),
00421 gsl_vector_get (s->x, 2),
00422 gsl_vector_get (s->x, 3),
00423 gsl_blas_dnrm2 (s->f));
00424 }
00425
00426
00427 void FittingFunctionSigmoidWithBackground::printResults()
00428 {
00429 printf ("pos = %g \t+/- %g\n", getParameter("Position") , _pParametersErrors[0]);
00430 printf ("width = %g \t+/- %g\n", getParameter("Width") , _pParametersErrors[1]);
00431 printf ("height = %g \t+/- %g\n", getParameter("Height") , _pParametersErrors[2]);
00432 printf ("background = %g \t+/- %g\n", getParameter("Background") , _pParametersErrors[3]);
00433 printf ("chisq/dof = %g\n", _dQualityFactor);
00434 printf ("Sigma = %g\n", _dSigma);
00435 printf ("HWHM = %g\n", _dHWHM);
00436 printf ("FHWM = %g\n", _dFHWM);
00437 }
00438
00439 double FittingFunctionSigmoidWithBackground::computeValue(double dX,double dPosition, double dWidth, double dHeight, double dBackground)
00440 {
00441
00442 double dXPos,dYi;
00443
00444 if (dWidth == 0) dWidth = 1E-20;
00445
00446 dXPos = dX - dPosition;
00447 dYi = dBackground + dHeight / (1 + exp(-dXPos/dWidth));
00448
00449 return dYi;
00450 }
00451
00452
00453 double FittingFunctionSigmoidWithBackground::computeDerivateValue(double dX,double dPosition, double dWidth, double dHeight)
00454 {
00455 double dXPos,dYi;
00456
00457 if (dWidth == 0) dWidth = 1E-20;
00458
00459 dXPos = dX - dPosition;
00460 dYi = dHeight * exp(-dXPos/dWidth) / ((1 + exp(-dXPos/dWidth)) * (1 + exp(-dXPos/dWidth)) * dWidth);
00461
00462 return dYi;
00463 }
00464
00465
00466 double FittingFunctionSigmoidWithBackground::estimateInitialPosition(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00467 {
00468 double dHeight = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00469 double dBackground = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00470 double dHalfMaximum = 0.5 * dHeight + dBackground;
00471
00472 size_t lVectorSize = vExperimentalDataX->size;;
00473 size_t iCurrentIndex = 0;
00474 double dYCurrentValue = gsl_vector_get(vExperimentalDataY, iCurrentIndex);
00475
00476 double dPosition;
00477
00478 while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue < dHalfMaximum)
00479 {
00480 dYCurrentValue = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00481 }
00482
00483 if (iCurrentIndex != lVectorSize)
00484 {
00485 dPosition = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00486 _lPositionIndex = (long)iCurrentIndex;
00487 }
00488 else
00489 {
00490
00491 dPosition = gsl_vector_get(vExperimentalDataX, lVectorSize / 2);
00492 _lPositionIndex = (long)(lVectorSize/2);
00493 }
00494
00495
00496 return dPosition;
00497 }
00498
00499 double FittingFunctionSigmoidWithBackground::estimateInitialHeight(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00500 {
00501 size_t max_index = gsl_vector_max_index (vExperimentalDataY);
00502
00503 double background = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00504 double maxValue = gsl_vector_get(vExperimentalDataY, max_index);
00505 double height = maxValue - background;
00506 return height;
00507 }
00508
00509
00510 double FittingFunctionSigmoidWithBackground::estimateInitialBackground(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00511 {
00512 double dBackground = 0.0;
00513
00514 size_t min_index = gsl_vector_min_index (vExperimentalDataY);
00515
00516 dBackground = gsl_vector_get(vExperimentalDataY, min_index);
00517
00518 return dBackground;
00519 }
00520
00521
00522 double FittingFunctionSigmoidWithBackground::estimateInitialSigma(gsl_vector *vExperimentalDataX, gsl_vector *vExperimentalDataY)
00523 {
00524 double dSigma;
00525 double dCovariance;
00526 double dHeight = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00527 double dBackground = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00528
00529 size_t iIndexMax = gsl_vector_max_index(vExperimentalDataY);
00530 double dYMaxValue = gsl_vector_get(vExperimentalDataY, iIndexMax);
00531
00532 double dYCurrentValue = dYMaxValue;
00533
00534 size_t iCurrentIndex = iIndexMax;
00535 size_t lVectorSize = vExperimentalDataX->size;
00536 double lHalfVectorSize = 0.5 * lVectorSize;
00537 double dHalfMaximum = 0.5 * dHeight + dBackground;
00538
00539 size_t iIndexLow = iCurrentIndex;
00540 size_t iIndexHigh = iCurrentIndex;
00541
00542 if (iIndexMax < lHalfVectorSize)
00543 {
00544
00545
00546 while(iCurrentIndex < lVectorSize - 1 && dYCurrentValue > dHalfMaximum)
00547 {
00548 dYCurrentValue = gsl_vector_get(vExperimentalDataY, ++iCurrentIndex);
00549 iIndexHigh = iCurrentIndex;
00550 iIndexLow = iIndexMax - iCurrentIndex;
00551 }
00552 }
00553 else
00554 {
00555
00556
00557 while( iCurrentIndex > 0 && dYCurrentValue > dHalfMaximum)
00558 {
00559 dYCurrentValue = gsl_vector_get(vExperimentalDataY, --iCurrentIndex);
00560 iIndexLow = iCurrentIndex;
00561 iIndexHigh = iIndexMax + iCurrentIndex;
00562 }
00563 }
00564
00565 if (iCurrentIndex != iIndexMax)
00566 {
00567 double dPosition = gsl_vector_get(vExperimentalDataX, iIndexMax);
00568 double dX2 = gsl_vector_get(vExperimentalDataX, iCurrentIndex);
00569 dSigma = ::fabs(dX2 - dPosition);
00570 }
00571 else
00572 {
00573
00574 dSigma = lVectorSize / 6;
00575 }
00576
00577 cout << "\n---------------------------------------------" << endl;
00578 cout << "Estimation Initial Width" << endl;
00579 double d34Value = 3.0/4*dHeight + dBackground;
00580 double d12Value = 1.0/2*dHeight + dBackground;
00581 double d14Value = 1.0/4*dHeight + dBackground;
00582
00583 int i34Index = getIndexValue(d34Value,vExperimentalDataY);
00584 int i12Index = getIndexValue(d12Value,vExperimentalDataY);
00585 int i14Index = getIndexValue(d14Value,vExperimentalDataY);
00586
00587 double d34XValue = gsl_vector_get(vExperimentalDataX,i34Index);
00588 double d34YValue = gsl_vector_get(vExperimentalDataY,i34Index);
00589
00590 double d12XValue = gsl_vector_get(vExperimentalDataX,i12Index);
00591 double d12YValue = gsl_vector_get(vExperimentalDataY,i12Index);
00592
00593 double d14XValue = gsl_vector_get(vExperimentalDataX,i14Index);
00594 double d14YValue = gsl_vector_get(vExperimentalDataY,i14Index);
00595
00596 cout << "3/4A=" << "(" << d34Value << ")-->" << i34Index << " (x,y)=(" << d34XValue << "," << d34YValue << ")" << endl;
00597 cout << "1/2A=" << "(" << d12Value << ")-->" << i12Index << " (x,y)=(" << d12XValue << "," << d12YValue << ")" << endl;
00598 cout << "1/4A=" << "(" << d14Value << ")-->" << i14Index << " (x,y)=(" << d14XValue << "," << d14YValue << ")" << endl;
00599
00600 double dWidthEstimated = ::fabs(d34XValue - d14XValue);
00601
00602 cout << "Width estimated = " << dWidthEstimated << endl;
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617 dCovariance = dSigma * dSigma;
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637 dSigma = dWidthEstimated/2.0;
00638
00639 return dSigma;
00640 }
00641
00642
00643
00644 void FittingFunctionSigmoidWithBackground::estimateInitialGuess(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY)
00645 {
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658 double dInitialPosition = estimateInitialPosition(vExperimentalDataX,vExperimentalDataY);
00659 double dInitialSigma = estimateInitialSigma(vExperimentalDataX,vExperimentalDataY);
00660 double dInitialHeight = estimateInitialHeight(vExperimentalDataX,vExperimentalDataY);
00661 double dInitialBackground = estimateInitialBackground(vExperimentalDataX,vExperimentalDataY);
00662
00663 gsl_vector_set(_vInitialGuess, 0, dInitialPosition);
00664 gsl_vector_set(_vInitialGuess, 1, dInitialSigma);
00665 gsl_vector_set(_vInitialGuess, 2, dInitialHeight);
00666 gsl_vector_set(_vInitialGuess, 3, dInitialBackground);
00667 }
00668
00669
00670 void FittingFunctionSigmoidWithBackground::initializeInitialsParameters(gsl_vector* vExperimentalDataX,gsl_vector* vExperimentalDataY,gsl_vector* vInitialGuess,bool bEstimate)
00671 {
00672 if (bEstimate)
00673 {
00674 estimateInitialGuess(vExperimentalDataX,vExperimentalDataY);
00675 }
00676 else
00677 {
00678 gsl_vector_set(_vInitialGuess,0,gsl_vector_get(vInitialGuess,0));
00679 gsl_vector_set(_vInitialGuess,1,gsl_vector_get(vInitialGuess,1));
00680 gsl_vector_set(_vInitialGuess,2,gsl_vector_get(vInitialGuess,2));
00681 gsl_vector_set(_vInitialGuess,3,gsl_vector_get(vInitialGuess,3));
00682 }
00683 }
00684
00685
00686 void FittingFunctionSigmoidWithBackground::estimateSigmoidIntersections()
00687 {
00688
00689 double xlow;
00690 double xhigh;
00691 double dBackground;
00692
00693
00694 dBackground = getParameter(3);
00695
00696
00697
00698 double ylow = dBackground;
00699 double yhigh = ylow + getParameter(2);
00700
00701 double dPosition = getParameter(0);
00702 double dWidth = getParameter(1);
00703 double dHeight = getParameter(2);
00704
00705
00706
00707
00708
00709
00710
00711
00712 double fa = computeValue(dPosition,dPosition,dWidth,dHeight,dBackground);
00713 double fprimea = computeDerivateValue(dPosition,dPosition,dWidth,dHeight);
00714
00715 xlow = (ylow - fa)/fprimea + dPosition;
00716 xhigh = (yhigh - fa)/fprimea + dPosition;
00717
00718 cout << "\n---------------------------------------------" << endl;
00719 cout << "Intersections" << endl;
00720 cout << "(xlow ,ylow )=(" << xlow << "," << ylow << ")" << endl;
00721 cout << "(xhigh,yhigh)=(" << xhigh << "," << yhigh << ")" << endl;
00722
00723 _dXLow = xlow;
00724 _dXHigh = xhigh;
00725 }
00726