checkmate is hosted by Hepforge, IPPP Durham
mt2_bisect.cc
Go to the documentation of this file.
1 /***********************************************************************/
2 /* */
3 /* Finding mt2 by Bisection */
4 /* */
5 /* Authors: Hsin-Chia Cheng, Zhenyu Han */
6 /* Dec 11, 2008, v1.01a */
7 /* */
8 /***********************************************************************/
9 
10 
11 /*******************************************************************************
12  Usage:
13 
14  1. Define an object of type "mt2":
15 
16  mt2_bisect::mt2 mt2_event;
17 
18  2. Set momenta and the mass of the invisible particle, mn:
19 
20  mt2_event.set_momenta( pa, pb, pmiss );
21  mt2_event.set_mass( mn );
22 
23  where array pa[0..2], pb[0..2], pmiss[0..2] contains (mass,px,py)
24  for the visible particles and the missing momentum. pmiss[0] is not used.
25  All quantities are given in double.
26 
27  3. Use mt2::get_mt2() to obtain the value of mt2:
28 
29  double mt2_value = mt2_event.get_mt2();
30 
31 *******************************************************************************/
32 
33 #include <iostream>
34 #include <math.h>
35 #include "mt2_bisect.h"
36 
37 using namespace std;
38 
39 namespace mt2_bisect
40 {
41 mt2::mt2()
42 {
43  solved = false;
44  momenta_set = false;
45  mt2_b = 0.;
46  scale = 1.;
47 }
48 
49 double mt2::get_mt2()
50 {
51  if (!momenta_set)
52  {
53  cout <<" Please set momenta first!" << endl;
54  return 0;
55  }
56 
57  if (!solved) mt2_bisect();
58  return mt2_b*scale;
59 }
60 
61 void mt2::set_momenta(double* pa0, double* pb0, double* pmiss0)
62 {
63  solved = false; //reset solved tag when momenta are changed.
64  momenta_set = true;
65 
66  ma = fabs(pa0[0]); // mass cannot be negative
67 
68  if (ma < ZERO_MASS) ma = ZERO_MASS;
69 
70  pax = pa0[1];
71  pay = pa0[2];
72  masq = ma*ma;
73  Easq = masq+pax*pax+pay*pay;
74  Ea = sqrt(Easq);
75 
76  mb = fabs(pb0[0]);
77 
78  if (mb < ZERO_MASS) mb = ZERO_MASS;
79 
80  pbx = pb0[1];
81  pby = pb0[2];
82  mbsq = mb*mb;
83  Ebsq = mbsq+pbx*pbx+pby*pby;
84  Eb = sqrt(Ebsq);
85 
86  pmissx = pmiss0[1]; pmissy = pmiss0[2];
87  pmissxsq = pmissx*pmissx;
88  pmissysq = pmissy*pmissy;
89 
90 // set ma>= mb
91  if(masq < mbsq)
92  {
93  double temp;
94  temp = pax; pax = pbx; pbx = temp;
95  temp = pay; pay = pby; pby = temp;
96  temp = Ea; Ea = Eb; Eb = temp;
97  temp = Easq; Easq = Ebsq; Ebsq = temp;
98  temp = masq; masq = mbsq; mbsq = temp;
99  temp = ma; ma = mb; mb = temp;
100  }
101 //normalize max{Ea, Eb} to 100
102  if (Ea > Eb) scale = Ea/100.;
103  else scale = Eb/100.;
104 
105  if (sqrt(pmissxsq+pmissysq)/100 > scale) scale = sqrt(pmissxsq+pmissysq)/100;
106  //scale = 1;
107  double scalesq = scale * scale;
108  ma = ma/scale;
109  mb = mb/scale;
110  masq = masq/scalesq;
111  mbsq = mbsq/scalesq;
112  pax = pax/scale; pay = pay/scale;
113  pbx = pbx/scale; pby = pby/scale;
114  Ea = Ea/scale; Eb = Eb/scale;
115 
116  Easq = Easq/scalesq;
117  Ebsq = Ebsq/scalesq;
118  pmissx = pmissx/scale;
119  pmissy = pmissy/scale;
120  pmissxsq = pmissxsq/scalesq;
121  pmissysq = pmissysq/scalesq;
122  mn = mn_unscale/scale;
123  mnsq = mn*mn;
124 
126  else precision = 100.*RELATIVE_PRECISION;
127 }
128 
129 void mt2::set_mn(double mn0)
130 {
131  solved = false; //reset solved tag when mn is changed.
132  mn_unscale = fabs(mn0); //mass cannot be negative
133  mn = mn_unscale/scale;
134  mnsq = mn*mn;
135 }
136 
137 void mt2::print()
138 {
139  cout << " pax = " << pax*scale << "; pay = " << pay*scale << "; ma = " << ma*scale <<";"<< endl;
140  cout << " pbx = " << pbx*scale << "; pby = " << pby*scale << "; mb = " << mb*scale <<";"<< endl;
141  cout << " pmissx = " << pmissx*scale << "; pmissy = " << pmissy*scale <<";"<< endl;
142  cout << " mn = " << mn_unscale<<";" << endl;
143 }
144 
145 //special case, the visible particle is massless
146 void mt2::mt2_massless()
147 {
148 
149 //rotate so that pay = 0
150  double theta,s,c;
151  theta = atan(pay/pax);
152  s = sin(theta);
153  c = cos(theta);
154 
155  double pxtemp,pytemp;
156  Easq = pax*pax+pay*pay;
157  Ebsq = pbx*pbx+pby*pby;
158  Ea = sqrt(Easq);
159  Eb = sqrt(Ebsq);
160 
161  pxtemp = pax*c+pay*s;
162  pax = pxtemp;
163  pay = 0;
164  pxtemp = pbx*c+pby*s;
165  pytemp = -s*pbx+c*pby;
166  pbx = pxtemp;
167  pby = pytemp;
168  pxtemp = pmissx*c+pmissy*s;
169  pytemp = -s*pmissx+c*pmissy;
170  pmissx = pxtemp;
171  pmissy = pytemp;
172 
173  a2 = 1-pbx*pbx/(Ebsq);
174  b2 = -pbx*pby/(Ebsq);
175  c2 = 1-pby*pby/(Ebsq);
176 
177  d21 = (Easq*pbx)/Ebsq;
178  d20 = - pmissx + (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
179  e21 = (Easq*pby)/Ebsq;
180  e20 = - pmissy + (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
181  f22 = -(Easq*Easq/Ebsq);
182  f21 = -2*Easq*(pbx*pmissx + pby*pmissy)/Ebsq;
183  f20 = mnsq + pmissxsq + pmissysq - (pbx*pmissx + pby*pmissy)*(pbx*pmissx + pby*pmissy)/Ebsq;
184 
185  double Deltasq0 = 0;
186  double Deltasq_low, Deltasq_high;
187  int nsols_high, nsols_low;
188 
189  Deltasq_low = Deltasq0 + precision;
190  nsols_low = nsols_massless(Deltasq_low);
191 
192  if(nsols_low > 1)
193  {
194  mt2_b = (double) sqrt(Deltasq0+mnsq);
195  return;
196  }
197 
198 /*
199  if( nsols_massless(Deltasq_high) > 0 )
200  {
201  mt2_b = (double) sqrt(mnsq+Deltasq0);
202  return;
203  }*/
204 
205 //look for when both parablos contain origin
206  double Deltasq_high1, Deltasq_high2;
207  Deltasq_high1 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy;
208  Deltasq_high2 = 2*Ea*mn;
209 
210  if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high2;
211  else Deltasq_high = Deltasq_high1;
212 
213  nsols_high=nsols_massless(Deltasq_high);
214 
215  int foundhigh;
216  if (nsols_high == nsols_low)
217  {
218 
219 
220  foundhigh=0;
221 
222  double minmass, maxmass;
223  minmass = mn ;
224  maxmass = sqrt(mnsq + Deltasq_high);
225  for(double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
226  {
227  Deltasq_high = mass*mass - mnsq;
228 
229  nsols_high = nsols_massless(Deltasq_high);
230  if(nsols_high>0)
231  {
232  foundhigh=1;
233  Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
234  break;
235  }
236  }
237  if(foundhigh==0)
238  {
239 
240  cout<<"Deltasq_high not found at event " << nevt <<endl;
241 
242 
243  mt2_b = (double)sqrt(Deltasq_low+mnsq);
244  return;
245  }
246  }
247 
248  if(nsols_high == nsols_low)
249  {
250  cout << "error: nsols_low=nsols_high=" << nsols_high << endl;
251  cout << "Deltasq_high=" << Deltasq_high << endl;
252  cout << "Deltasq_low= "<< Deltasq_low << endl;
253 
254  mt2_b = sqrt(mnsq + Deltasq_low);
255  return;
256  }
257  double minmass, maxmass;
258  minmass = sqrt(Deltasq_low+mnsq);
259  maxmass = sqrt(Deltasq_high+mnsq);
260  while(maxmass - minmass > precision)
261  {
262  double Delta_mid, midmass, nsols_mid;
263  midmass = (minmass+maxmass)/2.;
264  Delta_mid = midmass * midmass - mnsq;
265  nsols_mid = nsols_massless(Delta_mid);
266  if(nsols_mid != nsols_low) maxmass = midmass;
267  if(nsols_mid == nsols_low) minmass = midmass;
268  }
269  mt2_b = minmass;
270  return;
271 }
272 
273 int mt2::nsols_massless(double Dsq)
274 {
275  double delta;
276  delta = Dsq/(2*Easq);
277  d1 = d11*delta;
278  e1 = e11*delta;
279  f1 = f12*delta*delta+f10;
280  d2 = d21*delta+d20;
281  e2 = e21*delta+e20;
282  f2 = f22*delta*delta+f21*delta+f20;
283 
284  double a,b;
285  if (pax > 0) a = Ea/Dsq;
286  else a = -Ea/Dsq;
287  if (pax > 0) b = -Dsq/(4*Ea)+mnsq*Ea/Dsq;
288  else b = Dsq/(4*Ea)-mnsq*Ea/Dsq;
289 
290  double A4,A3,A2,A1,A0;
291 
292  A4 = a*a*a2;
293  A3 = 2*a*b2/Ea;
294  A2 = (2*a*a2*b+c2+2*a*d2)/(Easq);
295  A1 = (2*b*b2+2*e2)/(Easq*Ea);
296  A0 = (a2*b*b+2*b*d2+f2)/(Easq*Easq);
297 
298  long double A0sq, A1sq, A2sq, A3sq, A4sq;
299  A0sq = A0*A0;
300  A1sq = A1*A1;
301  A2sq = A2*A2;
302  A3sq = A3*A3;
303  A4sq = A4*A4;
304 
305  long double B3, B2, B1, B0;
306  B3 = 4*A4;
307  B2 = 3*A3;
308  B1 = 2*A2;
309  B0 = A1;
310  long double C2, C1, C0;
311  C2 = -(A2/2 - 3*A3sq/(16*A4));
312  C1 = -(3*A1/4. -A2*A3/(8*A4));
313  C0 = -A0 + A1*A3/(16*A4);
314  long double D1, D0;
315  D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
316  D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
317 
318  long double E0;
319  E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
320 
321  long double t1,t2,t3,t4,t5;
322 
323 //find the coefficients for the leading term in the Sturm sequence
324  t1 = A4;
325  t2 = A4;
326  t3 = C2;
327  t4 = D1;
328  t5 = E0;
329 
330  int nsol;
331  nsol = signchange_n(t1,t2,t3,t4,t5)-signchange_p(t1,t2,t3,t4,t5);
332  if( nsol < 0 ) nsol=0;
333 
334  return nsol;
335 
336 }
337 
338 void mt2::mt2_bisect()
339 {
340 
341 
342  solved = true;
343  cout.precision(11);
344 
345 //if masses are very small, use code for massless case.
346  if(masq < MIN_MASS && mbsq < MIN_MASS)
347  {
348  mt2_massless();
349  return;
350  }
351 
352 
353  double Deltasq0;
354  Deltasq0 = ma*(ma + 2*mn); //The minimum mass square to have two ellipses
355 
356 // find the coefficients for the two quadratic equations when Deltasq=Deltasq0.
357 
358  a1 = 1-pax*pax/(Easq);
359  b1 = -pax*pay/(Easq);
360  c1 = 1-pay*pay/(Easq);
361  d1 = -pax*(Deltasq0-masq)/(2*Easq);
362  e1 = -pay*(Deltasq0-masq)/(2*Easq);
363  a2 = 1-pbx*pbx/(Ebsq);
364  b2 = -pbx*pby/(Ebsq);
365  c2 = 1-pby*pby/(Ebsq);
366  d2 = -pmissx+pbx*(Deltasq0-mbsq)/(2*Ebsq)+pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
367  e2 = -pmissy+pby*(Deltasq0-mbsq)/(2*Ebsq)+pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
368  f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq0-mbsq)/(2*Eb)+
369  (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq0-mbsq)/(2*Eb)+
370  (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
371 
372 // find the center of the smaller ellipse
373  double x0,y0;
374  x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
375  y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
376 
377 
378 // Does the larger ellipse contain the smaller one?
379  double dis=a2*x0*x0+2*b2*x0*y0+c2*y0*y0+2*d2*x0+2*e2*y0+f2;
380 
381  if(dis<=0.01)
382  {
383  mt2_b = (double) sqrt(mnsq+Deltasq0);
384  return;
385  }
386 
387 
388 /* find the coefficients for the two quadratic equations */
389 /* coefficients for quadratic terms do not change */
390 /* coefficients for linear and constant terms are polynomials of */
391 /* delta=(Deltasq-m7sq)/(2 E7sq) */
392  d11 = -pax;
393  e11 = -pay;
394  f10 = mnsq;
395  f12 = -Easq;
396  d21 = (Easq*pbx)/Ebsq;
397  d20 = ((masq - mbsq)*pbx)/(2.*Ebsq) - pmissx +
398  (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
399  e21 = (Easq*pby)/Ebsq;
400  e20 = ((masq - mbsq)*pby)/(2.*Ebsq) - pmissy +
401  (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
402  f22 = -Easq*Easq/Ebsq;
403  f21 = (-2*Easq*((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb))/Eb;
404  f20 = mnsq + pmissx*pmissx + pmissy*pmissy -
405  ((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb)
406  *((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb);
407 
408 //Estimate upper bound of mT2
409 //when Deltasq > Deltasq_high1, the larger encloses the center of the smaller
410  double p2x0,p2y0;
411  double Deltasq_high1;
412  p2x0 = pmissx-x0;
413  p2y0 = pmissy-y0;
414  Deltasq_high1 = 2*Eb*sqrt(p2x0*p2x0+p2y0*p2y0+mnsq)-2*pbx*p2x0-2*pby*p2y0+mbsq;
415 
416 //Another estimate, if both ellipses enclose the origin, Deltasq > mT2
417 
418  double Deltasq_high2, Deltasq_high21, Deltasq_high22;
419  Deltasq_high21 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy+mbsq;
420  Deltasq_high22 = 2*Ea*mn+masq;
421 
422  if ( Deltasq_high21 < Deltasq_high22 ) Deltasq_high2 = Deltasq_high22;
423  else Deltasq_high2 = Deltasq_high21;
424 
425 //pick the smaller upper bound
426  double Deltasq_high;
427  if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high1;
428  else Deltasq_high = Deltasq_high2;
429 
430 
431  double Deltasq_low; //lower bound
432  Deltasq_low = Deltasq0;
433 
434 //number of solutions at Deltasq_low should not be larger than zero
435  if( nsols(Deltasq_low) > 0 )
436  {
437  //cout << "nsolutions(Deltasq_low) > 0"<<endl;
438  mt2_b = (double) sqrt(mnsq+Deltasq0);
439  return;
440  }
441 
442  int nsols_high, nsols_low;
443 
444  nsols_low = nsols(Deltasq_low);
445  int foundhigh;
446 
447 
448 //if nsols_high=nsols_low, we missed the region where the two ellipse overlap
449 //if nsols_high=4, also need a scan because we may find the wrong tangent point.
450 
451  nsols_high = nsols(Deltasq_high);
452 
453  if(nsols_high == nsols_low || nsols_high == 4)
454  {
455  //foundhigh = scan_high(Deltasq_high);
456  foundhigh = find_high(Deltasq_high);
457  if(foundhigh == 0)
458  {
459  cout << "Deltasq_high not found at event " << nevt << endl;
460  mt2_b = sqrt( Deltasq_low + mnsq );
461  return;
462  }
463 
464  }
465 
466  while(sqrt(Deltasq_high+mnsq) - sqrt(Deltasq_low+mnsq) > precision)
467  {
468  double Deltasq_mid,nsols_mid;
469  //bisect
470  Deltasq_mid = (Deltasq_high+Deltasq_low)/2.;
471  nsols_mid = nsols(Deltasq_mid);
472  // if nsols_mid = 4, rescan for Deltasq_high
473  if ( nsols_mid == 4 )
474  {
475  Deltasq_high = Deltasq_mid;
476  //scan_high(Deltasq_high);
477  find_high(Deltasq_high);
478  continue;
479  }
480 
481 
482  if(nsols_mid != nsols_low) Deltasq_high = Deltasq_mid;
483  if(nsols_mid == nsols_low) Deltasq_low = Deltasq_mid;
484  }
485  mt2_b = (double) sqrt( mnsq + Deltasq_high);
486  return;
487 }
488 
489 int mt2::find_high(double & Deltasq_high)
490 {
491  double x0,y0;
492  x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
493  y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
494  double Deltasq_low = (mn + ma)*(mn + ma) - mnsq;
495  do
496  {
497  double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.;
498  int nsols_mid = nsols(Deltasq_mid);
499  if ( nsols_mid == 2 )
500  {
501  Deltasq_high = Deltasq_mid;
502  return 1;
503  }
504  else if (nsols_mid == 4)
505  {
506  Deltasq_high = Deltasq_mid;
507  continue;
508  }
509  else if (nsols_mid ==0)
510  {
511  d1 = -pax*(Deltasq_mid-masq)/(2*Easq);
512  e1 = -pay*(Deltasq_mid-masq)/(2*Easq);
513  d2 = -pmissx + pbx*(Deltasq_mid - mbsq)/(2*Ebsq)
514  + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
515  e2 = -pmissy + pby*(Deltasq_mid - mbsq)/(2*Ebsq)
516  + pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
517  f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq)/(2*Eb)+
518  (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq)/(2*Eb)+
519  (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
520 // Does the larger ellipse contain the smaller one?
521  double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
522  if (dis < 0) Deltasq_high = Deltasq_mid;
523  else Deltasq_low = Deltasq_mid;
524  }
525 
526  } while ( Deltasq_high - Deltasq_low > 0.001);
527  return 0;
528 }
529 int mt2::scan_high(double & Deltasq_high)
530 {
531  int foundhigh = 0 ;
532  int nsols_high;
533 
534 
535  double Deltasq_low;
536  double tempmass, maxmass;
537  tempmass = mn + ma;
538  maxmass = sqrt(mnsq + Deltasq_high);
539  if (nevt == 32334) cout << "Deltasq_high = " << Deltasq_high << endl;
540  for(double mass = tempmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
541  {
542  Deltasq_high = mass*mass - mnsq;
543  nsols_high = nsols(Deltasq_high);
544 
545  if( nsols_high > 0)
546  {
547  Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
548  foundhigh = 1;
549  break;
550  }
551  }
552  return foundhigh;
553 }
554 int mt2::nsols( double Dsq)
555 {
556  double delta = (Dsq-masq)/(2*Easq);
557 
558 //calculate coefficients for the two quadratic equations
559  d1 = d11*delta;
560  e1 = e11*delta;
561  f1 = f12*delta*delta+f10;
562  d2 = d21*delta+d20;
563  e2 = e21*delta+e20;
564  f2 = f22*delta*delta+f21*delta+f20;
565 
566 //obtain the coefficients for the 4th order equation
567 //devided by Ea^n to make the variable dimensionless
568  long double A4, A3, A2, A1, A0;
569 
570  A4 =
571  -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 +
572  4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 +
573  a1*a1*c2*c2;
574 
575  A3 =
576  (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 +
577  8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 +
578  4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 -
579  4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea;
580 
581 
582  A2 =
583  (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 -
584  8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 +
585  4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 -
586  8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 +
587  4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 +
588  4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq;
589 
590  A1 =
591  (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 -
592  4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 -
593  4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 -
594  4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea);
595 
596  A0 =
597  (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 +
598  4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 +
599  a1*a1*f2*f2)/(Easq*Easq);
600 
601  long double A0sq, A1sq, A2sq, A3sq, A4sq;
602  A0sq = A0*A0;
603  A1sq = A1*A1;
604  A2sq = A2*A2;
605  A3sq = A3*A3;
606  A4sq = A4*A4;
607 
608  long double B3, B2, B1, B0;
609  B3 = 4*A4;
610  B2 = 3*A3;
611  B1 = 2*A2;
612  B0 = A1;
613 
614  long double C2, C1, C0;
615  C2 = -(A2/2 - 3*A3sq/(16*A4));
616  C1 = -(3*A1/4. -A2*A3/(8*A4));
617  C0 = -A0 + A1*A3/(16*A4);
618 
619  long double D1, D0;
620  D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
621  D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
622 
623  long double E0;
624  E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
625 
626  long double t1,t2,t3,t4,t5;
627 //find the coefficients for the leading term in the Sturm sequence
628  t1 = A4;
629  t2 = A4;
630  t3 = C2;
631  t4 = D1;
632  t5 = E0;
633 
634 
635 //The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
636  int nsol;
637  nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);
638 
639 //Cannot have negative number of solutions, must be roundoff effect
640  if (nsol < 0) nsol = 0;
641 
642  return nsol;
643 
644 }
645 
646 inline int mt2::signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5)
647 {
648  int nsc;
649  nsc=0;
650  if(t1*t2>0) nsc++;
651  if(t2*t3>0) nsc++;
652  if(t3*t4>0) nsc++;
653  if(t4*t5>0) nsc++;
654  return nsc;
655 }
656 inline int mt2::signchange_p( long double t1, long double t2, long double t3, long double t4, long double t5)
657 {
658  int nsc;
659  nsc=0;
660  if(t1*t2<0) nsc++;
661  if(t2*t3<0) nsc++;
662  if(t3*t4<0) nsc++;
663  if(t4*t5<0) nsc++;
664  return nsc;
665 }
666 
667 }//end namespace mt2_bisect
#define SCANSTEP
Definition: mt2_bisect.h:14
#define MIN_MASS
Definition: mt2_bisect.h:12
#define RELATIVE_PRECISION
Definition: mt2_bisect.h:7
#define ABSOLUTE_PRECISION
Definition: mt2_bisect.h:8
#define ZERO_MASS
Definition: mt2_bisect.h:13