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