mt2_bisect::mt2 Class Reference

#include <mt2_bisect.h>

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 ()
void	print ()

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 17 of file mt2_bisect.h.

Constructor & Destructor Documentation

mt2()

mt2_bisect::mt2::mt2

(

)

Definition at line 42 of file mt2_bisect.cpp.

{
   solved = false;
   momenta_set = false;
   mt2_b  = 0.;
   scale = 1.;
}

Member Function Documentation

find_high()

int mt2_bisect::mt2::find_high ( double &  Deltasq_high )

private

Definition at line 483 of file mt2_bisect.cpp.

{
   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;
}  

get_mt2()

double mt2_bisect::mt2::get_mt2

(

)

Definition at line 50 of file mt2_bisect.cpp.

Referenced by Gambit::ColliderBit::calcMT(), Gambit::ColliderBit::calcMT_1l(), Gambit::ColliderBit::Phi_mpi_pi(), Gambit::ColliderBit::sortByPT_2lep(), and Gambit::ColliderBit::sortByPT_l().

{
   if (!momenta_set)
   {
       cout <<" Please set momenta first!" << endl;
       return 0;
   }
        
   if (!solved) mt2_bisect();
   return mt2_b*scale;
}

Here is the caller graph for this function:

mt2_bisect()

void mt2_bisect::mt2::mt2_bisect

(

)

Definition at line 332 of file mt2_bisect.cpp.

References double, and MIN_MASS.

{
  
   
   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;
}

mt2_massless()

void mt2_bisect::mt2::mt2_massless

(

)

Definition at line 147 of file mt2_bisect.cpp.

References double, and SCANSTEP.

{
   
//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;
}

nsols()

int mt2_bisect::mt2::nsols ( double  Dsq )

private

Definition at line 545 of file mt2_bisect.cpp.

References Gambit::FlavBit::LoopFunctions::B0(), Gambit::FlavBit::LoopFunctions::B1(), Gambit::FlavBit::LoopFunctions::C0(), Gambit::FlavBit::LoopFunctions::D0(), and daFunk::delta().

{
   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;
  
}  

Here is the call graph for this function:

nsols_massless()

int mt2_bisect::mt2::nsols_massless ( double  Dsq )

private

Definition at line 274 of file mt2_bisect.cpp.

References b, Gambit::FlavBit::LoopFunctions::B0(), Gambit::FlavBit::LoopFunctions::B1(), Gambit::FlavBit::LoopFunctions::C0(), Gambit::FlavBit::LoopFunctions::D0(), and daFunk::delta().

{
  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;
  
}

Here is the call graph for this function:

print()

void mt2_bisect::mt2::print

(

)

Definition at line 138 of file mt2_bisect.cpp.

{
   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;
}

scan_high()

int mt2_bisect::mt2::scan_high ( double &  Deltasq_high )

private

Definition at line 523 of file mt2_bisect.cpp.

References SCANSTEP.

{
   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;
}

set_mn()

void mt2_bisect::mt2::set_mn

(

double

mn

)

Definition at line 130 of file mt2_bisect.cpp.

Referenced by Gambit::ColliderBit::calcMT(), Gambit::ColliderBit::calcMT_1l(), Gambit::ColliderBit::Phi_mpi_pi(), Gambit::ColliderBit::sortByPT_2lep(), and Gambit::ColliderBit::sortByPT_l().

{
   solved = false;    //reset solved tag when mn is changed.
   mn_unscale   = fabs(mn0);  //mass cannot be negative
   mn = mn_unscale/scale;
   mnsq = mn*mn;
}

Here is the caller graph for this function:

set_momenta()

void mt2_bisect::mt2::set_momenta	(	double *	pa0,
		double *	pb0,
		double *	pmiss0
	)

Definition at line 62 of file mt2_bisect.cpp.

References ABSOLUTE_PRECISION, RELATIVE_PRECISION, and ZERO_MASS.

Referenced by Gambit::ColliderBit::calcMT(), Gambit::ColliderBit::calcMT_1l(), Gambit::ColliderBit::Phi_mpi_pi(), Gambit::ColliderBit::sortByPT_2lep(), and Gambit::ColliderBit::sortByPT_l().

{
   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;
}

Here is the caller graph for this function:

signchange_n()

int mt2_bisect::mt2::signchange_n ( long double  t1, long double  t2, long double  t3, long double  t4, long double  t5  )

inlineprivate

Definition at line 630 of file mt2_bisect.cpp.

{
   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;
}

signchange_p()

int mt2_bisect::mt2::signchange_p ( long double  t1, long double  t2, long double  t3, long double  t4, long double  t5  )

inlineprivate

Definition at line 640 of file mt2_bisect.cpp.

{
   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;
}

Member Data Documentation

a1

double mt2_bisect::mt2::a1

private

Definition at line 54 of file mt2_bisect.h.

a2

double mt2_bisect::mt2::a2

private

Definition at line 54 of file mt2_bisect.h.

b1

double mt2_bisect::mt2::b1

private

Definition at line 54 of file mt2_bisect.h.

