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

Davismt2()

heppy::Davismt2::Davismt2

(

)

Definition at line 50 of file Davismt2.cc.

References L1EGammaClusterEmuProducer_cfi::scale, and verbose.

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

~Davismt2()

heppy::Davismt2::~Davismt2 ( )

virtual

Definition at line 58 of file Davismt2.cc.

{}

Member Function Documentation

find_high()

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

private

Definition at line 512 of file Davismt2.cc.

References testProducerWithPsetDescEmpty_cfi::a2, b1, b2, alignmentValidation::c1, d1, StorageManager_cfg::e1, and DeadROC_duringRun::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;
  }

get_mt2()

double heppy::Davismt2::get_mt2

(

)

Definition at line 60 of file Davismt2.cc.

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

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

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

mt2_bisect()

void heppy::Davismt2::mt2_bisect

(

)

Definition at line 371 of file Davismt2.cc.

References testProducerWithPsetDescEmpty_cfi::a2, b1, b2, alignmentValidation::c1, gather_cfg::cout, d1, StorageManager_cfg::e1, DeadROC_duringRun::f2, nevt, boostedTaus_cff::precision, mathSSE::sqrt(), and verbose.

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

mt2_massless()

void heppy::Davismt2::mt2_massless

(

)

Definition at line 177 of file Davismt2.cc.

References testProducerWithPsetDescEmpty_cfi::a2, b2, HltBtagPostValidation_cff::c, funct::cos(), gather_cfg::cout, EgHLTOffHistBins_cfi::mass, BPHMonitor_cfi::maxmass, BPHMonitor_cfi::minmass, nevt, boostedTaus_cff::precision, alignCSCRings::s, funct::sin(), mathSSE::sqrt(), theta(), and verbose.

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

nsols()

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

private

Definition at line 570 of file Davismt2.cc.

References testProducerWithPsetDescEmpty_cfi::a2, b1, b2, alignmentValidation::c1, d1, dumpMFGeometry_cfg::delta, StorageManager_cfg::e1, DeadROC_duringRun::f1, DeadROC_duringRun::f2, RandomServiceHelper::t1, RandomServiceHelper::t2, and RandomServiceHelper::t3.

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

nsols_massless()

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

private

Definition at line 299 of file Davismt2.cc.

References a, testProducerWithPsetDescEmpty_cfi::a2, b, b2, d1, dumpMFGeometry_cfg::delta, StorageManager_cfg::e1, DeadROC_duringRun::f1, DeadROC_duringRun::f2, RandomServiceHelper::t1, RandomServiceHelper::t2, and RandomServiceHelper::t3.

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

print()

void heppy::Davismt2::print

(

void

)

Definition at line 167 of file Davismt2.cc.

References gather_cfg::cout, and L1EGammaClusterEmuProducer_cfi::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;
  }

scan_high()

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

private

Definition at line 547 of file Davismt2.cc.

References gather_cfg::cout, EgHLTOffHistBins_cfi::mass, BPHMonitor_cfi::maxmass, nevt, 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;
  }

set_mn()

void heppy::Davismt2::set_mn

(

double

mn

)

Definition at line 160 of file Davismt2.cc.

References L1EGammaClusterEmuProducer_cfi::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;
  }

set_momenta()

void heppy::Davismt2::set_momenta	(	double *	pa0,
		double *	pb0,
		double *	pmiss0
	)

Definition at line 71 of file Davismt2.cc.

References boostedTaus_cff::precision, L1EGammaClusterEmuProducer_cfi::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;
  }

set_verbose()

void heppy::Davismt2::set_verbose ( int  vlevel )

inline

Definition at line 27 of file Davismt2.h.

References verbose.

{ verbose = vlevel; };

signchange_n()

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

private

Definition at line 660 of file Davismt2.cc.

References RandomServiceHelper::t1, RandomServiceHelper::t2, and RandomServiceHelper::t3.

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

signchange_p()

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

private

Definition at line 675 of file Davismt2.cc.

References RandomServiceHelper::t1, RandomServiceHelper::t2, and RandomServiceHelper::t3.

                                                                                                           {
    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

a1

double heppy::Davismt2::a1

private

Definition at line 59 of file Davismt2.h.

a2

double heppy::Davismt2::a2

private

Definition at line 59 of file Davismt2.h.

ABSOLUTE_PRECISION

const float heppy::Davismt2::ABSOLUTE_PRECISION = 0.0

static

Definition at line 16 of file Davismt2.h.

b1

double heppy::Davismt2::b1

private

Definition at line 59 of file Davismt2.h.

b2

double heppy::Davismt2::b2

private

Definition at line 59 of file Davismt2.h.

c1

double heppy::Davismt2::c1

private

Definition at line 59 of file Davismt2.h.

c2

double heppy::Davismt2::c2

private

Definition at line 59 of file Davismt2.h.

d1

double heppy::Davismt2::d1

private

Definition at line 59 of file Davismt2.h.

d11

double heppy::Davismt2::d11

private

Definition at line 60 of file Davismt2.h.

d2

double heppy::Davismt2::d2

private

Definition at line 59 of file Davismt2.h.

d20

double heppy::Davismt2::d20

private

Definition at line 60 of file Davismt2.h.

d21

double heppy::Davismt2::d21

private

Definition at line 60 of file Davismt2.h.

e1

double heppy::Davismt2::e1

private

Definition at line 59 of file Davismt2.h.

e11

double heppy::Davismt2::e11

private

Definition at line 60 of file Davismt2.h.

e2

double heppy::Davismt2::e2

private

Definition at line 59 of file Davismt2.h.

e20

double heppy::Davismt2::e20

private

Definition at line 60 of file Davismt2.h.

e21

double heppy::Davismt2::e21

private

Definition at line 60 of file Davismt2.h.

Ea

double heppy::Davismt2::Ea

private

Definition at line 47 of file Davismt2.h.

Easq

double heppy::Davismt2::Easq

private

Definition at line 53 of file Davismt2.h.

Eb

double heppy::Davismt2::Eb

private

Definition at line 49 of file Davismt2.h.

Ebsq

double heppy::Davismt2::Ebsq

private

Definition at line 54 of file Davismt2.h.

f1

double heppy::Davismt2::f1

private

Definition at line 59 of file Davismt2.h.

f10

double heppy::Davismt2::f10

private

Definition at line 60 of file Davismt2.h.

f12

double heppy::Davismt2::f12

private

Definition at line 60 of file Davismt2.h.

f2

double heppy::Davismt2::f2

private

Definition at line 59 of file Davismt2.h.

f20

double heppy::Davismt2::f20

private

Definition at line 60 of file Davismt2.h.

f21

double heppy::Davismt2::f21

private

Definition at line 60 of file Davismt2.h.

f22

double heppy::Davismt2::f22

private

Definition at line 60 of file Davismt2.h.

ma

double heppy::Davismt2::ma

private

Definition at line 47 of file Davismt2.h.

masq

double heppy::Davismt2::masq

private

Definition at line 53 of file Davismt2.h.

mb

double heppy::Davismt2::mb

private

Definition at line 49 of file Davismt2.h.

mbsq

double heppy::Davismt2::mbsq

private

Definition at line 54 of file Davismt2.h.

MIN_MASS

const float heppy::Davismt2::MIN_MASS = 0.1

static

Definition at line 17 of file Davismt2.h.

mn

double heppy::Davismt2::mn

private

Definition at line 50 of file Davismt2.h.

mn_unscale

double heppy::Davismt2::mn_unscale

private

Definition at line 50 of file Davismt2.h.

mnsq

double heppy::Davismt2::mnsq

private

Definition at line 56 of file Davismt2.h.

momenta_set

bool heppy::Davismt2::momenta_set

private

Definition at line 35 of file Davismt2.h.

mt2_b

double heppy::Davismt2::mt2_b

private

Definition at line 36 of file Davismt2.h.

nevt

int heppy::Davismt2::nevt

Definition at line 30 of file Davismt2.h.

pax

double heppy::Davismt2::pax

private

Definition at line 47 of file Davismt2.h.

pay

double heppy::Davismt2::pay

private

Definition at line 47 of file Davismt2.h.

pbx

double heppy::Davismt2::pbx

private

Definition at line 49 of file Davismt2.h.

pby

double heppy::Davismt2::pby

private

Definition at line 49 of file Davismt2.h.

pmissx

double heppy::Davismt2::pmissx

private

Definition at line 48 of file Davismt2.h.

pmissxsq

double heppy::Davismt2::pmissxsq

private

Definition at line 55 of file Davismt2.h.

pmissy

double heppy::Davismt2::pmissy

private

Definition at line 48 of file Davismt2.h.

pmissysq

double heppy::Davismt2::pmissysq

private

Definition at line 55 of file Davismt2.h.

precision

double heppy::Davismt2::precision

private

Definition at line 63 of file Davismt2.h.

Referenced by output.OutputBranch::fill().

RELATIVE_PRECISION

const float heppy::Davismt2::RELATIVE_PRECISION = 0.00001

static

Definition at line 15 of file Davismt2.h.

scale

double heppy::Davismt2::scale

private

Definition at line 62 of file Davismt2.h.

Referenced by python.rootplot.rootmath.Target::__repr__().

SCANSTEP

const float heppy::Davismt2::SCANSTEP = 0.1

static

Definition at line 19 of file Davismt2.h.

solved

bool heppy::Davismt2::solved

private

Definition at line 34 of file Davismt2.h.

verbose

int heppy::Davismt2::verbose

private

Definition at line 33 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(), objects.VertexAnalyzer.VertexAnalyzer::process(), and set_verbose().

ZERO_MASS

const float heppy::Davismt2::ZERO_MASS = 0.0

static

Definition at line 18 of file Davismt2.h.