Davismt2.cc

Go to the documentation of this file.

// #ifndef MT2_BISECT_C
// #define MT2_BISECT_C
/***********************************************************************/
/*                                                                     */
/*              Finding mt2 by Bisection                               */
/*                                                                     */
/*              Authors: Hsin-Chia Cheng, Zhenyu Han                   */ 
/*              Dec 11, 2008, v1.01a                                   */
/*                                                                     */  
/***********************************************************************/


/*******************************************************************************
Usage: 

1. Define an object of type "mt2":

mt2_bisect::mt2 mt2_event;

2. Set momenta and the mass of the invisible particle, mn:

mt2_event.set_momenta( pa, pb, pmiss );
mt2_event.set_mass( mn );

where array pa[0..2], pb[0..2], pmiss[0..2] contains (mass,px,py) 
    for the visible particles and the missing momentum. pmiss[0] is not used. 
    All quantities are given in double.    

    3. Use Davismt2::get_mt2() to obtain the value of mt2:

double mt2_value = mt2_event.get_mt2();       

*******************************************************************************/ 

#include "PhysicsTools/Heppy/interface/Davismt2.h"
// ClassImp(Davismt2);

using namespace std;

namespace heppy {

/*The user can change the desired precision below, the larger one of the following two definitions is used. Relative precision less than 0.00001 is not guaranteed to be 
  achievable--use with caution*/ 

const float Davismt2::RELATIVE_PRECISION = 0.00001;
const float Davismt2::ABSOLUTE_PRECISION = 0.0;
const float Davismt2::ZERO_MASS = 0.0; 
const float Davismt2::MIN_MASS = 0.1; 
const float Davismt2::SCANSTEP = 0.1;

Davismt2::Davismt2(){
    solved = false;
    momenta_set = false;
    mt2_b  = 0.;
    scale = 1.;
    verbose = 1;
}

Davismt2::~Davismt2(){}

double Davismt2::get_mt2(){
    if (!momenta_set)
    {
        cout << "Davismt2::get_mt2() ==> Please set momenta first!" << endl;
        return 0;
    }

    if (!solved) mt2_bisect();
    return mt2_b*scale;
}

void Davismt2::set_momenta(double* pa0, double* pb0, double* pmiss0){
    solved = false;     //reset solved tag when momenta are changed.
    momenta_set = true;

    ma = fabs(pa0[0]);  // mass cannot be negative

    if (ma < ZERO_MASS) ma = ZERO_MASS;

    pax  = pa0[1]; 
    pay  = pa0[2];
    masq = ma*ma;
    Easq = masq+pax*pax+pay*pay;
    Ea   = sqrt(Easq);

    mb = fabs(pb0[0]);

    if (mb < ZERO_MASS) mb = ZERO_MASS;

    pbx  = pb0[1]; 
    pby  = pb0[2];
    mbsq = mb*mb;
    Ebsq = mbsq+pbx*pbx+pby*pby;
    Eb   = sqrt(Ebsq);

    pmissx   = pmiss0[1]; pmissy = pmiss0[2];
    pmissxsq = pmissx*pmissx;
    pmissysq = pmissy*pmissy;

// set ma>= mb
    if(masq < mbsq)
    {
        double temp;
        temp = pax;  pax  = pbx;  pbx  = temp;
        temp = pay;  pay  = pby;  pby  = temp;
        temp = Ea;   Ea   = Eb;   Eb   = temp;
        temp = Easq; Easq = Ebsq; Ebsq = temp;
        temp = masq; masq = mbsq; mbsq = temp;
        temp = ma;   ma   = mb;   mb   = temp;   
    }
//normalize max{Ea, Eb} to 100
    if (Ea > Eb) scale = Ea/100.;
    else scale = Eb/100.;

    if (sqrt(pmissxsq+pmissysq)/100 > scale) scale = sqrt(pmissxsq+pmissysq)/100;
    //scale = 1;
    double scalesq = scale * scale;
    ma  = ma/scale;
    mb  = mb/scale;
    masq = masq/scalesq;
    mbsq = mbsq/scalesq;
    pax = pax/scale; pay = pay/scale;
    pbx = pbx/scale; pby = pby/scale;
    Ea  = Ea/scale;  Eb = Eb/scale;

    Easq = Easq/scalesq;
    Ebsq = Ebsq/scalesq;
    pmissx = pmissx/scale;
    pmissy = pmissy/scale;
    pmissxsq = pmissxsq/scalesq;
    pmissysq = pmissysq/scalesq;
    mn   = mn_unscale/scale; 
    mnsq = mn*mn;

    if (ABSOLUTE_PRECISION > 100.*RELATIVE_PRECISION) precision = ABSOLUTE_PRECISION;
    else precision = 100.*RELATIVE_PRECISION;
}

void Davismt2::set_mn(double mn0){
    solved = false;    //reset solved tag when mn is changed.
    mn_unscale   = fabs(mn0);  //mass cannot be negative
    mn = mn_unscale/scale;
    mnsq = mn*mn;
}

void Davismt2::print(){
    cout << " Davismt2::print() ==> pax = " << pax*scale << ";   pay = " << pay*scale << ";   ma = " << ma*scale <<";"<< endl;
    cout << " Davismt2::print() ==> pbx = " << pbx*scale << ";   pby = " << pby*scale << ";   mb = " << mb*scale <<";"<< endl;
    cout << " Davismt2::print() ==> pmissx = " << pmissx*scale << ";   pmissy = " << pmissy*scale <<";"<< endl;
    cout << " Davismt2::print() ==> mn = " << mn_unscale<<";" << endl;
}

//special case, the visible particle is massless
void Davismt2::mt2_massless(){

//rotate so that pay = 0 
    double theta,s,c;
    theta = atan(pay/pax);
    s = sin(theta);
    c = cos(theta);

    double pxtemp,pytemp;
    Easq   = pax*pax+pay*pay;
    Ebsq   = pbx*pbx+pby*pby;
    Ea     = sqrt(Easq);
    Eb     = sqrt(Ebsq);

    pxtemp = pax*c+pay*s;
    pax    = pxtemp;
    pay    = 0;
    pxtemp = pbx*c+pby*s;
    pytemp = -s*pbx+c*pby;
    pbx    = pxtemp;
    pby    = pytemp;
    pxtemp = pmissx*c+pmissy*s;
    pytemp = -s*pmissx+c*pmissy;
    pmissx = pxtemp;
    pmissy = pytemp;

    a2  = 1-pbx*pbx/(Ebsq);
    b2  = -pbx*pby/(Ebsq);
    c2  = 1-pby*pby/(Ebsq);

    d21 = (Easq*pbx)/Ebsq;
    d20 = - pmissx +  (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
    e21 = (Easq*pby)/Ebsq;
    e20 = - pmissy +  (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
    f22 = -(Easq*Easq/Ebsq);
    f21 = -2*Easq*(pbx*pmissx + pby*pmissy)/Ebsq;
    f20 = mnsq + pmissxsq + pmissysq - (pbx*pmissx + pby*pmissy)*(pbx*pmissx + pby*pmissy)/Ebsq;

    double Deltasq0    = 0; 
    double Deltasq_low, Deltasq_high;
    int    nsols_high, nsols_low;

    Deltasq_low = Deltasq0 + precision;
    nsols_low = nsols_massless(Deltasq_low);

    if(nsols_low > 1) 
    { 
        mt2_b = (double) sqrt(Deltasq0+mnsq);
        return;
    }

/*   
    if( nsols_massless(Deltasq_high) > 0 )
    {
        mt2_b = (double) sqrt(mnsq+Deltasq0);
        return;
        }*/

//look for when both parablos contain origin  
    double Deltasq_high1, Deltasq_high2;
    Deltasq_high1 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy;
    Deltasq_high2 = 2*Ea*mn;

    if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high2;
    else Deltasq_high = Deltasq_high1;

    nsols_high=nsols_massless(Deltasq_high);

    int foundhigh;
    if (nsols_high == nsols_low)
    {


        foundhigh=0;

        double minmass, maxmass;
        minmass  = mn ;
        maxmass  = sqrt(mnsq + Deltasq_high);
        for(double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
        {
            Deltasq_high = mass*mass - mnsq;

            nsols_high = nsols_massless(Deltasq_high);
            if(nsols_high>0)
            {
                foundhigh=1;
                Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
                break;
            }
        }
        if(foundhigh==0) 
        {

            if(verbose > 0) cout << "Davismt2::mt2_massless() ==> Deltasq_high not found at event " << nevt <<endl;


            mt2_b = (double)sqrt(Deltasq_low+mnsq);
            return;
        }
    }

    if(nsols_high == nsols_low)
    { 
        if(verbose > 0) cout << "Davismt2::mt2_massless() ==> error: nsols_low=nsols_high=" << nsols_high << endl;
        if(verbose > 0) cout << "Davismt2::mt2_massless() ==> Deltasq_high=" << Deltasq_high << endl;
        if(verbose > 0) cout << "Davismt2::mt2_massless() ==> Deltasq_low= "<< Deltasq_low << endl;

        mt2_b = sqrt(mnsq + Deltasq_low);
        return;
    }
    double minmass, maxmass;
    minmass = sqrt(Deltasq_low+mnsq);
    maxmass = sqrt(Deltasq_high+mnsq);
    while(maxmass - minmass > precision)
    {
        double Delta_mid, midmass, nsols_mid;
        midmass   = (minmass+maxmass)/2.;
        Delta_mid = midmass * midmass - mnsq;
        nsols_mid = nsols_massless(Delta_mid);
        if(nsols_mid != nsols_low) maxmass = midmass;
        if(nsols_mid == nsols_low) minmass = midmass;
    }
    mt2_b = minmass;
    return;
}

int Davismt2::nsols_massless(double Dsq){
    double delta;
    delta = Dsq/(2*Easq);
    d1    = d11*delta;
    e1    = e11*delta;
    f1    = f12*delta*delta+f10;
    d2    = d21*delta+d20;
    e2    = e21*delta+e20;
    f2    = f22*delta*delta+f21*delta+f20;

    double a,b;
    if (pax > 0) a = Ea/Dsq;
    else         a = -Ea/Dsq;
    if (pax > 0) b = -Dsq/(4*Ea)+mnsq*Ea/Dsq;
    else         b = Dsq/(4*Ea)-mnsq*Ea/Dsq;

    double A4,A3,A2,A1,A0;

    A4 = a*a*a2;
    A3 = 2*a*b2/Ea;
    A2 = (2*a*a2*b+c2+2*a*d2)/(Easq);
    A1 = (2*b*b2+2*e2)/(Easq*Ea);
    A0 = (a2*b*b+2*b*d2+f2)/(Easq*Easq);

        long double A3sq;
        A3sq = A3*A3;
        /*  
        long  double A0sq, A1sq, A2sq, A3sq, A4sq;
    A0sq = A0*A0;
    A1sq = A1*A1;
    A2sq = A2*A2;
    A3sq = A3*A3;
    A4sq = A4*A4;
        */

    long double B3, B2, B1, B0;
    B3 = 4*A4;
    B2 = 3*A3;
    B1 = 2*A2;
    B0 = A1;
    long double C2, C1, C0;
    C2 = -(A2/2 - 3*A3sq/(16*A4));
    C1 = -(3*A1/4. -A2*A3/(8*A4));
    C0 = -A0 + A1*A3/(16*A4);
    long double  D1, D0;
    D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
    D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;

    long double E0;
    E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;

    long  double t1,t2,t3,t4,t5;

//find the coefficients for the leading term in the Sturm sequence  
    t1 = A4;
    t2 = A4;
    t3 = C2;
    t4 = D1;
    t5 = E0;

    int nsol;
    nsol = signchange_n(t1,t2,t3,t4,t5)-signchange_p(t1,t2,t3,t4,t5);
    if( nsol < 0 ) nsol=0;

    return nsol;
}

void Davismt2::mt2_bisect(){

    solved = true;
    cout.precision(11);

//if masses are very small, use code for massless case.  
    if(masq < MIN_MASS && mbsq < MIN_MASS){ 
        mt2_massless();
        return;
    }


    double Deltasq0;     
    Deltasq0 = ma*(ma + 2*mn); //The minimum mass square to have two ellipses 

// find the coefficients for the two quadratic equations when Deltasq=Deltasq0.

    a1 = 1-pax*pax/(Easq);
    b1 = -pax*pay/(Easq);
    c1 = 1-pay*pay/(Easq);
    d1 = -pax*(Deltasq0-masq)/(2*Easq);
    e1 = -pay*(Deltasq0-masq)/(2*Easq);
    a2 = 1-pbx*pbx/(Ebsq);
    b2 = -pbx*pby/(Ebsq);
    c2 = 1-pby*pby/(Ebsq);
    d2 = -pmissx+pbx*(Deltasq0-mbsq)/(2*Ebsq)+pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
    e2 = -pmissy+pby*(Deltasq0-mbsq)/(2*Ebsq)+pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
    f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq0-mbsq)/(2*Eb)+
        (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq0-mbsq)/(2*Eb)+
        (pbx*pmissx+pby*pmissy)/Eb)+mnsq;

// find the center of the smaller ellipse 
    double x0,y0;
    x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
    y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);


// Does the larger ellipse contain the smaller one? 
    double dis=a2*x0*x0+2*b2*x0*y0+c2*y0*y0+2*d2*x0+2*e2*y0+f2;

    if(dis<=0.01)
    { 
        mt2_b  = (double) sqrt(mnsq+Deltasq0);
        return;
    }


/* find the coefficients for the two quadratic equations           */
/* coefficients for quadratic terms do not change                  */
/* coefficients for linear and constant terms are polynomials of   */
/*       delta=(Deltasq-m7sq)/(2 E7sq)                             */  
    d11 = -pax;
    e11 = -pay;
    f10 = mnsq;
    f12 = -Easq;
    d21 = (Easq*pbx)/Ebsq;
    d20 = ((masq - mbsq)*pbx)/(2.*Ebsq) - pmissx +
        (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
    e21 = (Easq*pby)/Ebsq;
    e20 = ((masq - mbsq)*pby)/(2.*Ebsq) - pmissy +
        (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
    f22 = -Easq*Easq/Ebsq;
    f21 = (-2*Easq*((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb))/Eb;
    f20 = mnsq + pmissx*pmissx + pmissy*pmissy - 
        ((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb)
        *((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb);

//Estimate upper bound of mT2
//when Deltasq > Deltasq_high1, the larger encloses the center of the smaller 
    double p2x0,p2y0;
    double Deltasq_high1;
    p2x0 = pmissx-x0;
    p2y0 = pmissy-y0;
    Deltasq_high1 = 2*Eb*sqrt(p2x0*p2x0+p2y0*p2y0+mnsq)-2*pbx*p2x0-2*pby*p2y0+mbsq;

//Another estimate, if both ellipses enclose the origin, Deltasq > mT2

    double Deltasq_high2, Deltasq_high21, Deltasq_high22;
    Deltasq_high21 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy+mbsq;
    Deltasq_high22 = 2*Ea*mn+masq;

    if ( Deltasq_high21 < Deltasq_high22 ) Deltasq_high2 = Deltasq_high22;
    else Deltasq_high2 = Deltasq_high21;

//pick the smaller upper bound   
    double Deltasq_high;
    if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high1;
    else Deltasq_high = Deltasq_high2;


    double Deltasq_low; //lower bound
    Deltasq_low = Deltasq0;

//number of solutions at Deltasq_low should not be larger than zero
    if( nsols(Deltasq_low) > 0 )
    {
//cout << "Davismt2::mt2_bisect() ==> nsolutions(Deltasq_low) > 0"<<endl;
        mt2_b = (double) sqrt(mnsq+Deltasq0);
        return;
    }

    int nsols_high, nsols_low;

    nsols_low  = nsols(Deltasq_low);
    int foundhigh;


//if nsols_high=nsols_low, we missed the region where the two ellipse overlap 
//if nsols_high=4, also need a scan because we may find the wrong tangent point.

    nsols_high = nsols(Deltasq_high);

    if(nsols_high == nsols_low || nsols_high == 4)
    {
    //foundhigh = scan_high(Deltasq_high);
        foundhigh = find_high(Deltasq_high);
        if(foundhigh == 0) 
        {
            if(verbose > 0) cout << "Davismt2::mt2_bisect() ==> Deltasq_high not found at event " << nevt << endl;
            mt2_b = sqrt( Deltasq_low + mnsq );
            return;
        }

    }

    while(sqrt(Deltasq_high+mnsq) - sqrt(Deltasq_low+mnsq) > precision)
    {
        double Deltasq_mid,nsols_mid;
    //bisect
        Deltasq_mid = (Deltasq_high+Deltasq_low)/2.;
        nsols_mid = nsols(Deltasq_mid);
    // if nsols_mid = 4, rescan for Deltasq_high
        if ( nsols_mid == 4 ) 
        {
            Deltasq_high = Deltasq_mid;
        //scan_high(Deltasq_high);
            find_high(Deltasq_high);
            continue;
        } 


        if(nsols_mid != nsols_low) Deltasq_high = Deltasq_mid;
        if(nsols_mid == nsols_low) Deltasq_low  = Deltasq_mid;
    }
    mt2_b = (double) sqrt( mnsq + Deltasq_high);
    return;
}

int Davismt2::find_high(double & Deltasq_high){
    double x0,y0;
    x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
    y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
    double Deltasq_low = (mn + ma)*(mn + ma) - mnsq;
    do 
    {
        double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.;
        int nsols_mid = nsols(Deltasq_mid);
        if ( nsols_mid == 2 )
        {
            Deltasq_high = Deltasq_mid;
            return 1;
        }
        else if (nsols_mid == 4)
        {
            Deltasq_high = Deltasq_mid;
            continue;
        }
        else if (nsols_mid ==0)
        {
            d1 = -pax*(Deltasq_mid-masq)/(2*Easq);
            e1 = -pay*(Deltasq_mid-masq)/(2*Easq);
            d2 = -pmissx + pbx*(Deltasq_mid - mbsq)/(2*Ebsq)
                + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
            e2 = -pmissy + pby*(Deltasq_mid - mbsq)/(2*Ebsq)
                + pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
            f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq)/(2*Eb)+
                (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq)/(2*Eb)+
                (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
// Does the larger ellipse contain the smaller one? 
            double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
            if (dis < 0) Deltasq_high = Deltasq_mid;
            else Deltasq_low = Deltasq_mid;
        }

    } while ( Deltasq_high - Deltasq_low > 0.001);
    return 0;
}

int Davismt2::scan_high(double & Deltasq_high){
    int foundhigh = 0 ;
    int nsols_high;


        //        double Deltasq_low;
    double tempmass, maxmass;
    tempmass = mn + ma;
    maxmass  = sqrt(mnsq + Deltasq_high);
    if (nevt == 32334) cout << "Davismt2::scan_high() ==> Deltasq_high = " << Deltasq_high << endl; // ???
    for(double mass = tempmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
    {
        Deltasq_high = mass*mass - mnsq;
        nsols_high   = nsols(Deltasq_high);

        if( nsols_high > 0)
        {
                  // Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
            foundhigh   = 1;
            break;
        }
    }
    return foundhigh;
}

int Davismt2::nsols(double Dsq){
    double delta = (Dsq-masq)/(2*Easq);

//calculate coefficients for the two quadratic equations
    d1 = d11*delta;
    e1 = e11*delta;
    f1 = f12*delta*delta+f10;
    d2 = d21*delta+d20;
    e2 = e21*delta+e20;
    f2 = f22*delta*delta+f21*delta+f20;

//obtain the coefficients for the 4th order equation 
//devided by Ea^n to make the variable dimensionless
    long double A4, A3, A2, A1, A0;

    A4 = 
        -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 + 
        4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 + 
        a1*a1*c2*c2;  

    A3 =
        (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 + 
        8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 + 
        4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 - 
        4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea;


    A2 =
        (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 - 
        8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 + 
        4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 - 
        8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 + 
        4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 + 
        4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq;

    A1 =
        (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 - 
        4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 - 
        4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 - 
        4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea);

    A0 =
        (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 + 
        4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 + 
        a1*a1*f2*f2)/(Easq*Easq);

        long double A3sq;
        /*
    long  double A0sq, A1sq, A2sq, A3sq, A4sq;
    A0sq = A0*A0;
    A1sq = A1*A1;
    A2sq = A2*A2;
    A4sq = A4*A4;
        */
    A3sq = A3*A3;

    long double B3, B2, B1, B0;
    B3 = 4*A4;
    B2 = 3*A3;
    B1 = 2*A2;
    B0 = A1;

    long double C2, C1, C0;
    C2 = -(A2/2 - 3*A3sq/(16*A4));
    C1 = -(3*A1/4. -A2*A3/(8*A4));
    C0 = -A0 + A1*A3/(16*A4);

    long double D1, D0;
    D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
    D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;

    long double E0;
    E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;

    long  double t1,t2,t3,t4,t5;
//find the coefficients for the leading term in the Sturm sequence  
    t1 = A4;
    t2 = A4;
    t3 = C2;
    t4 = D1;
    t5 = E0;


//The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
    int nsol;
    nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);

//Cannot have negative number of solutions, must be roundoff effect
    if (nsol < 0) nsol = 0;

    return nsol;

}  

//inline
int Davismt2::signchange_n(long double t1, long double t2, long double t3, long double t4, long double t5){
    int nsc;
    nsc=0;
    if(t1*t2>0) nsc++;
    if(t2*t3>0) nsc++;
    if(t3*t4>0) nsc++;
    if(t4*t5>0) nsc++;
    return nsc;
}

//inline
int Davismt2::signchange_p(long double t1, long double t2, long double t3, long double t4, long double t5){
    int nsc;
    nsc=0;
    if(t1*t2<0) nsc++;
    if(t2*t3<0) nsc++;
    if(t3*t4<0) nsc++;
    if(t4*t5<0) nsc++;
    return nsc;
}

}