file src/mt2w_bisect.cpp
[No description available]
Namespaces
| Name |
|---|
| mt2w_bisect |
Source code
/***********************************************************************/
/* */
/* 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 <math.h>
#include "gambit/ColliderBit/mt2w_bisect.h"
using namespace std;
namespace mt2w_bisect
{
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 B3, B2, B1, B0;
B3 = 4*A4;
B2 = 3*A3;
B1 = 2*A2;
B0 = A1;
long double C2, C1, C0;
C2 = -(A2/2 - 3*A3*A3/(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
Updated on 2023-06-26 at 21:36:56 +0000