b2

double mt2_bisect::mt2::b2

private

Definition at line 54 of file mt2_bisect.h.

c1

double mt2_bisect::mt2::c1

private

Definition at line 54 of file mt2_bisect.h.

c2

double mt2_bisect::mt2::c2

private

Definition at line 54 of file mt2_bisect.h.

d1

double mt2_bisect::mt2::d1

private

Definition at line 54 of file mt2_bisect.h.

d11

double mt2_bisect::mt2::d11

private

Definition at line 55 of file mt2_bisect.h.

d2

double mt2_bisect::mt2::d2

private

Definition at line 54 of file mt2_bisect.h.

d20

double mt2_bisect::mt2::d20

private

Definition at line 55 of file mt2_bisect.h.

d21

double mt2_bisect::mt2::d21

private

Definition at line 55 of file mt2_bisect.h.

e1

double mt2_bisect::mt2::e1

private

Definition at line 54 of file mt2_bisect.h.

e11

double mt2_bisect::mt2::e11

private

Definition at line 55 of file mt2_bisect.h.

e2

double mt2_bisect::mt2::e2

private

Definition at line 54 of file mt2_bisect.h.

e20

double mt2_bisect::mt2::e20

private

Definition at line 55 of file mt2_bisect.h.

e21

double mt2_bisect::mt2::e21

private

Definition at line 55 of file mt2_bisect.h.

Ea

double mt2_bisect::mt2::Ea

private

Definition at line 42 of file mt2_bisect.h.

Easq

double mt2_bisect::mt2::Easq

private

Definition at line 48 of file mt2_bisect.h.

Eb

double mt2_bisect::mt2::Eb

private

Definition at line 44 of file mt2_bisect.h.

Ebsq

double mt2_bisect::mt2::Ebsq

private

Definition at line 49 of file mt2_bisect.h.

f1

double mt2_bisect::mt2::f1

private

Definition at line 54 of file mt2_bisect.h.

f10

double mt2_bisect::mt2::f10

private

Definition at line 55 of file mt2_bisect.h.

f12

double mt2_bisect::mt2::f12

private

Definition at line 55 of file mt2_bisect.h.

f2

double mt2_bisect::mt2::f2

private

Definition at line 54 of file mt2_bisect.h.

f20

double mt2_bisect::mt2::f20

private

Definition at line 55 of file mt2_bisect.h.

f21

double mt2_bisect::mt2::f21

private

Definition at line 55 of file mt2_bisect.h.

f22

double mt2_bisect::mt2::f22

private

Definition at line 55 of file mt2_bisect.h.

ma

double mt2_bisect::mt2::ma

private

Definition at line 42 of file mt2_bisect.h.

masq

double mt2_bisect::mt2::masq

private

Definition at line 48 of file mt2_bisect.h.

mb

double mt2_bisect::mt2::mb

private

Definition at line 44 of file mt2_bisect.h.

mbsq

double mt2_bisect::mt2::mbsq

private

Definition at line 49 of file mt2_bisect.h.

mn

double mt2_bisect::mt2::mn

private

Definition at line 45 of file mt2_bisect.h.

mn_unscale

double mt2_bisect::mt2::mn_unscale

private

Definition at line 45 of file mt2_bisect.h.

mnsq

double mt2_bisect::mt2::mnsq

private

Definition at line 51 of file mt2_bisect.h.

momenta_set

bool mt2_bisect::mt2::momenta_set

private

Definition at line 32 of file mt2_bisect.h.

mt2_b

double mt2_bisect::mt2::mt2_b

private

Definition at line 33 of file mt2_bisect.h.

nevt

int mt2_bisect::mt2::nevt

Definition at line 28 of file mt2_bisect.h.

pax

double mt2_bisect::mt2::pax

private

Definition at line 42 of file mt2_bisect.h.

pay

double mt2_bisect::mt2::pay

private

Definition at line 42 of file mt2_bisect.h.

pbx

double mt2_bisect::mt2::pbx

private

Definition at line 44 of file mt2_bisect.h.

pby

double mt2_bisect::mt2::pby

private

Definition at line 44 of file mt2_bisect.h.

pmissx

double mt2_bisect::mt2::pmissx

private

Definition at line 43 of file mt2_bisect.h.

pmissxsq

double mt2_bisect::mt2::pmissxsq

private

Definition at line 50 of file mt2_bisect.h.

pmissy

double mt2_bisect::mt2::pmissy

private

Definition at line 43 of file mt2_bisect.h.

pmissysq

double mt2_bisect::mt2::pmissysq

private

Definition at line 50 of file mt2_bisect.h.

precision

double mt2_bisect::mt2::precision

private

Definition at line 58 of file mt2_bisect.h.

scale

double mt2_bisect::mt2::scale

private

Definition at line 57 of file mt2_bisect.h.

solved

bool mt2_bisect::mt2::solved

private

Definition at line 31 of file mt2_bisect.h.

---

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

• ColliderBit/include/gambit/ColliderBit/mt2_bisect.h
• ColliderBit/src/mt2_bisect.cpp