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