|
|
[Rivet-svn] r3054 - in trunk: include/Rivet include/Rivet/Tools src/Tools
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Wed Apr 13 15:24:22 BST 2011
Author: dgrell
Date: Wed Apr 13 15:24:21 2011
New Revision: 3054
Log:
mT2 support: Cheng/Han code and wrapper
Added:
trunk/include/Rivet/Tools/RivetMT2.hh
trunk/src/Tools/RivetMT2.cc
trunk/src/Tools/mt2_bisect.cc (contents, props changed)
trunk/src/Tools/mt2_bisect.hh
Modified:
trunk/include/Rivet/Makefile.am
trunk/src/Tools/Makefile.am
Modified: trunk/include/Rivet/Makefile.am
==============================================================================
--- trunk/include/Rivet/Makefile.am Fri Apr 8 18:03:59 2011 (r3053)
+++ trunk/include/Rivet/Makefile.am Wed Apr 13 15:24:21 2011 (r3054)
@@ -90,7 +90,8 @@
Tools/RivetPaths.hh \
Tools/BinnedHistogram.hh \
Tools/ParticleIdUtils.hh \
- Tools/TypeTraits.hh
+ Tools/TypeTraits.hh \
+ Tools/RivetMT2.hh
nobase_dist_noinst_HEADERS += \
Tools/osdir.hh
Added: trunk/include/Rivet/Tools/RivetMT2.hh
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/include/Rivet/Tools/RivetMT2.hh Wed Apr 13 15:24:21 2011 (r3054)
@@ -0,0 +1,23 @@
+// -*- C++ -*-
+#ifndef RIVET_MT2_HH
+#define RIVET_MT2_HH
+
+// Convenience wrapper for the m_T2 calculator of Cheng/Han
+// (arXiv:0810.5178). Could be adapted for other backends in future
+
+namespace Rivet {
+
+ class FourMomentum;
+
+ namespace mT2 {
+
+ double mT2(const FourMomentum & a,
+ const FourMomentum & b,
+ const FourMomentum & pTmiss,
+ double invisiblesMass);
+
+ }
+}
+
+
+#endif
Modified: trunk/src/Tools/Makefile.am
==============================================================================
--- trunk/src/Tools/Makefile.am Fri Apr 8 18:03:59 2011 (r3053)
+++ trunk/src/Tools/Makefile.am Wed Apr 13 15:24:21 2011 (r3054)
@@ -14,9 +14,11 @@
TinyXML/tinyxml.cpp \
TinyXML/tinyxmlerror.cpp \
TinyXML/tinyxmlparser.cpp \
- BinnedHistogram.cc
+ BinnedHistogram.cc \
+ mt2_bisect.cc \
+ RivetMT2.cc
-dist_noinst_HEADERS = binreloc.h
+dist_noinst_HEADERS = binreloc.h mt2_bisect.hh
libRivetTools_la_CPPFLAGS = \
$(AM_CPPFLAGS) \
Added: trunk/src/Tools/RivetMT2.cc
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Tools/RivetMT2.cc Wed Apr 13 15:24:21 2011 (r3054)
@@ -0,0 +1,29 @@
+#include "Rivet/Tools/RivetMT2.hh"
+#include "Rivet/Math/Vector4.hh"
+#include "mt2_bisect.hh"
+
+
+namespace Rivet {
+ namespace mT2 {
+
+ double mT2(const FourMomentum & a,
+ const FourMomentum & b,
+ const FourMomentum & pTmiss,
+ double invisiblesMass)
+ {
+ mt2_bisect::mt2 mt2_event;
+
+ double unused = -999.999;
+
+ double pa[3] = { 0.0, a.x(), a.y() };
+ double pb[3] = { 0.0, b.x(), b.y() };
+ double pmiss[3] = { unused, pTmiss.x(), pTmiss.y() };
+
+ mt2_event.set_momenta( pa, pb, pmiss );
+ mt2_event.set_mn( invisiblesMass );
+
+ return mt2_event.get_mt2();
+ }
+
+ }
+}
Added: trunk/src/Tools/mt2_bisect.cc
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Tools/mt2_bisect.cc Wed Apr 13 15:24:21 2011 (r3054)
@@ -0,0 +1,676 @@
+/***********************************************************************/
+/* */
+/* Finding mt2 by Bisection */
+/* */
+/* Authors: Hsin-Chia Cheng, Zhenyu Han */
+/* Dec 11, 2008, v1.01a */
+/* */
+/* see arXiv:0810.5178 */
+/* */
+/* Wrapped for Rivet by D. Grellscheid */
+/* Apr 13, 2011 */
+/* */
+/***********************************************************************/
+
+
+/*******************************************************************************
+ Usage:
+
+ 1. Define an object of type "mt2":
+
+ mt2_bisect::mt2 mt2_event;
+
+ 2. Set momenta and the mass of the invisible particle, mn:
+
+ mt2_event.set_momenta( pa, pb, pmiss );
+ mt2_event.set_mass( mn );
+
+ where array pa[0..2], pb[0..2], pmiss[0..2] contains (mass,px,py)
+ for the visible particles and the missing momentum. pmiss[0] is not used.
+ All quantities are given in double.
+
+ 3. Use mt2::get_mt2() to obtain the value of mt2:
+
+ double mt2_value = mt2_event.get_mt2();
+
+*******************************************************************************/
+#include "mt2_bisect.hh"
+#include "Rivet/Tools/Logging.hh"
+#include <cmath>
+#include <cassert>
+
+namespace Rivet {
+
+namespace mt2_bisect
+{
+using namespace std;
+
+mt2::mt2()
+{
+ solved = false;
+ momenta_set = false;
+ mt2_b = 0.;
+ scale = 1.;
+}
+
+double mt2::get_mt2()
+{
+ assert(momenta_set);
+ if (!solved) mt2_bisect();
+ return mt2_b*scale;
+}
+
+void mt2::set_momenta(double* pa0, double* pb0, double* pmiss0)
+{
+ solved = false; //reset solved tag when momenta are changed.
+ momenta_set = true;
+
+ ma = fabs(pa0[0]); // mass cannot be negative
+
+ if (ma < ZERO_MASS) ma = ZERO_MASS;
+
+ pax = pa0[1];
+ pay = pa0[2];
+ masq = ma*ma;
+ Easq = masq+pax*pax+pay*pay;
+ Ea = sqrt(Easq);
+
+ mb = fabs(pb0[0]);
+
+ if (mb < ZERO_MASS) mb = ZERO_MASS;
+
+ pbx = pb0[1];
+ pby = pb0[2];
+ mbsq = mb*mb;
+ Ebsq = mbsq+pbx*pbx+pby*pby;
+ Eb = sqrt(Ebsq);
+
+ pmissx = pmiss0[1]; pmissy = pmiss0[2];
+ pmissxsq = pmissx*pmissx;
+ pmissysq = pmissy*pmissy;
+
+// set ma>= mb
+ if(masq < mbsq)
+ {
+ double temp;
+ temp = pax; pax = pbx; pbx = temp;
+ temp = pay; pay = pby; pby = temp;
+ temp = Ea; Ea = Eb; Eb = temp;
+ temp = Easq; Easq = Ebsq; Ebsq = temp;
+ temp = masq; masq = mbsq; mbsq = temp;
+ temp = ma; ma = mb; mb = temp;
+ }
+//normalize max{Ea, Eb} to 100
+ if (Ea > Eb) scale = Ea/100.;
+ else scale = Eb/100.;
+
+ if (sqrt(pmissxsq+pmissysq)/100 > scale) scale = sqrt(pmissxsq+pmissysq)/100;
+ //scale = 1;
+ double scalesq = scale * scale;
+ ma = ma/scale;
+ mb = mb/scale;
+ masq = masq/scalesq;
+ mbsq = mbsq/scalesq;
+ pax = pax/scale; pay = pay/scale;
+ pbx = pbx/scale; pby = pby/scale;
+ Ea = Ea/scale; Eb = Eb/scale;
+
+ Easq = Easq/scalesq;
+ Ebsq = Ebsq/scalesq;
+ pmissx = pmissx/scale;
+ pmissy = pmissy/scale;
+ pmissxsq = pmissxsq/scalesq;
+ pmissysq = pmissysq/scalesq;
+ mn = mn_unscale/scale;
+ mnsq = mn*mn;
+
+ if (ABSOLUTE_PRECISION > 100.*RELATIVE_PRECISION) precision = ABSOLUTE_PRECISION;
+ else precision = 100.*RELATIVE_PRECISION;
+}
+
+void mt2::set_mn(double mn0)
+{
+ solved = false; //reset solved tag when mn is changed.
+ mn_unscale = fabs(mn0); //mass cannot be negative
+ mn = mn_unscale/scale;
+ mnsq = mn*mn;
+}
+
+// void mt2::print()
+// {
+// cout << " pax = " << pax*scale << "; pay = " << pay*scale << "; ma = " << ma*scale <<";"<< endl;
+// cout << " pbx = " << pbx*scale << "; pby = " << pby*scale << "; mb = " << mb*scale <<";"<< endl;
+// cout << " pmissx = " << pmissx*scale << "; pmissy = " << pmissy*scale <<";"<< endl;
+// cout << " mn = " << mn_unscale<<";" << endl;
+// }
+
+//special case, the visible particle is massless
+void mt2::mt2_massless()
+{
+
+//rotate so that pay = 0
+ double theta,s,c;
+ theta = atan(pay/pax);
+ s = sin(theta);
+ c = cos(theta);
+
+ double pxtemp,pytemp;
+ Easq = pax*pax+pay*pay;
+ Ebsq = pbx*pbx+pby*pby;
+ Ea = sqrt(Easq);
+ Eb = sqrt(Ebsq);
+
+ pxtemp = pax*c+pay*s;
+ pax = pxtemp;
+ pay = 0;
+ pxtemp = pbx*c+pby*s;
+ pytemp = -s*pbx+c*pby;
+ pbx = pxtemp;
+ pby = pytemp;
+ pxtemp = pmissx*c+pmissy*s;
+ pytemp = -s*pmissx+c*pmissy;
+ pmissx = pxtemp;
+ pmissy = pytemp;
+
+ a2 = 1-pbx*pbx/(Ebsq);
+ b2 = -pbx*pby/(Ebsq);
+ c2 = 1-pby*pby/(Ebsq);
+
+ d21 = (Easq*pbx)/Ebsq;
+ d20 = - pmissx + (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
+ e21 = (Easq*pby)/Ebsq;
+ e20 = - pmissy + (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
+ f22 = -(Easq*Easq/Ebsq);
+ f21 = -2*Easq*(pbx*pmissx + pby*pmissy)/Ebsq;
+ f20 = mnsq + pmissxsq + pmissysq - (pbx*pmissx + pby*pmissy)*(pbx*pmissx + pby*pmissy)/Ebsq;
+
+ double Deltasq0 = 0;
+ double Deltasq_low, Deltasq_high;
+ int nsols_high, nsols_low;
+
+ Deltasq_low = Deltasq0 + precision;
+ nsols_low = nsols_massless(Deltasq_low);
+
+ if(nsols_low > 1)
+ {
+ mt2_b = (double) sqrt(Deltasq0+mnsq);
+ return;
+ }
+
+/*
+ if( nsols_massless(Deltasq_high) > 0 )
+ {
+ mt2_b = (double) sqrt(mnsq+Deltasq0);
+ return;
+ }*/
+
+//look for when both parablos contain origin
+ double Deltasq_high1, Deltasq_high2;
+ Deltasq_high1 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy;
+ Deltasq_high2 = 2*Ea*mn;
+
+ if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high2;
+ else Deltasq_high = Deltasq_high1;
+
+ nsols_high=nsols_massless(Deltasq_high);
+
+ int foundhigh;
+ if (nsols_high == nsols_low)
+ {
+
+
+ foundhigh=0;
+
+ double minmass, maxmass;
+ minmass = mn ;
+ maxmass = sqrt(mnsq + Deltasq_high);
+ for(double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
+ {
+ Deltasq_high = mass*mass - mnsq;
+
+ nsols_high = nsols_massless(Deltasq_high);
+ if(nsols_high>0)
+ {
+ foundhigh=1;
+ Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
+ break;
+ }
+ }
+ if(foundhigh==0)
+ {
+
+ 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 mt2::nsols_massless(double Dsq)
+{
+ double delta;
+ delta = Dsq/(2*Easq);
+ d1 = d11*delta;
+ e1 = e11*delta;
+ f1 = f12*delta*delta+f10;
+ d2 = d21*delta+d20;
+ e2 = e21*delta+e20;
+ f2 = f22*delta*delta+f21*delta+f20;
+
+ double a,b;
+ if (pax > 0) a = Ea/Dsq;
+ else a = -Ea/Dsq;
+ if (pax > 0) b = -Dsq/(4*Ea)+mnsq*Ea/Dsq;
+ else b = Dsq/(4*Ea)-mnsq*Ea/Dsq;
+
+ double A4,A3,A2,A1,A0;
+
+ A4 = a*a*a2;
+ A3 = 2*a*b2/Ea;
+ A2 = (2*a*a2*b+c2+2*a*d2)/(Easq);
+ A1 = (2*b*b2+2*e2)/(Easq*Ea);
+ A0 = (a2*b*b+2*b*d2+f2)/(Easq*Easq);
+
+ long double A0sq, A1sq, A2sq, A3sq, A4sq;
+ A0sq = A0*A0;
+ A1sq = A1*A1;
+ A2sq = A2*A2;
+ A3sq = A3*A3;
+ A4sq = A4*A4;
+
+ 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*A3sq/(16*A4));
+ C1 = -(3*A1/4. -A2*A3/(8*A4));
+ C0 = -A0 + A1*A3/(16*A4);
+ long double D1, D0;
+ D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
+ D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
+
+ long double E0;
+ E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
+
+ long double t1,t2,t3,t4,t5;
+
+//find the coefficients for the leading term in the Sturm sequence
+ t1 = A4;
+ t2 = A4;
+ t3 = C2;
+ t4 = D1;
+ t5 = E0;
+
+ int nsol;
+ nsol = signchange_n(t1,t2,t3,t4,t5)-signchange_p(t1,t2,t3,t4,t5);
+ if( nsol < 0 ) nsol=0;
+
+ return nsol;
+
+}
+
+void mt2::mt2_bisect()
+{
+ solved = true;
+
+//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;
+}
+
+int mt2::find_high(double & Deltasq_high)
+{
+ double x0,y0;
+ x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
+ y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
+ double Deltasq_low = (mn + ma)*(mn + ma) - mnsq;
+ do
+ {
+ double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.;
+ int nsols_mid = nsols(Deltasq_mid);
+ if ( nsols_mid == 2 )
+ {
+ Deltasq_high = Deltasq_mid;
+ return 1;
+ }
+ else if (nsols_mid == 4)
+ {
+ Deltasq_high = Deltasq_mid;
+ continue;
+ }
+ else if (nsols_mid ==0)
+ {
+ d1 = -pax*(Deltasq_mid-masq)/(2*Easq);
+ e1 = -pay*(Deltasq_mid-masq)/(2*Easq);
+ d2 = -pmissx + pbx*(Deltasq_mid - mbsq)/(2*Ebsq)
+ + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
+ e2 = -pmissy + pby*(Deltasq_mid - mbsq)/(2*Ebsq)
+ + pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
+ f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq)/(2*Eb)+
+ (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq)/(2*Eb)+
+ (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
+// Does the larger ellipse contain the smaller one?
+ double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
+ if (dis < 0) Deltasq_high = Deltasq_mid;
+ else Deltasq_low = Deltasq_mid;
+ }
+
+ } while ( Deltasq_high - Deltasq_low > 0.001);
+ return 0;
+}
+int mt2::scan_high(double & Deltasq_high)
+{
+ int foundhigh = 0 ;
+ int nsols_high;
+
+
+ double Deltasq_low;
+ 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)
+ {
+ Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
+ foundhigh = 1;
+ break;
+ }
+ }
+ return foundhigh;
+}
+int mt2::nsols( double Dsq)
+{
+ double delta = (Dsq-masq)/(2*Easq);
+
+//calculate coefficients for the two quadratic equations
+ d1 = d11*delta;
+ e1 = e11*delta;
+ f1 = f12*delta*delta+f10;
+ d2 = d21*delta+d20;
+ e2 = e21*delta+e20;
+ f2 = f22*delta*delta+f21*delta+f20;
+
+//obtain the coefficients for the 4th order equation
+//devided by Ea^n to make the variable dimensionless
+ long double A4, A3, A2, A1, A0;
+
+ A4 =
+ -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 +
+ 4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 +
+ a1*a1*c2*c2;
+
+ A3 =
+ (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 +
+ 8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 +
+ 4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 -
+ 4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea;
+
+
+ A2 =
+ (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 -
+ 8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 +
+ 4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 -
+ 8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 +
+ 4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 +
+ 4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq;
+
+ A1 =
+ (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 -
+ 4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 -
+ 4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 -
+ 4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea);
+
+ A0 =
+ (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 +
+ 4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 +
+ a1*a1*f2*f2)/(Easq*Easq);
+
+ long double A0sq, A1sq, A2sq, A3sq, A4sq;
+ A0sq = A0*A0;
+ A1sq = A1*A1;
+ A2sq = A2*A2;
+ A3sq = A3*A3;
+ A4sq = A4*A4;
+
+ 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*A3sq/(16*A4));
+ C1 = -(3*A1/4. -A2*A3/(8*A4));
+ C0 = -A0 + A1*A3/(16*A4);
+
+ long double D1, D0;
+ D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
+ D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;
+
+ long double E0;
+ E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;
+
+ long double t1,t2,t3,t4,t5;
+//find the coefficients for the leading term in the Sturm sequence
+ t1 = A4;
+ t2 = A4;
+ t3 = C2;
+ t4 = D1;
+ t5 = E0;
+
+
+//The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
+ int nsol;
+ nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);
+
+//Cannot have negative number of solutions, must be roundoff effect
+ if (nsol < 0) nsol = 0;
+
+ return nsol;
+
+}
+
+inline int mt2::signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5)
+{
+ int nsc;
+ nsc=0;
+ if(t1*t2>0) nsc++;
+ if(t2*t3>0) nsc++;
+ if(t3*t4>0) nsc++;
+ if(t4*t5>0) nsc++;
+ return nsc;
+}
+inline int mt2::signchange_p( long double t1, long double t2, long double t3, long double t4, long double t5)
+{
+ int nsc;
+ nsc=0;
+ if(t1*t2<0) nsc++;
+ if(t2*t3<0) nsc++;
+ if(t3*t4<0) nsc++;
+ if(t4*t5<0) nsc++;
+ return nsc;
+}
+
+}//end namespace mt2_bisect
+
+}//end namespace Rivet
Added: trunk/src/Tools/mt2_bisect.hh
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Tools/mt2_bisect.hh Wed Apr 13 15:24:21 2011 (r3054)
@@ -0,0 +1,85 @@
+#ifndef MT2_BISECT_H
+#define MT2_BISECT_H
+
+/***********************************************************************/
+/* */
+/* Finding mt2 by Bisection */
+/* */
+/* Authors: Hsin-Chia Cheng, Zhenyu Han */
+/* Dec 11, 2008, v1.01a */
+/* */
+/* see arXiv:0810.5178 */
+/* */
+/* Wrapped for Rivet by D. Grellscheid */
+/* Apr 13, 2011 */
+/* */
+/***********************************************************************/
+
+namespace Rivet {
+
+namespace mt2_bisect
+{
+
+/*The user can change the desired precision below, the larger one of the following two definitions is used. Relative precision less than 0.00001 is not guaranteed to be achievable--use with caution*/
+
+ const double RELATIVE_PRECISION = 0.00001;
+ //defined as precision = RELATIVE_PRECISION * scale, where scale = max{Ea, Eb}
+
+ const double ABSOLUTE_PRECISION = 0.0;
+ //absolute precision for mt2, unused by default
+
+
+//Reserved for expert
+ const double MIN_MASS = 0.1; //if ma<MINMASS and mb<MINMASS, use massless code
+ const double ZERO_MASS = 0.000; //give massless particles a small mass
+ const double SCANSTEP = 0.1;
+
+class mt2
+{
+ public:
+
+ 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();
+ int nevt;
+ private:
+
+ bool solved;
+ bool momenta_set;
+ double mt2_b;
+
+ int nsols(double Dsq);
+ int nsols_massless(double Dsq);
+ inline int signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5);
+ inline 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);
+ //data members
+ double pax, pay, ma, Ea;
+ double pmissx, pmissy;
+ double pbx, pby, mb, Eb;
+ double mn, mn_unscale;
+
+ //auxiliary definitions
+ double masq, Easq;
+ double mbsq, Ebsq;
+ double pmissxsq, pmissysq;
+ double mnsq;
+
+ //auxiliary coefficients
+ double a1, b1, c1, a2, b2, c2, d1, e1, f1, d2, e2, f2;
+ double d11, e11, f12, f10, d21, d20, e21, e20, f22, f21, f20;
+
+ double scale;
+ double precision;
+};
+
+}//end namespace mt2_bisect
+
+}//end namespace Rivet
+
+#endif
More information about the Rivet-svn
mailing list
|