67struct gsl_function_struct
73#define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params)
77 double epsabs,
double epsrel,
78 double *result,
double *abserr,
88 static void registerIntegrator()
95struct Roo_reg_GKInteg1D {
96 Roo_reg_GKInteg1D() { Roo_internal_GKInteg1D::registerIntegrator(); }
130 _epsAbs(config.epsRel()),
131 _epsRel(config.epsAbs())
145 _epsAbs(config.epsRel()),
146 _epsRel(config.epsAbs()),
199 oocoutE((
TObject*)0,
Eval) <<
"RooGaussKronrodIntegrator1D::setLimits: cannot override integrand's limits" << endl;
253 double result, error;
287#define GSL_EBADTOL 13
289#define GSL_ERROR(a,b) oocoutE((TObject*)0,Eval) << "RooGaussKronrodIntegrator1D::integral() ERROR: " << a << endl ; return b ;
290#define GSL_DBL_MIN 2.2250738585072014e-308
291#define GSL_DBL_EPSILON 2.2204460492503131e-16
298 double epsabs,
double epsrel,
299 double *result,
double *abserr,
305static double rescale_error (
double err,
const double result_abs,
const double result_asc) ;
312 if (result_asc != 0 && err != 0)
314 double scale =
TMath::Power((200 * err / result_asc), 1.5) ;
318 err = result_asc * scale ;
346static const double x1[5] = {
347 0.973906528517171720077964012084452,
348 0.865063366688984510732096688423493,
349 0.679409568299024406234327365114874,
350 0.433395394129247190799265943165784,
351 0.148874338981631210884826001129720
355static const double w10[5] = {
356 0.066671344308688137593568809893332,
357 0.149451349150580593145776339657697,
358 0.219086362515982043995534934228163,
359 0.269266719309996355091226921569469,
360 0.295524224714752870173892994651338
364static const double x2[5] = {
365 0.995657163025808080735527280689003,
366 0.930157491355708226001207180059508,
367 0.780817726586416897063717578345042,
368 0.562757134668604683339000099272694,
369 0.294392862701460198131126603103866
374 0.032558162307964727478818972459390,
375 0.075039674810919952767043140916190,
376 0.109387158802297641899210590325805,
377 0.134709217311473325928054001771707,
378 0.147739104901338491374841515972068
383 0.011694638867371874278064396062192,
384 0.054755896574351996031381300244580,
385 0.093125454583697605535065465083366,
386 0.123491976262065851077958109831074,
387 0.142775938577060080797094273138717,
388 0.149445554002916905664936468389821
392static const double x3[11] = {
393 0.999333360901932081394099323919911,
394 0.987433402908088869795961478381209,
395 0.954807934814266299257919200290473,
396 0.900148695748328293625099494069092,
397 0.825198314983114150847066732588520,
398 0.732148388989304982612354848755461,
399 0.622847970537725238641159120344323,
400 0.499479574071056499952214885499755,
401 0.364901661346580768043989548502644,
402 0.222254919776601296498260928066212,
403 0.074650617461383322043914435796506
407static const double w43a[10] = {
408 0.016296734289666564924281974617663,
409 0.037522876120869501461613795898115,
410 0.054694902058255442147212685465005,
411 0.067355414609478086075553166302174,
412 0.073870199632393953432140695251367,
413 0.005768556059769796184184327908655,
414 0.027371890593248842081276069289151,
415 0.046560826910428830743339154433824,
416 0.061744995201442564496240336030883,
417 0.071387267268693397768559114425516
421static const double w43b[12] = {
422 0.001844477640212414100389106552965,
423 0.010798689585891651740465406741293,
424 0.021895363867795428102523123075149,
425 0.032597463975345689443882222526137,
426 0.042163137935191811847627924327955,
427 0.050741939600184577780189020092084,
428 0.058379395542619248375475369330206,
429 0.064746404951445885544689259517511,
430 0.069566197912356484528633315038405,
431 0.072824441471833208150939535192842,
432 0.074507751014175118273571813842889,
433 0.074722147517403005594425168280423
437static const double x4[22] = {
438 0.999902977262729234490529830591582,
439 0.997989895986678745427496322365960,
440 0.992175497860687222808523352251425,
441 0.981358163572712773571916941623894,
442 0.965057623858384619128284110607926,
443 0.943167613133670596816416634507426,
444 0.915806414685507209591826430720050,
445 0.883221657771316501372117548744163,
446 0.845710748462415666605902011504855,
447 0.803557658035230982788739474980964,
448 0.757005730685495558328942793432020,
449 0.706273209787321819824094274740840,
450 0.651589466501177922534422205016736,
451 0.593223374057961088875273770349144,
452 0.531493605970831932285268948562671,
453 0.466763623042022844871966781659270,
454 0.399424847859218804732101665817923,
455 0.329874877106188288265053371824597,
456 0.258503559202161551802280975429025,
457 0.185695396568346652015917141167606,
458 0.111842213179907468172398359241362,
459 0.037352123394619870814998165437704
463static const double w87a[21] = {
464 0.008148377384149172900002878448190,
465 0.018761438201562822243935059003794,
466 0.027347451050052286161582829741283,
467 0.033677707311637930046581056957588,
468 0.036935099820427907614589586742499,
469 0.002884872430211530501334156248695,
470 0.013685946022712701888950035273128,
471 0.023280413502888311123409291030404,
472 0.030872497611713358675466394126442,
473 0.035693633639418770719351355457044,
474 0.000915283345202241360843392549948,
475 0.005399280219300471367738743391053,
476 0.010947679601118931134327826856808,
477 0.016298731696787335262665703223280,
478 0.021081568889203835112433060188190,
479 0.025370969769253827243467999831710,
480 0.029189697756475752501446154084920,
481 0.032373202467202789685788194889595,
482 0.034783098950365142750781997949596,
483 0.036412220731351787562801163687577,
484 0.037253875503047708539592001191226
488static const double w87b[23] = {
489 0.000274145563762072350016527092881,
490 0.001807124155057942948341311753254,
491 0.004096869282759164864458070683480,
492 0.006758290051847378699816577897424,
493 0.009549957672201646536053581325377,
494 0.012329447652244853694626639963780,
495 0.015010447346388952376697286041943,
496 0.017548967986243191099665352925900,
497 0.019938037786440888202278192730714,
498 0.022194935961012286796332102959499,
499 0.024339147126000805470360647041454,
500 0.026374505414839207241503786552615,
501 0.028286910788771200659968002987960,
502 0.030052581128092695322521110347341,
503 0.031646751371439929404586051078883,
504 0.033050413419978503290785944862689,
505 0.034255099704226061787082821046821,
506 0.035262412660156681033782717998428,
507 0.036076989622888701185500318003895,
508 0.036698604498456094498018047441094,
509 0.037120549269832576114119958413599,
510 0.037334228751935040321235449094698,
511 0.037361073762679023410321241766599
518 double epsabs,
double epsrel,
519 double * result,
double * abserr,
size_t * neval)
521 double fv1[5], fv2[5], fv3[5], fv4[5];
523 double res10, res21, res43, res87;
524 double result_kronrod, err ;
528 const double half_length = 0.5 * (
b -
a);
529 const double abs_half_length =
fabs (half_length);
530 const double center = 0.5 * (
b +
a);
535 if (epsabs <= 0 && (epsrel < 50 *
GSL_DBL_EPSILON || epsrel < 0.5e-28))
540 GSL_ERROR (
"tolerance cannot be acheived with given epsabs and epsrel",
547 res21 =
w21b[5] * f_center;
550 for (k = 0; k < 5; k++)
552 const double abscissa = half_length *
x1[k];
555 const double fval = fval1 + fval2;
556 res10 +=
w10[k] * fval;
557 res21 +=
w21a[k] * fval;
564 for (k = 0; k < 5; k++)
566 const double abscissa = half_length *
x2[k];
569 const double fval = fval1 + fval2;
570 res21 +=
w21b[k] * fval;
572 savfun[k + 5] = fval;
577 resabs *= abs_half_length ;
580 const double mean = 0.5 * res21;
582 resasc =
w21b[5] *
fabs (f_center - mean);
584 for (k = 0; k < 5; k++)
587 (
w21a[k] * (
fabs (fv1[k] - mean) +
fabs (fv2[k] - mean))
588 +
w21b[k] * (
fabs (fv3[k] - mean) +
fabs (fv4[k] - mean)));
590 resasc *= abs_half_length ;
593 result_kronrod = res21 * half_length;
595 err =
rescale_error ((res21 - res10) * half_length, resabs, resasc) ;
599 if (err < epsabs || err < epsrel *
fabs (result_kronrod))
601 * result = result_kronrod ;
609 res43 =
w43b[11] * f_center;
611 for (k = 0; k < 10; k++)
613 res43 += savfun[k] *
w43a[k];
616 for (k = 0; k < 11; k++)
618 const double abscissa = half_length *
x3[k];
621 res43 += fval *
w43b[k];
622 savfun[k + 10] = fval;
627 result_kronrod = res43 * half_length;
628 err =
rescale_error ((res43 - res21) * half_length, resabs, resasc);
630 if (err < epsabs || err < epsrel *
fabs (result_kronrod))
632 * result = result_kronrod ;
640 res87 =
w87b[22] * f_center;
642 for (k = 0; k < 21; k++)
644 res87 += savfun[k] *
w87a[k];
647 for (k = 0; k < 22; k++)
649 const double abscissa = half_length *
x4[k];
656 result_kronrod = res87 * half_length ;
658 err =
rescale_error ((res87 - res43) * half_length, resabs, resasc);
660 if (err < epsabs || err < epsrel *
fabs (result_kronrod))
662 * result = result_kronrod ;
670 * result = result_kronrod ;
struct gsl_function_struct gsl_function
static const double x2[5]
static double rescale_error(double err, const double result_abs, const double result_asc)
static const double w10[5]
static const double w43a[10]
static const double x4[22]
static const double w21b[6]
static const double x1[5]
double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
static const double w87a[21]
static const double w21a[5]
static const double w43b[12]
static const double x3[11]
#define GSL_FN_EVAL(F, x)
static const double w87b[23]
int gsl_integration_qng(const gsl_function *f, double a, double b, double epsabs, double epsrel, double *result, double *abserr, size_t *neval)
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
virtual Double_t getMinLimit(UInt_t dimension) const =0
virtual Double_t getMaxLimit(UInt_t dimension) const =0
UInt_t getDimension() const
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
const RooAbsFunc * _function
const RooAbsFunc * integrand() const
RooArgSet is a container object that can hold multiple RooAbsArg objects.
RooGaussKronrodIntegrator1D implements the Gauss-Kronrod integration algorithm.
Double_t _epsAbs
do not persist
friend double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
Bool_t initialize()
Perform one-time initialization of integrator.
RooGaussKronrodIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
Double_t _xmax
Lower integration bound.
virtual ~RooGaussKronrodIntegrator1D()
Destructor.
virtual RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const
Clone integrator with given function and configuration. Needed for RooNumIntFactory.
virtual Bool_t checkLimits() const
Check that our integration range is finite and otherwise return kFALSE.
Bool_t _useIntegrandLimits
virtual Double_t integral(const Double_t *yvec=0)
Calculate and return integral.
Bool_t setLimits(Double_t *xmin, Double_t *xmax)
Change our integration limits.
static void registerIntegrator(RooNumIntFactory &fact)
Register RooGaussKronrodIntegrator1D, its parameters and capabilities with RooNumIntConfig.
RooNumIntConfig holds the configuration parameters of the various numeric integrators used by RooReal...
RooNumIntFactory is a factory to instantiate numeric integrators from a given function binding and a ...
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 RooNumIntFactory & instance()
Static method returning reference to singleton instance of factory.
Mother of all ROOT objects.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
static Roo_reg_GKInteg1D instance
static Roo_reg_AGKInteg1D instance
LongDouble_t Power(LongDouble_t x, LongDouble_t y)