Heppy/src/Davismt2.cc

0001 // #ifndef MT2_BISECT_C
0002 // #define MT2_BISECT_C
0003 /***********************************************************************/
0004 /*                                                                     */
0005 /*              Finding mt2 by Bisection                               */
0006 /*                                                                     */
0007 /*              Authors: Hsin-Chia Cheng, Zhenyu Han                   */
0008 /*              Dec 11, 2008, v1.01a                                   */
0009 /*                                                                     */
0010 /***********************************************************************/
0011
0012 /*******************************************************************************
0013 Usage:
0014
0015 1. Define an object of type "mt2":
0016
0017 mt2_bisect::mt2 mt2_event;
0018
0019 2. Set momenta and the mass of the invisible particle, mn:
0020
0021 mt2_event.set_momenta( pa, pb, pmiss );
0022 mt2_event.set_mass( mn );
0023
0024 where array pa[0..2], pb[0..2], pmiss[0..2] contains (mass,px,py)
0025     for the visible particles and the missing momentum. pmiss[0] is not used.
0026     All quantities are given in double.
0027
0028     3. Use Davismt2::get_mt2() to obtain the value of mt2:
0029
0030 double mt2_value = mt2_event.get_mt2();
0031
0032 *******************************************************************************/
0033
0034 #include "PhysicsTools/Heppy/interface/Davismt2.h"
0035 // ClassImp(Davismt2);
0036
0037 using namespace std;
0038
0039 namespace heppy {
0040
0041   /*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
0042   achievable--use with caution*/
0043
0044   const float Davismt2::RELATIVE_PRECISION = 0.00001;
0045   const float Davismt2::ABSOLUTE_PRECISION = 0.0;
0046   const float Davismt2::ZERO_MASS = 0.0;
0047   const float Davismt2::MIN_MASS = 0.1;
0048   const float Davismt2::SCANSTEP = 0.1;
0049
0050   Davismt2::Davismt2() {
0051     solved = false;
0052     momenta_set = false;
0053     mt2_b = 0.;
0054     scale = 1.;
0055     verbose = 1;
0056   }
0057
0058   Davismt2::~Davismt2() {}
0059
0060   double Davismt2::get_mt2() {
0061     if (!momenta_set) {
0062       cout << "Davismt2::get_mt2() ==> Please set momenta first!" << endl;
0063       return 0;
0064     }
0065
0066     if (!solved)
0067       mt2_bisect();
0068     return mt2_b * scale;
0069   }
0070
0071   void Davismt2::set_momenta(double* pa0, double* pb0, double* pmiss0) {
0072     solved = false;  //reset solved tag when momenta are changed.
0073     momenta_set = true;
0074
0075     ma = fabs(pa0[0]);  // mass cannot be negative
0076
0077     if (ma < ZERO_MASS)
0078       ma = ZERO_MASS;
0079
0080     pax = pa0[1];
0081     pay = pa0[2];
0082     masq = ma * ma;
0083     Easq = masq + pax * pax + pay * pay;
0084     Ea = sqrt(Easq);
0085
0086     mb = fabs(pb0[0]);
0087
0088     if (mb < ZERO_MASS)
0089       mb = ZERO_MASS;
0090
0091     pbx = pb0[1];
0092     pby = pb0[2];
0093     mbsq = mb * mb;
0094     Ebsq = mbsq + pbx * pbx + pby * pby;
0095     Eb = sqrt(Ebsq);
0096
0097     pmissx = pmiss0[1];
0098     pmissy = pmiss0[2];
0099     pmissxsq = pmissx * pmissx;
0100     pmissysq = pmissy * pmissy;
0101
0102     // set ma>= mb
0103     if (masq < mbsq) {
0104       double temp;
0105       temp = pax;
0106       pax = pbx;
0107       pbx = temp;
0108       temp = pay;
0109       pay = pby;
0110       pby = temp;
0111       temp = Ea;
0112       Ea = Eb;
0113       Eb = temp;
0114       temp = Easq;
0115       Easq = Ebsq;
0116       Ebsq = temp;
0117       temp = masq;
0118       masq = mbsq;
0119       mbsq = temp;
0120       temp = ma;
0121       ma = mb;
0122       mb = temp;
0123     }
0124     //normalize max{Ea, Eb} to 100
0125     if (Ea > Eb)
0126       scale = Ea / 100.;
0127     else
0128       scale = Eb / 100.;
0129
0130     if (sqrt(pmissxsq + pmissysq) / 100 > scale)
0131       scale = sqrt(pmissxsq + pmissysq) / 100;
0132     //scale = 1;
0133     double scalesq = scale * scale;
0134     ma = ma / scale;
0135     mb = mb / scale;
0136     masq = masq / scalesq;
0137     mbsq = mbsq / scalesq;
0138     pax = pax / scale;
0139     pay = pay / scale;
0140     pbx = pbx / scale;
0141     pby = pby / scale;
0142     Ea = Ea / scale;
0143     Eb = Eb / scale;
0144
0145     Easq = Easq / scalesq;
0146     Ebsq = Ebsq / scalesq;
0147     pmissx = pmissx / scale;
0148     pmissy = pmissy / scale;
0149     pmissxsq = pmissxsq / scalesq;
0150     pmissysq = pmissysq / scalesq;
0151     mn = mn_unscale / scale;
0152     mnsq = mn * mn;
0153
0154     if (ABSOLUTE_PRECISION > 100. * RELATIVE_PRECISION)
0155       precision = ABSOLUTE_PRECISION;
0156     else
0157       precision = 100. * RELATIVE_PRECISION;
0158   }
0159
0160   void Davismt2::set_mn(double mn0) {
0161     solved = false;          //reset solved tag when mn is changed.
0162     mn_unscale = fabs(mn0);  //mass cannot be negative
0163     mn = mn_unscale / scale;
0164     mnsq = mn * mn;
0165   }
0166
0167   void Davismt2::print() {
0168     cout << " Davismt2::print() ==> pax = " << pax * scale << ";   pay = " << pay * scale << ";   ma = " << ma * scale
0169          << ";" << endl;
0170     cout << " Davismt2::print() ==> pbx = " << pbx * scale << ";   pby = " << pby * scale << ";   mb = " << mb * scale
0171          << ";" << endl;
0172     cout << " Davismt2::print() ==> pmissx = " << pmissx * scale << ";   pmissy = " << pmissy * scale << ";" << endl;
0173     cout << " Davismt2::print() ==> mn = " << mn_unscale << ";" << endl;
0174   }
0175
0176   //special case, the visible particle is massless
0177   void Davismt2::mt2_massless() {
0178     //rotate so that pay = 0
0179     double theta, s, c;
0180     theta = atan(pay / pax);
0181     s = sin(theta);
0182     c = cos(theta);
0183
0184     double pxtemp, pytemp;
0185     Easq = pax * pax + pay * pay;
0186     Ebsq = pbx * pbx + pby * pby;
0187     Ea = sqrt(Easq);
0188     Eb = sqrt(Ebsq);
0189
0190     pxtemp = pax * c + pay * s;
0191     pax = pxtemp;
0192     pay = 0;
0193     pxtemp = pbx * c + pby * s;
0194     pytemp = -s * pbx + c * pby;
0195     pbx = pxtemp;
0196     pby = pytemp;
0197     pxtemp = pmissx * c + pmissy * s;
0198     pytemp = -s * pmissx + c * pmissy;
0199     pmissx = pxtemp;
0200     pmissy = pytemp;
0201
0202     a2 = 1 - pbx * pbx / (Ebsq);
0203     b2 = -pbx * pby / (Ebsq);
0204     c2 = 1 - pby * pby / (Ebsq);
0205
0206     d21 = (Easq * pbx) / Ebsq;
0207     d20 = -pmissx + (pbx * (pbx * pmissx + pby * pmissy)) / Ebsq;
0208     e21 = (Easq * pby) / Ebsq;
0209     e20 = -pmissy + (pby * (pbx * pmissx + pby * pmissy)) / Ebsq;
0210     f22 = -(Easq * Easq / Ebsq);
0211     f21 = -2 * Easq * (pbx * pmissx + pby * pmissy) / Ebsq;
0212     f20 = mnsq + pmissxsq + pmissysq - (pbx * pmissx + pby * pmissy) * (pbx * pmissx + pby * pmissy) / Ebsq;
0213
0214     double Deltasq0 = 0;
0215     double Deltasq_low, Deltasq_high;
0216     int nsols_high, nsols_low;
0217
0218     Deltasq_low = Deltasq0 + precision;
0219     nsols_low = nsols_massless(Deltasq_low);
0220
0221     if (nsols_low > 1) {
0222       mt2_b = (double)sqrt(Deltasq0 + mnsq);
0223       return;
0224     }
0225
0226     /*
0227     if( nsols_massless(Deltasq_high) > 0 )
0228     {
0229         mt2_b = (double) sqrt(mnsq+Deltasq0);
0230         return;
0231         }*/
0232
0233     //look for when both parablos contain origin
0234     double Deltasq_high1, Deltasq_high2;
0235     Deltasq_high1 = 2 * Eb * sqrt(pmissx * pmissx + pmissy * pmissy + mnsq) - 2 * pbx * pmissx - 2 * pby * pmissy;
0236     Deltasq_high2 = 2 * Ea * mn;
0237
0238     if (Deltasq_high1 < Deltasq_high2)
0239       Deltasq_high = Deltasq_high2;
0240     else
0241       Deltasq_high = Deltasq_high1;
0242
0243     nsols_high = nsols_massless(Deltasq_high);
0244
0245     int foundhigh;
0246     if (nsols_high == nsols_low) {
0247       foundhigh = 0;
0248
0249       double minmass, maxmass;
0250       minmass = mn;
0251       maxmass = sqrt(mnsq + Deltasq_high);
0252       for (double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP) {
0253         Deltasq_high = mass * mass - mnsq;
0254
0255         nsols_high = nsols_massless(Deltasq_high);
0256         if (nsols_high > 0) {
0257           foundhigh = 1;
0258           Deltasq_low = (mass - SCANSTEP) * (mass - SCANSTEP) - mnsq;
0259           break;
0260         }
0261       }
0262       if (foundhigh == 0) {
0263         if (verbose > 0)
0264           cout << "Davismt2::mt2_massless() ==> Deltasq_high not found at event " << nevt << endl;
0265
0266         mt2_b = (double)sqrt(Deltasq_low + mnsq);
0267         return;
0268       }
0269     }
0270
0271     if (nsols_high == nsols_low) {
0272       if (verbose > 0)
0273         cout << "Davismt2::mt2_massless() ==> error: nsols_low=nsols_high=" << nsols_high << endl;
0274       if (verbose > 0)
0275         cout << "Davismt2::mt2_massless() ==> Deltasq_high=" << Deltasq_high << endl;
0276       if (verbose > 0)
0277         cout << "Davismt2::mt2_massless() ==> Deltasq_low= " << Deltasq_low << endl;
0278
0279       mt2_b = sqrt(mnsq + Deltasq_low);
0280       return;
0281     }
0282     double minmass, maxmass;
0283     minmass = sqrt(Deltasq_low + mnsq);
0284     maxmass = sqrt(Deltasq_high + mnsq);
0285     while (maxmass - minmass > precision) {
0286       double Delta_mid, midmass, nsols_mid;
0287       midmass = (minmass + maxmass) / 2.;
0288       Delta_mid = midmass * midmass - mnsq;
0289       nsols_mid = nsols_massless(Delta_mid);
0290       if (nsols_mid != nsols_low)
0291         maxmass = midmass;
0292       if (nsols_mid == nsols_low)
0293         minmass = midmass;
0294     }
0295     mt2_b = minmass;
0296     return;
0297   }
0298
0299   int Davismt2::nsols_massless(double Dsq) {
0300     double delta;
0301     delta = Dsq / (2 * Easq);
0302     d1 = d11 * delta;
0303     e1 = e11 * delta;
0304     f1 = f12 * delta * delta + f10;
0305     d2 = d21 * delta + d20;
0306     e2 = e21 * delta + e20;
0307     f2 = f22 * delta * delta + f21 * delta + f20;
0308
0309     double a, b;
0310     if (pax > 0)
0311       a = Ea / Dsq;
0312     else
0313       a = -Ea / Dsq;
0314     if (pax > 0)
0315       b = -Dsq / (4 * Ea) + mnsq * Ea / Dsq;
0316     else
0317       b = Dsq / (4 * Ea) - mnsq * Ea / Dsq;
0318
0319     double A4, A3, A2, A1, A0;
0320
0321     A4 = a * a * a2;
0322     A3 = 2 * a * b2 / Ea;
0323     A2 = (2 * a * a2 * b + c2 + 2 * a * d2) / (Easq);
0324     A1 = (2 * b * b2 + 2 * e2) / (Easq * Ea);
0325     A0 = (a2 * b * b + 2 * b * d2 + f2) / (Easq * Easq);
0326
0327     long double A3sq;
0328     A3sq = A3 * A3;
0329     /*
0330         long  double A0sq, A1sq, A2sq, A3sq, A4sq;
0331     A0sq = A0*A0;
0332     A1sq = A1*A1;
0333     A2sq = A2*A2;
0334     A3sq = A3*A3;
0335     A4sq = A4*A4;
0336         */
0337
0338     long double B3, B2, B1, B0;
0339     B3 = 4 * A4;
0340     B2 = 3 * A3;
0341     B1 = 2 * A2;
0342     B0 = A1;
0343     long double C2, C1, C0;
0344     C2 = -(A2 / 2 - 3 * A3sq / (16 * A4));
0345     C1 = -(3 * A1 / 4. - A2 * A3 / (8 * A4));
0346     C0 = -A0 + A1 * A3 / (16 * A4);
0347     long double D1, D0;
0348     D1 = -B1 - (B3 * C1 * C1 / C2 - B3 * C0 - B2 * C1) / C2;
0349     D0 = -B0 - B3 * C0 * C1 / (C2 * C2) + B2 * C0 / C2;
0350
0351     long double E0;
0352     E0 = -C0 - C2 * D0 * D0 / (D1 * D1) + C1 * D0 / D1;
0353
0354     long double t1, t2, t3, t4, t5;
0355
0356     //find the coefficients for the leading term in the Sturm sequence
0357     t1 = A4;
0358     t2 = A4;
0359     t3 = C2;
0360     t4 = D1;
0361     t5 = E0;
0362
0363     int nsol;
0364     nsol = signchange_n(t1, t2, t3, t4, t5) - signchange_p(t1, t2, t3, t4, t5);
0365     if (nsol < 0)
0366       nsol = 0;
0367
0368     return nsol;
0369   }
0370
0371   void Davismt2::mt2_bisect() {
0372     solved = true;
0373     cout.precision(11);
0374
0375     //if masses are very small, use code for massless case.
0376     if (masq < MIN_MASS && mbsq < MIN_MASS) {
0377       mt2_massless();
0378       return;
0379     }
0380
0381     double Deltasq0;
0382     Deltasq0 = ma * (ma + 2 * mn);  //The minimum mass square to have two ellipses
0383
0384     // find the coefficients for the two quadratic equations when Deltasq=Deltasq0.
0385
0386     a1 = 1 - pax * pax / (Easq);
0387     b1 = -pax * pay / (Easq);
0388     c1 = 1 - pay * pay / (Easq);
0389     d1 = -pax * (Deltasq0 - masq) / (2 * Easq);
0390     e1 = -pay * (Deltasq0 - masq) / (2 * Easq);
0391     a2 = 1 - pbx * pbx / (Ebsq);
0392     b2 = -pbx * pby / (Ebsq);
0393     c2 = 1 - pby * pby / (Ebsq);
0394     d2 = -pmissx + pbx * (Deltasq0 - mbsq) / (2 * Ebsq) + pbx * (pbx * pmissx + pby * pmissy) / (Ebsq);
0395     e2 = -pmissy + pby * (Deltasq0 - mbsq) / (2 * Ebsq) + pby * (pbx * pmissx + pby * pmissy) / (Ebsq);
0396     f2 = pmissx * pmissx + pmissy * pmissy -
0397          ((Deltasq0 - mbsq) / (2 * Eb) + (pbx * pmissx + pby * pmissy) / Eb) *
0398              ((Deltasq0 - mbsq) / (2 * Eb) + (pbx * pmissx + pby * pmissy) / Eb) +
0399          mnsq;
0400
0401     // find the center of the smaller ellipse
0402     double x0, y0;
0403     x0 = (c1 * d1 - b1 * e1) / (b1 * b1 - a1 * c1);
0404     y0 = (a1 * e1 - b1 * d1) / (b1 * b1 - a1 * c1);
0405
0406     // Does the larger ellipse contain the smaller one?
0407     double dis = a2 * x0 * x0 + 2 * b2 * x0 * y0 + c2 * y0 * y0 + 2 * d2 * x0 + 2 * e2 * y0 + f2;
0408
0409     if (dis <= 0.01) {
0410       mt2_b = (double)sqrt(mnsq + Deltasq0);
0411       return;
0412     }
0413
0414     /* find the coefficients for the two quadratic equations           */
0415     /* coefficients for quadratic terms do not change                  */
0416     /* coefficients for linear and constant terms are polynomials of   */
0417     /*       delta=(Deltasq-m7sq)/(2 E7sq)                             */
0418     d11 = -pax;
0419     e11 = -pay;
0420     f10 = mnsq;
0421     f12 = -Easq;
0422     d21 = (Easq * pbx) / Ebsq;
0423     d20 = ((masq - mbsq) * pbx) / (2. * Ebsq) - pmissx + (pbx * (pbx * pmissx + pby * pmissy)) / Ebsq;
0424     e21 = (Easq * pby) / Ebsq;
0425     e20 = ((masq - mbsq) * pby) / (2. * Ebsq) - pmissy + (pby * (pbx * pmissx + pby * pmissy)) / Ebsq;
0426     f22 = -Easq * Easq / Ebsq;
0427     f21 = (-2 * Easq * ((masq - mbsq) / (2. * Eb) + (pbx * pmissx + pby * pmissy) / Eb)) / Eb;
0428     f20 = mnsq + pmissx * pmissx + pmissy * pmissy -
0429           ((masq - mbsq) / (2. * Eb) + (pbx * pmissx + pby * pmissy) / Eb) *
0430               ((masq - mbsq) / (2. * Eb) + (pbx * pmissx + pby * pmissy) / Eb);
0431
0432     //Estimate upper bound of mT2
0433     //when Deltasq > Deltasq_high1, the larger encloses the center of the smaller
0434     double p2x0, p2y0;
0435     double Deltasq_high1;
0436     p2x0 = pmissx - x0;
0437     p2y0 = pmissy - y0;
0438     Deltasq_high1 = 2 * Eb * sqrt(p2x0 * p2x0 + p2y0 * p2y0 + mnsq) - 2 * pbx * p2x0 - 2 * pby * p2y0 + mbsq;
0439
0440     //Another estimate, if both ellipses enclose the origin, Deltasq > mT2
0441
0442     double Deltasq_high2, Deltasq_high21, Deltasq_high22;
0443     Deltasq_high21 =
0444         2 * Eb * sqrt(pmissx * pmissx + pmissy * pmissy + mnsq) - 2 * pbx * pmissx - 2 * pby * pmissy + mbsq;
0445     Deltasq_high22 = 2 * Ea * mn + masq;
0446
0447     if (Deltasq_high21 < Deltasq_high22)
0448       Deltasq_high2 = Deltasq_high22;
0449     else
0450       Deltasq_high2 = Deltasq_high21;
0451
0452     //pick the smaller upper bound
0453     double Deltasq_high;
0454     if (Deltasq_high1 < Deltasq_high2)
0455       Deltasq_high = Deltasq_high1;
0456     else
0457       Deltasq_high = Deltasq_high2;
0458
0459     double Deltasq_low;  //lower bound
0460     Deltasq_low = Deltasq0;
0461
0462     //number of solutions at Deltasq_low should not be larger than zero
0463     if (nsols(Deltasq_low) > 0) {
0464       //cout << "Davismt2::mt2_bisect() ==> nsolutions(Deltasq_low) > 0"<<endl;
0465       mt2_b = (double)sqrt(mnsq + Deltasq0);
0466       return;
0467     }
0468
0469     int nsols_high, nsols_low;
0470
0471     nsols_low = nsols(Deltasq_low);
0472     int foundhigh;
0473
0474     //if nsols_high=nsols_low, we missed the region where the two ellipse overlap
0475     //if nsols_high=4, also need a scan because we may find the wrong tangent point.
0476
0477     nsols_high = nsols(Deltasq_high);
0478
0479     if (nsols_high == nsols_low || nsols_high == 4) {
0480       //foundhigh = scan_high(Deltasq_high);
0481       foundhigh = find_high(Deltasq_high);
0482       if (foundhigh == 0) {
0483         if (verbose > 0)
0484           cout << "Davismt2::mt2_bisect() ==> Deltasq_high not found at event " << nevt << endl;
0485         mt2_b = sqrt(Deltasq_low + mnsq);
0486         return;
0487       }
0488     }
0489
0490     while (sqrt(Deltasq_high + mnsq) - sqrt(Deltasq_low + mnsq) > precision) {
0491       double Deltasq_mid, nsols_mid;
0492       //bisect
0493       Deltasq_mid = (Deltasq_high + Deltasq_low) / 2.;
0494       nsols_mid = nsols(Deltasq_mid);
0495       // if nsols_mid = 4, rescan for Deltasq_high
0496       if (nsols_mid == 4) {
0497         Deltasq_high = Deltasq_mid;
0498         //scan_high(Deltasq_high);
0499         find_high(Deltasq_high);
0500         continue;
0501       }
0502
0503       if (nsols_mid != nsols_low)
0504         Deltasq_high = Deltasq_mid;
0505       if (nsols_mid == nsols_low)
0506         Deltasq_low = Deltasq_mid;
0507     }
0508     mt2_b = (double)sqrt(mnsq + Deltasq_high);
0509     return;
0510   }
0511
0512   int Davismt2::find_high(double& Deltasq_high) {
0513     double x0, y0;
0514     x0 = (c1 * d1 - b1 * e1) / (b1 * b1 - a1 * c1);
0515     y0 = (a1 * e1 - b1 * d1) / (b1 * b1 - a1 * c1);
0516     double Deltasq_low = (mn + ma) * (mn + ma) - mnsq;
0517     do {
0518       double Deltasq_mid = (Deltasq_high + Deltasq_low) / 2.;
0519       int nsols_mid = nsols(Deltasq_mid);
0520       if (nsols_mid == 2) {
0521         Deltasq_high = Deltasq_mid;
0522         return 1;
0523       } else if (nsols_mid == 4) {
0524         Deltasq_high = Deltasq_mid;
0525         continue;
0526       } else if (nsols_mid == 0) {
0527         d1 = -pax * (Deltasq_mid - masq) / (2 * Easq);
0528         e1 = -pay * (Deltasq_mid - masq) / (2 * Easq);
0529         d2 = -pmissx + pbx * (Deltasq_mid - mbsq) / (2 * Ebsq) + pbx * (pbx * pmissx + pby * pmissy) / (Ebsq);
0530         e2 = -pmissy + pby * (Deltasq_mid - mbsq) / (2 * Ebsq) + pby * (pbx * pmissx + pby * pmissy) / (Ebsq);
0531         f2 = pmissx * pmissx + pmissy * pmissy -
0532              ((Deltasq_mid - mbsq) / (2 * Eb) + (pbx * pmissx + pby * pmissy) / Eb) *
0533                  ((Deltasq_mid - mbsq) / (2 * Eb) + (pbx * pmissx + pby * pmissy) / Eb) +
0534              mnsq;
0535         // Does the larger ellipse contain the smaller one?
0536         double dis = a2 * x0 * x0 + 2 * b2 * x0 * y0 + c2 * y0 * y0 + 2 * d2 * x0 + 2 * e2 * y0 + f2;
0537         if (dis < 0)
0538           Deltasq_high = Deltasq_mid;
0539         else
0540           Deltasq_low = Deltasq_mid;
0541       }
0542
0543     } while (Deltasq_high - Deltasq_low > 0.001);
0544     return 0;
0545   }
0546
0547   int Davismt2::scan_high(double& Deltasq_high) {
0548     int foundhigh = 0;
0549     int nsols_high;
0550
0551     //        double Deltasq_low;
0552     double tempmass, maxmass;
0553     tempmass = mn + ma;
0554     maxmass = sqrt(mnsq + Deltasq_high);
0555     if (nevt == 32334)
0556       cout << "Davismt2::scan_high() ==> Deltasq_high = " << Deltasq_high << endl;  // ???
0557     for (double mass = tempmass + SCANSTEP; mass < maxmass; mass += SCANSTEP) {
0558       Deltasq_high = mass * mass - mnsq;
0559       nsols_high = nsols(Deltasq_high);
0560
0561       if (nsols_high > 0) {
0562         // Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
0563         foundhigh = 1;
0564         break;
0565       }
0566     }
0567     return foundhigh;
0568   }
0569
0570   int Davismt2::nsols(double Dsq) {
0571     double delta = (Dsq - masq) / (2 * Easq);
0572
0573     //calculate coefficients for the two quadratic equations
0574     d1 = d11 * delta;
0575     e1 = e11 * delta;
0576     f1 = f12 * delta * delta + f10;
0577     d2 = d21 * delta + d20;
0578     e2 = e21 * delta + e20;
0579     f2 = f22 * delta * delta + f21 * delta + f20;
0580
0581     //obtain the coefficients for the 4th order equation
0582     //devided by Ea^n to make the variable dimensionless
0583     long double A4, A3, A2, A1, A0;
0584
0585     A4 = -4 * a2 * b1 * b2 * c1 + 4 * a1 * b2 * b2 * c1 + a2 * a2 * c1 * c1 + 4 * a2 * b1 * b1 * c2 -
0586          4 * a1 * b1 * b2 * c2 - 2 * a1 * a2 * c1 * c2 + a1 * a1 * c2 * c2;
0587
0588     A3 = (-4 * a2 * b2 * c1 * d1 + 8 * a2 * b1 * c2 * d1 - 4 * a1 * b2 * c2 * d1 - 4 * a2 * b1 * c1 * d2 +
0589           8 * a1 * b2 * c1 * d2 - 4 * a1 * b1 * c2 * d2 - 8 * a2 * b1 * b2 * e1 + 8 * a1 * b2 * b2 * e1 +
0590           4 * a2 * a2 * c1 * e1 - 4 * a1 * a2 * c2 * e1 + 8 * a2 * b1 * b1 * e2 - 8 * a1 * b1 * b2 * e2 -
0591           4 * a1 * a2 * c1 * e2 + 4 * a1 * a1 * c2 * e2) /
0592          Ea;
0593
0594     A2 = (4 * a2 * c2 * d1 * d1 - 4 * a2 * c1 * d1 * d2 - 4 * a1 * c2 * d1 * d2 + 4 * a1 * c1 * d2 * d2 -
0595           8 * a2 * b2 * d1 * e1 - 8 * a2 * b1 * d2 * e1 + 16 * a1 * b2 * d2 * e1 + 4 * a2 * a2 * e1 * e1 +
0596           16 * a2 * b1 * d1 * e2 - 8 * a1 * b2 * d1 * e2 - 8 * a1 * b1 * d2 * e2 - 8 * a1 * a2 * e1 * e2 +
0597           4 * a1 * a1 * e2 * e2 - 4 * a2 * b1 * b2 * f1 + 4 * a1 * b2 * b2 * f1 + 2 * a2 * a2 * c1 * f1 -
0598           2 * a1 * a2 * c2 * f1 + 4 * a2 * b1 * b1 * f2 - 4 * a1 * b1 * b2 * f2 - 2 * a1 * a2 * c1 * f2 +
0599           2 * a1 * a1 * c2 * f2) /
0600          Easq;
0601
0602     A1 = (-8 * a2 * d1 * d2 * e1 + 8 * a1 * d2 * d2 * e1 + 8 * a2 * d1 * d1 * e2 - 8 * a1 * d1 * d2 * e2 -
0603           4 * a2 * b2 * d1 * f1 - 4 * a2 * b1 * d2 * f1 + 8 * a1 * b2 * d2 * f1 + 4 * a2 * a2 * e1 * f1 -
0604           4 * a1 * a2 * e2 * f1 + 8 * a2 * b1 * d1 * f2 - 4 * a1 * b2 * d1 * f2 - 4 * a1 * b1 * d2 * f2 -
0605           4 * a1 * a2 * e1 * f2 + 4 * a1 * a1 * e2 * f2) /
0606          (Easq * Ea);
0607
0608     A0 = (-4 * a2 * d1 * d2 * f1 + 4 * a1 * d2 * d2 * f1 + a2 * a2 * f1 * f1 + 4 * a2 * d1 * d1 * f2 -
0609           4 * a1 * d1 * d2 * f2 - 2 * a1 * a2 * f1 * f2 + a1 * a1 * f2 * f2) /
0610          (Easq * Easq);
0611
0612     long double A3sq;
0613     /*
0614     long  double A0sq, A1sq, A2sq, A3sq, A4sq;
0615     A0sq = A0*A0;
0616     A1sq = A1*A1;
0617     A2sq = A2*A2;
0618     A4sq = A4*A4;
0619         */
0620     A3sq = A3 * A3;
0621
0622     long double B3, B2, B1, B0;
0623     B3 = 4 * A4;
0624     B2 = 3 * A3;
0625     B1 = 2 * A2;
0626     B0 = A1;
0627
0628     long double C2, C1, C0;
0629     C2 = -(A2 / 2 - 3 * A3sq / (16 * A4));
0630     C1 = -(3 * A1 / 4. - A2 * A3 / (8 * A4));
0631     C0 = -A0 + A1 * A3 / (16 * A4);
0632
0633     long double D1, D0;
0634     D1 = -B1 - (B3 * C1 * C1 / C2 - B3 * C0 - B2 * C1) / C2;
0635     D0 = -B0 - B3 * C0 * C1 / (C2 * C2) + B2 * C0 / C2;
0636
0637     long double E0;
0638     E0 = -C0 - C2 * D0 * D0 / (D1 * D1) + C1 * D0 / D1;
0639
0640     long double t1, t2, t3, t4, t5;
0641     //find the coefficients for the leading term in the Sturm sequence
0642     t1 = A4;
0643     t2 = A4;
0644     t3 = C2;
0645     t4 = D1;
0646     t5 = E0;
0647
0648     //The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
0649     int nsol;
0650     nsol = signchange_n(t1, t2, t3, t4, t5) - signchange_p(t1, t2, t3, t4, t5);
0651
0652     //Cannot have negative number of solutions, must be roundoff effect
0653     if (nsol < 0)
0654       nsol = 0;
0655
0656     return nsol;
0657   }
0658
0659   //inline
0660   int Davismt2::signchange_n(long double t1, long double t2, long double t3, long double t4, long double t5) {
0661     int nsc;
0662     nsc = 0;
0663     if (t1 * t2 > 0)
0664       nsc++;
0665     if (t2 * t3 > 0)
0666       nsc++;
0667     if (t3 * t4 > 0)
0668       nsc++;
0669     if (t4 * t5 > 0)
0670       nsc++;
0671     return nsc;
0672   }
0673
0674   //inline
0675   int Davismt2::signchange_p(long double t1, long double t2, long double t3, long double t4, long double t5) {
0676     int nsc;
0677     nsc = 0;
0678     if (t1 * t2 < 0)
0679       nsc++;
0680     if (t2 * t3 < 0)
0681       nsc++;
0682     if (t3 * t4 < 0)
0683       nsc++;
0684     if (t4 * t5 < 0)
0685       nsc++;
0686     return nsc;
0687   }
0688
0689 }  // namespace heppy