heppy::mt2w_bisect::mt2w Class Reference

#include <mt2w_bisect.h>

Public Member Functions
double	get_mt2w ()
	mt2w (double upper_bound=500.0, double error_value=499.0, double scan_step=0.5)
void	set_momenta (double *pl0, double *pb10, double *pb20, double *pmiss0)
void	set_momenta (double El, double plx, double ply, double plz, double Eb1, double pb1x, double pb1y, double pb1z, double Eb2, double pb2x, double pb2y, double pb2z, double pmissx, double pmissy)

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

Protected Member Functions
void	mt2w_bisect ()

Private Member Functions
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)
int	teco (double mtop)

Private Attributes
double	a1
double	a2
double	b1
double	b2
double	c1
double	c2
double	d1
double	d2
double	d2o
double	e1
double	e2
double	e2o
double	Eb1
double	Eb1sq
double	Eb2
double	Eb2sq
double	El
double	Elsq
double	error_value
double	f1
double	f2
double	f2o
double	mb1
double	mb1sq
double	mb2
double	mb2sq
double	ml
double	mlsq
bool	momenta_set
double	mt2w_b
double	mv
double	mw
double	pb1x
double	pb1y
double	pb1z
double	pb2x
double	pb2y
double	pb2z
double	plx
double	ply
double	plz
double	pmissx
double	pmissy
double	precision
double	scan_step
bool	solved
double	upper_bound

Detailed Description

Definition at line 18 of file mt2w_bisect.h.

Constructor & Destructor Documentation

heppy::mt2w_bisect::mt2w::mt2w	(	double	upper_bound = 500.0,
		double	error_value = 499.0,
		double	scan_step = 0.5
	)

Definition at line 55 of file mt2w_bisect.cc.

  {
    solved = false;
    momenta_set = false;
    mt2w_b  = 0.;  // The result field.  Start it off at zero.
    this->upper_bound = upper_bound;  // the upper bound of search for MT2W, default value is 500 GeV 
    this->error_value = error_value;  // if we couldn't find any compatible region below the upper_bound, output mt2w = error_value;
    this->scan_step = scan_step;    // if we need to scan to find the compatible region, this is the step of the scan
  }

Member Function Documentation

double heppy::mt2w_bisect::mt2w::get_mt2w

(

)

Definition at line 65 of file mt2w_bisect.cc.

References gather_cfg::cout.

  {
    if (!momenta_set)
      {
    cout <<" Please set momenta first!" << endl;
    return error_value;
      }
        
    if (!solved) mt2w_bisect();
    return mt2w_b;
  }

void heppy::mt2w_bisect::mt2w::mt2w_bisect ( )

protected

Definition at line 162 of file mt2w_bisect.cc.

References gather_cfg::cout, MuonStationSelectors_cff::mb1, and MuonStationSelectors_cff::mb2.

  {
  
   
    solved = true;
    cout.precision(11);

    // In normal running, mtop_high WILL be compatible, and mtop_low will NOT.
    double mtop_high = upper_bound; //set the upper bound of the search region
    double mtop_low;                //the lower bound of the search region is best chosen as m_W + m_b

    if (mb1 >= mb2) {mtop_low = mw + mb1;}
    else {mtop_low = mw + mb2;}

    // The following if and while deal with the case where there might be a compatable region
    // between mtop_low and 500 GeV, but it doesn't extend all the way up to 500.
    // 

    // If our starting high guess is not compatible, start the high guess from the low guess...
    if (teco(mtop_high)==0) {mtop_high = mtop_low;}

    // .. and scan up until a compatible high bound is found.
    //We can also raise the lower bound since we scaned over a region that is not compatible
    while (teco(mtop_high)==0 && mtop_high < upper_bound + 2.*scan_step) {

      mtop_low=mtop_high;
      mtop_high = mtop_high + scan_step;
    }

    // if we can not find a compatible region under the upper bound, output the error value
    if (mtop_high > upper_bound) {
      mt2w_b = error_value;
      return;
    }

    // Once we have an compatible mtop_high, we can find mt2w using bisection method
    while(mtop_high - mtop_low > precision)
      {
    double mtop_mid,teco_mid;
    //bisect
    mtop_mid = (mtop_high+mtop_low)/2.;
    teco_mid = teco(mtop_mid);
      
    if(teco_mid == 0) {mtop_low  = mtop_mid;}
    else {mtop_high  = mtop_mid;}

      }
    mt2w_b = mtop_high;   //output the value of mt2w
    return;
  }

void heppy::mt2w_bisect::mt2w::set_momenta	(	double *	pl0,
		double *	pb10,
		double *	pb20,
		double *	pmiss0
	)

Definition at line 78 of file mt2w_bisect.cc.

  {
    // Pass in pointers to 4-vectors {E, px, py, px} of doubles.  
    // and pmiss must have [1] and [2] components for x and y.  The [0] component is ignored.
    set_momenta(pl[0],  pl[1],  pl[2],  pl[3],
        pb1[0], pb1[1], pb1[2], pb1[3],
        pb2[0], pb2[1], pb2[2], pb2[3],
        pmiss[1], pmiss[2]);
  }

void heppy::mt2w_bisect::mt2w::set_momenta	(	double	El,
		double	plx,
		double	ply,
		double	plz,
		double	Eb1,
		double	pb1x,
		double	pb1y,
		double	pb1z,
		double	Eb2,
		double	pb2x,
		double	pb2y,
		double	pb2z,
		double	pmissx,
		double	pmissy
	)

Definition at line 90 of file mt2w_bisect.cc.

References MuonStationSelectors_cff::mb1, MuonStationSelectors_cff::mb2, and mathSSE::sqrt().

  {
    solved = false;     //reset solved tag when momenta are changed.
    momenta_set = true;

    double msqtemp;   //used for saving the mass squared temporarily

    //l is the visible lepton

    this->El  = El;
    this->plx = plx;
    this->ply = ply;
    this->plz = plz;

    Elsq = El*El;

    msqtemp = El*El-plx*plx-ply*ply-plz*plz;
    if (msqtemp > 0.0) {mlsq = msqtemp;}
    else {mlsq = 0.0;}                           //mass squared can not be negative
    ml = sqrt(mlsq);                             // all the input masses are calculated from sqrt(p^2)

    //b1 is the bottom on the same side as the visible lepton

    this->Eb1  = Eb1;
    this->pb1x = pb1x;
    this->pb1y = pb1y;
    this->pb1z = pb1z;

    Eb1sq = Eb1*Eb1;

    msqtemp = Eb1*Eb1-pb1x*pb1x-pb1y*pb1y-pb1z*pb1z;
    if (msqtemp > 0.0) {mb1sq = msqtemp;}
    else {mb1sq = 0.0;}                          //mass squared can not be negative
    mb1 = sqrt(mb1sq);                           // all the input masses are calculated from sqrt(p^2)

    //b2 is the other bottom (paired with the invisible W)

    this->Eb2  = Eb2;
    this->pb2x = pb2x;
    this->pb2y = pb2y;
    this->pb2z = pb2z;

    Eb2sq = Eb2*Eb2;

    msqtemp = Eb2*Eb2-pb2x*pb2x-pb2y*pb2y-pb2z*pb2z;
    if (msqtemp > 0.0) {mb2sq = msqtemp;}
    else {mb2sq = 0.0;}                          //mass squared can not be negative
    mb2 = sqrt(mb2sq);                           // all the input masses are calculated from sqrt(p^2)


    //missing pt


    this->pmissx = pmissx; 
    this->pmissy = pmissy;

    //set the values of masses

    mv = 0.0;   //mass of neutrino
    mw = 80.4;  //mass of W-boson


    //precision?

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

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

inlineprivate

Definition at line 375 of file mt2w_bisect.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::mt2w_bisect::mt2w::signchange_p ( long double  t1, long double  t2, long double  t3, long double  t4, long double  t5  )

inlineprivate

Definition at line 385 of file mt2w_bisect.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::mt2w_bisect::mt2w::teco ( double  mtop )

private

Definition at line 216 of file mt2w_bisect.cc.

References alignmentValidation::c1, delta, python.connectstrParser::f1, python.connectstrParser::f2, MuonStationSelectors_cff::mb1, MuonStationSelectors_cff::mb2, dbtoconf::out, and mathSSE::sqrt().

  {

    //first test if mtop is larger than mb+mw

    if (mtop < mb1+mw || mtop < mb2+mw) {return 0;}

    //define delta for convenience, note the definition is different from the one in mathematica code by 2*E^2_{b2}

    double ETb2sq = Eb2sq - pb2z*pb2z;  //transverse energy of b2
    double delta = (mtop*mtop-mw*mw-mb2sq)/(2.*ETb2sq);


    //del1 and del2 are \Delta'_1 and \Delta'_2 in the notes eq. 10,11

    double del1 = mw*mw - mv*mv - mlsq;
    double del2 = mtop*mtop - mw*mw - mb1sq - 2*(El*Eb1-plx*pb1x-ply*pb1y-plz*pb1z);

    // aa bb cc are A B C in the notes eq.15

    double aa = (El*pb1x-Eb1*plx)/(Eb1*plz-El*pb1z);
    double bb = (El*pb1y-Eb1*ply)/(Eb1*plz-El*pb1z);
    double cc = (El*del2-Eb1*del1)/(2.*Eb1*plz-2.*El*pb1z);

  
    //calculate coefficients for the two quadratic equations (ellipses), which are
    //
    //  a1 x^2 + 2 b1 x y + c1 y^2 + 2 d1 x + 2 e1 y + f1 = 0 ,  from the 2 steps decay chain (with visible lepton)
    //
    //  a2 x^2 + 2 b2 x y + c2 y^2 + 2 d2 x + 2 e2 y + f2 <= 0 , from the 1 stop decay chain (with W missing)
    //
    //  where x and y are px and py of the neutrino on the visible lepton chain

    a1 = Eb1sq*(1.+aa*aa)-(pb1x+pb1z*aa)*(pb1x+pb1z*aa);
    b1 = Eb1sq*aa*bb - (pb1x+pb1z*aa)*(pb1y+pb1z*bb);
    c1 = Eb1sq*(1.+bb*bb)-(pb1y+pb1z*bb)*(pb1y+pb1z*bb);
    d1 = Eb1sq*aa*cc - (pb1x+pb1z*aa)*(pb1z*cc+del2/2.0);
    e1 = Eb1sq*bb*cc - (pb1y+pb1z*bb)*(pb1z*cc+del2/2.0);
    f1 = Eb1sq*(mv*mv+cc*cc) - (pb1z*cc+del2/2.0)*(pb1z*cc+del2/2.0);

    //  First check if ellipse 1 is real (don't need to do this for ellipse 2, ellipse 2 is always real for mtop > mw+mb)

    double det1 = (a1*(c1*f1 - e1*e1) - b1*(b1*f1 - d1*e1) + d1*(b1*e1-c1*d1))/(a1+c1);

    if (det1 > 0.0) {return 0;}

    //coefficients of the ellptical region

    a2 = 1-pb2x*pb2x/(ETb2sq);
    b2 = -pb2x*pb2y/(ETb2sq);
    c2 = 1-pb2y*pb2y/(ETb2sq);

    // d2o e2o f2o are coefficients in the p2x p2y plane (p2 is the momentum of the missing W-boson)
    // it is convenient to calculate them first and transfer the ellipse to the p1x p1y plane
    d2o = -delta*pb2x;
    e2o = -delta*pb2y;
    f2o = mw*mw - delta*delta*ETb2sq;

    d2 = -d2o -a2*pmissx -b2*pmissy;
    e2 = -e2o -c2*pmissy -b2*pmissx;
    f2 = a2*pmissx*pmissx + 2*b2*pmissx*pmissy + c2*pmissy*pmissy + 2*d2o*pmissx + 2*e2o*pmissy + f2o;

    //find a point in ellipse 1 and see if it's within the ellipse 2, define h0 for convenience
    double x0, h0, y0, r0;
    x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
    h0 = (b1*x0 + e1)*(b1*x0 + e1) - c1*(a1*x0*x0 + 2*d1*x0 + f1);
    if (h0 < 0.0) {return 0;}  // if h0 < 0, y0 is not real and ellipse 1 is not real, this is a redundant check.
    y0 = (-b1*x0 -e1 + sqrt(h0))/c1;
    r0 = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
    if (r0 < 0.0) {return 1;}  // if the point is within the 2nd ellipse, mtop is compatible


    //obtain the coefficients for the 4th order equation 
    //devided by Eb1^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)/Eb1;


    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)/Eb1sq;

    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)/(Eb1sq*Eb1);

    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)/(Eb1sq*Eb1sq);

    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;

    int out;
    if (nsol == 0) {out = 0;}  //output 0 if there is no solution, 1 if there is solution
    else {out = 1;}

    return out;
  
  }  

Member Data Documentation

double heppy::mt2w_bisect::mt2w::a1

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::a2

private

Definition at line 75 of file mt2w_bisect.h.

const float heppy::mt2w_bisect::mt2w::ABSOLUTE_PRECISION = 0.0

static

Definition at line 22 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::b1

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::b2

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::c1

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::c2

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::d1

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::d2

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::d2o

private

Definition at line 76 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::e1

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::e2

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::e2o

private

Definition at line 76 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::Eb1

private

Definition at line 64 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::Eb1sq

private

Definition at line 71 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::Eb2

private

Definition at line 65 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::Eb2sq

private

Definition at line 72 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::El

private

Definition at line 63 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::Elsq

private

Definition at line 70 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::error_value

private

Definition at line 54 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::f1

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::f2

private

Definition at line 75 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::f2o

private

Definition at line 76 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::mb1

private

Definition at line 64 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::mb1sq

private

Definition at line 71 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::mb2

private

Definition at line 65 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::mb2sq

private

Definition at line 72 of file mt2w_bisect.h.

const float heppy::mt2w_bisect::mt2w::MIN_MASS = 0.1

static

Definition at line 23 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::ml

private

Definition at line 63 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::mlsq

private

Definition at line 70 of file mt2w_bisect.h.

bool heppy::mt2w_bisect::mt2w::momenta_set

