checkmate is hosted by Hepforge, IPPP Durham
mt2bl_bisect.cc
Go to the documentation of this file.
1 /***********************************************************************/
2 /* */
3 /* Finding MT2bl */
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 /* Nov 9, 2012, v1.00a */
9 /* */
10 /***********************************************************************/
11 
12 
13 /*******************************************************************************
14  Usage:
15 
16  1. Define an object of type "mt2bl":
17 
18  mt2blw_bisect::mt2bl mt2bl_event;
19 
20  2. Set momenta:
21 
22  mt2bl_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 mt2bl::get_mt2bl() to obtain the value of mt2bl:
31 
32  double mt2bl_value = mt2bl_event.get_mt2bl();
33 
34  *******************************************************************************/
35 
36 
37 
38 //note: For a given trial mother mass my (for top in this case), Deltasq is defined as my^2-mn1^2,
39 // where mn1 is the mass of one of the daughter particle (the neutrino on the visible lepton
40 // chain or the invisible W-boson, see the code for details.) We try to find the minimum
41 // compatible Deltasq and transfer it to minimum my.
42 //
43 // There is also variables "delta" and "delta2" which are NOT the sqrt of Delasq.
44 //
45 //another note: This code can be easily modified to calculate mT2 with different daughter masses
46 // in general.
47 
48 
49 
50 #include <iostream>
51 #include <math.h>
52 #include "mt2bl_bisect.h"
53 
54 using namespace std;
55 
56 namespace mt2bl_bisect
57 {
58 mt2bl::mt2bl()
59 {
60  solved = false;
61  momenta_set = false;
62  mt2bl_b = 0.;
63  scan_step = 0.1;
64 }
65 
66 double mt2bl::get_mt2bl()
67 {
68  if (!momenta_set)
69  {
70  cout <<" Please set momenta first!" << endl;
71  return 0;
72  }
73 
74  if (!solved) mt2bl_bisect();
75  return mt2bl_b;
76 }
77 
78 void mt2bl::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 mt2bl::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  //combine pl and pb1 and call it pa, pb2 is renamed as pb.
99 
100  pax = plx + pb1x;
101  pay = ply + pb1y;
102  Easq = (El+Eb1)*(El+Eb1) - (plz+pb1z)*(plz+pb1z); //Ea is actually the ET of a
103  Ea = sqrt(Easq);
104  masq = max(0., Easq - pax*pax - pay*pay);
105  ma = sqrt(masq);
106 
107 
108  pbx = pb2x;
109  pby = pb2y;
110  Ebsq = Eb2*Eb2 - pb2z*pb2z; //Eb is actually the ET of b
111  Eb = sqrt(Ebsq);
112  mbsq = max(0., Ebsq - pbx*pbx - pby*pby);
113  mb = sqrt(mbsq);
114 
115  this->pmissx = pmissx;
116  this->pmissy = pmissy;
117 
118  pmissxsq = pmissx*pmissx;
119  pmissysq = pmissy*pmissy;
120 
121  mn1 = 0.0; //mn1 is the mass of neutrino by default
122  mn2 = 80.4; //mass of W-boson
123 
124 
125 // set ma + mn1 >= mb +mn2 (i.e. if ma + mn1 < mb + mn2, just switch the names of the two sides)
126  if(ma + mn1 < mb + mn2)
127  {
128  double temp;
129  temp = pax; pax = pbx; pbx = temp;
130  temp = pay; pay = pby; pby = temp;
131  temp = Ea; Ea = Eb; Eb = temp;
132  temp = Easq; Easq = Ebsq; Ebsq = temp;
133  temp = masq; masq = mbsq; mbsq = temp;
134  temp = ma; ma = mb; mb = temp;
135  temp = mn1; mn1 = mn2; mn2 = temp;
136  }
137 
138 // cout << endl << "ma+mn1 = " << ma + mn1 << endl;
139 
140  masq = ma*ma;
141  mbsq = mb*mb;
142 
143  mn1sq = mn1*mn1;
144  mn2sq = mn2*mn2;
145 
146 
148  else precision = 100.*RELATIVE_PRECISION;
149 }
150 
151 //void mt2bl::print()
152 //{
153 // cout << " pax = " << pax*scale << "; pay = " << pay*scale << "; ma = " << ma*scale <<";"<< endl;
154 // cout << " pbx = " << pbx*scale << "; pby = " << pby*scale << "; mb = " << mb*scale <<";"<< endl;
155 // cout << " pmissx = " << pmissx*scale << "; pmissy = " << pmissy*scale <<";"<< endl;
156 // cout << " mn = " << mn_unscale<<";" << endl;
157 //}
158 
159 
160 void mt2bl::mt2bl_bisect()
161 {
162 
163 
164  solved = true;
165  cout.precision(11);
166 
167  double Deltasq0 = ma*(ma + 2*mn1); //The minimum mass square to have two ellipses (in general Deltasq = my^2 - mn1^2)
168  double Delta2sq0 = Deltasq0 + mn1sq - mn2sq;
169 
170 // find the coefficients for the two quadratic equations when Deltasq=Deltasq0.
171 
172  a1 = 1-pax*pax/(Easq);
173  b1 = -pax*pay/(Easq);
174  c1 = 1-pay*pay/(Easq);
175  d1 = -pax*(Deltasq0-masq)/(2*Easq);
176  e1 = -pay*(Deltasq0-masq)/(2*Easq);
177  a2 = 1-pbx*pbx/(Ebsq);
178  b2 = -pbx*pby/(Ebsq);
179  c2 = 1-pby*pby/(Ebsq);
180  d2 = -pmissx+pbx*(Delta2sq0-mbsq)/(2*Ebsq)+pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
181  e2 = -pmissy+pby*(Delta2sq0-mbsq)/(2*Ebsq)+pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
182  f2 = pmissx*pmissx+pmissy*pmissy-((Delta2sq0-mbsq)/(2*Eb)+
183  (pbx*pmissx+pby*pmissy)/Eb)*((Delta2sq0-mbsq)/(2*Eb)+
184  (pbx*pmissx+pby*pmissy)/Eb)+mn2sq;
185 
186 
187 //cout << "pmissx = " << pmissx << " , pmissy = " << pmissy << endl;
188 //cout << "ma = " << ma << " , mb = " << mb <<endl;
189 //cout << "a2 = " << a2 << " , d2=" << d2 << " , e2=" << e2 << " , f2=" << f2 << endl;
190 
191 // find the center of the smaller ellipse
192  double x0,y0;
193  x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
194  y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
195 
196 
197 // Does the larger ellipse contain the smaller one?
198  double dis=a2*x0*x0+2*b2*x0*y0+c2*y0*y0+2*d2*x0+2*e2*y0+f2;
199 
200  if(dis<=0.01)
201  {
202  mt2bl_b = (double) sqrt(mn1sq+Deltasq0);
203  return;
204  }
205 
206 
207 // find the coefficients for the two quadratic equations
208 // coefficients for quadratic terms do not change
209 // Coefficients for linear and constant terms are are polynomials of delta=(my^2 - mn1^2 -mb^2)/(2 Eb^2) for ellipse 1
210 // and polynomials of delta2 = (my^2 - mn2^2 -mb^2)/(2 Eb^2) for ellipse 2.
211 
212  d11 = -pax;
213  e11 = -pay;
214  f10 = mn1sq;
215  f12 = -Easq;
216  d21 = pbx;
217  d20 = -a2*pmissx - b2*pmissy;
218  e21 = pby;
219  e20 = -c2*pmissy - b2*pmissx;
220  f22 = - Ebsq;
221  f21 = -2.*(pmissx*pbx+pmissy*pby);
222  f20 = a2*pmissxsq + 2.*b2*pmissx*pmissy + c2*pmissysq + mn2sq;
223 
224 //Estimate upper bound of mt2bl
225 //when Deltasq > Deltasq_high1, the larger encloses the center of the smaller
226  double p2x0,p2y0;
227  double Deltasq_high1;
228  p2x0 = pmissx-x0;
229  p2y0 = pmissy-y0;
230  Deltasq_high1 = 2*Eb*sqrt(p2x0*p2x0+p2y0*p2y0+mn2sq)-2*pbx*p2x0-2*pby*p2y0+mbsq + mn2sq - mn1sq; // need mn2sq - mn1sq for ellipse 2
231 
232 //Another estimate, if both ellipses enclose the origin, Deltasq > mt2bl
233 
234  double Deltasq_high2, Deltasq_high21, Deltasq_high22;
235  Deltasq_high21 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mn2sq)-2*pbx*pmissx-2*pby*pmissy+mbsq + mn2sq - mn1sq;
236  Deltasq_high22 = 2*Ea*mn1+masq;
237 
238  if ( Deltasq_high21 < Deltasq_high22 ) Deltasq_high2 = Deltasq_high22;
239  else Deltasq_high2 = Deltasq_high21;
240 
241 //pick the smaller upper bound
242  double Deltasq_high;
243  if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high1;
244  else Deltasq_high = Deltasq_high2;
245 
246  double Deltasq_low; //lower bound
247  Deltasq_low = Deltasq0;
248 
249 //number of solutions at Deltasq_low should not be larger than zero
250  if( nsols(Deltasq_low) > 0 )
251  {
252  //cout << "nsolutions(Deltasq_low) > 0"<<endl;
253  mt2bl_b = (double) sqrt(mn1sq+Deltasq0);
254  return;
255  }
256 
257  int nsols_high, nsols_low;
258 
259  nsols_low = nsols(Deltasq_low);
260  int foundhigh;
261 
262 
263 //if nsols_high=nsols_low, we missed the region where the two ellipse overlap
264 //if nsols_high=4, also need a scan because we may find the wrong tangent point.
265 
266  nsols_high = nsols(Deltasq_high);
267 
268  if(nsols_high == nsols_low || nsols_high == 4)
269  {
270  foundhigh = scan_high(Deltasq_high);
271  // foundhigh = find_high(Deltasq_high);
272  if(foundhigh == 0)
273  {
274  cout << "Deltasq_high not found, output ma+mn1" << endl;
275 // cout << "myhigh1 = " << sqrt(Deltasq_high1 + mn1sq) << endl;
276 // cout << "myhigh2 = " << sqrt(Deltasq_high2 + mn1sq) << endl;
277 // cout << "nsol(high1) = " << nsols(Deltasq_high1) << endl;
278 // cout << "nsol(high2) = " << nsols(Deltasq_high2) << endl;
279 
280 
281 
282  mt2bl_b = sqrt( Deltasq_low + mn1sq );
283  return;
284  }
285 
286  }
287 
288  while(sqrt(Deltasq_high+mn1sq) - sqrt(Deltasq_low+mn1sq) > precision)
289  {
290  double Deltasq_mid,nsols_mid;
291  //bisect
292  Deltasq_mid = (Deltasq_high+Deltasq_low)/2.;
293  nsols_mid = nsols(Deltasq_mid);
294  // if nsols_mid = 4, rescan for Deltasq_high
295  if ( nsols_mid == 4 )
296  {
297  Deltasq_high = Deltasq_mid;
298  // scan_high(Deltasq_high);
299  find_high(Deltasq_high);
300  continue;
301  }
302 
303 
304  if(nsols_mid != nsols_low) Deltasq_high = Deltasq_mid;
305  if(nsols_mid == nsols_low) Deltasq_low = Deltasq_mid;
306  }
307  mt2bl_b = (double) sqrt( mn1sq + Deltasq_high);
308  return;
309 }
310 
311 
312 //For some reason this sophiscicated way of finding the upper bound does not seem to work for mt2bl...
313 int mt2bl::find_high(double & Deltasq_high)
314 {
315  double x0,y0;
316  x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
317  y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
318  double Deltasq_low = ma*(ma + 2*mn1);
319  do
320  {
321  double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.;
322  int nsols_mid = nsols(Deltasq_mid);
323  if ( nsols_mid == 2 )
324  {
325  Deltasq_high = Deltasq_mid;
326  return 1;
327  }
328  else if (nsols_mid == 4)
329  {
330  Deltasq_high = Deltasq_mid;
331  continue;
332  }
333  else if (nsols_mid ==0)
334  {
335  d1 = -pax*(Deltasq_mid-masq)/(2*Easq);
336  e1 = -pay*(Deltasq_mid-masq)/(2*Easq);
337  d2 = -pmissx + pbx*(Deltasq_mid - mbsq - mn2sq + mn1sq )/(2*Ebsq)
338  + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
339  e2 = -pmissy + pby*(Deltasq_mid - mbsq - mn2sq + mn1sq)/(2*Ebsq)
340  + pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
341  f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq- mn2sq + mn1sq)/(2*Eb)+
342  (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq- mn2sq + mn1sq)/(2*Eb)+
343  (pbx*pmissx+pby*pmissy)/Eb)+mn2sq;
344 // Does the larger ellipse contain the smaller one?
345  double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
346  if (dis < 0) Deltasq_high = Deltasq_mid;
347  else Deltasq_low = Deltasq_mid;
348 
349  }
350 
351  } while ( Deltasq_high - Deltasq_low > 0.001);
352  return 0;
353 }
354 
355 //...so we just scan to find the upper bound if our initial guess does not work
356 int mt2bl::scan_high(double & Deltasq_high)
357 {
358  int foundhigh = 0 ;
359  int nsols_high;
360 
361 
362  double Deltasq_low;
363  double tempmass, maxmass;
364  tempmass = mn1 + ma;
365  maxmass = sqrt(mn1sq + Deltasq_high);
366  //if (nevt == 32334) cout << "Deltasq_high = " << Deltasq_high << endl;
367  for(double mass = tempmass + scan_step; mass < maxmass; mass += scan_step)
368  {
369  Deltasq_high = mass*mass - mn1sq;
370  nsols_high = nsols(Deltasq_high);
371 
372  if( nsols_high > 0)
373  {
374  Deltasq_low = (mass-scan_step)*(mass-scan_step) - mn1sq;
375  foundhigh = 1;
376  break;
377  }
378  }
379  return foundhigh;
380 }
381 
382 
383 
384 int mt2bl::nsols( double Dsq)
385 {
386  double delta = (Dsq-masq)/(2*Easq);
387  double delta2 = (Dsq-mbsq - mn2sq + mn1sq )/(2*Ebsq);
388 
389 
390 //calculate coefficients for the two quadratic equations
391  d1 = d11*delta;
392  e1 = e11*delta;
393  f1 = f12*delta*delta+f10;
394  d2 = d21*delta2+d20;
395  e2 = e21*delta2+e20;
396  f2 = f22*delta2*delta2+f21*delta2+f20;
397 
398 //obtain the coefficients for the 4th order equation
399 //devided by Ea^n to make the variable dimensionless
400  long double A4, A3, A2, A1, A0;
401 
402  A4 =
403  -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 +
404  4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 +
405  a1*a1*c2*c2;
406 
407  A3 =
408  (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 +
409  8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 +
410  4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 -
411  4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea;
412 
413 
414  A2 =
415  (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 -
416  8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 +
417  4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 -
418  8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 +
419  4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 +
420  4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq;
421 
422  A1 =
423  (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 -
424  4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 -
425  4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 -
426  4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea);
427 
428  A0 =
429  (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 +
430  4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 +
431  a1*a1*f2*f2)/(Easq*Easq);
432 
433  long double A0sq, A1sq, A2sq, A3sq, A4sq;
434  A0sq = A0*A0;
435  A1sq = A1*A1;
436  A2sq = A2*A2;
437  A3sq = A3*A3;
438  A4sq = A4*A4;
439 
440  long double B3, B2, B1, B0;
441  B3 = 4*A4;
442  B2 = 3*A3;
443  B1 = 2*A2;
444  B0 = A1;
445 
446  long double C2, C1, C0;
447  C2 = -(A2/2 - 3*A3sq/(16*A4));
448  C1 = -(3*A1/4. -A2*A3/(8*A4));
449  C0 = -A0 + A1*A3/(16*A4);
450 
451  long double D1, D0;
452  D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
453  D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
454 
455  long double E0;
456  E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
457 
458  long double t1,t2,t3,t4,t5;
459 //find the coefficients for the leading term in the Sturm sequence
460  t1 = A4;
461  t2 = A4;
462  t3 = C2;
463  t4 = D1;
464  t5 = E0;
465 
466 
467 //The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
468  int nsol;
469  nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);
470 
471 //Cannot have negative number of solutions, must be roundoff effect
472  if (nsol < 0) nsol = 0;
473 
474  return nsol;
475 
476 }
477 
478 inline int mt2bl::signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5)
479 {
480  int nsc;
481  nsc=0;
482  if(t1*t2>0) nsc++;
483  if(t2*t3>0) nsc++;
484  if(t3*t4>0) nsc++;
485  if(t4*t5>0) nsc++;
486  return nsc;
487 }
488 inline int mt2bl::signchange_p( long double t1, long double t2, long double t3, long double t4, long double t5)
489 {
490  int nsc;
491  nsc=0;
492  if(t1*t2<0) nsc++;
493  if(t2*t3<0) nsc++;
494  if(t3*t4<0) nsc++;
495  if(t4*t5<0) nsc++;
496  return nsc;
497 }
498 
499 }//end namespace mt2bl_bisect
#define RELATIVE_PRECISION
Definition: mt2_bisect.h:7
#define ABSOLUTE_PRECISION
Definition: mt2_bisect.h:8