69#define GSL_FN_EVAL(F,x) (*((F)->function))(x,(F)->params)
73 double epsabs,
double epsrel,
74 double *
result,
double *abserr,
105 oocoutI(
nullptr,
Integration) <<
"RooGaussKronrodIntegrator1D has been registered" << std::endl;
126 _epsAbs(config.epsRel()),
127 _epsRel(config.epsAbs())
141 _epsAbs(config.epsRel()),
142 _epsRel(config.epsAbs()),
195 oocoutE(
nullptr,
Eval) <<
"RooGaussKronrodIntegrator1D::setLimits: cannot override integrand's limits" << endl;
283#define GSL_EBADTOL 13
285#define GSL_ERROR(a,b) oocoutE(nullptr,Eval) << "RooGaussKronrodIntegrator1D::integral() ERROR: " << a << endl ; return b ;
286#define GSL_DBL_MIN 2.2250738585072014e-308
287#define GSL_DBL_EPSILON 2.2204460492503131e-16
294 double epsabs,
double epsrel,
295 double *
result,
double *abserr,
301static double rescale_error (
double err,
const double result_abs,
const double result_asc) ;
308 if (result_asc != 0 && err != 0)
310 double scale =
TMath::Power((200 * err / result_asc), 1.5) ;
314 err = result_asc * scale ;
342static const double x1[5] = {
343 0.973906528517171720077964012084452,
344 0.865063366688984510732096688423493,
345 0.679409568299024406234327365114874,
346 0.433395394129247190799265943165784,
347 0.148874338981631210884826001129720
351static const double w10[5] = {
352 0.066671344308688137593568809893332,
353 0.149451349150580593145776339657697,
354 0.219086362515982043995534934228163,
355 0.269266719309996355091226921569469,
356 0.295524224714752870173892994651338
360static const double x2[5] = {
361 0.995657163025808080735527280689003,
362 0.930157491355708226001207180059508,
363 0.780817726586416897063717578345042,
364 0.562757134668604683339000099272694,
365 0.294392862701460198131126603103866
370 0.032558162307964727478818972459390,
371 0.075039674810919952767043140916190,
372 0.109387158802297641899210590325805,
373 0.134709217311473325928054001771707,
374 0.147739104901338491374841515972068
379 0.011694638867371874278064396062192,
380 0.054755896574351996031381300244580,
381 0.093125454583697605535065465083366,
382 0.123491976262065851077958109831074,
383 0.142775938577060080797094273138717,
384 0.149445554002916905664936468389821
388static const double x3[11] = {
389 0.999333360901932081394099323919911,
390 0.987433402908088869795961478381209,
391 0.954807934814266299257919200290473,
392 0.900148695748328293625099494069092,
393 0.825198314983114150847066732588520,
394 0.732148388989304982612354848755461,
395 0.622847970537725238641159120344323,
396 0.499479574071056499952214885499755,
397 0.364901661346580768043989548502644,
398 0.222254919776601296498260928066212,
399 0.074650617461383322043914435796506
403static const double w43a[10] = {
404 0.016296734289666564924281974617663,
405 0.037522876120869501461613795898115,
406 0.054694902058255442147212685465005,
407 0.067355414609478086075553166302174,
408 0.073870199632393953432140695251367,
409 0.005768556059769796184184327908655,
410 0.027371890593248842081276069289151,
411 0.046560826910428830743339154433824,
412 0.061744995201442564496240336030883,
413 0.071387267268693397768559114425516
417static const double w43b[12] = {
418 0.001844477640212414100389106552965,
419 0.010798689585891651740465406741293,
420 0.021895363867795428102523123075149,
421 0.032597463975345689443882222526137,
422 0.042163137935191811847627924327955,
423 0.050741939600184577780189020092084,
424 0.058379395542619248375475369330206,
425 0.064746404951445885544689259517511,
426 0.069566197912356484528633315038405,
427 0.072824441471833208150939535192842,
428 0.074507751014175118273571813842889,
429 0.074722147517403005594425168280423
433static const double x4[22] = {
434 0.999902977262729234490529830591582,
435 0.997989895986678745427496322365960,
436 0.992175497860687222808523352251425,
437 0.981358163572712773571916941623894,
438 0.965057623858384619128284110607926,
439 0.943167613133670596816416634507426,
440 0.915806414685507209591826430720050,
441 0.883221657771316501372117548744163,
442 0.845710748462415666605902011504855,
443 0.803557658035230982788739474980964,
444 0.757005730685495558328942793432020,
445 0.706273209787321819824094274740840,
446 0.651589466501177922534422205016736,
447 0.593223374057961088875273770349144,
448 0.531493605970831932285268948562671,
449 0.466763623042022844871966781659270,
450 0.399424847859218804732101665817923,
451 0.329874877106188288265053371824597,
452 0.258503559202161551802280975429025,
453 0.185695396568346652015917141167606,
454 0.111842213179907468172398359241362,
455 0.037352123394619870814998165437704
459static const double w87a[21] = {
460 0.008148377384149172900002878448190,
461 0.018761438201562822243935059003794,
462 0.027347451050052286161582829741283,
463 0.033677707311637930046581056957588,
464 0.036935099820427907614589586742499,
465 0.002884872430211530501334156248695,
466 0.013685946022712701888950035273128,
467 0.023280413502888311123409291030404,
468 0.030872497611713358675466394126442,
469 0.035693633639418770719351355457044,
470 0.000915283345202241360843392549948,
471 0.005399280219300471367738743391053,
472 0.010947679601118931134327826856808,
473 0.016298731696787335262665703223280,
474 0.021081568889203835112433060188190,
475 0.025370969769253827243467999831710,
476 0.029189697756475752501446154084920,
477 0.032373202467202789685788194889595,
478 0.034783098950365142750781997949596,
479 0.036412220731351787562801163687577,
480 0.037253875503047708539592001191226
484static const double w87b[23] = {
485 0.000274145563762072350016527092881,
486 0.001807124155057942948341311753254,
487 0.004096869282759164864458070683480,
488 0.006758290051847378699816577897424,
489 0.009549957672201646536053581325377,
490 0.012329447652244853694626639963780,
491 0.015010447346388952376697286041943,
492 0.017548967986243191099665352925900,
493 0.019938037786440888202278192730714,
494 0.022194935961012286796332102959499,
495 0.024339147126000805470360647041454,
496 0.026374505414839207241503786552615,
497 0.028286910788771200659968002987960,
498 0.030052581128092695322521110347341,
499 0.031646751371439929404586051078883,
500 0.033050413419978503290785944862689,
501 0.034255099704226061787082821046821,
502 0.035262412660156681033782717998428,
503 0.036076989622888701185500318003895,
504 0.036698604498456094498018047441094,
505 0.037120549269832576114119958413599,
506 0.037334228751935040321235449094698,
507 0.037361073762679023410321241766599
514 double epsabs,
double epsrel,
515 double *
result,
double * abserr,
size_t * neval)
517 double fv1[5], fv2[5], fv3[5], fv4[5];
519 double res10, res21, res43, res87;
520 double result_kronrod, err ;
524 const double half_length = 0.5 * (
b -
a);
525 const double abs_half_length =
fabs (half_length);
526 const double center = 0.5 * (
b +
a);
531 if (epsabs <= 0 && (epsrel < 50 *
GSL_DBL_EPSILON || epsrel < 0.5e-28))
536 GSL_ERROR (
"tolerance cannot be acheived with given epsabs and epsrel",
543 res21 =
w21b[5] * f_center;
546 for (k = 0; k < 5; k++)
548 const double abscissa = half_length *
x1[k];
551 const double fval = fval1 + fval2;
552 res10 +=
w10[k] * fval;
553 res21 +=
w21a[k] * fval;
560 for (k = 0; k < 5; k++)
562 const double abscissa = half_length *
x2[k];
565 const double fval = fval1 + fval2;
566 res21 +=
w21b[k] * fval;
568 savfun[k + 5] = fval;
573 resabs *= abs_half_length ;
576 const double mean = 0.5 * res21;
578 resasc =
w21b[5] *
fabs (f_center - mean);
580 for (k = 0; k < 5; k++)
583 (
w21a[k] * (
fabs (fv1[k] - mean) +
fabs (fv2[k] - mean))
584 +
w21b[k] * (
fabs (fv3[k] - mean) +
fabs (fv4[k] - mean)));
586 resasc *= abs_half_length ;
589 result_kronrod = res21 * half_length;
591 err =
rescale_error ((res21 - res10) * half_length, resabs, resasc) ;
595 if (err < epsabs || err < epsrel *
fabs (result_kronrod))
597 *
result = result_kronrod ;
605 res43 =
w43b[11] * f_center;
607 for (k = 0; k < 10; k++)
609 res43 += savfun[k] *
w43a[k];
612 for (k = 0; k < 11; k++)
614 const double abscissa = half_length *
x3[k];
617 res43 += fval *
w43b[k];
618 savfun[k + 10] = fval;
623 result_kronrod = res43 * half_length;
624 err =
rescale_error ((res43 - res21) * half_length, resabs, resasc);
626 if (err < epsabs || err < epsrel *
fabs (result_kronrod))
628 *
result = result_kronrod ;
636 res87 =
w87b[22] * f_center;
638 for (k = 0; k < 21; k++)
640 res87 += savfun[k] *
w87a[k];
643 for (k = 0; k < 22; k++)
645 const double abscissa = half_length *
x4[k];
652 result_kronrod = res87 * half_length ;
654 err =
rescale_error ((res87 - res43) * half_length, resabs, resasc);
656 if (err < epsabs || err < epsrel *
fabs (result_kronrod))
658 *
result = result_kronrod ;
666 *
result = result_kronrod ;
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)
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
virtual double getMaxLimit(UInt_t dimension) const =0
virtual double getMinLimit(UInt_t dimension) const =0
UInt_t getDimension() const
RooAbsIntegrator is the abstract interface for integrators of real-valued functions that implement th...
bool isValid() const
Is integrator in valid state.
const RooAbsFunc * _function
Pointer to function binding of integrand.
const RooAbsFunc * integrand() const
Return integrand function binding.
bool _valid
Is integrator in valid state?
RooArgSet is a container object that can hold multiple RooAbsArg objects.
RooGaussKronrodIntegrator1D implements the Gauss-Kronrod integration algorithm.
double integral(const double *yvec=0) override
Calculate and return integral.
double _epsAbs
do not persist
double _xmax
Lower integration bound.
~RooGaussKronrodIntegrator1D() override
Destructor.
friend double RooGaussKronrodIntegrator1D_GSL_GlueFunction(double x, void *data)
bool initialize()
Perform one-time initialization of integrator.
RooGaussKronrodIntegrator1D()
coverity[UNINIT_CTOR] Default constructor
bool checkLimits() const override
Check that our integration range is finite and otherwise return false.
RooAbsIntegrator * clone(const RooAbsFunc &function, const RooNumIntConfig &config) const override
Clone integrator with given function and configuration. Needed for RooNumIntFactory.
bool setLimits(double *xmin, double *xmax) override
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 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.
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)
static void registerIntegrator()
double(* function)(double x, void *params)