file src/mt2_bisect.cpp
[No description available]
Namespaces
| Name |
|---|
| mt2_bisect |
Source code
/***********************************************************************/
/* */
/* 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 <iostream>
#include <math.h>
#include "gambit/ColliderBit/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 <<endl;
mt2_b = (double)sqrt(Deltasq_low+mnsq);
return;
}
}
if(nsols_high == nsols_low)
{
cout << "error: nsols_low=nsols_high=" << nsols_high << endl;
cout << "Deltasq_high=" << Deltasq_high << endl;
cout << "Deltasq_low= "<< Deltasq_low << endl;
mt2_b = sqrt(mnsq + Deltasq_low);
return;
}
double minmass, maxmass;
minmass = sqrt(Deltasq_low+mnsq);
maxmass = sqrt(Deltasq_high+mnsq);
while(maxmass - minmass > precision)
{
double Delta_mid, midmass, nsols_mid;
midmass = (minmass+maxmass)/2.;
Delta_mid = midmass * midmass - mnsq;
nsols_mid = nsols_massless(Delta_mid);
if(nsols_mid != nsols_low) maxmass = midmass;
if(nsols_mid == nsols_low) minmass = midmass;
}
mt2_b = minmass;
return;
}
int mt2::nsols_massless(double Dsq)
{
double delta;
delta = Dsq/(2*Easq);
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;
double a,b;
if (pax > 0) a = Ea/Dsq;
else a = -Ea/Dsq;
if (pax > 0) b = -Dsq/(4*Ea)+mnsq*Ea/Dsq;
else b = Dsq/(4*Ea)-mnsq*Ea/Dsq;
double A4,A3,A2,A1,A0;
A4 = a*a*a2;
A3 = 2*a*b2/Ea;
A2 = (2*a*a2*b+c2+2*a*d2)/(Easq);
A1 = (2*b*b2+2*e2)/(Easq*Ea);
A0 = (a2*b*b+2*b*d2+f2)/(Easq*Easq);
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*A3*A3/(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;
int nsol;
nsol = signchange_n(t1,t2,t3,t4,t5)-signchange_p(t1,t2,t3,t4,t5);
if( nsol < 0 ) nsol=0;
return nsol;
}
void mt2::mt2_bisect()
{
solved = true;
cout.precision(11);
//if masses are very small, use code for massless case.
if(masq < MIN_MASS && mbsq < MIN_MASS)
{
mt2_massless();
return;
}
double Deltasq0;
Deltasq0 = ma*(ma + 2*mn); //The minimum mass square to have two ellipses
// find the coefficients for the two quadratic equations when Deltasq=Deltasq0.
a1 = 1-pax*pax/(Easq);
b1 = -pax*pay/(Easq);
c1 = 1-pay*pay/(Easq);
d1 = -pax*(Deltasq0-masq)/(2*Easq);
e1 = -pay*(Deltasq0-masq)/(2*Easq);
a2 = 1-pbx*pbx/(Ebsq);
b2 = -pbx*pby/(Ebsq);
c2 = 1-pby*pby/(Ebsq);
d2 = -pmissx+pbx*(Deltasq0-mbsq)/(2*Ebsq)+pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
e2 = -pmissy+pby*(Deltasq0-mbsq)/(2*Ebsq)+pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq0-mbsq)/(2*Eb)+
(pbx*pmissx+pby*pmissy)/Eb)*((Deltasq0-mbsq)/(2*Eb)+
(pbx*pmissx+pby*pmissy)/Eb)+mnsq;
// find the center of the smaller ellipse
double x0,y0;
x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
// 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.01)
{
mt2_b = (double) sqrt(mnsq+Deltasq0);
return;
}
/* find the coefficients for the two quadratic equations */
/* coefficients for quadratic terms do not change */
/* coefficients for linear and constant terms are polynomials of */
/* delta=(Deltasq-m7sq)/(2 E7sq) */
d11 = -pax;
e11 = -pay;
f10 = mnsq;
f12 = -Easq;
d21 = (Easq*pbx)/Ebsq;
d20 = ((masq - mbsq)*pbx)/(2.*Ebsq) - pmissx +
(pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
e21 = (Easq*pby)/Ebsq;
e20 = ((masq - mbsq)*pby)/(2.*Ebsq) - pmissy +
(pby*(pbx*pmissx + pby*pmissy))/Ebsq;
f22 = -Easq*Easq/Ebsq;
f21 = (-2*Easq*((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb))/Eb;
f20 = mnsq + pmissx*pmissx + pmissy*pmissy -
((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb)
*((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb);
//Estimate upper bound of mT2
//when Deltasq > Deltasq_high1, the larger encloses the center of the smaller
double p2x0,p2y0;
double Deltasq_high1;
p2x0 = pmissx-x0;
p2y0 = pmissy-y0;
Deltasq_high1 = 2*Eb*sqrt(p2x0*p2x0+p2y0*p2y0+mnsq)-2*pbx*p2x0-2*pby*p2y0+mbsq;
//Another estimate, if both ellipses enclose the origin, Deltasq > mT2
double Deltasq_high2, Deltasq_high21, Deltasq_high22;
Deltasq_high21 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy+mbsq;
Deltasq_high22 = 2*Ea*mn+masq;
if ( Deltasq_high21 < Deltasq_high22 ) Deltasq_high2 = Deltasq_high22;
else Deltasq_high2 = Deltasq_high21;
//pick the smaller upper bound
double Deltasq_high;
if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high1;
else Deltasq_high = Deltasq_high2;
double Deltasq_low; //lower bound
Deltasq_low = Deltasq0;
//number of solutions at Deltasq_low should not be larger than zero
if( nsols(Deltasq_low) > 0 )
{
//cout << "nsolutions(Deltasq_low) > 0"<<endl;
mt2_b = (double) sqrt(mnsq+Deltasq0);
return;
}
int nsols_high, nsols_low;
nsols_low = nsols(Deltasq_low);
int foundhigh;
//if nsols_high=nsols_low, we missed the region where the two ellipse overlap
//if nsols_high=4, also need a scan because we may find the wrong tangent point.
nsols_high = nsols(Deltasq_high);
if(nsols_high == nsols_low || nsols_high == 4)
{
//foundhigh = scan_high(Deltasq_high);
foundhigh = find_high(Deltasq_high);
if(foundhigh == 0)
{
cout << "Deltasq_high not found at event " << nevt << endl;
mt2_b = sqrt( Deltasq_low + mnsq );
return;
}
}
while(sqrt(Deltasq_high+mnsq) - sqrt(Deltasq_low+mnsq) > 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 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)
{
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 B3, B2, B1, B0;
B3 = 4*A4;
B2 = 3*A3;
B1 = 2*A2;
B0 = A1;
long double C2, C1, C0;
C2 = -(A2/2 - 3*A3*A3/(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
Updated on 2023-06-26 at 21:36:56 +0000