rivet is hosted by Hepforge, IPPP Durham

#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 48 of file mt2_bisect.cc.

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

Member Function Documentation

int find_high ( double &  Deltasq_high) [private]

Definition at line 486 of file mt2_bisect.cc.

{
   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;
}  
double get_mt2 ( )

Definition at line 56 of file mt2_bisect.cc.

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

{
   assert(momenta_set);
   if (!solved) mt2_bisect();
   return mt2_b*scale;
}
void mt2_bisect ( )

Definition at line 336 of file mt2_bisect.cc.

References Log::getLog(), Rivet::mt2_bisect::MIN_MASS, and Log::WARN.

{   
  solved = true;

//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) 
      {
     Log::getLog("Rivet.Tools.mt2") 
       << Log::WARN  
       << "Deltasq_high not found at event " << nevt << '\n';
         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;
}
void mt2_massless ( )

Definition at line 148 of file mt2_bisect.cc.

References Log::ERROR, Log::getLog(), Rivet::mass(), Rivet::s, Rivet::mt2_bisect::SCANSTEP, Rivet::theta(), and Log::WARN.

{
   
//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) 
      {
       
    Log::getLog("Rivet.Tools.mt2") 
      << Log::WARN  
      << "Deltasq_high not found at event " << nevt << '\n';
       
         mt2_b = (double)sqrt(Deltasq_low+mnsq);
         return;
      }
   }

   if(nsols_high == nsols_low)
   { 
     Log::getLog("Rivet.Tools.mt2") 
       << Log::ERROR  
       << "error: nsols_low=nsols_high=" << nsols_high << '\n'
       << "Deltasq_high=" << Deltasq_high << '\n'
       << "Deltasq_low= "<< Deltasq_low << '\n';
    
      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 nsols ( double  Dsq) [private]

Definition at line 552 of file mt2_bisect.cc.

References Rivet::PID::B0, and Rivet::PID::D0.

{
   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;
  
}  
int nsols_massless ( double  Dsq) [private]

Definition at line 278 of file mt2_bisect.cc.

References Rivet::PID::B0, and Rivet::PID::D0.

{
  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;
  
}
int scan_high ( double &  Deltasq_high) [private]

Definition at line 526 of file mt2_bisect.cc.

References Rivet::mass(), and Rivet::mt2_bisect::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;
}
void set_mn ( double  mn)

Definition at line 131 of file mt2_bisect.cc.

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

{
   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 set_momenta ( double *  pa0,
double *  pb0,
double *  pmiss0 
)

Definition at line 63 of file mt2_bisect.cc.

References Rivet::mt2_bisect::ABSOLUTE_PRECISION, Rivet::mt2_bisect::RELATIVE_PRECISION, and Rivet::mt2_bisect::ZERO_MASS.

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

{
   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;
}
int signchange_n ( long double  t1,
long double  t2,
long double  t3,
long double  t4,
long double  t5 
) [inline, private]

Definition at line 637 of file mt2_bisect.cc.

{
   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;
}
int signchange_p ( long double  t1,
long double  t2,
long double  t3,
long double  t4,
long double  t5 
) [inline, private]

Definition at line 647 of file mt2_bisect.cc.

{
   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

double a1 [private]

Definition at line 74 of file mt2_bisect.hh.

double a2 [private]

Definition at line 74 of file mt2_bisect.hh.

double b1 [private]

Definition at line 74 of file mt2_bisect.hh.

double b2 [private]

Definition at line 74 of file mt2_bisect.hh.

double c1 [private]

Definition at line 74 of file mt2_bisect.hh.

double c2 [private]

Definition at line 74 of file mt2_bisect.hh.

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.

double d21 [private]

Definition at line 75 of file mt2_bisect.hh.

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.

double e21 [private]

Definition at line 75 of file mt2_bisect.hh.

double Ea [private]

Definition at line 62 of file mt2_bisect.hh.

double Easq [private]

Definition at line 68 of file mt2_bisect.hh.

double Eb [private]

Definition at line 64 of file mt2_bisect.hh.

double Ebsq [private]

Definition at line 69 of file mt2_bisect.hh.

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.

double f21 [private]

Definition at line 75 of file mt2_bisect.hh.

double f22 [private]

Definition at line 75 of file mt2_bisect.hh.

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.

bool momenta_set [private]

Definition at line 52 of file mt2_bisect.hh.

double mt2_b [private]

Definition at line 53 of file mt2_bisect.hh.

int nevt

Definition at line 48 of file mt2_bisect.hh.

double pax [private]

Definition at line 62 of file mt2_bisect.hh.

double pay [private]

Definition at line 62 of file mt2_bisect.hh.

double pbx [private]

Definition at line 64 of file mt2_bisect.hh.

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.

double pmissy [private]

Definition at line 63 of file mt2_bisect.hh.

double pmissysq [private]

Definition at line 70 of file mt2_bisect.hh.

double precision [private]

Definition at line 78 of file mt2_bisect.hh.

double scale [private]

Definition at line 77 of file mt2_bisect.hh.

bool solved [private]

Definition at line 51 of file mt2_bisect.hh.


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