mt2 Class Reference

#include <mt2_bisect.hh>

List of all members.

Public Member Functions

 mt2 ()
void mt2_bisect ()
void mt2_massless ()
void set_momenta (double *pa0, double *pb0, double *pmiss0)
void set_mn (double mn)
double get_mt2 ()

Public Attributes

int nevt

Private Member Functions

int nsols (double Dsq)
int nsols_massless (double Dsq)
int signchange_n (long double t1, long double t2, long double t3, long double t4, long double t5)
int signchange_p (long double t1, long double t2, long double t3, long double t4, long double t5)
int scan_high (double &Deltasq_high)
int find_high (double &Deltasq_high)

Private Attributes

bool solved
bool momenta_set
double mt2_b
double pax
double pay
double ma
double Ea
double pmissx
double pmissy
double pbx
double pby
double mb
double Eb
double mn
double mn_unscale
double masq
double Easq
double mbsq
double Ebsq
double pmissxsq
double pmissysq
double mnsq
double a1
double b1
double c1
double a2
double b2
double c2
double d1
double e1
double f1
double d2
double e2
double f2
double d11
double e11
double f12
double f10
double d21
double d20
double e21
double e20
double f22
double f21
double f20
double scale
double precision

Detailed Description

Definition at line 37 of file mt2_bisect.hh.

Constructor & Destructor Documentation

mt2 (  ) 

Definition at line 41 of file mt2_bisect.cc.

