heppy::Davismt2 Class Reference

#include <Davismt2.h>

Public Member Functions
	Davismt2 ()
double	get_mt2 ()
void	mt2_bisect ()
void	mt2_massless ()
void	print ()
void	set_mn (double mn)
void	set_momenta (double *pa0, double *pb0, double *pmiss0)
void	set_verbose (int vlevel)
virtual	~Davismt2 ()

Public Attributes
int	nevt

Static Public Attributes
static const float	ABSOLUTE_PRECISION = 0.0
static const float	MIN_MASS = 0.1
static const float	RELATIVE_PRECISION = 0.00001
static const float	SCANSTEP = 0.1
static const float	ZERO_MASS = 0.0

Private Member Functions
int	find_high (double &Deltasq_high)
int	nsols (double Dsq)
int	nsols_massless (double Dsq)
int	scan_high (double &Deltasq_high)
int	signchange_n (long double t1, long double t2, long double t3, long double t4, long double t5)
int	signchange_p (long double t1, long double t2, long double t3, long double t4, long double t5)

Private Attributes
double	a1
double	a2
double	b1
double	b2
double	c1
double	c2
double	d1
double	d11
double	d2
double	d20
double	d21
double	e1
double	e11
double	e2
double	e20
double	e21
double	Ea
double	Easq
double	Eb
double	Ebsq
double	f1
double	f10
double	f12
double	f2
double	f20
double	f21
double	f22
double	ma
double	masq
double	mb
double	mbsq
double	mn
double	mn_unscale
double	mnsq
bool	momenta_set
double	mt2_b
double	pax
double	pay
double	pbx
double	pby
double	pmissx
double	pmissxsq
double	pmissy
double	pmissysq
double	precision
double	scale
bool	solved
int	verbose

Detailed Description

Definition at line 12 of file Davismt2.h.

Constructor & Destructor Documentation

heppy::Davismt2::Davismt2

(

)

Definition at line 51 of file Davismt2.cc.

References pileupReCalc_HLTpaths::scale.

                  {
    solved = false;
    momenta_set = false;
    mt2_b  = 0.;
    scale = 1.;
    verbose = 1;
}

heppy::Davismt2::~Davismt2 ( )

virtual

Definition at line 59 of file Davismt2.cc.

{}

Member Function Documentation

int heppy::Davismt2::find_high ( double &  Deltasq_high )

private

Definition at line 495 of file Davismt2.cc.

References alignmentValidation::c1, reco::e1, reco::e2, and validate-o2o-wbm::f2.

                                            {
    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;
}

double heppy::Davismt2::get_mt2

(

)

Definition at line 61 of file Davismt2.cc.

References gather_cfg::cout, and pileupReCalc_HLTpaths::scale.

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

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

void heppy::Davismt2::mt2_bisect

(

)

Definition at line 347 of file Davismt2.cc.

References alignmentValidation::c1, gather_cfg::cout, reco::e1, reco::e2, validate-o2o-wbm::f2, and mathSSE::sqrt().

                         {

    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;
}

void heppy::Davismt2::mt2_massless

(

)

Definition at line 154 of file Davismt2.cc.

References EnergyCorrector::c, funct::cos(), gather_cfg::cout, ResonanceBuilder::mass, alignCSCRings::s, funct::sin(), mathSSE::sqrt(), and theta().

                           {

//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 heppy::Davismt2::nsols ( double  Dsq )

private

Definition at line 560 of file Davismt2.cc.

References alignmentValidation::c1, delta, reco::e1, reco::e2, validate-o2o-wbm::f1, validate-o2o-wbm::f2, and reco::t2.

                             {
    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;

}  

int heppy::Davismt2::nsols_massless ( double  Dsq )

private

Definition at line 280 of file Davismt2.cc.

References a, b, delta, reco::e1, reco::e2, validate-o2o-wbm::f1, validate-o2o-wbm::f2, and reco::t2.

                                      {
    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 heppy::Davismt2::print

(

void

)

Definition at line 146 of file Davismt2.cc.

References gather_cfg::cout, and pileupReCalc_HLTpaths::scale.

                    {
    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;
}

int heppy::Davismt2::scan_high ( double &  Deltasq_high )

private

Definition at line 535 of file Davismt2.cc.

References gather_cfg::cout, ResonanceBuilder::mass, and mathSSE::sqrt().

                                            {
    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;
}

void heppy::Davismt2::set_mn

(

double

mn

)

Definition at line 139 of file Davismt2.cc.

References pileupReCalc_HLTpaths::scale.

                               {
    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 heppy::Davismt2::set_momenta	(	double *	pa0,
		double *	pb0,
		double *	pmiss0
	)

Definition at line 72 of file Davismt2.cc.

References pileupReCalc_HLTpaths::scale, mathSSE::sqrt(), and groupFilesInBlocks::temp.

                                                                  {
    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 heppy::Davismt2::set_verbose ( int  vlevel )

inline

Definition at line 27 of file Davismt2.h.

{verbose = vlevel;};

int heppy::Davismt2::signchange_n ( long double  t1, long double  t2, long double  t3, long double  t4, long double  t5  )

private

Definition at line 655 of file Davismt2.cc.

                                                                                                        {
    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;
}

int heppy::Davismt2::signchange_p ( long double  t1, long double  t2, long double  t3, long double  t4, long double  t5  )

private

Definition at line 666 of file Davismt2.cc.

                                                                                                        {
    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;
}

Member Data Documentation

double heppy::Davismt2::a1

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::a2

private

Definition at line 60 of file Davismt2.h.

const float heppy::Davismt2::ABSOLUTE_PRECISION = 0.0

static

Definition at line 16 of file Davismt2.h.

double heppy::Davismt2::b1

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::b2

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::c1

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::c2

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::d1

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::d11

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::d2

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::d20

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::d21

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::e1

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::e11

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::e2

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::e20

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::e21

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::Ea

private

Definition at line 48 of file Davismt2.h.

double heppy::Davismt2::Easq

private

Definition at line 54 of file Davismt2.h.

double heppy::Davismt2::Eb

private

Definition at line 50 of file Davismt2.h.

double heppy::Davismt2::Ebsq

private

Definition at line 55 of file Davismt2.h.

double heppy::Davismt2::f1

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::f10

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::f12

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::f2

private

Definition at line 60 of file Davismt2.h.

double heppy::Davismt2::f20

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::f21

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::f22

private

Definition at line 61 of file Davismt2.h.

double heppy::Davismt2::ma

private

Definition at line 48 of file Davismt2.h.

double heppy::Davismt2::masq

private

Definition at line 54 of file Davismt2.h.

double heppy::Davismt2::mb

private

Definition at line 50 of file Davismt2.h.

double heppy::Davismt2::mbsq

private

Definition at line 55 of file Davismt2.h.

const float heppy::Davismt2::MIN_MASS = 0.1

static

Definition at line 17 of file Davismt2.h.

double heppy::Davismt2::mn

private

Definition at line 51 of file Davismt2.h.

double heppy::Davismt2::mn_unscale

private

Definition at line 51 of file Davismt2.h.

double heppy::Davismt2::mnsq

private

Definition at line 57 of file Davismt2.h.

bool heppy::Davismt2::momenta_set

private

Definition at line 36 of file Davismt2.h.

double heppy::Davismt2::mt2_b

private

Definition at line 37 of file Davismt2.h.

int heppy::Davismt2::nevt

Definition at line 30 of file Davismt2.h.

double heppy::Davismt2::pax

private

Definition at line 48 of file Davismt2.h.

double heppy::Davismt2::pay

private

Definition at line 48 of file Davismt2.h.

double heppy::Davismt2::pbx

private

Definition at line 50 of file Davismt2.h.

double heppy::Davismt2::pby

private

Definition at line 50 of file Davismt2.h.

double heppy::Davismt2::pmissx

private

Definition at line 49 of file Davismt2.h.

double heppy::Davismt2::pmissxsq

private

Definition at line 56 of file Davismt2.h.

double heppy::Davismt2::pmissy

private

Definition at line 49 of file Davismt2.h.

double heppy::Davismt2::pmissysq

private

Definition at line 56 of file Davismt2.h.

double heppy::Davismt2::precision

private

Definition at line 64 of file Davismt2.h.

const float heppy::Davismt2::RELATIVE_PRECISION = 0.00001

static

Definition at line 15 of file Davismt2.h.

double heppy::Davismt2::scale

private

Definition at line 63 of file Davismt2.h.

Referenced by python.rootplot.rootmath.Target::__repr__(), Vispa.Gui.PortWidget.PortWidget::connectionPoint(), Vispa.Gui.VispaWidget.VispaWidget::drawHeaderBackground(), Vispa.Gui.ConnectableWidget.ConnectableWidget::drawPortNames(), Vispa.Gui.ConnectableWidget.ConnectableWidget::getEffectivePortHeight(), Vispa.Gui.VispaWidget.VispaWidget::imageRectF(), Vispa.Gui.VispaWidget.VispaWidget::move(), Vispa.Gui.VispaWidget.VispaWidget::setZoom(), and Vispa.Gui.ConnectableWidget.ConnectableWidget::sizeHint().

const float heppy::Davismt2::SCANSTEP = 0.1

static

Definition at line 19 of file Davismt2.h.

bool heppy::Davismt2::solved

private

Definition at line 35 of file Davismt2.h.

int heppy::Davismt2::verbose

private

Definition at line 34 of file Davismt2.h.

Referenced by python.diff_provenance.difference::list_diff(), python.diffProv.difference::list_diff(), python.diff_provenance.difference::module_diff(), python.diffProv.difference::module_diff(), python.diff_provenance.difference::onefilemodules(), python.diffProv.difference::onefilemodules(), core.TriggerMatchAnalyzer.TriggerMatchAnalyzer::process(), core.SkimAnalyzerCount.SkimAnalyzerCount::process(), and objects.VertexAnalyzer.VertexAnalyzer::process().

const float heppy::Davismt2::ZERO_MASS = 0.0

static

Definition at line 18 of file Davismt2.h.