gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool

#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.

43 {
44  solved = false;
45  momenta_set = false;
46  mt2_b = 0.;
47  scale = 1.;
48 }

Member Function Documentation

◆ find_high()

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

Definition at line 483 of file mt2_bisect.cpp.

484 {
485  double x0,y0;
486  x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
487  y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
488  double Deltasq_low = (mn + ma)*(mn + ma) - mnsq;
489  do
490  {
491  double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.;
492  int nsols_mid = nsols(Deltasq_mid);
493  if ( nsols_mid == 2 )
494  {
495  Deltasq_high = Deltasq_mid;
496  return 1;
497  }
498  else if (nsols_mid == 4)
499  {
500  Deltasq_high = Deltasq_mid;
501  continue;
502  }
503  else if (nsols_mid ==0)
504  {
505  d1 = -pax*(Deltasq_mid-masq)/(2*Easq);
506  e1 = -pay*(Deltasq_mid-masq)/(2*Easq);
507  d2 = -pmissx + pbx*(Deltasq_mid - mbsq)/(2*Ebsq)
508  + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
509  e2 = -pmissy + pby*(Deltasq_mid - mbsq)/(2*Ebsq)
510  + pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
511  f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq)/(2*Eb)+
512  (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq)/(2*Eb)+
513  (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
514 // Does the larger ellipse contain the smaller one?
515  double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
516  if (dis < 0) Deltasq_high = Deltasq_mid;
517  else Deltasq_low = Deltasq_mid;
518  }
519 
520  } while ( Deltasq_high - Deltasq_low > 0.001);
521  return 0;
522 }
int nsols(double Dsq)
Definition: mt2_bisect.cpp:545

◆ 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().

51 {
52  if (!momenta_set)
53  {
54  cout <<" Please set momenta first!" << endl;
55  return 0;
56  }
57 
58  if (!solved) mt2_bisect();
59  return mt2_b*scale;
60 }
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.

333 {
334 
335 
336  solved = true;
337  cout.precision(11);
338 
339 //if masses are very small, use code for massless case.
340  if(masq < MIN_MASS && mbsq < MIN_MASS)
341  {
342  mt2_massless();
343  return;
344  }
345 
346 
347  double Deltasq0;
348  Deltasq0 = ma*(ma + 2*mn); //The minimum mass square to have two ellipses
349 
350 // find the coefficients for the two quadratic equations when Deltasq=Deltasq0.
351 
352  a1 = 1-pax*pax/(Easq);
353  b1 = -pax*pay/(Easq);
354  c1 = 1-pay*pay/(Easq);
355  d1 = -pax*(Deltasq0-masq)/(2*Easq);
356  e1 = -pay*(Deltasq0-masq)/(2*Easq);
357  a2 = 1-pbx*pbx/(Ebsq);
358  b2 = -pbx*pby/(Ebsq);
359  c2 = 1-pby*pby/(Ebsq);
360  d2 = -pmissx+pbx*(Deltasq0-mbsq)/(2*Ebsq)+pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
361  e2 = -pmissy+pby*(Deltasq0-mbsq)/(2*Ebsq)+pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
362  f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq0-mbsq)/(2*Eb)+
363  (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq0-mbsq)/(2*Eb)+
364  (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
365 
366 // find the center of the smaller ellipse
367  double x0,y0;
368  x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
369  y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
370 
371 
372 // Does the larger ellipse contain the smaller one?
373  double dis=a2*x0*x0+2*b2*x0*y0+c2*y0*y0+2*d2*x0+2*e2*y0+f2;
374 
375  if(dis<=0.01)
376  {
377  mt2_b = (double) sqrt(mnsq+Deltasq0);
378  return;
379  }
380 
381 
382 /* find the coefficients for the two quadratic equations */
383 /* coefficients for quadratic terms do not change */
384 /* coefficients for linear and constant terms are polynomials of */
385 /* delta=(Deltasq-m7sq)/(2 E7sq) */
386  d11 = -pax;
387  e11 = -pay;
388  f10 = mnsq;
389  f12 = -Easq;
390  d21 = (Easq*pbx)/Ebsq;
391  d20 = ((masq - mbsq)*pbx)/(2.*Ebsq) - pmissx +
392  (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
393  e21 = (Easq*pby)/Ebsq;
394  e20 = ((masq - mbsq)*pby)/(2.*Ebsq) - pmissy +
395  (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
396  f22 = -Easq*Easq/Ebsq;
397  f21 = (-2*Easq*((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb))/Eb;
398  f20 = mnsq + pmissx*pmissx + pmissy*pmissy -
399  ((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb)
400  *((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb);
401 
402 //Estimate upper bound of mT2
403 //when Deltasq > Deltasq_high1, the larger encloses the center of the smaller
404  double p2x0,p2y0;
405  double Deltasq_high1;
406  p2x0 = pmissx-x0;
407  p2y0 = pmissy-y0;
408  Deltasq_high1 = 2*Eb*sqrt(p2x0*p2x0+p2y0*p2y0+mnsq)-2*pbx*p2x0-2*pby*p2y0+mbsq;
409 
410 //Another estimate, if both ellipses enclose the origin, Deltasq > mT2
411 
412  double Deltasq_high2, Deltasq_high21, Deltasq_high22;
413  Deltasq_high21 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy+mbsq;
414  Deltasq_high22 = 2*Ea*mn+masq;
415 
416  if ( Deltasq_high21 < Deltasq_high22 ) Deltasq_high2 = Deltasq_high22;
417  else Deltasq_high2 = Deltasq_high21;
418 
419 //pick the smaller upper bound
420  double Deltasq_high;
421  if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high1;
422  else Deltasq_high = Deltasq_high2;
423 
424 
425  double Deltasq_low; //lower bound
426  Deltasq_low = Deltasq0;
427 
428 //number of solutions at Deltasq_low should not be larger than zero
429  if( nsols(Deltasq_low) > 0 )
430  {
431  //cout << "nsolutions(Deltasq_low) > 0"<<endl;
432  mt2_b = (double) sqrt(mnsq+Deltasq0);
433  return;
434  }
435 
436  int nsols_high, nsols_low;
437 
438  nsols_low = nsols(Deltasq_low);
439  int foundhigh;
440 
441 
442 //if nsols_high=nsols_low, we missed the region where the two ellipse overlap
443 //if nsols_high=4, also need a scan because we may find the wrong tangent point.
444 
445  nsols_high = nsols(Deltasq_high);
446 
447  if(nsols_high == nsols_low || nsols_high == 4)
448  {
449  //foundhigh = scan_high(Deltasq_high);
450  foundhigh = find_high(Deltasq_high);
451  if(foundhigh == 0)
452  {
453  cout << "Deltasq_high not found at event " << nevt << endl;
454  mt2_b = sqrt( Deltasq_low + mnsq );
455  return;
456  }
457 
458  }
459 
460  while(sqrt(Deltasq_high+mnsq) - sqrt(Deltasq_low+mnsq) > precision)
461  {
462  double Deltasq_mid,nsols_mid;
463  //bisect
464  Deltasq_mid = (Deltasq_high+Deltasq_low)/2.;
465  nsols_mid = nsols(Deltasq_mid);
466  // if nsols_mid = 4, rescan for Deltasq_high
467  if ( nsols_mid == 4 )
468  {
469  Deltasq_high = Deltasq_mid;
470  //scan_high(Deltasq_high);
471  find_high(Deltasq_high);
472  continue;
473  }
474 
475 
476  if(nsols_mid != nsols_low) Deltasq_high = Deltasq_mid;
477  if(nsols_mid == nsols_low) Deltasq_low = Deltasq_mid;
478  }
479  mt2_b = (double) sqrt( mnsq + Deltasq_high);
480  return;
481 }
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
#define MIN_MASS
Definition: mt2_bisect.h:12
int nsols(double Dsq)
Definition: mt2_bisect.cpp:545
double precision
Definition: mt2_bisect.h:58
int find_high(double &Deltasq_high)
Definition: mt2_bisect.cpp:483

◆ mt2_massless()

void mt2_bisect::mt2::mt2_massless ( )

Definition at line 147 of file mt2_bisect.cpp.

References double, and SCANSTEP.

148 {
149 
150 //rotate so that pay = 0
151  double theta,s,c;
152  theta = atan(pay/pax);
153  s = sin(theta);
154  c = cos(theta);
155 
156  double pxtemp,pytemp;
157  Easq = pax*pax+pay*pay;
158  Ebsq = pbx*pbx+pby*pby;
159  Ea = sqrt(Easq);
160  Eb = sqrt(Ebsq);
161 
162  pxtemp = pax*c+pay*s;
163  pax = pxtemp;
164  pay = 0;
165  pxtemp = pbx*c+pby*s;
166  pytemp = -s*pbx+c*pby;
167  pbx = pxtemp;
168  pby = pytemp;
169  pxtemp = pmissx*c+pmissy*s;
170  pytemp = -s*pmissx+c*pmissy;
171  pmissx = pxtemp;
172  pmissy = pytemp;
173 
174  a2 = 1-pbx*pbx/(Ebsq);
175  b2 = -pbx*pby/(Ebsq);
176  c2 = 1-pby*pby/(Ebsq);
177 
178  d21 = (Easq*pbx)/Ebsq;
179  d20 = - pmissx + (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
180  e21 = (Easq*pby)/Ebsq;
181  e20 = - pmissy + (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
182  f22 = -(Easq*Easq/Ebsq);
183  f21 = -2*Easq*(pbx*pmissx + pby*pmissy)/Ebsq;
185 
186  double Deltasq0 = 0;
187  double Deltasq_low, Deltasq_high;
188  int nsols_high, nsols_low;
189 
190  Deltasq_low = Deltasq0 + precision;
191  nsols_low = nsols_massless(Deltasq_low);
192 
193  if(nsols_low > 1)
194  {
195  mt2_b = (double) sqrt(Deltasq0+mnsq);
196  return;
197  }
198 
199 /*
200  if( nsols_massless(Deltasq_high) > 0 )
201  {
202  mt2_b = (double) sqrt(mnsq+Deltasq0);
203  return;
204  }*/
205 
206 //look for when both parablos contain origin
207  double Deltasq_high1, Deltasq_high2;
208  Deltasq_high1 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy;
209  Deltasq_high2 = 2*Ea*mn;
210 
211  if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high2;
212  else Deltasq_high = Deltasq_high1;
213 
214  nsols_high=nsols_massless(Deltasq_high);
215 
216  int foundhigh;
217  if (nsols_high == nsols_low)
218  {
219 
220 
221  foundhigh=0;
222 
223  double minmass, maxmass;
224  minmass = mn ;
225  maxmass = sqrt(mnsq + Deltasq_high);
226  for(double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
227  {
228  Deltasq_high = mass*mass - mnsq;
229 
230  nsols_high = nsols_massless(Deltasq_high);
231  if(nsols_high>0)
232  {
233  foundhigh=1;
234  Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
235  break;
236  }
237  }
238  if(foundhigh==0)
239  {
240 
241  //cout<<"Deltasq_high not found at event " << nevt <<endl;
242 
243 
244  mt2_b = (double)sqrt(Deltasq_low+mnsq);
245  return;
246  }
247  }
248 
249  if(nsols_high == nsols_low)
250  {
251  cout << "error: nsols_low=nsols_high=" << nsols_high << endl;
252  cout << "Deltasq_high=" << Deltasq_high << endl;
253  cout << "Deltasq_low= "<< Deltasq_low << endl;
254 
255  mt2_b = sqrt(mnsq + Deltasq_low);
256  return;
257  }
258  double minmass, maxmass;
259  minmass = sqrt(Deltasq_low+mnsq);
260  maxmass = sqrt(Deltasq_high+mnsq);
261  while(maxmass - minmass > precision)
262  {
263  double Delta_mid, midmass, nsols_mid;
264  midmass = (minmass+maxmass)/2.;
265  Delta_mid = midmass * midmass - mnsq;
266  nsols_mid = nsols_massless(Delta_mid);
267  if(nsols_mid != nsols_low) maxmass = midmass;
268  if(nsols_mid == nsols_low) minmass = midmass;
269  }
270  mt2_b = minmass;
271  return;
272 }
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
#define SCANSTEP
Definition: mt2_bisect.h:14
double pmissysq
Definition: mt2_bisect.h:50
double pmissxsq
Definition: mt2_bisect.h:50
double precision
Definition: mt2_bisect.h:58
int nsols_massless(double Dsq)
Definition: mt2_bisect.cpp:274

◆ 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().

546 {
547  double delta = (Dsq-masq)/(2*Easq);
548 
549 //calculate coefficients for the two quadratic equations
550  d1 = d11*delta;
551  e1 = e11*delta;
552  f1 = f12*delta*delta+f10;
553  d2 = d21*delta+d20;
554  e2 = e21*delta+e20;
555  f2 = f22*delta*delta+f21*delta+f20;
556 
557 //obtain the coefficients for the 4th order equation
558 //devided by Ea^n to make the variable dimensionless
559  long double A4, A3, A2, A1, A0;
560 
561  A4 =
562  -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 +
563  4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 +
564  a1*a1*c2*c2;
565 
566  A3 =
567  (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 +
568  8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 +
569  4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 -
570  4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea;
571 
572 
573  A2 =
574  (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 -
575  8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 +
576  4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 -
577  8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 +
578  4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 +
579  4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq;
580 
581  A1 =
582  (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 -
583  4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 -
584  4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 -
585  4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea);
586 
587  A0 =
588  (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 +
589  4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 +
590  a1*a1*f2*f2)/(Easq*Easq);
591 
592  long double B3, B2, B1, B0;
593  B3 = 4*A4;
594  B2 = 3*A3;
595  B1 = 2*A2;
596  B0 = A1;
597 
598  long double C2, C1, C0;
599  C2 = -(A2/2 - 3*A3*A3/(16*A4));
600  C1 = -(3*A1/4. -A2*A3/(8*A4));
601  C0 = -A0 + A1*A3/(16*A4);
602 
603  long double D1, D0;
604  D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
605  D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
606 
607  long double E0;
608  E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
609 
610  long double t1,t2,t3,t4,t5;
611 //find the coefficients for the leading term in the Sturm sequence
612  t1 = A4;
613  t2 = A4;
614  t3 = C2;
615  t4 = D1;
616  t5 = E0;
617 
618 
619 //The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
620  int nsol;
621  nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);
622 
623 //Cannot have negative number of solutions, must be roundoff effect
624  if (nsol < 0) nsol = 0;
625 
626  return nsol;
627 
628 }
Funk delta(std::string arg, double pos, double width)
Definition: daFunk.hpp:902
int signchange_p(long double t1, long double t2, long double t3, long double t4, long double t5)
Definition: mt2_bisect.cpp:640
double D0(const double a, const double b, const double c, const double d)
double C0(const double a, const double b, const double c)
double B0(const double a, const double b, const double Q)
int signchange_n(long double t1, long double t2, long double t3, long double t4, long double t5)
Definition: mt2_bisect.cpp:630
double B1(const double a, const double b, const double Q)
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().

275 {
276  double delta;
277  delta = Dsq/(2*Easq);
278  d1 = d11*delta;
279  e1 = e11*delta;
280  f1 = f12*delta*delta+f10;
281  d2 = d21*delta+d20;
282  e2 = e21*delta+e20;
283  f2 = f22*delta*delta+f21*delta+f20;
284 
285  double a,b;
286  if (pax > 0) a = Ea/Dsq;
287  else a = -Ea/Dsq;
288  if (pax > 0) b = -Dsq/(4*Ea)+mnsq*Ea/Dsq;
289  else b = Dsq/(4*Ea)-mnsq*Ea/Dsq;
290 
291  double A4,A3,A2,A1,A0;
292 
293  A4 = a*a*a2;
294  A3 = 2*a*b2/Ea;
295  A2 = (2*a*a2*b+c2+2*a*d2)/(Easq);
296  A1 = (2*b*b2+2*e2)/(Easq*Ea);
297  A0 = (a2*b*b+2*b*d2+f2)/(Easq*Easq);
298 
299  long double B3, B2, B1, B0;
300  B3 = 4*A4;
301  B2 = 3*A3;
302  B1 = 2*A2;
303  B0 = A1;
304  long double C2, C1, C0;
305  C2 = -(A2/2 - 3*A3*A3/(16*A4));
306  C1 = -(3*A1/4. -A2*A3/(8*A4));
307  C0 = -A0 + A1*A3/(16*A4);
308  long double D1, D0;
309  D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
310  D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
311 
312  long double E0;
313  E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
314 
315  long double t1,t2,t3,t4,t5;
316 
317 //find the coefficients for the leading term in the Sturm sequence
318  t1 = A4;
319  t2 = A4;
320  t3 = C2;
321  t4 = D1;
322  t5 = E0;
323 
324  int nsol;
325  nsol = signchange_n(t1,t2,t3,t4,t5)-signchange_p(t1,t2,t3,t4,t5);
326  if( nsol < 0 ) nsol=0;
327 
328  return nsol;
329 
330 }
Funk delta(std::string arg, double pos, double width)
Definition: daFunk.hpp:902
int signchange_p(long double t1, long double t2, long double t3, long double t4, long double t5)
Definition: mt2_bisect.cpp:640
double D0(const double a, const double b, const double c, const double d)
START_MODEL b
Definition: demo.hpp:270
double C0(const double a, const double b, const double c)
double B0(const double a, const double b, const double Q)
int signchange_n(long double t1, long double t2, long double t3, long double t4, long double t5)
Definition: mt2_bisect.cpp:630
double B1(const double a, const double b, const double Q)
Here is the call graph for this function:

◆ print()

void mt2_bisect::mt2::print ( )

Definition at line 138 of file mt2_bisect.cpp.

139 {
140  cout << " pax = " << pax*scale << "; pay = " << pay*scale << "; ma = " << ma*scale <<";"<< endl;
141  cout << " pbx = " << pbx*scale << "; pby = " << pby*scale << "; mb = " << mb*scale <<";"<< endl;
142  cout << " pmissx = " << pmissx*scale << "; pmissy = " << pmissy*scale <<";"<< endl;
143  cout << " mn = " << mn_unscale<<";" << endl;
144 }
double mn_unscale
Definition: mt2_bisect.h:45

◆ scan_high()

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

Definition at line 523 of file mt2_bisect.cpp.

References SCANSTEP.

524 {
525  int foundhigh = 0 ;
526  int nsols_high;
527 
528  double tempmass, maxmass;
529  tempmass = mn + ma;
530  maxmass = sqrt(mnsq + Deltasq_high);
531  if (nevt == 32334) cout << "Deltasq_high = " << Deltasq_high << endl;
532  for(double mass = tempmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
533  {
534  Deltasq_high = mass*mass - mnsq;
535  nsols_high = nsols(Deltasq_high);
536 
537  if( nsols_high > 0)
538  {
539  foundhigh = 1;
540  break;
541  }
542  }
543  return foundhigh;
544 }
#define SCANSTEP
Definition: mt2_bisect.h:14
int nsols(double Dsq)
Definition: mt2_bisect.cpp:545

◆ 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().

131 {
132  solved = false; //reset solved tag when mn is changed.
133  mn_unscale = fabs(mn0); //mass cannot be negative
134  mn = mn_unscale/scale;
135  mnsq = mn*mn;
136 }
double mn_unscale
Definition: mt2_bisect.h:45
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().

63 {
64  solved = false; //reset solved tag when momenta are changed.
65  momenta_set = true;
66 
67  ma = fabs(pa0[0]); // mass cannot be negative
68 
69  if (ma < ZERO_MASS) ma = ZERO_MASS;
70 
71  pax = pa0[1];
72  pay = pa0[2];
73  masq = ma*ma;
74  Easq = masq+pax*pax+pay*pay;
75  Ea = sqrt(Easq);
76 
77  mb = fabs(pb0[0]);
78 
79  if (mb < ZERO_MASS) mb = ZERO_MASS;
80 
81  pbx = pb0[1];
82  pby = pb0[2];
83  mbsq = mb*mb;
84  Ebsq = mbsq+pbx*pbx+pby*pby;
85  Eb = sqrt(Ebsq);
86 
87  pmissx = pmiss0[1]; pmissy = pmiss0[2];
90 
91 // set ma>= mb
92  if(masq < mbsq)
93  {
94  double temp;
95  temp = pax; pax = pbx; pbx = temp;
96  temp = pay; pay = pby; pby = temp;
97  temp = Ea; Ea = Eb; Eb = temp;
98  temp = Easq; Easq = Ebsq; Ebsq = temp;
99  temp = masq; masq = mbsq; mbsq = temp;
100  temp = ma; ma = mb; mb = temp;
101  }
102 //normalize max{Ea, Eb} to 100
103  if (Ea > Eb) scale = Ea/100.;
104  else scale = Eb/100.;
105 
106  if (sqrt(pmissxsq+pmissysq)/100 > scale) scale = sqrt(pmissxsq+pmissysq)/100;
107  //scale = 1;
108  double scalesq = scale * scale;
109  ma = ma/scale;
110  mb = mb/scale;
111  masq = masq/scalesq;
112  mbsq = mbsq/scalesq;
113  pax = pax/scale; pay = pay/scale;
114  pbx = pbx/scale; pby = pby/scale;
115  Ea = Ea/scale; Eb = Eb/scale;
116 
117  Easq = Easq/scalesq;
118  Ebsq = Ebsq/scalesq;
119  pmissx = pmissx/scale;
120  pmissy = pmissy/scale;
121  pmissxsq = pmissxsq/scalesq;
122  pmissysq = pmissysq/scalesq;
123  mn = mn_unscale/scale;
124  mnsq = mn*mn;
125 
127  else precision = 100.*RELATIVE_PRECISION;
128 }
double pmissysq
Definition: mt2_bisect.h:50
double pmissxsq
Definition: mt2_bisect.h:50
double mn_unscale
Definition: mt2_bisect.h:45
double precision
Definition: mt2_bisect.h:58
#define RELATIVE_PRECISION
Definition: mt2_bisect.h:7
#define ABSOLUTE_PRECISION
Definition: mt2_bisect.h:8
#define ZERO_MASS
Definition: mt2_bisect.h:13
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.

631 {
632  int nsc;
633  nsc=0;
634  if(t1*t2>0) nsc++;
635  if(t2*t3>0) nsc++;
636  if(t3*t4>0) nsc++;
637  if(t4*t5>0) nsc++;
638  return nsc;
639 }

◆ 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.

641 {
642  int nsc;
643  nsc=0;
644  if(t1*t2<0) nsc++;
645  if(t2*t3<0) nsc++;
646  if(t3*t4<0) nsc++;
647  if(t4*t5<0) nsc++;
648  return nsc;
649 }

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: