Heppy/src/mt2w_bisect.cc

0001 /***********************************************************************/
0002 /*                                                                     */
0003 /*              Finding MT2W                                           */
0004 /*              Reference:  arXiv:1203.4813 [hep-ph]                   */
0005 /*              Authors: Yang Bai, Hsin-Chia Cheng,                    */
0006 /*                       Jason Gallicchio, Jiayin Gu                   */
0007 /*              Based on MT2 by: Hsin-Chia Cheng, Zhenyu Han           */
0008 /*              May 8, 2012, v1.00a                                    */
0009 /*                                                                     */
0010 /***********************************************************************/
0011
0012 /*******************************************************************************
0013   Usage:
0014
0015   1. Define an object of type "mt2w":
0016
0017      mt2w_bisect::mt2w mt2w_event;
0018
0019   2. Set momenta:
0020
0021      mt2w_event.set_momenta(pl,pb1,pb2,pmiss);
0022
0023      where array pl[0..3], pb1[0..3], pb2[0..3] contains (E,px,py,pz), pmiss[0..2] contains (0,px,py)
0024      for the visible particles and the missing momentum. pmiss[0] is not used.
0025      All quantities are given in double.
0026
0027      (Or use non-pointer method to do the same.)
0028
0029   3. Use mt2w::get_mt2w() to obtain the value of mt2w:
0030
0031      double mt2w_value = mt2w_event.get_mt2w();
0032
0033 *******************************************************************************/
0034
0035 #include <iostream>
0036 #include <cmath>
0037 #include "PhysicsTools/Heppy/interface/mt2w_bisect.h"
0038
0039 using namespace std;
0040
0041 namespace heppy {
0042
0043   namespace mt2w_bisect {
0044
0045     /*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*/
0046
0047     const float mt2w::RELATIVE_PRECISION = 0.00001;
0048     const float mt2w::ABSOLUTE_PRECISION = 0.0;
0049     const float mt2w::MIN_MASS = 0.1;
0050     const float mt2w::ZERO_MASS = 0.0;
0051     const float mt2w::SCANSTEP = 0.1;
0052
0053     mt2w::mt2w(double upper_bound, double error_value, double scan_step) {
0054       solved = false;
0055       momenta_set = false;
0056       mt2w_b = 0.;                      // The result field.  Start it off at zero.
0057       this->upper_bound = upper_bound;  // the upper bound of search for MT2W, default value is 500 GeV
0058       this->error_value =
0059           error_value;  // if we couldn't find any compatible region below the upper_bound, output mt2w = error_value;
0060       this->scan_step = scan_step;  // if we need to scan to find the compatible region, this is the step of the scan
0061     }
0062
0063     double mt2w::get_mt2w() {
0064       if (!momenta_set) {
0065         cout << " Please set momenta first!" << endl;
0066         return error_value;
0067       }
0068
0069       if (!solved)
0070         mt2w_bisect();
0071       return mt2w_b;
0072     }
0073
0074     void mt2w::set_momenta(double *pl, double *pb1, double *pb2, double *pmiss) {
0075       // Pass in pointers to 4-vectors {E, px, py, px} of doubles.
0076       // and pmiss must have [1] and [2] components for x and y.  The [0] component is ignored.
0077       set_momenta(
0078           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]);
0079     }
0080
0081     void mt2w::set_momenta(double El,
0082                            double plx,
0083                            double ply,
0084                            double plz,
0085                            double Eb1,
0086                            double pb1x,
0087                            double pb1y,
0088                            double pb1z,
0089                            double Eb2,
0090                            double pb2x,
0091                            double pb2y,
0092                            double pb2z,
0093                            double pmissx,
0094                            double pmissy) {
0095       solved = false;  //reset solved tag when momenta are changed.
0096       momenta_set = true;
0097
0098       double msqtemp;  //used for saving the mass squared temporarily
0099
0100       //l is the visible lepton
0101
0102       this->El = El;
0103       this->plx = plx;
0104       this->ply = ply;
0105       this->plz = plz;
0106
0107       Elsq = El * El;
0108
0109       msqtemp = El * El - plx * plx - ply * ply - plz * plz;
0110       if (msqtemp > 0.0) {
0111         mlsq = msqtemp;
0112       } else {
0113         mlsq = 0.0;
0114       }                 //mass squared can not be negative
0115       ml = sqrt(mlsq);  // all the input masses are calculated from sqrt(p^2)
0116
0117       //b1 is the bottom on the same side as the visible lepton
0118
0119       this->Eb1 = Eb1;
0120       this->pb1x = pb1x;
0121       this->pb1y = pb1y;
0122       this->pb1z = pb1z;
0123
0124       Eb1sq = Eb1 * Eb1;
0125
0126       msqtemp = Eb1 * Eb1 - pb1x * pb1x - pb1y * pb1y - pb1z * pb1z;
0127       if (msqtemp > 0.0) {
0128         mb1sq = msqtemp;
0129       } else {
0130         mb1sq = 0.0;
0131       }                   //mass squared can not be negative
0132       mb1 = sqrt(mb1sq);  // all the input masses are calculated from sqrt(p^2)
0133
0134       //b2 is the other bottom (paired with the invisible W)
0135
0136       this->Eb2 = Eb2;
0137       this->pb2x = pb2x;
0138       this->pb2y = pb2y;
0139       this->pb2z = pb2z;
0140
0141       Eb2sq = Eb2 * Eb2;
0142
0143       msqtemp = Eb2 * Eb2 - pb2x * pb2x - pb2y * pb2y - pb2z * pb2z;
0144       if (msqtemp > 0.0) {
0145         mb2sq = msqtemp;
0146       } else {
0147         mb2sq = 0.0;
0148       }                   //mass squared can not be negative
0149       mb2 = sqrt(mb2sq);  // all the input masses are calculated from sqrt(p^2)
0150
0151       //missing pt
0152
0153       this->pmissx = pmissx;
0154       this->pmissy = pmissy;
0155
0156       //set the values of masses
0157
0158       mv = 0.0;   //mass of neutrino
0159       mw = 80.4;  //mass of W-boson
0160
0161       //precision?
0162
0163       if (ABSOLUTE_PRECISION > 100. * RELATIVE_PRECISION)
0164         precision = ABSOLUTE_PRECISION;
0165       else
0166         precision = 100. * RELATIVE_PRECISION;
0167     }
0168
0169     void mt2w::mt2w_bisect() {
0170       solved = true;
0171       cout.precision(11);
0172
0173       // In normal running, mtop_high WILL be compatible, and mtop_low will NOT.
0174       double mtop_high = upper_bound;  //set the upper bound of the search region
0175       double mtop_low;                 //the lower bound of the search region is best chosen as m_W + m_b
0176
0177       if (mb1 >= mb2) {
0178         mtop_low = mw + mb1;
0179       } else {
0180         mtop_low = mw + mb2;
0181       }
0182
0183       // The following if and while deal with the case where there might be a compatable region
0184       // between mtop_low and 500 GeV, but it doesn't extend all the way up to 500.
0185       //
0186
0187       // If our starting high guess is not compatible, start the high guess from the low guess...
0188       if (teco(mtop_high) == 0) {
0189         mtop_high = mtop_low;
0190       }
0191
0192       // .. and scan up until a compatible high bound is found.
0193       //We can also raise the lower bound since we scaned over a region that is not compatible
0194       while (teco(mtop_high) == 0 && mtop_high < upper_bound + 2. * scan_step) {
0195         mtop_low = mtop_high;
0196         mtop_high = mtop_high + scan_step;
0197       }
0198
0199       // if we can not find a compatible region under the upper bound, output the error value
0200       if (mtop_high > upper_bound) {
0201         mt2w_b = error_value;
0202         return;
0203       }
0204
0205       // Once we have an compatible mtop_high, we can find mt2w using bisection method
0206       while (mtop_high - mtop_low > precision) {
0207         double mtop_mid, teco_mid;
0208         //bisect
0209         mtop_mid = (mtop_high + mtop_low) / 2.;
0210         teco_mid = teco(mtop_mid);
0211
0212         if (teco_mid == 0) {
0213           mtop_low = mtop_mid;
0214         } else {
0215           mtop_high = mtop_mid;
0216         }
0217       }
0218       mt2w_b = mtop_high;  //output the value of mt2w
0219       return;
0220     }
0221
0222     // for a given event, teco ( mtop ) gives 1 if trial top mass mtop is compatible, 0 if mtop is not.
0223
0224     int mt2w::teco(double mtop) {
0225       //first test if mtop is larger than mb+mw
0226
0227       if (mtop < mb1 + mw || mtop < mb2 + mw) {
0228         return 0;
0229       }
0230
0231       //define delta for convenience, note the definition is different from the one in mathematica code by 2*E^2_{b2}
0232
0233       double ETb2sq = Eb2sq - pb2z * pb2z;  //transverse energy of b2
0234       double delta = (mtop * mtop - mw * mw - mb2sq) / (2. * ETb2sq);
0235
0236       //del1 and del2 are \Delta'_1 and \Delta'_2 in the notes eq. 10,11
0237
0238       double del1 = mw * mw - mv * mv - mlsq;
0239       double del2 = mtop * mtop - mw * mw - mb1sq - 2 * (El * Eb1 - plx * pb1x - ply * pb1y - plz * pb1z);
0240
0241       // aa bb cc are A B C in the notes eq.15
0242
0243       double aa = (El * pb1x - Eb1 * plx) / (Eb1 * plz - El * pb1z);
0244       double bb = (El * pb1y - Eb1 * ply) / (Eb1 * plz - El * pb1z);
0245       double cc = (El * del2 - Eb1 * del1) / (2. * Eb1 * plz - 2. * El * pb1z);
0246
0247       //calculate coefficients for the two quadratic equations (ellipses), which are
0248       //
0249       //  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)
0250       //
0251       //  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)
0252       //
0253       //  where x and y are px and py of the neutrino on the visible lepton chain
0254
0255       a1 = Eb1sq * (1. + aa * aa) - (pb1x + pb1z * aa) * (pb1x + pb1z * aa);
0256       b1 = Eb1sq * aa * bb - (pb1x + pb1z * aa) * (pb1y + pb1z * bb);
0257       c1 = Eb1sq * (1. + bb * bb) - (pb1y + pb1z * bb) * (pb1y + pb1z * bb);
0258       d1 = Eb1sq * aa * cc - (pb1x + pb1z * aa) * (pb1z * cc + del2 / 2.0);
0259       e1 = Eb1sq * bb * cc - (pb1y + pb1z * bb) * (pb1z * cc + del2 / 2.0);
0260       f1 = Eb1sq * (mv * mv + cc * cc) - (pb1z * cc + del2 / 2.0) * (pb1z * cc + del2 / 2.0);
0261
0262       //  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)
0263
0264       double det1 = (a1 * (c1 * f1 - e1 * e1) - b1 * (b1 * f1 - d1 * e1) + d1 * (b1 * e1 - c1 * d1)) / (a1 + c1);
0265
0266       if (det1 > 0.0) {
0267         return 0;
0268       }
0269
0270       //coefficients of the ellptical region
0271
0272       a2 = 1 - pb2x * pb2x / (ETb2sq);
0273       b2 = -pb2x * pb2y / (ETb2sq);
0274       c2 = 1 - pb2y * pb2y / (ETb2sq);
0275
0276       // d2o e2o f2o are coefficients in the p2x p2y plane (p2 is the momentum of the missing W-boson)
0277       // it is convenient to calculate them first and transfer the ellipse to the p1x p1y plane
0278       d2o = -delta * pb2x;
0279       e2o = -delta * pb2y;
0280       f2o = mw * mw - delta * delta * ETb2sq;
0281
0282       d2 = -d2o - a2 * pmissx - b2 * pmissy;
0283       e2 = -e2o - c2 * pmissy - b2 * pmissx;
0284       f2 = a2 * pmissx * pmissx + 2 * b2 * pmissx * pmissy + c2 * pmissy * pmissy + 2 * d2o * pmissx +
0285            2 * e2o * pmissy + f2o;
0286
0287       //find a point in ellipse 1 and see if it's within the ellipse 2, define h0 for convenience
0288       double x0, h0, y0, r0;
0289       x0 = (c1 * d1 - b1 * e1) / (b1 * b1 - a1 * c1);
0290       h0 = (b1 * x0 + e1) * (b1 * x0 + e1) - c1 * (a1 * x0 * x0 + 2 * d1 * x0 + f1);
0291       if (h0 < 0.0) {
0292         return 0;
0293       }  // if h0 < 0, y0 is not real and ellipse 1 is not real, this is a redundant check.
0294       y0 = (-b1 * x0 - e1 + sqrt(h0)) / c1;
0295       r0 = a2 * x0 * x0 + 2 * b2 * x0 * y0 + c2 * y0 * y0 + 2 * d2 * x0 + 2 * e2 * y0 + f2;
0296       if (r0 < 0.0) {
0297         return 1;
0298       }  // if the point is within the 2nd ellipse, mtop is compatible
0299
0300       //obtain the coefficients for the 4th order equation
0301       //devided by Eb1^n to make the variable dimensionless
0302       long double A4, A3, A2, A1, A0;
0303
0304       A4 = -4 * a2 * b1 * b2 * c1 + 4 * a1 * b2 * b2 * c1 + a2 * a2 * c1 * c1 + 4 * a2 * b1 * b1 * c2 -
0305            4 * a1 * b1 * b2 * c2 - 2 * a1 * a2 * c1 * c2 + a1 * a1 * c2 * c2;
0306
0307       A3 = (-4 * a2 * b2 * c1 * d1 + 8 * a2 * b1 * c2 * d1 - 4 * a1 * b2 * c2 * d1 - 4 * a2 * b1 * c1 * d2 +
0308             8 * a1 * b2 * c1 * d2 - 4 * a1 * b1 * c2 * d2 - 8 * a2 * b1 * b2 * e1 + 8 * a1 * b2 * b2 * e1 +
0309             4 * a2 * a2 * c1 * e1 - 4 * a1 * a2 * c2 * e1 + 8 * a2 * b1 * b1 * e2 - 8 * a1 * b1 * b2 * e2 -
0310             4 * a1 * a2 * c1 * e2 + 4 * a1 * a1 * c2 * e2) /
0311            Eb1;
0312
0313       A2 = (4 * a2 * c2 * d1 * d1 - 4 * a2 * c1 * d1 * d2 - 4 * a1 * c2 * d1 * d2 + 4 * a1 * c1 * d2 * d2 -
0314             8 * a2 * b2 * d1 * e1 - 8 * a2 * b1 * d2 * e1 + 16 * a1 * b2 * d2 * e1 + 4 * a2 * a2 * e1 * e1 +
0315             16 * a2 * b1 * d1 * e2 - 8 * a1 * b2 * d1 * e2 - 8 * a1 * b1 * d2 * e2 - 8 * a1 * a2 * e1 * e2 +
0316             4 * a1 * a1 * e2 * e2 - 4 * a2 * b1 * b2 * f1 + 4 * a1 * b2 * b2 * f1 + 2 * a2 * a2 * c1 * f1 -
0317             2 * a1 * a2 * c2 * f1 + 4 * a2 * b1 * b1 * f2 - 4 * a1 * b1 * b2 * f2 - 2 * a1 * a2 * c1 * f2 +
0318             2 * a1 * a1 * c2 * f2) /
0319            Eb1sq;
0320
0321       A1 = (-8 * a2 * d1 * d2 * e1 + 8 * a1 * d2 * d2 * e1 + 8 * a2 * d1 * d1 * e2 - 8 * a1 * d1 * d2 * e2 -
0322             4 * a2 * b2 * d1 * f1 - 4 * a2 * b1 * d2 * f1 + 8 * a1 * b2 * d2 * f1 + 4 * a2 * a2 * e1 * f1 -
0323             4 * a1 * a2 * e2 * f1 + 8 * a2 * b1 * d1 * f2 - 4 * a1 * b2 * d1 * f2 - 4 * a1 * b1 * d2 * f2 -
0324             4 * a1 * a2 * e1 * f2 + 4 * a1 * a1 * e2 * f2) /
0325            (Eb1sq * Eb1);
0326
0327       A0 = (-4 * a2 * d1 * d2 * f1 + 4 * a1 * d2 * d2 * f1 + a2 * a2 * f1 * f1 + 4 * a2 * d1 * d1 * f2 -
0328             4 * a1 * d1 * d2 * f2 - 2 * a1 * a2 * f1 * f2 + a1 * a1 * f2 * f2) /
0329            (Eb1sq * Eb1sq);
0330
0331       long double A3sq;
0332       /*
0333     long  double A0sq, A1sq, A2sq, A3sq, A4sq;
0334     A0sq = A0*A0;
0335     A1sq = A1*A1;
0336     A2sq = A2*A2;
0337     A4sq = A4*A4;
0338     */
0339       A3sq = A3 * A3;
0340
0341       long double B3, B2, B1, B0;
0342       B3 = 4 * A4;
0343       B2 = 3 * A3;
0344       B1 = 2 * A2;
0345       B0 = A1;
0346
0347       long double C2, C1, C0;
0348       C2 = -(A2 / 2 - 3 * A3sq / (16 * A4));
0349       C1 = -(3 * A1 / 4. - A2 * A3 / (8 * A4));
0350       C0 = -A0 + A1 * A3 / (16 * A4);
0351
0352       long double D1, D0;
0353       D1 = -B1 - (B3 * C1 * C1 / C2 - B3 * C0 - B2 * C1) / C2;
0354       D0 = -B0 - B3 * C0 * C1 / (C2 * C2) + B2 * C0 / C2;
0355
0356       long double E0;
0357       E0 = -C0 - C2 * D0 * D0 / (D1 * D1) + C1 * D0 / D1;
0358
0359       long double t1, t2, t3, t4, t5;
0360       //find the coefficients for the leading term in the Sturm sequence
0361       t1 = A4;
0362       t2 = A4;
0363       t3 = C2;
0364       t4 = D1;
0365       t5 = E0;
0366
0367       //The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
0368       int nsol;
0369       nsol = signchange_n(t1, t2, t3, t4, t5) - signchange_p(t1, t2, t3, t4, t5);
0370
0371       //Cannot have negative number of solutions, must be roundoff effect
0372       if (nsol < 0)
0373         nsol = 0;
0374
0375       int out;
0376       if (nsol == 0) {
0377         out = 0;
0378       }  //output 0 if there is no solution, 1 if there is solution
0379       else {
0380         out = 1;
0381       }
0382
0383       return out;
0384     }
0385
0386     inline int mt2w::signchange_n(long double t1, long double t2, long double t3, long double t4, long double t5) {
0387       int nsc;
0388       nsc = 0;
0389       if (t1 * t2 > 0)
0390         nsc++;
0391       if (t2 * t3 > 0)
0392         nsc++;
0393       if (t3 * t4 > 0)
0394         nsc++;
0395       if (t4 * t5 > 0)
0396         nsc++;
0397       return nsc;
0398     }
0399     inline int mt2w::signchange_p(long double t1, long double t2, long double t3, long double t4, long double t5) {
0400       int nsc;
0401       nsc = 0;
0402       if (t1 * t2 < 0)
0403         nsc++;
0404       if (t2 * t3 < 0)
0405         nsc++;
0406       if (t3 * t4 < 0)
0407         nsc++;
0408       if (t4 * t5 < 0)
0409         nsc++;
0410       return nsc;
0411     }
0412
0413   }  //end namespace mt2w_bisect
0414
0415 }  // namespace heppy