67 struct gsl_function_struct
69 double (*
function) (
double x,
void * params);
73 #define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params)
77 double epsabs,
double epsrel,
78 double *
result,
double *abserr,
110 _epsAbs(config.epsRel()),
111 _epsRel(config.epsAbs())
125 _epsAbs(config.epsRel()),
126 _epsRel(config.epsAbs()),
179 oocoutE((
TObject*)0,
Eval) <<
"RooGaussKronrodIntegrator1D::setLimits: cannot override integrand's limits" << endl;
266 #define GSL_SUCCESS 0
267 #define GSL_EBADTOL 13
269 #define GSL_ERROR(a,b) oocoutE((TObject*)0,Eval) << "RooGaussKronrodIntegrator1D::integral() ERROR: " << a << endl ; return b ;
270 #define GSL_DBL_MIN 2.2250738585072014e-308
271 #define GSL_DBL_EPSILON 2.2204460492503131e-16
278 double epsabs,
double epsrel,
279 double *
result,
double *abserr,
285 static double rescale_error (
double err,
const double result_abs,
const double result_asc) ;
288 rescale_error (
double err,
const double result_abs,
const double result_asc)
292 if (result_asc != 0 && err != 0)
294 double scale =
TMath::Power((200 * err / result_asc), 1.5) ;
298 err = result_asc * scale ;
326 static const double x1[5] = {
327 0.973906528517171720077964012084452,
328 0.865063366688984510732096688423493,
329 0.679409568299024406234327365114874,
330 0.433395394129247190799265943165784,
331 0.148874338981631210884826001129720
335 static const double w10[5] = {
336 0.066671344308688137593568809893332,
337 0.149451349150580593145776339657697,
338 0.219086362515982043995534934228163,
339 0.269266719309996355091226921569469,
340 0.295524224714752870173892994651338
344 static const double x2[5] = {
345 0.995657163025808080735527280689003,
346 0.930157491355708226001207180059508,
347 0.780817726586416897063717578345042,
348 0.562757134668604683339000099272694,
349 0.294392862701460198131126603103866
353 static const double w21a[5] = {
354 0.032558162307964727478818972459390,
355 0.075039674810919952767043140916190,
356 0.109387158802297641899210590325805,
357 0.134709217311473325928054001771707,
358 0.147739104901338491374841515972068
362 static const double w21b[6] = {
363 0.011694638867371874278064396062192,
364 0.054755896574351996031381300244580,
365 0.093125454583697605535065465083366,
366 0.123491976262065851077958109831074,
367 0.142775938577060080797094273138717,
368 0.149445554002916905664936468389821
372 static const double x3[11] = {
373 0.999333360901932081394099323919911,
374 0.987433402908088869795961478381209,
375 0.954807934814266299257919200290473,
376 0.900148695748328293625099494069092,
377 0.825198314983114150847066732588520,
378 0.732148388989304982612354848755461,
379 0.622847970537725238641159120344323,
380 0.499479574071056499952214885499755,
381 0.364901661346580768043989548502644,
382 0.222254919776601296498260928066212,
383 0.074650617461383322043914435796506
387 static const double w43a[10] = {
388 0.016296734289666564924281974617663,
389 0.037522876120869501461613795898115,
390 0.054694902058255442147212685465005,
391 0.067355414609478086075553166302174,
392 0.073870199632393953432140695251367,
393 0.005768556059769796184184327908655,
394 0.027371890593248842081276069289151,
395 0.046560826910428830743339154433824,
396 0.061744995201442564496240336030883,
397 0.071387267268693397768559114425516
401 static const double w43b[12] = {
402 0.001844477640212414100389106552965,
403 0.010798689585891651740465406741293,
404 0.021895363867795428102523123075149,
405 0.032597463975345689443882222526137,
406 0.042163137935191811847627924327955,
407 0.050741939600184577780189020092084,
408 0.058379395542619248375475369330206,
409 0.064746404951445885544689259517511,
410 0.069566197912356484528633315038405,
411 0.072824441471833208150939535192842,
412 0.074507751014175118273571813842889,
413 0.074722147517403005594425168280423
417 static const double x4[22] = {
418 0.999902977262729234490529830591582,
419 0.997989895986678745427496322365960,
420 0.992175497860687222808523352251425,
421 0.981358163572712773571916941623894,
422 0.965057623858384619128284110607926,
423 0.943167613133670596816416634507426,
424 0.915806414685507209591826430720050,
425 0.883221657771316501372117548744163,
426 0.845710748462415666605902011504855,
427 0.803557658035230982788739474980964,
428 0.757005730685495558328942793432020,
429 0.706273209787321819824094274740840,
430 0.651589466501177922534422205016736,
431 0.593223374057961088875273770349144,
432 0.531493605970831932285268948562671,
433 0.466763623042022844871966781659270,
434 0.399424847859218804732101665817923,
435 0.329874877106188288265053371824597,
436 0.258503559202161551802280975429025,
437 0.185695396568346652015917141167606,
438 0.111842213179907468172398359241362,
439 0.037352123394619870814998165437704
443 static const double w87a[21] = {
444 0.008148377384149172900002878448190,
445 0.018761438201562822243935059003794,
446 0.027347451050052286161582829741283,
447 0.033677707311637930046581056957588,
448 0.036935099820427907614589586742499,
449 0.002884872430211530501334156248695,
450 0.013685946022712701888950035273128,
451 0.023280413502888311123409291030404,
452 0.030872497611713358675466394126442,
453 0.035693633639418770719351355457044,
454 0.000915283345202241360843392549948,
455 0.005399280219300471367738743391053,
456 0.010947679601118931134327826856808,
457 0.016298731696787335262665703223280,
458 0.021081568889203835112433060188190,
459 0.025370969769253827243467999831710,
460 0.029189697756475752501446154084920,
461 0.032373202467202789685788194889595,
462 0.034783098950365142750781997949596,
463 0.036412220731351787562801163687577,
464 0.037253875503047708539592001191226
468 static const double w87b[23] = {
469 0.000274145563762072350016527092881,
470 0.001807124155057942948341311753254,
471 0.004096869282759164864458070683480,
472 0.006758290051847378699816577897424,
473 0.009549957672201646536053581325377,
474 0.012329447652244853694626639963780,
475 0.015010447346388952376697286041943,
476 0.017548967986243191099665352925900,
477 0.019938037786440888202278192730714,
478 0.022194935961012286796332102959499,
479 0.024339147126000805470360647041454,
480 0.026374505414839207241503786552615,
481 0.028286910788771200659968002987960,
482 0.030052581128092695322521110347341,
483 0.031646751371439929404586051078883,
484 0.033050413419978503290785944862689,
485 0.034255099704226061787082821046821,
486 0.035262412660156681033782717998428,
487 0.036076989622888701185500318003895,
488 0.036698604498456094498018047441094,
489 0.037120549269832576114119958413599,
490 0.037334228751935040321235449094698,
491 0.037361073762679023410321241766599
498 double epsabs,
double epsrel,
499 double *
result,
double * abserr,
size_t * neval)
501 double fv1[5], fv2[5], fv3[5], fv4[5];
503 double res10, res21, res43, res87;
504 double result_kronrod, err ;
508 const double half_length = 0.5 * (b -
a);
509 const double abs_half_length =
fabs (half_length);
510 const double center = 0.5 * (b +
a);
515 if (epsabs <= 0 && (epsrel < 50 *
GSL_DBL_EPSILON || epsrel < 0.5e-28))
520 GSL_ERROR (
"tolerance cannot be acheived with given epsabs and epsrel",
527 res21 =
w21b[5] * f_center;
530 for (k = 0; k < 5; k++)
532 const double abscissa = half_length *
x1[k];
533 const double fval1 =
GSL_FN_EVAL(f, center + abscissa);
534 const double fval2 =
GSL_FN_EVAL(f, center - abscissa);
535 const double fval = fval1 + fval2;
536 res10 +=
w10[k] * fval;
537 res21 +=
w21a[k] * fval;
544 for (k = 0; k < 5; k++)
546 const double abscissa = half_length *
x2[k];
547 const double fval1 =
GSL_FN_EVAL(f, center + abscissa);
548 const double fval2 =
GSL_FN_EVAL(f, center - abscissa);
549 const double fval = fval1 + fval2;
550 res21 +=
w21b[k] * fval;
552 savfun[k + 5] = fval;
557 resabs *= abs_half_length ;
560 const double mean = 0.5 * res21;
562 resasc =
w21b[5] *
fabs (f_center - mean);
564 for (k = 0; k < 5; k++)
567 (
w21a[k] * (
fabs (fv1[k] - mean) +
fabs (fv2[k] - mean))
568 +
w21b[k] * (fabs (fv3[k] - mean) +
fabs (fv4[k] - mean)));
570 resasc *= abs_half_length ;
573 result_kronrod = res21 * half_length;
575 err =
rescale_error ((res21 - res10) * half_length, resabs, resasc) ;
579 if (err < epsabs || err < epsrel * fabs (result_kronrod))
581 * result = result_kronrod ;
589 res43 =
w43b[11] * f_center;
591 for (k = 0; k < 10; k++)
593 res43 += savfun[k] *
w43a[k];
596 for (k = 0; k < 11; k++)
598 const double abscissa = half_length *
x3[k];
599 const double fval = (
GSL_FN_EVAL(f, center + abscissa)
601 res43 += fval *
w43b[k];
602 savfun[k + 10] = fval;
607 result_kronrod = res43 * half_length;
608 err =
rescale_error ((res43 - res21) * half_length, resabs, resasc);
610 if (err < epsabs || err < epsrel * fabs (result_kronrod))
612 * result = result_kronrod ;
620 res87 =
w87b[22] * f_center;
622 for (k = 0; k < 21; k++)
624 res87 += savfun[k] *
w87a[k];
627 for (k = 0; k < 22; k++)
629 const double abscissa = half_length *
x4[k];
636 result_kronrod = res87 * half_length ;
638 err =
rescale_error ((res87 - res43) * half_length, resabs, resasc);
640 if (err < epsabs || err < epsrel * fabs (result_kronrod))
642 * result = result_kronrod ;
650 * result = result_kronrod ;
static double rescale_error(double err, const double result_abs, const double result_asc)
Bool_t _useIntegrandLimits
int gsl_integration_qng(const gsl_function *f, double a, double b, double epsabs, double epsrel, double *result, double *abserr, size_t *neval)
static const double w87a[21]
static const double w43b[12]
static const double w21b[6]
RooGaussKronrodIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
virtual Bool_t checkLimits() const
Check that our integration range is finite and otherwise return kFALSE.
Double_t _xmax
Lower integration bound.
UInt_t getDimension() const
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Clone integrator with given function and configuration. Needed for RooNumIntFactory.
static const double x2[5]
double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
Double_t _epsAbs
do not persist
static const double x4[22]
Bool_t initialize()
Perform one-time initialization of integrator.
static const double w10[5]
static const double w43a[10]
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static void registerIntegrator(RooNumIntFactory &fact)
Register RooGaussKronrodIntegrator1D, its parameters and capabilities with RooNumIntConfig.
void function(const char *name_, T fun, const char *docstring=0)
virtual Double_t integral(const Double_t *yvec=0)
Calculate and return integral.
virtual ~RooGaussKronrodIntegrator1D()
Destructor.
const RooAbsFunc * _function
virtual Double_t getMinLimit(UInt_t dimension) const =0
Double_t integrand(const Double_t x[]) const
static const double w21a[5]
static const double w87b[23]
static const double x1[5]
struct gsl_function_struct gsl_function
virtual Double_t getMaxLimit(UInt_t dimension) const =0
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
Mother of all ROOT objects.
const RooAbsFunc * integrand() const
Double_t * xvec(Double_t &xx)
friend double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
ClassImp(RooGaussKronrodIntegrator1D)
#define GSL_FN_EVAL(F, x)
Bool_t storeProtoIntegrator(RooAbsIntegrator *proto, const RooArgSet &defConfig, const char *depName="")
Method accepting registration of a prototype numeric integrator along with a RooArgSet of its default...
static const double x3[11]