00042                 {
00043 
00044 namespace mt2_bisect
00045 {
00046 using namespace std;
00047

Member Function Documentation

int find_high ( double &  Deltasq_high  )  [private]

Definition at line 479 of file mt2_bisect.cc.

00487 {
00488    double x0,y0;
00489    x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
00490    y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
00491    double Deltasq_low = (mn + ma)*(mn + ma) - mnsq;
00492    do 
00493    {
00494       double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.;
00495       int nsols_mid = nsols(Deltasq_mid);
00496       if ( nsols_mid == 2 )
00497       {
00498          Deltasq_high = Deltasq_mid;
00499          return 1;
00500       }
00501       else if (nsols_mid == 4)
00502       {
00503          Deltasq_high = Deltasq_mid;
00504          continue;
00505       }
00506       else if (nsols_mid ==0)
00507       {
00508          d1 = -pax*(Deltasq_mid-masq)/(2*Easq);
00509          e1 = -pay*(Deltasq_mid-masq)/(2*Easq);
00510          d2 = -pmissx + pbx*(Deltasq_mid - mbsq)/(2*Ebsq)
00511               + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
00512          e2 = -pmissy + pby*(Deltasq_mid - mbsq)/(2*Ebsq)
00513               + pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
00514          f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq)/(2*Eb)+
00515               (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq)/(2*Eb)+
00516               (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
00517 // Does the larger ellipse contain the smaller one? 
00518          double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;

double get_mt2 (  ) 

Definition at line 49 of file mt2_bisect.cc.

References mt2::momenta_set, mt2::mt2_b, mt2::scale, and mt2::solved.

Referenced by Rivet::mT2::mT2().

00049 {
00050    solved = false;
00051    momenta_set = false;
00052    mt2_b  = 0.;
00053    scale = 1.;
00054 }

void mt2_bisect (  ) 

Definition at line 329 of file mt2_bisect.cc.

Referenced by mt2::set_momenta().

00337 {   
00338   solved = true;
00339 
00340 //if masses are very small, use code for massless case.  
00341    if(masq < MIN_MASS && mbsq < MIN_MASS) 
00342    { 
00343       mt2_massless();
00344       return;
00345    }
00346  
00347 
00348    double Deltasq0;     
00349    Deltasq0 = ma*(ma + 2*mn); //The minimum mass square to have two ellipses 
00350  
00351 // find the coefficients for the two quadratic equations when Deltasq=Deltasq0.
00352   
00353    a1 = 1-pax*pax/(Easq);
00354    b1 = -pax*pay/(Easq);
00355    c1 = 1-pay*pay/(Easq);
00356    d1 = -pax*(Deltasq0-masq)/(2*Easq);
00357    e1 = -pay*(Deltasq0-masq)/(2*Easq);
00358    a2 = 1-pbx*pbx/(Ebsq);
00359    b2 = -pbx*pby/(Ebsq);
00360    c2 = 1-pby*pby/(Ebsq);
00361    d2 = -pmissx+pbx*(Deltasq0-mbsq)/(2*Ebsq)+pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
00362    e2 = -pmissy+pby*(Deltasq0-mbsq)/(2*Ebsq)+pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
00363    f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq0-mbsq)/(2*Eb)+
00364         (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq0-mbsq)/(2*Eb)+
00365         (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
00366    
00367 // find the center of the smaller ellipse 
00368    double x0,y0;
00369    x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
00370    y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
00371 
00372    
00373 // Does the larger ellipse contain the smaller one? 
00374    double dis=a2*x0*x0+2*b2*x0*y0+c2*y0*y0+2*d2*x0+2*e2*y0+f2;
00375 
00376    if(dis<=0.01)
00377    { 
00378       mt2_b  = (double) sqrt(mnsq+Deltasq0);
00379       return;
00380    }
00381    
00382 
00383 /* find the coefficients for the two quadratic equations           */
00384 /* coefficients for quadratic terms do not change                  */
00385 /* coefficients for linear and constant terms are polynomials of   */
00386 /*       delta=(Deltasq-m7sq)/(2 E7sq)                             */  
00387    d11 = -pax;
00388    e11 = -pay;
00389    f10 = mnsq;
00390    f12 = -Easq;
00391    d21 = (Easq*pbx)/Ebsq;
00392    d20 = ((masq - mbsq)*pbx)/(2.*Ebsq) - pmissx +
00393          (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
00394    e21 = (Easq*pby)/Ebsq;
00395    e20 = ((masq - mbsq)*pby)/(2.*Ebsq) - pmissy +
00396          (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
00397    f22 = -Easq*Easq/Ebsq;
00398    f21 = (-2*Easq*((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb))/Eb;
00399    f20 = mnsq + pmissx*pmissx + pmissy*pmissy - 
00400          ((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb)
00401          *((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb);
00402 
00403 //Estimate upper bound of mT2
00404 //when Deltasq > Deltasq_high1, the larger encloses the center of the smaller 
00405    double p2x0,p2y0;
00406    double Deltasq_high1;
00407    p2x0 = pmissx-x0;
00408    p2y0 = pmissy-y0;
00409    Deltasq_high1 = 2*Eb*sqrt(p2x0*p2x0+p2y0*p2y0+mnsq)-2*pbx*p2x0-2*pby*p2y0+mbsq;
00410    
00411 //Another estimate, if both ellipses enclose the origin, Deltasq > mT2
00412 
00413    double Deltasq_high2, Deltasq_high21, Deltasq_high22;
00414    Deltasq_high21 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy+mbsq;
00415    Deltasq_high22 = 2*Ea*mn+masq;
00416   
00417    if ( Deltasq_high21 < Deltasq_high22 ) Deltasq_high2 = Deltasq_high22;
00418    else Deltasq_high2 = Deltasq_high21;
00419 
00420 //pick the smaller upper bound   
00421    double Deltasq_high;
00422    if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high1;
00423    else Deltasq_high = Deltasq_high2;
00424    
00425   
00426    double Deltasq_low; //lower bound
00427    Deltasq_low = Deltasq0;
00428 
00429 //number of solutions at Deltasq_low should not be larger than zero
00430    if( nsols(Deltasq_low) > 0 )
00431    {
00432      //cout << "nsolutions(Deltasq_low) > 0"<<endl;
00433      mt2_b = (double) sqrt(mnsq+Deltasq0);
00434      return;
00435    }
00436   
00437    int nsols_high, nsols_low;
00438 
00439    nsols_low  = nsols(Deltasq_low);
00440    int foundhigh;
00441   
00442 
00443 //if nsols_high=nsols_low, we missed the region where the two ellipse overlap 
00444 //if nsols_high=4, also need a scan because we may find the wrong tangent point.
00445 
00446    nsols_high = nsols(Deltasq_high);
00447   
00448    if(nsols_high == nsols_low || nsols_high == 4)
00449    {
00450       //foundhigh = scan_high(Deltasq_high);
00451       foundhigh = find_high(Deltasq_high);
00452       if(foundhigh == 0) 
00453       {
00454      Log::getLog("Rivet.Tools.mt2") 
00455        << Log::WARN  
00456        << "Deltasq_high not found at event " << nevt << '\n';
00457          mt2_b = sqrt( Deltasq_low + mnsq );
00458          return;
00459       }
00460       
00461    }
00462 
00463    while(sqrt(Deltasq_high+mnsq) - sqrt(Deltasq_low+mnsq) > precision)
00464    {
00465       double Deltasq_mid,nsols_mid;
00466       //bisect
00467       Deltasq_mid = (Deltasq_high+Deltasq_low)/2.;
00468       nsols_mid = nsols(Deltasq_mid);
00469       // if nsols_mid = 4, rescan for Deltasq_high
00470       if ( nsols_mid == 4 ) 
00471       {
00472          Deltasq_high = Deltasq_mid;
00473          //scan_high(Deltasq_high);
00474          find_high(Deltasq_high);
00475          continue;
00476       } 
00477

void mt2_massless (  ) 

Definition at line 141 of file mt2_bisect.cc.

References mt2::a2, mt2::b2, mt2::c2, mt2::d20, mt2::d21, mt2::e20, mt2::e21, mt2::Ea, mt2::Easq, mt2::Eb, mt2::Ebsq, Log::ERROR, mt2::f20, mt2::f21, mt2::f22, Log::getLog(), Rivet::mass(), mt2::mnsq, mt2::mt2_b, mt2::nevt, mt2::nsols_massless(), mt2::pax, mt2::pbx, mt2::pmissxsq, mt2::pmissysq, mt2::precision, Rivet::s, Rivet::mt2_bisect::SCANSTEP, Rivet::theta(), and Log::WARN.

00149 {
00150    
00151 //rotate so that pay = 0 
00152    double theta,s,c;
00153    theta = atan(pay/pax);
00154    s = sin(theta);
00155    c = cos(theta);
00156 
00157    double pxtemp,pytemp;
00158    Easq   = pax*pax+pay*pay;
00159    Ebsq   = pbx*pbx+pby*pby;
00160    Ea     = sqrt(Easq);
00161    Eb     = sqrt(Ebsq);
00162   
00163    pxtemp = pax*c+pay*s;
00164    pax    = pxtemp;
00165    pay    = 0;
00166    pxtemp = pbx*c+pby*s;
00167    pytemp = -s*pbx+c*pby;
00168    pbx    = pxtemp;
00169    pby    = pytemp;
00170    pxtemp = pmissx*c+pmissy*s;
00171    pytemp = -s*pmissx+c*pmissy;
00172    pmissx = pxtemp;
00173    pmissy = pytemp;
00174 
00175    a2  = 1-pbx*pbx/(Ebsq);
00176    b2  = -pbx*pby/(Ebsq);
00177    c2  = 1-pby*pby/(Ebsq);
00178 
00179    d21 = (Easq*pbx)/Ebsq;
00180    d20 = - pmissx +  (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
00181    e21 = (Easq*pby)/Ebsq;
00182    e20 = - pmissy +  (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
00183    f22 = -(Easq*Easq/Ebsq);
00184    f21 = -2*Easq*(pbx*pmissx + pby*pmissy)/Ebsq;
00185    f20 = mnsq + pmissxsq + pmissysq - (pbx*pmissx + pby*pmissy)*(pbx*pmissx + pby*pmissy)/Ebsq;
00186 
00187    double Deltasq0    = 0; 
00188    double Deltasq_low, Deltasq_high;
00189    int    nsols_high, nsols_low;
00190 
00191    Deltasq_low = Deltasq0 + precision;
00192    nsols_low = nsols_massless(Deltasq_low);
00193    
00194    if(nsols_low > 1) 
00195    { 
00196       mt2_b = (double) sqrt(Deltasq0+mnsq);
00197       return;
00198    }
00199 
00200 /*   
00201    if( nsols_massless(Deltasq_high) > 0 )
00202    {
00203       mt2_b = (double) sqrt(mnsq+Deltasq0);
00204       return;
00205    }*/
00206 
00207 //look for when both parablos contain origin  
00208    double Deltasq_high1, Deltasq_high2;
00209    Deltasq_high1 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy;
00210    Deltasq_high2 = 2*Ea*mn;
00211  
00212    if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high2;
00213      else Deltasq_high = Deltasq_high1;
00214 
00215    nsols_high=nsols_massless(Deltasq_high);
00216   
00217    int foundhigh;
00218    if (nsols_high == nsols_low)
00219    {
00220       
00221       
00222       foundhigh=0;
00223       
00224       double minmass, maxmass;
00225       minmass  = mn ;
00226       maxmass  = sqrt(mnsq + Deltasq_high);
00227       for(double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
00228       {
00229      Deltasq_high = mass*mass - mnsq;
00230       
00231          nsols_high = nsols_massless(Deltasq_high);
00232          if(nsols_high>0)
00233          {
00234             foundhigh=1;
00235             Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
00236             break;
00237          }
00238       }
00239       if(foundhigh==0) 
00240       {
00241        
00242     Log::getLog("Rivet.Tools.mt2") 
00243       << Log::WARN  
00244       << "Deltasq_high not found at event " << nevt << '\n';
00245        
00246          mt2_b = (double)sqrt(Deltasq_low+mnsq);
00247          return;
00248       }
00249    }
00250 
00251    if(nsols_high == nsols_low)
00252    { 
00253      Log::getLog("Rivet.Tools.mt2") 
00254        << Log::ERROR  
00255        << "error: nsols_low=nsols_high=" << nsols_high << '\n'
00256        << "Deltasq_high=" << Deltasq_high << '\n'
00257        << "Deltasq_low= "<< Deltasq_low << '\n';
00258     
00259       mt2_b = sqrt(mnsq + Deltasq_low);
00260       return;
00261    }
00262    double minmass, maxmass;
00263    minmass = sqrt(Deltasq_low+mnsq);
00264    maxmass = sqrt(Deltasq_high+mnsq);
00265    while(maxmass - minmass > precision)
00266    {
00267       double Delta_mid, midmass, nsols_mid;
00268       midmass   = (minmass+maxmass)/2.;
00269       Delta_mid = midmass * midmass - mnsq;

int nsols ( double  Dsq  )  [private]

Definition at line 545 of file mt2_bisect.cc.

00545       {
00546          foundhigh   = 1;
00547          break;
00548       }
00549     }
00550     return foundhigh;
00551 }
00552 int mt2::nsols(  double Dsq)
00553 {
00554    double delta = (Dsq-masq)/(2*Easq);
00555   
00556 //calculate coefficients for the two quadratic equations
00557    d1 = d11*delta;
00558    e1 = e11*delta;
00559    f1 = f12*delta*delta+f10;
00560    d2 = d21*delta+d20;
00561    e2 = e21*delta+e20;
00562    f2 = f22*delta*delta+f21*delta+f20;
00563 
00564 //obtain the coefficients for the 4th order equation 
00565 //devided by Ea^n to make the variable dimensionless
00566    long double A4, A3, A2, A1, A0;
00567 
00568    A4 = 
00569    -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 + 
00570    4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 + 
00571    a1*a1*c2*c2;  
00572 
00573    A3 =
00574      (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 + 
00575    8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 + 
00576    4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 - 
00577      4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea;
00578 
00579    
00580    A2 =
00581      (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 - 
00582    8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 + 
00583    4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 - 
00584    8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 + 
00585    4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 + 
00586      4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq;
00587   
00588    A1 =
00589      (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 - 
00590    4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 - 
00591    4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 - 
00592      4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea);
00593   
00594    A0 =
00595      (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 + 
00596    4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 + 
00597      a1*a1*f2*f2)/(Easq*Easq);
00598    
00599    long double B3, B2, B1, B0;
00600    B3 = 4*A4;
00601    B2 = 3*A3;
00602    B1 = 2*A2;
00603    B0 = A1;
00604    
00605    long double C2, C1, C0;
00606    C2 = -(A2/2 - 3*A3*A3/(16*A4));
00607    C1 = -(3*A1/4. -A2*A3/(8*A4));
00608    C0 = -A0 + A1*A3/(16*A4);
00609    
00610    long double D1, D0;
00611    D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
00612    D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
00613    
00614    long double E0;
00615    E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
00616    
00617    long  double t1,t2,t3,t4,t5;
00618 //find the coefficients for the leading term in the Sturm sequence  
00619    t1 = A4;
00620    t2 = A4;
00621    t3 = C2;
00622    t4 = D1;
00623    t5 = E0;
00624  
00625 
00626 //The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
00627    int nsol;
00628    nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);

int nsols_massless ( double  Dsq  )  [private]

Definition at line 271 of file mt2_bisect.cc.

Referenced by mt2::mt2_massless().

00279 {
00280   double delta;
00281   delta = Dsq/(2*Easq);
00282   d1    = d11*delta;
00283   e1    = e11*delta;
00284   f1    = f12*delta*delta+f10;
00285   d2    = d21*delta+d20;
00286   e2    = e21*delta+e20;
00287   f2    = f22*delta*delta+f21*delta+f20;
00288   
00289   double a,b;
00290   if (pax > 0) a = Ea/Dsq;
00291   else         a = -Ea/Dsq;
00292   if (pax > 0) b = -Dsq/(4*Ea)+mnsq*Ea/Dsq;
00293   else         b = Dsq/(4*Ea)-mnsq*Ea/Dsq;
00294   
00295   double A4,A3,A2,A1,A0;
00296 
00297   A4 = a*a*a2;
00298   A3 = 2*a*b2/Ea;
00299   A2 = (2*a*a2*b+c2+2*a*d2)/(Easq);
00300   A1 = (2*b*b2+2*e2)/(Easq*Ea);
00301   A0 = (a2*b*b+2*b*d2+f2)/(Easq*Easq);
00302   
00303   long double B3, B2, B1, B0;
00304   B3 = 4*A4;
00305   B2 = 3*A3;
00306   B1 = 2*A2;
00307   B0 = A1;
00308   long double C2, C1, C0;
00309   C2 = -(A2/2 - 3*A3*A3/(16*A4));
00310   C1 = -(3*A1/4. -A2*A3/(8*A4));
00311   C0 = -A0 + A1*A3/(16*A4);
00312   long double  D1, D0;
00313   D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
00314   D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
00315  
00316   long double E0;
00317   E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
00318  
00319   long  double t1,t2,t3,t4,t5;
00320    
00321 //find the coefficients for the leading term in the Sturm sequence  
00322    t1 = A4;
00323    t2 = A4;
00324    t3 = C2;
00325    t4 = D1;
00326    t5 = E0;
00327

int scan_high ( double &  Deltasq_high  )  [private]

Definition at line 519 of file mt2_bisect.cc.

00527 {
00528    int foundhigh = 0 ;
00529    int nsols_high;
00530 
00531    
00532    double tempmass, maxmass;
00533    tempmass = mn + ma;
00534    maxmass  = sqrt(mnsq + Deltasq_high);
00535    // if (nevt == 32334) {
00536      
00537    //   cout << "Deltasq_high = " << Deltasq_high << endl;
00538    // }
00539    for(double mass = tempmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
00540    {
00541       Deltasq_high = mass*mass - mnsq;
00542       nsols_high   = nsols(Deltasq_high);
00543       
00544       if( nsols_high > 0)

void set_mn ( double  mn  ) 

Definition at line 124 of file mt2_bisect.cc.

Referenced by Rivet::mT2::mT2().

void set_momenta ( double *  pa0,
double *  pb0,
double *  pmiss0 
)

Definition at line 56 of file mt2_bisect.cc.

References mt2::momenta_set, mt2::mt2_b, mt2::mt2_bisect(), mt2::scale, and mt2::solved.

Referenced by Rivet::mT2::mT2().

00057 {
00058    assert(momenta_set);
00059    if (!solved) mt2_bisect();
00060    return mt2_b*scale;
00061 }
00062 
00063 void mt2::set_momenta(double* pa0, double* pb0, double* pmiss0)
00064 {
00065    solved = false;     //reset solved tag when momenta are changed.
00066    momenta_set = true;
00067 
00068    ma = fabs(pa0[0]);  // mass cannot be negative
00069 
00070    if (ma < ZERO_MASS) ma = ZERO_MASS;
00071 
00072    pax  = pa0[1]; 
00073    pay  = pa0[2];
00074    masq = ma*ma;
00075    Easq = masq+pax*pax+pay*pay;
00076    Ea   = sqrt(Easq);
00077    
00078    mb = fabs(pb0[0]);
00079 
00080    if (mb < ZERO_MASS) mb = ZERO_MASS;
00081 
00082    pbx  = pb0[1]; 
00083    pby  = pb0[2];
00084    mbsq = mb*mb;
00085    Ebsq = mbsq+pbx*pbx+pby*pby;
00086    Eb   = sqrt(Ebsq);
00087      
00088    pmissx   = pmiss0[1]; pmissy = pmiss0[2];
00089    pmissxsq = pmissx*pmissx;
00090    pmissysq = pmissy*pmissy;
00091 
00092 // set ma>= mb
00093    if(masq < mbsq)
00094    {
00095       double temp;
00096       temp = pax;  pax  = pbx;  pbx  = temp;
00097       temp = pay;  pay  = pby;  pby  = temp;
00098       temp = Ea;   Ea   = Eb;   Eb   = temp;
00099       temp = Easq; Easq = Ebsq; Ebsq = temp;
00100       temp = masq; masq = mbsq; mbsq = temp;
00101       temp = ma;   ma   = mb;   mb   = temp;   
00102    }
00103 //normalize max{Ea, Eb} to 100
00104    if (Ea > Eb) scale = Ea/100.;
00105    else scale = Eb/100.;
00106    
00107    if (sqrt(pmissxsq+pmissysq)/100 > scale) scale = sqrt(pmissxsq+pmissysq)/100;
00108    //scale = 1;
00109    double scalesq = scale * scale;
00110    ma  = ma/scale;
00111    mb  = mb/scale;
00112    masq = masq/scalesq;
00113    mbsq = mbsq/scalesq;
00114    pax = pax/scale; pay = pay/scale;
00115    pbx = pbx/scale; pby = pby/scale;
00116    Ea  = Ea/scale;  Eb = Eb/scale;
00117    
00118    Easq = Easq/scalesq;
00119    Ebsq = Ebsq/scalesq;
00120    pmissx = pmissx/scale;
00121    pmissy = pmissy/scale;
00122    pmissxsq = pmissxsq/scalesq;

int signchange_n ( long double  t1,
long double  t2,
long double  t3,
long double  t4,
long double  t5 
) [inline, private]

Definition at line 630 of file mt2_bisect.cc.

00638 {
00639    int nsc;

int signchange_p ( long double  t1,
long double  t2,
long double  t3,
long double  t4,
long double  t5 
) [inline, private]

Definition at line 640 of file mt2_bisect.cc.

00648 {
00649    int nsc;

Member Data Documentation

double a1 [private]

Definition at line 74 of file mt2_bisect.hh.

double a2 [private]

Definition at line 74 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double b1 [private]

Definition at line 74 of file mt2_bisect.hh.

double b2 [private]

Definition at line 74 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double c1 [private]

Definition at line 74 of file mt2_bisect.hh.

double c2 [private]

Definition at line 74 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double d1 [private]

Definition at line 74 of file mt2_bisect.hh.

double d11 [private]

Definition at line 75 of file mt2_bisect.hh.

double d2 [private]

Definition at line 74 of file mt2_bisect.hh.

double d20 [private]

Definition at line 75 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double d21 [private]

Definition at line 75 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double e1 [private]

Definition at line 74 of file mt2_bisect.hh.

double e11 [private]

Definition at line 75 of file mt2_bisect.hh.

double e2 [private]

Definition at line 74 of file mt2_bisect.hh.

double e20 [private]

Definition at line 75 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double e21 [private]

Definition at line 75 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double Ea [private]

Definition at line 62 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double Easq [private]

Definition at line 68 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double Eb [private]

Definition at line 64 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double Ebsq [private]

Definition at line 69 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double f1 [private]

Definition at line 74 of file mt2_bisect.hh.

double f10 [private]

Definition at line 75 of file mt2_bisect.hh.

double f12 [private]

Definition at line 75 of file mt2_bisect.hh.

double f2 [private]

Definition at line 74 of file mt2_bisect.hh.

double f20 [private]

Definition at line 75 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double f21 [private]

Definition at line 75 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double f22 [private]

Definition at line 75 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double ma [private]

Definition at line 62 of file mt2_bisect.hh.

double masq [private]

Definition at line 68 of file mt2_bisect.hh.

double mb [private]

Definition at line 64 of file mt2_bisect.hh.

double mbsq [private]

Definition at line 69 of file mt2_bisect.hh.

double mn [private]

Definition at line 65 of file mt2_bisect.hh.

double mn_unscale [private]

Definition at line 65 of file mt2_bisect.hh.

double mnsq [private]

Definition at line 71 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

bool momenta_set [private]

Definition at line 52 of file mt2_bisect.hh.

Referenced by mt2::get_mt2(), and mt2::set_momenta().

double mt2_b [private]

Definition at line 53 of file mt2_bisect.hh.

Referenced by mt2::get_mt2(), mt2::mt2_massless(), and mt2::set_momenta().

int nevt

Definition at line 48 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double pax [private]

Definition at line 62 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double pay [private]

Definition at line 62 of file mt2_bisect.hh.

double pbx [private]

Definition at line 64 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double pby [private]

Definition at line 64 of file mt2_bisect.hh.

double pmissx [private]

Definition at line 63 of file mt2_bisect.hh.

double pmissxsq [private]

Definition at line 70 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double pmissy [private]

Definition at line 63 of file mt2_bisect.hh.

double pmissysq [private]

Definition at line 70 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double precision [private]

Definition at line 78 of file mt2_bisect.hh.

Referenced by mt2::mt2_massless().

double scale [private]

Definition at line 77 of file mt2_bisect.hh.

Referenced by mt2::get_mt2(), and mt2::set_momenta().

bool solved [private]

Definition at line 51 of file mt2_bisect.hh.

Referenced by mt2::get_mt2(), and mt2::set_momenta().

The documentation for this class was generated from the following files: