mt2w_bisect.cc

Go to the documentation of this file.

/***********************************************************************/
/*                                                                     */
/*              Finding MT2W                                           */
/*              Reference:  arXiv:1203.4813 [hep-ph]                   */
/*              Authors: Yang Bai, Hsin-Chia Cheng,                    */
/*                       Jason Gallicchio, Jiayin Gu                   */
/*              Based on MT2 by: Hsin-Chia Cheng, Zhenyu Han           */ 
/*              May 8, 2012, v1.00a                                    */
/*                                                                     */  
/***********************************************************************/


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

  1. Define an object of type "mt2w":
     
     mt2w_bisect::mt2w mt2w_event;
 
  2. Set momenta:
 
     mt2w_event.set_momenta(pl,pb1,pb2,pmiss);
 
     where array pl[0..3], pb1[0..3], pb2[0..3] contains (E,px,py,pz), pmiss[0..2] contains (0,px,py) 
     for the visible particles and the missing momentum. pmiss[0] is not used. 
     All quantities are given in double.    
     
     (Or use non-pointer method to do the same.)

  3. Use mt2w::get_mt2w() to obtain the value of mt2w:

     double mt2w_value = mt2w_event.get_mt2w();       
          
*******************************************************************************/ 
              
#include <iostream>
#include <cmath>
#include "PhysicsTools/Heppy/interface/mt2w_bisect.h"

using namespace std;

namespace heppy {
  
namespace mt2w_bisect
{

/*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 mt2w::RELATIVE_PRECISION = 0.00001;
  const float mt2w::ABSOLUTE_PRECISION = 0.0;
  const float mt2w::MIN_MASS = 0.1; 
  const float mt2w::ZERO_MASS = 0.0; 
  const float mt2w::SCANSTEP = 0.1;

  mt2w::mt2w(double upper_bound, double error_value, double scan_step)
  {
    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
  }

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


  void mt2w::set_momenta(double *pl, double *pb1, double *pb2, double* pmiss)
  {
    // 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 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)
  {
    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;
  }


  void mt2w::mt2w_bisect()
  {
  
   
    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;
  }


  // for a given event, teco ( mtop ) gives 1 if trial top mass mtop is compatible, 0 if mtop is not.

  int mt2w::teco(  double mtop)
  {

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

  inline int mt2w::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 mt2w::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;
  }

}//end namespace mt2w_bisect

}