RooHypatia2 GSL errors

ROOT
1

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

Please read tips for efficient and successful posting and posting code

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

2

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

3

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

4

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

5

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