private

Definition at line 52 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::mt2w_b

private

Definition at line 56 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::mv

private

Definition at line 67 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::mw

private

Definition at line 67 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::pb1x

private

Definition at line 64 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::pb1y

private

Definition at line 64 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::pb1z

private

Definition at line 64 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::pb2x

private

Definition at line 65 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::pb2y

private

Definition at line 65 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::pb2z

private

Definition at line 65 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::plx

private

Definition at line 63 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::ply

private

Definition at line 63 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::plz

private

Definition at line 63 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::pmissx

private

Definition at line 66 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::pmissy

private

Definition at line 66 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::precision

private

Definition at line 78 of file mt2w_bisect.h.

const float heppy::mt2w_bisect::mt2w::RELATIVE_PRECISION = 0.00001

static

Definition at line 21 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::scan_step

private

Definition at line 55 of file mt2w_bisect.h.

const float heppy::mt2w_bisect::mt2w::SCANSTEP = 0.1

static

Definition at line 25 of file mt2w_bisect.h.

bool heppy::mt2w_bisect::mt2w::solved

private

Definition at line 51 of file mt2w_bisect.h.

double heppy::mt2w_bisect::mt2w::upper_bound

private

Definition at line 53 of file mt2w_bisect.h.

const float heppy::mt2w_bisect::mt2w::ZERO_MASS = 0.0

static

Definition at line 24 of file mt2w_bisect.h.