CMS 3D CMS Logo

mt2w_bisect.cc
Go to the documentation of this file.
1 /***********************************************************************/
2 /* */
3 /* Finding MT2W */
4 /* Reference: arXiv:1203.4813 [hep-ph] */
5 /* Authors: Yang Bai, Hsin-Chia Cheng, */
6 /* Jason Gallicchio, Jiayin Gu */
7 /* Based on MT2 by: Hsin-Chia Cheng, Zhenyu Han */
8 /* May 8, 2012, v1.00a */
9 /* */
10 /***********************************************************************/
11 
12 
13 /*******************************************************************************
14  Usage:
15 
16  1. Define an object of type "mt2w":
17 
18  mt2w_bisect::mt2w mt2w_event;
19 
20  2. Set momenta:
21 
22  mt2w_event.set_momenta(pl,pb1,pb2,pmiss);
23 
24  where array pl[0..3], pb1[0..3], pb2[0..3] contains (E,px,py,pz), pmiss[0..2] contains (0,px,py)
25  for the visible particles and the missing momentum. pmiss[0] is not used.
26  All quantities are given in double.
27 
28  (Or use non-pointer method to do the same.)
29 
30  3. Use mt2w::get_mt2w() to obtain the value of mt2w:
31 
32  double mt2w_value = mt2w_event.get_mt2w();
33 
34 *******************************************************************************/
35 
36 #include <iostream>
37 #include <math.h>
39 
40 using namespace std;
41 
42 namespace heppy {
43 
44 namespace mt2w_bisect
45 {
46 
47 /*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*/
48 
49  const float mt2w::RELATIVE_PRECISION = 0.00001;
50  const float mt2w::ABSOLUTE_PRECISION = 0.0;
51  const float mt2w::MIN_MASS = 0.1;
52  const float mt2w::ZERO_MASS = 0.0;
53  const float mt2w::SCANSTEP = 0.1;
54 
55  mt2w::mt2w(double upper_bound, double error_value, double scan_step)
56  {
57  solved = false;
58  momenta_set = false;
59  mt2w_b = 0.; // The result field. Start it off at zero.
60  this->upper_bound = upper_bound; // the upper bound of search for MT2W, default value is 500 GeV
61  this->error_value = error_value; // if we couldn't find any compatible region below the upper_bound, output mt2w = error_value;
62  this->scan_step = scan_step; // if we need to scan to find the compatible region, this is the step of the scan
63  }
64 
65  double mt2w::get_mt2w()
66  {
67  if (!momenta_set)
68  {
69  cout <<" Please set momenta first!" << endl;
70  return error_value;
71  }
72 
73  if (!solved) mt2w_bisect();
74  return mt2w_b;
75  }
76 
77 
78  void mt2w::set_momenta(double *pl, double *pb1, double *pb2, double* pmiss)
79  {
80  // Pass in pointers to 4-vectors {E, px, py, px} of doubles.
81  // and pmiss must have [1] and [2] components for x and y. The [0] component is ignored.
82  set_momenta(pl[0], pl[1], pl[2], pl[3],
83  pb1[0], pb1[1], pb1[2], pb1[3],
84  pb2[0], pb2[1], pb2[2], pb2[3],
85  pmiss[1], pmiss[2]);
86  }
87 
88 
89 
90  void mt2w::set_momenta(double El, double plx, double ply, double plz,
91  double Eb1, double pb1x, double pb1y, double pb1z,
92  double Eb2, double pb2x, double pb2y, double pb2z,
93  double pmissx, double pmissy)
94  {
95  solved = false; //reset solved tag when momenta are changed.
96  momenta_set = true;
97 
98  double msqtemp; //used for saving the mass squared temporarily
99 
100  //l is the visible lepton
101 
102  this->El = El;
103  this->plx = plx;
104  this->ply = ply;
105  this->plz = plz;
106 
107  Elsq = El*El;
108 
109  msqtemp = El*El-plx*plx-ply*ply-plz*plz;
110  if (msqtemp > 0.0) {mlsq = msqtemp;}
111  else {mlsq = 0.0;} //mass squared can not be negative
112  ml = sqrt(mlsq); // all the input masses are calculated from sqrt(p^2)
113 
114  //b1 is the bottom on the same side as the visible lepton
115 
116  this->Eb1 = Eb1;
117  this->pb1x = pb1x;
118  this->pb1y = pb1y;
119  this->pb1z = pb1z;
120 
121  Eb1sq = Eb1*Eb1;
122 
123  msqtemp = Eb1*Eb1-pb1x*pb1x-pb1y*pb1y-pb1z*pb1z;
124  if (msqtemp > 0.0) {mb1sq = msqtemp;}
125  else {mb1sq = 0.0;} //mass squared can not be negative
126  mb1 = sqrt(mb1sq); // all the input masses are calculated from sqrt(p^2)
127 
128  //b2 is the other bottom (paired with the invisible W)
129 
130  this->Eb2 = Eb2;
131  this->pb2x = pb2x;
132  this->pb2y = pb2y;
133  this->pb2z = pb2z;
134 
135  Eb2sq = Eb2*Eb2;
136 
137  msqtemp = Eb2*Eb2-pb2x*pb2x-pb2y*pb2y-pb2z*pb2z;
138  if (msqtemp > 0.0) {mb2sq = msqtemp;}
139  else {mb2sq = 0.0;} //mass squared can not be negative
140  mb2 = sqrt(mb2sq); // all the input masses are calculated from sqrt(p^2)
141 
142 
143  //missing pt
144 
145 
146  this->pmissx = pmissx;
147  this->pmissy = pmissy;
148 
149  //set the values of masses
150 
151  mv = 0.0; //mass of neutrino
152  mw = 80.4; //mass of W-boson
153 
154 
155  //precision?
156 
157  if (ABSOLUTE_PRECISION > 100.*RELATIVE_PRECISION) precision = ABSOLUTE_PRECISION;
158  else precision = 100.*RELATIVE_PRECISION;
159  }
160 
161 
162  void mt2w::mt2w_bisect()
163  {
164 
165 
166  solved = true;
167  cout.precision(11);
168 
169  // In normal running, mtop_high WILL be compatible, and mtop_low will NOT.
170  double mtop_high = upper_bound; //set the upper bound of the search region
171  double mtop_low; //the lower bound of the search region is best chosen as m_W + m_b
172 
173  if (mb1 >= mb2) {mtop_low = mw + mb1;}
174  else {mtop_low = mw + mb2;}
175 
176  // The following if and while deal with the case where there might be a compatable region
177  // between mtop_low and 500 GeV, but it doesn't extend all the way up to 500.
178  //
179 
180  // If our starting high guess is not compatible, start the high guess from the low guess...
181  if (teco(mtop_high)==0) {mtop_high = mtop_low;}
182 
183  // .. and scan up until a compatible high bound is found.
184  //We can also raise the lower bound since we scaned over a region that is not compatible
185  while (teco(mtop_high)==0 && mtop_high < upper_bound + 2.*scan_step) {
186 
187  mtop_low=mtop_high;
188  mtop_high = mtop_high + scan_step;
189  }
190 
191  // if we can not find a compatible region under the upper bound, output the error value
192  if (mtop_high > upper_bound) {
193  mt2w_b = error_value;
194  return;
195  }
196 
197  // Once we have an compatible mtop_high, we can find mt2w using bisection method
198  while(mtop_high - mtop_low > precision)
199  {
200  double mtop_mid,teco_mid;
201  //bisect
202  mtop_mid = (mtop_high+mtop_low)/2.;
203  teco_mid = teco(mtop_mid);
204 
205  if(teco_mid == 0) {mtop_low = mtop_mid;}
206  else {mtop_high = mtop_mid;}
207 
208  }
209  mt2w_b = mtop_high; //output the value of mt2w
210  return;
211  }
212 
213 
214  // for a given event, teco ( mtop ) gives 1 if trial top mass mtop is compatible, 0 if mtop is not.
215 
216  int mt2w::teco( double mtop)
217  {
218 
219  //first test if mtop is larger than mb+mw
220 
221  if (mtop < mb1+mw || mtop < mb2+mw) {return 0;}
222 
223  //define delta for convenience, note the definition is different from the one in mathematica code by 2*E^2_{b2}
224 
225  double ETb2sq = Eb2sq - pb2z*pb2z; //transverse energy of b2
226  double delta = (mtop*mtop-mw*mw-mb2sq)/(2.*ETb2sq);
227 
228 
229  //del1 and del2 are \Delta'_1 and \Delta'_2 in the notes eq. 10,11
230 
231  double del1 = mw*mw - mv*mv - mlsq;
232  double del2 = mtop*mtop - mw*mw - mb1sq - 2*(El*Eb1-plx*pb1x-ply*pb1y-plz*pb1z);
233 
234  // aa bb cc are A B C in the notes eq.15
235 
236  double aa = (El*pb1x-Eb1*plx)/(Eb1*plz-El*pb1z);
237  double bb = (El*pb1y-Eb1*ply)/(Eb1*plz-El*pb1z);
238  double cc = (El*del2-Eb1*del1)/(2.*Eb1*plz-2.*El*pb1z);
239 
240 
241  //calculate coefficients for the two quadratic equations (ellipses), which are
242  //
243  // 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)
244  //
245  // 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)
246  //
247  // where x and y are px and py of the neutrino on the visible lepton chain
248 
249  a1 = Eb1sq*(1.+aa*aa)-(pb1x+pb1z*aa)*(pb1x+pb1z*aa);
250  b1 = Eb1sq*aa*bb - (pb1x+pb1z*aa)*(pb1y+pb1z*bb);
251  c1 = Eb1sq*(1.+bb*bb)-(pb1y+pb1z*bb)*(pb1y+pb1z*bb);
252  d1 = Eb1sq*aa*cc - (pb1x+pb1z*aa)*(pb1z*cc+del2/2.0);
253  e1 = Eb1sq*bb*cc - (pb1y+pb1z*bb)*(pb1z*cc+del2/2.0);
254  f1 = Eb1sq*(mv*mv+cc*cc) - (pb1z*cc+del2/2.0)*(pb1z*cc+del2/2.0);
255 
256  // 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)
257 
258  double det1 = (a1*(c1*f1 - e1*e1) - b1*(b1*f1 - d1*e1) + d1*(b1*e1-c1*d1))/(a1+c1);
259 
260  if (det1 > 0.0) {return 0;}
261 
262  //coefficients of the ellptical region
263 
264  a2 = 1-pb2x*pb2x/(ETb2sq);
265  b2 = -pb2x*pb2y/(ETb2sq);
266  c2 = 1-pb2y*pb2y/(ETb2sq);
267 
268  // d2o e2o f2o are coefficients in the p2x p2y plane (p2 is the momentum of the missing W-boson)
269  // it is convenient to calculate them first and transfer the ellipse to the p1x p1y plane
270  d2o = -delta*pb2x;
271  e2o = -delta*pb2y;
272  f2o = mw*mw - delta*delta*ETb2sq;
273 
274  d2 = -d2o -a2*pmissx -b2*pmissy;
275  e2 = -e2o -c2*pmissy -b2*pmissx;
276  f2 = a2*pmissx*pmissx + 2*b2*pmissx*pmissy + c2*pmissy*pmissy + 2*d2o*pmissx + 2*e2o*pmissy + f2o;
277 
278  //find a point in ellipse 1 and see if it's within the ellipse 2, define h0 for convenience
279  double x0, h0, y0, r0;
280  x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
281  h0 = (b1*x0 + e1)*(b1*x0 + e1) - c1*(a1*x0*x0 + 2*d1*x0 + f1);
282  if (h0 < 0.0) {return 0;} // if h0 < 0, y0 is not real and ellipse 1 is not real, this is a redundant check.
283  y0 = (-b1*x0 -e1 + sqrt(h0))/c1;
284  r0 = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
285  if (r0 < 0.0) {return 1;} // if the point is within the 2nd ellipse, mtop is compatible
286 
287 
288  //obtain the coefficients for the 4th order equation
289  //devided by Eb1^n to make the variable dimensionless
290  long double A4, A3, A2, A1, A0;
291 
292  A4 =
293  -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 +
294  4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 +
295  a1*a1*c2*c2;
296 
297  A3 =
298  (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 +
299  8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 +
300  4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 -
301  4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Eb1;
302 
303 
304  A2 =
305  (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 -
306  8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 +
307  4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 -
308  8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 +
309  4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 +
310  4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Eb1sq;
311 
312  A1 =
313  (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 -
314  4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 -
315  4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 -
316  4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Eb1sq*Eb1);
317 
318  A0 =
319  (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 +
320  4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 +
321  a1*a1*f2*f2)/(Eb1sq*Eb1sq);
322 
323  long double A3sq;
324  /*
325  long double A0sq, A1sq, A2sq, A3sq, A4sq;
326  A0sq = A0*A0;
327  A1sq = A1*A1;
328  A2sq = A2*A2;
329  A4sq = A4*A4;
330  */
331  A3sq = A3*A3;
332 
333  long double B3, B2, B1, B0;
334  B3 = 4*A4;
335  B2 = 3*A3;
336  B1 = 2*A2;
337  B0 = A1;
338 
339  long double C2, C1, C0;
340  C2 = -(A2/2 - 3*A3sq/(16*A4));
341  C1 = -(3*A1/4. -A2*A3/(8*A4));
342  C0 = -A0 + A1*A3/(16*A4);
343 
344  long double D1, D0;
345  D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
346  D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
347 
348  long double E0;
349  E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
350 
351  long double t1,t2,t3,t4,t5;
352  //find the coefficients for the leading term in the Sturm sequence
353  t1 = A4;
354  t2 = A4;
355  t3 = C2;
356  t4 = D1;
357  t5 = E0;
358 
359 
360  //The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
361  int nsol;
362  nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);
363 
364  //Cannot have negative number of solutions, must be roundoff effect
365  if (nsol < 0) nsol = 0;
366 
367  int out;
368  if (nsol == 0) {out = 0;} //output 0 if there is no solution, 1 if there is solution
369  else {out = 1;}
370 
371  return out;
372 
373  }
374 
375  inline int mt2w::signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5)
376  {
377  int nsc;
378  nsc=0;
379  if(t1*t2>0) nsc++;
380  if(t2*t3>0) nsc++;
381  if(t3*t4>0) nsc++;
382  if(t4*t5>0) nsc++;
383  return nsc;
384  }
385  inline int mt2w::signchange_p( long double t1, long double t2, long double t3, long double t4, long double t5)
386  {
387  int nsc;
388  nsc=0;
389  if(t1*t2<0) nsc++;
390  if(t2*t3<0) nsc++;
391  if(t3*t4<0) nsc++;
392  if(t4*t5<0) nsc++;
393  return nsc;
394  }
395 
396 }//end namespace mt2w_bisect
397 
398 }
dbl * delta
Definition: mlp_gen.cc:36
Divides< arg, void > D0
Definition: Factorize.h:143
auto const T2 &decltype(t1.eta()) t2
Definition: deltaR.h:16
T sqrt(T t)
Definition: SSEVec.h:18
Divides< A, C > D1
Definition: Factorize.h:144
TAKEN FROM http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/ElectroWeakAnalysis/Utilities/src/PdfWeig...
Definition: AlphaT.h:17
Float e1
Definition: deltaR.h:20
Float e2
Definition: deltaR.h:21