# RooHypatia2 GSL errors

Dear ROoFit Experts,

I am running some fits using the Hypatia2 function ( ported to an older release of ROOT ( 6/18.04)
I copied most of the methods and I have enabled the analytical integral .
When running the fit i get

Error in <GSLError>: Error 24 in hyperg_2F1.c at 773 : error

Which boils down to come from this call

double Gauss2F1(double a, double b, double c, double x){
if (fabs(x) <= 1.) {
return ROOT::Math::hyperg(a, b, c, x);
} else {
return ROOT::Math::hyperg(c-a, b, c, 1-1/(1-x))/std::pow(1-x, b);
}
}

Is there any workaround i can do?
Thanks
Renato

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided

I realized that the code implemented (Commented out) in the ROOT/master is not compatible to a version of RooIpatia2 i have locally, In particular the files i am using look like this :

Double_t RooIpatia2::evaluate() const {
Double_t d     = x - mu;
Double_t cons0 = TMath::Sqrt(zeta);
Double_t alpha, beta, delta, delta2, cons1, phi, A, B, k1, k2;
Double_t asigma  = a * sigma;
Double_t a2sigma = a2 * sigma;
Double_t out     = 0.;
if (zeta != 0.) {
phi   = RooIpatia2__BK(l + 1., zeta) / RooIpatia2__BK(l, zeta);   // careful if zeta -> 0. You can implement a function for the ratio, but carefull again that |nu + 1 | != |nu| + 1 so you jave to deal wiht the signs
cons1 = sigma / TMath::Sqrt(phi);
alpha = cons0 / cons1;   //*TMath::Sqrt((1 - fb*fb));
beta  = fb;              //*alpha;
delta = cons0 * cons1;

if (d < -asigma) {
// printf("-_-\n");
// printf("alpha %e\n",alpha);
// printf("beta %e\n",beta);
// printf("delta %e\n",delta);

k1 = RooIpatia2__LogEval(-asigma, l, alpha, beta, delta);
k2 = RooIpatia2__diff_eval(-asigma, l, alpha, beta, delta);
B  = -asigma + n * k1 / k2;
A  = k1 * TMath::Power(B + asigma, n);
// printf("k1 is %e\n",k1);
// printf("k2 is %e\n",k2);
// printf("A is%e\n",A);
// printf("B is%e\n",B);
out = A * TMath::Power(B - d, -n);
} else if (d > a2sigma) {
// printf("uoeo\n");
k1 = RooIpatia2__LogEval(a2sigma, l, alpha, beta, delta);
k2 = RooIpatia2__diff_eval(a2sigma, l, alpha, beta, delta);

B = -a2sigma - n2 * k1 / k2;

A = k1 * TMath::Power(B + a2sigma, n2);

out = A * TMath::Power(B + d, -n2);
} else {
// printf("HERE\n");
out = RooIpatia2__LogEval(d, l, alpha, beta, delta);
}
} else if (l < 0.) {
beta  = fb;
cons1 = -2. * l;
// delta = sigma;
if (l <= -1.0) {
delta = sigma * sqrt(-2 + cons1);
} else {
printf("WARNING: zeta ==0 and l > -1 ==> not defined rms. Changing the meaning of sigma, but I keep fitting anyway\n");
delta = sigma;
}
delta2 = delta * delta;
if (d < -asigma) {
cons1 = TMath::Exp(-beta * asigma);
phi   = 1. + asigma * asigma / delta2;
k1    = cons1 * TMath::Power(phi, l - 0.5);
k2    = beta * k1 - cons1 * (l - 0.5) * TMath::Power(phi, l - 1.5) * 2 * asigma / delta2;
B     = -asigma + n * k1 / k2;
A     = k1 * TMath::Power(B + asigma, n);
out   = A * TMath::Power(B - d, -n);
} else if (d > a2sigma) {
cons1 = TMath::Exp(beta * a2sigma);
phi   = 1. + a2sigma * a2sigma / delta2;
k1    = cons1 * TMath::Power(phi, l - 0.5);
k2    = beta * k1 + cons1 * (l - 0.5) * TMath::Power(phi, l - 1.5) * 2. * a2sigma / delta2;
B     = -a2sigma - n2 * k1 / k2;
A     = k1 * TMath::Power(B + a2sigma, n2);
out   = A * TMath::Power(B + d, -n2);
} else {
out = TMath::Exp(beta * d) * TMath::Power(1. + d * d / delta2, l - 0.5);
}
} else {
printf("zeta = 0 only suported for l < 0, while l = %e\n",0);
}
// printf("result is %e\n",out);
return out;
}

Double_t RooIpatia2::analyticalIntegral(Int_t code, const char* rangeName) const
{
// Double_t d = x-mu;
//Double_t alpha, delta, delta2, cons1, phi, A, B, k1, k2;
Double_t delta, delta2, cons1, phi, A, B, k1, k2;
Double_t asigma = a*sigma;
Double_t a2sigma = a2*sigma;
Double_t out = 0.;
Double_t I0,I1;//,I2;
Double_t I1a = 0;
Double_t I1b = 0;
Double_t beta = 0;

Double_t d0 = x.min(rangeName)-mu;
Double_t d1 = x.max(rangeName)-mu;
if (code==1) {
cons1 = -2.*l;
if (l<=-1.0) { delta = sigma *sqrt(-2+cons1);}
else {delta = sigma;}
delta2 = delta*delta;
//printf("Here ^_^ \n");

if ((d0 > -asigma) && (d1 < a2sigma)){ return  stIntegral(d1,delta, l) - stIntegral(d0,delta, l);}

if (d0 > a2sigma) {
cons1 = 1.;
phi = 1. + a2sigma*a2sigma/delta2;
k1 = cons1*TMath::Power(phi,l-0.5);
k2 = beta*k1+ cons1*(l-0.5)*TMath::Power(phi,l-1.5)*2.*a2sigma/delta2;
B = -a2sigma - n2*k1/k2;
A = k1*TMath::Power(B+a2sigma,n2);
return A*(TMath::Power(B+d1,1-n2)/(1-n2) -TMath::Power(B+d0,1-n2)/(1-n2) ) ;
}

if (d1 < -asigma) {
cons1 = 1.;
phi = 1. + asigma*asigma/delta2;
k1 = cons1*TMath::Power(phi,l-0.5);
k2 = beta*k1- cons1*(l-0.5)*TMath::Power(phi,l-1.5)*2*asigma/delta2;
B = -asigma + n*k1/k2;
A = k1*TMath::Power(B+asigma,n);
I0 = A*TMath::Power(B-d0,1-n)/(n-1);
I1 = A*TMath::Power(B-d1,1-n)/(n-1);
return I1 - I0;
}

if (d0 <-asigma) {
cons1 = 1.;
phi = 1. + asigma*asigma/delta2;
k1 = cons1*TMath::Power(phi,l-0.5);
k2 = beta*k1- cons1*(l-0.5)*TMath::Power(phi,l-1.5)*2*asigma/delta2;
B = -asigma + n*k1/k2;
A = k1*TMath::Power(B+asigma,n);
I0 = A*TMath::Power(B-d0,1-n)/(n-1);
I1a = A*TMath::Power(B+asigma,1-n)/(n-1) - stIntegral(-asigma,delta, l);
}
else {  I0 = stIntegral(d0,delta, l);}

if (d1 > a2sigma) {
cons1 = 1.;
phi = 1. + a2sigma*a2sigma/delta2;
k1 = cons1*TMath::Power(phi,l-0.5);
k2 = beta*k1+ cons1*(l-0.5)*TMath::Power(phi,l-1.5)*2.*a2sigma/delta2;
B = -a2sigma - n2*k1/k2;
A = k1*TMath::Power(B+a2sigma,n2);
I1b = A*(TMath::Power(B+d1,1-n2)/(1-n2) -TMath::Power(B+a2sigma,1-n2)/(1-n2) ) - stIntegral(d1,delta, l) +  stIntegral(a2sigma,delta, l) ;
}
I1 = stIntegral(d1,delta, l) + I1a + I1b;
return I1 - I0;
// return (x.max(rangeName)-x.min(rangeName))
}
return out;
}

By comparison to the root master i noticed that the commented out code in evaluate() and getAnalyticalIntegral are not compatible in the part of code

else if (l < 0.) {
beta  = fb;
cons1 = -2. * l;
// delta = sigma;
if (l <= -1.0) {
delta = sigma * sqrt(-2 + cons1);
} else {
printf("WARNING: zeta ==0 and l > -1 ==> not defined rms. Changing the meaning of sigma, but I keep fitting anyway\n");
delta = sigma;
}

@moneta do you want me to send you the Ipatia implementation i am using for further debugging?
I believe the RooHyPatia2 which has been ported to ROOT/master migth be an outdated one

Hi,

Can you please open a GitHub issue with this and attach the code (or a link to the code).
This will be great so we will investigate and not forget about it

Thank you

Lorenzo

@moneta here the issue and the zip file with the version of the code i use.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.