/***********************************************************************/ /* */ /* Finding mt2 by Bisection */ /* */ /* Authors: Hsin-Chia Cheng, Zhenyu Han */ /* Dec 11, 2008, v1.01a */ /* */ /***********************************************************************/ /******************************************************************************* Usage: 1. Define an object of type "mt2": mt2_bisect::mt2 mt2_event; 2. Set momenta and the mass of the invisible particle, mn: mt2_event.set_momenta( pa, pb, pmiss ); mt2_event.set_mass( mn ); where array pa[0..2], pb[0..2], pmiss[0..2] contains (mass,px,py) for the visible particles and the missing momentum. pmiss[0] is not used. All quantities are given in double. 3. Use mt2::get_mt2() to obtain the value of mt2: double mt2_value = mt2_event.get_mt2(); *******************************************************************************/ #include #include #include "mt2_bisect.h" using namespace std; namespace mt2_bisect { mt2::mt2() { solved = false; momenta_set = false; mt2_b = 0.; scale = 1.; } double mt2::get_mt2() { if (!momenta_set) { cout <<" Please set momenta first!" << endl; return 0; } if (!solved) mt2_bisect(); return mt2_b*scale; } void mt2::set_momenta(double* pa0, double* pb0, double* pmiss0) { solved = false; //reset solved tag when momenta are changed. momenta_set = true; ma = fabs(pa0[0]); // mass cannot be negative if (ma < ZERO_MASS) ma = ZERO_MASS; pax = pa0[1]; pay = pa0[2]; masq = ma*ma; Easq = masq+pax*pax+pay*pay; Ea = sqrt(Easq); mb = fabs(pb0[0]); if (mb < ZERO_MASS) mb = ZERO_MASS; pbx = pb0[1]; pby = pb0[2]; mbsq = mb*mb; Ebsq = mbsq+pbx*pbx+pby*pby; Eb = sqrt(Ebsq); pmissx = pmiss0[1]; pmissy = pmiss0[2]; pmissxsq = pmissx*pmissx; pmissysq = pmissy*pmissy; // set ma>= mb if(masq < mbsq) { double temp; temp = pax; pax = pbx; pbx = temp; temp = pay; pay = pby; pby = temp; temp = Ea; Ea = Eb; Eb = temp; temp = Easq; Easq = Ebsq; Ebsq = temp; temp = masq; masq = mbsq; mbsq = temp; temp = ma; ma = mb; mb = temp; } //normalize max{Ea, Eb} to 100 if (Ea > Eb) scale = Ea/100.; else scale = Eb/100.; if (sqrt(pmissxsq+pmissysq)/100 > scale) scale = sqrt(pmissxsq+pmissysq)/100; //scale = 1; double scalesq = scale * scale; ma = ma/scale; mb = mb/scale; masq = masq/scalesq; mbsq = mbsq/scalesq; pax = pax/scale; pay = pay/scale; pbx = pbx/scale; pby = pby/scale; Ea = Ea/scale; Eb = Eb/scale; Easq = Easq/scalesq; Ebsq = Ebsq/scalesq; pmissx = pmissx/scale; pmissy = pmissy/scale; pmissxsq = pmissxsq/scalesq; pmissysq = pmissysq/scalesq; mn = mn_unscale/scale; mnsq = mn*mn; if (ABSOLUTE_PRECISION > 100.*RELATIVE_PRECISION) precision = ABSOLUTE_PRECISION; else precision = 100.*RELATIVE_PRECISION; } void mt2::set_mn(double mn0) { solved = false; //reset solved tag when mn is changed. mn_unscale = fabs(mn0); //mass cannot be negative mn = mn_unscale/scale; mnsq = mn*mn; } void mt2::print() { cout << " pax = " << pax*scale << "; pay = " << pay*scale << "; ma = " << ma*scale <<";"<< endl; cout << " pbx = " << pbx*scale << "; pby = " << pby*scale << "; mb = " << mb*scale <<";"<< endl; cout << " pmissx = " << pmissx*scale << "; pmissy = " << pmissy*scale <<";"<< endl; cout << " mn = " << mn_unscale<<";" << endl; } //special case, the visible particle is massless void mt2::mt2_massless() { //rotate so that pay = 0 double theta,s,c; theta = atan(pay/pax); s = sin(theta); c = cos(theta); double pxtemp,pytemp; Easq = pax*pax+pay*pay; Ebsq = pbx*pbx+pby*pby; Ea = sqrt(Easq); Eb = sqrt(Ebsq); pxtemp = pax*c+pay*s; pax = pxtemp; pay = 0; pxtemp = pbx*c+pby*s; pytemp = -s*pbx+c*pby; pbx = pxtemp; pby = pytemp; pxtemp = pmissx*c+pmissy*s; pytemp = -s*pmissx+c*pmissy; pmissx = pxtemp; pmissy = pytemp; a2 = 1-pbx*pbx/(Ebsq); b2 = -pbx*pby/(Ebsq); c2 = 1-pby*pby/(Ebsq); d21 = (Easq*pbx)/Ebsq; d20 = - pmissx + (pbx*(pbx*pmissx + pby*pmissy))/Ebsq; e21 = (Easq*pby)/Ebsq; e20 = - pmissy + (pby*(pbx*pmissx + pby*pmissy))/Ebsq; f22 = -(Easq*Easq/Ebsq); f21 = -2*Easq*(pbx*pmissx + pby*pmissy)/Ebsq; f20 = mnsq + pmissxsq + pmissysq - (pbx*pmissx + pby*pmissy)*(pbx*pmissx + pby*pmissy)/Ebsq; double Deltasq0 = 0; double Deltasq_low, Deltasq_high; int nsols_high, nsols_low; Deltasq_low = Deltasq0 + precision; nsols_low = nsols_massless(Deltasq_low); if(nsols_low > 1) { mt2_b = (double) sqrt(Deltasq0+mnsq); return; } /* if( nsols_massless(Deltasq_high) > 0 ) { mt2_b = (double) sqrt(mnsq+Deltasq0); return; }*/ //look for when both parablos contain origin double Deltasq_high1, Deltasq_high2; Deltasq_high1 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy; Deltasq_high2 = 2*Ea*mn; if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high2; else Deltasq_high = Deltasq_high1; nsols_high=nsols_massless(Deltasq_high); int foundhigh; if (nsols_high == nsols_low) { foundhigh=0; double minmass, maxmass; minmass = mn ; maxmass = sqrt(mnsq + Deltasq_high); for(double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP) { Deltasq_high = mass*mass - mnsq; nsols_high = nsols_massless(Deltasq_high); if(nsols_high>0) { foundhigh=1; Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq; break; } } if(foundhigh==0) { cout<<"Deltasq_high not found at event " << nevt < 0"< precision) { double Deltasq_mid,nsols_mid; //bisect Deltasq_mid = (Deltasq_high+Deltasq_low)/2.; nsols_mid = nsols(Deltasq_mid); // if nsols_mid = 4, rescan for Deltasq_high if ( nsols_mid == 4 ) { Deltasq_high = Deltasq_mid; //scan_high(Deltasq_high); find_high(Deltasq_high); continue; } if(nsols_mid != nsols_low) Deltasq_high = Deltasq_mid; if(nsols_mid == nsols_low) Deltasq_low = Deltasq_mid; } mt2_b = (double) sqrt( mnsq + Deltasq_high); return; } int mt2::find_high(double & Deltasq_high) { double x0,y0; x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1); y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1); double Deltasq_low = (mn + ma)*(mn + ma) - mnsq; do { double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.; int nsols_mid = nsols(Deltasq_mid); if ( nsols_mid == 2 ) { Deltasq_high = Deltasq_mid; return 1; } else if (nsols_mid == 4) { Deltasq_high = Deltasq_mid; continue; } else if (nsols_mid ==0) { d1 = -pax*(Deltasq_mid-masq)/(2*Easq); e1 = -pay*(Deltasq_mid-masq)/(2*Easq); d2 = -pmissx + pbx*(Deltasq_mid - mbsq)/(2*Ebsq) + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq); e2 = -pmissy + pby*(Deltasq_mid - mbsq)/(2*Ebsq) + pby*(pbx*pmissx+pby*pmissy)/(Ebsq); f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq)/(2*Eb)+ (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq)/(2*Eb)+ (pbx*pmissx+pby*pmissy)/Eb)+mnsq; // Does the larger ellipse contain the smaller one? double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2; if (dis < 0) Deltasq_high = Deltasq_mid; else Deltasq_low = Deltasq_mid; } } while ( Deltasq_high - Deltasq_low > 0.001); return 0; } int mt2::scan_high(double & Deltasq_high) { int foundhigh = 0 ; int nsols_high; double Deltasq_low; double tempmass, maxmass; tempmass = mn + ma; maxmass = sqrt(mnsq + Deltasq_high); if (nevt == 32334) cout << "Deltasq_high = " << Deltasq_high << endl; for(double mass = tempmass + SCANSTEP; mass < maxmass; mass += SCANSTEP) { Deltasq_high = mass*mass - mnsq; nsols_high = nsols(Deltasq_high); if( nsols_high > 0) { Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq; foundhigh = 1; break; } } return foundhigh; } int mt2::nsols( double Dsq) { double delta = (Dsq-masq)/(2*Easq); //calculate coefficients for the two quadratic equations d1 = d11*delta; e1 = e11*delta; f1 = f12*delta*delta+f10; d2 = d21*delta+d20; e2 = e21*delta+e20; f2 = f22*delta*delta+f21*delta+f20; //obtain the coefficients for the 4th order equation //devided by Ea^n to make the variable dimensionless long double A4, A3, A2, A1, A0; A4 = -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 + 4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 + a1*a1*c2*c2; A3 = (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 + 8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 + 4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 - 4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea; A2 = (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 - 8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 + 4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 - 8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 + 4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 + 4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq; A1 = (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 - 4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 - 4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 - 4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea); A0 = (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 + 4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 + a1*a1*f2*f2)/(Easq*Easq); long double A0sq, A1sq, A2sq, A3sq, A4sq; A0sq = A0*A0; A1sq = A1*A1; A2sq = A2*A2; A3sq = A3*A3; A4sq = A4*A4; long double B3, B2, B1, B0; B3 = 4*A4; B2 = 3*A3; B1 = 2*A2; B0 = A1; long double C2, C1, C0; C2 = -(A2/2 - 3*A3sq/(16*A4)); C1 = -(3*A1/4. -A2*A3/(8*A4)); C0 = -A0 + A1*A3/(16*A4); long double D1, D0; D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2; D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2; long double E0; E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1; long double t1,t2,t3,t4,t5; //find the coefficients for the leading term in the Sturm sequence t1 = A4; t2 = A4; t3 = C2; t4 = D1; t5 = E0; //The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf int nsol; nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5); //Cannot have negative number of solutions, must be roundoff effect if (nsol < 0) nsol = 0; return nsol; } inline int mt2::signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5) { int nsc; nsc=0; if(t1*t2>0) nsc++; if(t2*t3>0) nsc++; if(t3*t4>0) nsc++; if(t4*t5>0) nsc++; return nsc; } inline int mt2::signchange_p( long double t1, long double t2, long double t3, long double t4, long double t5) { int nsc; nsc=0; if(t1*t2<0) nsc++; if(t2*t3<0) nsc++; if(t3*t4<0) nsc++; if(t4*t5<0) nsc++; return nsc; } }//end namespace mt2_bisect