The results are compared and check at the numerical precision levels. Some small discrepancy can appear when the macro is executed on different architectures where it has been calibrated (Power PC G5) The macro is divided in 4 parts:
int ntest = 0;
int nfail = 0;
int ok = 0;
int compare( double v1, double v2, const char* name, double Scale = 1.0) {
ntest = ntest + 1;
double eps = Scale* 2.22044604925031308e-16;
int iret = 0;
double delta = v2 - v1;
double d = 0;
if (delta < 0 ) delta = - delta;
if (v1 == 0 || v2 == 0) {
if (delta > eps ) {
iret = 1;
}
}
else {
d = v1;
if ( v1 < 0) d = -d;
if ( delta/d > eps && delta > eps )
iret = 1;
}
if (iret == 0)
std::cout << ".";
else {
int pr = std::cout.precision (18);
int discr;
if (d != 0)
discr = int(delta/d/eps);
else
discr = int(delta/eps);
std::cout << "\nDiscrepancy in " << name << "() : " << v1 << " != " << v2 << " discr = " << discr
<< " (Allowed discrepancy is " << eps << ")\n";
std::cout.precision (pr);
nfail = nfail + 1;
}
return iret;
}
int testVector3D() {
std::cout << "\n************************************************************************\n "
<< " Vector 3D Test"
<< "\n************************************************************************\n";
std::cout << "Test Cartesian-Polar : " ;
ok = 0;
ok+= compare(v1.X(), v2.X(), "x");
ok+= compare(v1.Y(), v2.Y(), "y");
ok+= compare(v1.Z(), v2.Z(), "z");
ok+= compare(v1.Phi(), v2.Phi(), "phi");
ok+= compare(v1.Theta(), v2.Theta(), "theta");
ok+= compare(v1.R(), v2.R(), "r");
ok+= compare(v1.Eta(), v2.Eta(), "eta");
ok+= compare(v1.Rho(), v2.Rho(), "rho");
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test Cartesian-CylindricalEta : ";
ok = 0;
ok+= compare(v1.X(), v3.X(), "x");
ok+= compare(v1.Y(), v3.Y(), "y");
ok+= compare(v1.Z(), v3.Z(), "z");
ok+= compare(v1.Phi(), v3.Phi(), "phi");
ok+= compare(v1.Theta(), v3.Theta(), "theta");
ok+= compare(v1.R(), v3.R(), "r");
ok+= compare(v1.Eta(), v3.Eta(), "eta");
ok+= compare(v1.Rho(), v3.Rho(), "rho");
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test Cartesian-Cylindrical : ";
ok = 0;
ok+= compare(v1.X(), v4.X(), "x");
ok+= compare(v1.Y(), v4.Y(), "y");
ok+= compare(v1.Z(), v4.Z(), "z");
ok+= compare(v1.Phi(), v4.Phi(), "phi");
ok+= compare(v1.Theta(), v4.Theta(), "theta");
ok+= compare(v1.R(), v4.R(), "r");
ok+= compare(v1.Eta(), v4.Eta(), "eta");
ok+= compare(v1.Rho(), v4.Rho(), "rho");
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test Operations : " ;
ok = 0;
ok+= compare( Dot, v1.Mag2(),"dot" );
ok+= compare( vcross.R(), 0,"cross" );
ok+= compare( v1.R(), vscale2.
R(),
"scale");
ok+= compare(v2.Phi(),vu.Phi(),"unit Phi");
ok+= compare(v2.Theta(),vu.Theta(),"unit Theta");
ok+= compare(1.0,vu.R(),"unit ");
ok+= compare( q4.
X(), q1.X(),
"op X" );
ok+= compare( q4.
Y(), q1.Y(),
"op Y" );
ok+= compare( q4.
Z(), q1.Z(),
"op Z" );
ok+= compare( w1 == v1, static_cast<double>(true), "== XYZ");
ok+= compare( w2 == v2, static_cast<double>(true), "== Polar");
ok+= compare( w3 == v3, static_cast<double>(true), "== RhoEtaPhi");
ok+= compare( w4 == v4, static_cast<double>(true), "== RhoZPhi");
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test Setters : " ;
q2.SetXYZ(q1.X(), q1.Y(), q1.Z() );
ok+= compare( q2.X(), q1.X(), "setXYZ X" );
ok+= compare( q2.Y(), q1.Y(), "setXYZ Y" );
ok+= compare( q2.Z(), q1.Z(), "setXYZ Z" );
q2.SetCoordinates( 2.0*q1.Rho(), q1.Eta(), q1.Phi() );
ok+= compare( q2.X(), q1s.
X(),
"set X" );
ok+= compare( q2.Y(), q1s.
Y(),
"set Y" );
ok+= compare( q2.Z(), q1s.
Z(),
"set Z" );
if (ok == 0) std::cout << "\t\t OK " << std::endl;
std::cout << "Test Linear Algebra conversion: " ;
vxyz1.Coordinates().GetCoordinates(vla1.GetMatrixArray() );
vla2[0] = 1.; vla2[1] = -2.; vla2[2] = 1.;
ok = 0;
double prod1 = vxyz1.Dot(vxyz2);
double prod2 = vla1*vla2;
ok+= compare( prod1, prod2, "la test" );
if (ok == 0) std::cout << "\t\t OK " << std::endl;
return ok;
}
int testPoint3D() {
std::cout << "\n************************************************************************\n "
<< " Point 3D Tests"
<< "\n************************************************************************\n";
std::cout << "Test Cartesian-Polar : ";
ok = 0;
ok+= compare(p1.x(), p2.X(), "x");
ok+= compare(p1.y(), p2.Y(), "y");
ok+= compare(p1.z(), p2.Z(), "z");
ok+= compare(p1.phi(), p2.Phi(), "phi");
ok+= compare(p1.theta(), p2.Theta(), "theta");
ok+= compare(p1.r(), p2.R(), "r");
ok+= compare(p1.eta(), p2.Eta(), "eta");
ok+= compare(p1.rho(), p2.Rho(), "rho");
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test Polar-CylindricalEta : ";
ok = 0;
ok+= compare(p2.X(), p3.X(), "x");
ok+= compare(p2.Y(), p3.Y(), "y");
ok+= compare(p2.Z(), p3.Z(), "z",3);
ok+= compare(p2.Phi(), p3.Phi(), "phi");
ok+= compare(p2.Theta(), p3.Theta(), "theta");
ok+= compare(p2.R(), p3.R(), "r");
ok+= compare(p2.Eta(), p3.Eta(), "eta");
ok+= compare(p2.Rho(), p3.Rho(), "rho");
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test operations : ";
double Dot = p1.Dot(vperp);
ok+= compare( Dot, 0.0,"dot", 10 );
ok+= compare( vcross.
R(), p1.R(),
"cross mag" );
ok+= compare( vcross.
Dot(vperp), 0.0,
"cross dir" );
ok+= compare( p1.R(), pscale2.
R(),
"scale");
ok+= compare( p1 == pscale2, static_cast<double>(true), "== Point");
if (ok == 0) std::cout << "\t OK " << std::endl;
return ok;
}
int testLorentzVector() {
std::cout << "\n************************************************************************\n "
<< " Lorentz Vector Tests"
<< "\n************************************************************************\n";
std::cout << "Test XYZT - PtEtaPhiE Vectors: ";
ok = 0;
ok+= compare(v1.Px(), v2.X(), "x");
ok+= compare(v1.Py(), v2.Y(), "y");
ok+= compare(v1.Pz(), v2.Z(), "z", 2);
ok+= compare(v1.E(), v2.T(), "e");
ok+= compare(v1.Phi(), v2.Phi(), "phi");
ok+= compare(v1.Theta(), v2.Theta(), "theta");
ok+= compare(v1.Pt(), v2.Pt(), "pt");
ok+= compare(v1.M(), v2.M(), "mass", 5);
ok+= compare(v1.Et(), v2.Et(), "et");
ok+= compare(v1.Mt(), v2.Mt(), "mt", 3);
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test XYZT - PtEtaPhiM Vectors: ";
ok = 0;
ok+= compare(v1.Px(), v3.X(), "x");
ok+= compare(v1.Py(), v3.Y(), "y");
ok+= compare(v1.Pz(), v3.Z(), "z", 2);
ok+= compare(v1.E(), v3.T(), "e");
ok+= compare(v1.Phi(), v3.Phi(), "phi");
ok+= compare(v1.Theta(), v3.Theta(), "theta");
ok+= compare(v1.Pt(), v3.Pt(), "pt");
ok+= compare(v1.M(), v3.M(), "mass", 5);
ok+= compare(v1.Et(), v3.Et(), "et");
ok+= compare(v1.Mt(), v3.Mt(), "mt", 3);
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test PtEtaPhiE - PxPyPzM Vect.: ";
ok = 0;
ok+= compare(v4.Px(), v3.X(), "x");
ok+= compare(v4.Py(), v3.Y(), "y");
ok+= compare(v4.Pz(), v3.Z(), "z",2);
ok+= compare(v4.E(), v3.T(), "e");
ok+= compare(v4.Phi(), v3.Phi(), "phi");
ok+= compare(v4.Theta(), v3.Theta(), "theta");
ok+= compare(v4.Pt(), v3.Pt(), "pt");
ok+= compare(v4.M(), v3.M(), "mass",5);
ok+= compare(v4.Et(), v3.Et(), "et");
ok+= compare(v4.Mt(), v3.Mt(), "mt",3);
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test operations : ";
ok = 0;
ok+= compare( Dot, v1.M2(),"dot" , 10 );
ok+= compare( v1.M(), vscale2.
M(),
"scale");
ok+= compare( q4.
x(), q1.X(),
"op X" );
ok+= compare( q4.
y(), q1.Y(),
"op Y" );
ok+= compare( q4.
z(), q1.Z(),
"op Z" );
ok+= compare( q4.
t(), q1.E(),
"op E" );
ok+= compare( w1 == v1, static_cast<double>(true), "== PxPyPzE");
ok+= compare( w2 == v2, static_cast<double>(true), "== PtEtaPhiE");
ok+= compare( w3 == v3, static_cast<double>(true), "== PtEtaPhiM");
ok+= compare( w4 == v4, static_cast<double>(true), "== PxPyPzM");
double gamma = q1.Gamma();
ok += compare( b.
R(),
beta,
"beta" );
ok += compare( gamma, 1./
sqrt( 1 - beta*beta ),
"gamma");
if (ok == 0) std::cout << "\t OK " << std::endl;
std::cout << "Test Setters : " ;
q2.SetXYZT(q1.Px(), q1.Py(), q1.Pz(), q1.E() );
ok+= compare( q2.X(), q1.X(), "setXYZT X" );
ok+= compare( q2.Y(), q1.Y(), "setXYZT Y" );
ok+= compare( q2.Z(), q1.Z(), "setXYZT Z" ,2);
ok+= compare( q2.T(), q1.E(), "setXYZT E" );
q2.SetCoordinates( 2.0*q1.Rho(), q1.Eta(), q1.Phi(), 2.0*q1.E() );
ok+= compare( q2.X(), q1s.
X(),
"set X" );
ok+= compare( q2.Y(), q1s.
Y(),
"set Y" );
ok+= compare( q2.Z(), q1s.
Z(),
"set Z" ,2);
ok+= compare( q2.T(), q1s.
T(),
"set E" );
if (ok == 0) std::cout << "\t OK " << std::endl;
return ok;
}
int testVectorUtil() {
std::cout << "\n************************************************************************\n "
<< " Utility Function Tests"
<< "\n************************************************************************\n";
std::cout << "Test Vector utility functions : ";
ok = 0;
v2 = v2cyl;
if (ok == 0) std::cout << "\t\t OK " << std::endl;
std::cout << "Test Point utility functions : ";
ok = 0;
p2 = p2cyl;
if (ok == 0) std::cout << "\t\t OK " << std::endl;
std::cout << "LorentzVector utility funct.: ";
ok = 0;
if (ok == 0) std::cout << "\t\t OK " << std::endl;
return ok;
}
int testRotation() {
std::cout << "\n************************************************************************\n "
<< " Rotation and Transformation Tests"
<< "\n************************************************************************\n";
std::cout << "Test Vector Rotations : ";
ok = 0;
ok+= compare(v1.
X(), v2.
X(),
"x",2);
ok+= compare(v1.
Y(), v2.
Y(),
"y",2);
ok+= compare(v1.
Z(), v2.
Z(),
"z",2);
ok+= compare(v1.
X(), v3.
X(),
"x",2);
ok+= compare(v1.
Y(), v3.
Y(),
"y",2);
ok+= compare(v1.
Z(), v3.
Z(),
"z",2);
ok+= compare(v1.
X(), v4.
X(),
"x",5);
ok+= compare(v1.
Y(), v4.
Y(),
"y",5);
ok+= compare(v1.
Z(), v4.
Z(),
"z",5);
ok+= compare(v1.
X(), v5.
X(),
"x",2);
ok+= compare(v1.
Y(), v5.
Y(),
"y",2);
ok+= compare(v1.
Z(), v5.
Z(),
"z",2);
double rdata[9];
r2.GetComponents(rdata, rdata+9);
double vdata[3];
v.GetCoordinates(vdata);
ok+= compare(v1.
X(), v6.
X(),
"x");
ok+= compare(v1.
Y(), v6.
Y(),
"y");
ok+= compare(v1.
Z(), v6.
Z(),
"z");
if (ok == 0) std::cout << "\t OK " << std::endl;
else std::cout << std::endl;
std::cout << "Test Axial Rotations : ";
ok = 0;
RotationZYX rzyx( rz.Angle(), ry.Angle(), rx.Angle() );
ok+= compare(vrot1.X(), vrot2.
X(),
"x");
ok+= compare(vrot1.Y(), vrot2.
Y(),
"y");
ok+= compare(vrot1.Z(), vrot2.
Z(),
"z");
vrot2 = qx * qy * qz *
v;
ok+= compare(vrot1.X(), vrot2.X(), "x",2);
ok+= compare(vrot1.Y(), vrot2.Y(), "y",2);
ok+= compare(vrot1.Z(), vrot2.Z(), "z",2);
ok+= compare(vrot1.X(), vrot2.X(), "x");
ok+= compare(vrot1.Y(), vrot2.Y(), "y");
ok+= compare(vrot1.Z(), vrot2.Z(), "z");
vrot1 = rz * ry * rx *
v;
vrot2 = r3z * r3y * r3x *
v;
ok+= compare(vrot1.X(), vrot2.X(), "x");
ok+= compare(vrot1.Y(), vrot2.Y(), "y");
ok+= compare(vrot1.Z(), vrot2.Z(), "z");
XYZPoint vinv1 = rx.Inverse()*ry.Inverse()*rz.Inverse()*vrot1;
ok+= compare(vinv1.X(), v.X(), "x",2);
ok+= compare(vinv1.Y(), v.Y(), "y");
ok+= compare(vinv1.Z(), v.Z(), "z");
if (ok == 0) std::cout << "\t OK " << std::endl;
else std::cout << std::endl;
std::cout << "Test Rotations by a PI angle : ";
ok = 0;
double b[4] = { 6,8,10,3.14159265358979323 };
ok+= compare(arPi.Axis().X(), a1.
Axis().
X(),
"x");
ok+= compare(arPi.Axis().Y(), a1.
Axis().
Y(),
"y");
ok+= compare(arPi.Axis().Z(), a1.
Axis().
Z(),
"z");
ok+= compare(arPi.Angle(), a1.
Angle(),
"angle");
ok+= compare(ePi.Phi(), e1.
Phi(),
"phi");
ok+= compare(ePi.Theta(), e1.
Theta(),
"theta");
ok+= compare(ePi.Psi(), e1.
Psi(),
"ps1");
if (ok == 0) std::cout << "\t\t OK " << std::endl;
else std::cout << std::endl;
std::cout << "Test Inversions : ";
ok = 0;
ok+= compare(p.
X(), v.X(),
"x",10);
ok+= compare(p.
Y(), v.Y(),
"y",10);
ok+= compare(p.
Z(), v.Z(),
"z",10);
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
ok+= compare(p.X(), v.X(), "x",1E9);
ok+= compare(p.Y(), v.Y(), "y",1E9);
ok+= compare(p.Z(), v.Z(), "z",1E9);
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
if (ok == 0) std::cout << "\t OK " << std::endl;
else std::cout << std::endl;
std::cout << "Test rectify : ";
ok = 0;
XYZVector u1(0.999498,-0.00118212,-0.0316611);
ok+= compare(v.R(), vrr.
R(),
"R",1.E9);
if (ok == 0) std::cout << "\t\t OK " << std::endl;
else std::cout << std::endl;
std::cout << "Test Transform3D : ";
ok = 0;
ok+= compare(pd.
X(), vd.
X(),
"x");
ok+= compare(pd.
Y(), vd.
Y(),
"y");
ok+= compare(pd.
Z(), vd.
Z(),
"z");
double tdata[12];
t.GetComponents(tdata);
double vData[4];
v.GetCoordinates(vData);
vData[3] = 1;
ok+= compare(pd.
X(), qt(0),
"x");
ok+= compare(pd.
Y(), qt(1),
"y");
ok+= compare(pd.
Z(), qt(2),
"z");
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
Transform3D tinv2 ( r2.Inverse(), r2.Inverse() *( -d) ) ;
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
ok+= compare( tc == tb*ta, static_cast<double>(true), "== Rot*Tra");
ok+= compare( td == ta*tb, static_cast<double>(true), "== Rot*Tra");
if (ok == 0) std::cout << "\t OK " << std::endl;
else std::cout << std::endl;
std::cout << "Test Plane3D : ";
ok = 0;
ok+= compare(n.
Dot(p2-p1), 0.0,
"n.v12",10);
ok+= compare(n.
Dot(p3-p1), 0.0,
"n.v13",10);
ok+= compare(n.
Dot(p3-p2), 0.0,
"n.v23",10);
ok+= compare(n1.
X(), n2.
X(),
"a",10);
ok+= compare(n1.
Y(), n2.
Y(),
"b",10);
ok+= compare(n1.
Z(), n2.
Z(),
"c",10);
ok+= compare(plane1.HesseDistance(), plane2.HesseDistance(), "d",10);
ok += compare(plane1.Distance(pt1), 0.0, "distance",10);
if (ok == 0) std::cout << "\t OK " << std::endl;
else std::cout << std::endl;
std::cout << "Test LorentzRotation : ";
ok = 0;
ok+= compare(lv1== lv2,true,"V0==V2");
ok+= compare(lv1== lv2,true,"V1==V2");
double rlData[16];
double lvData[4];
lv.GetCoordinates(lvData);
ok+= compare(lv1.
X(), qlr(0),
"x");
ok+= compare(lv1.
Y(), qlr(1),
"y");
ok+= compare(lv1.
Z(), qlr(2),
"z");
ok+= compare(lv1.
E(), qlr(3),
"t");
ok+= compare(lv0.X(), lv.X(), "x");
ok+= compare(lv0.Y(), lv.Y(), "y");
ok+= compare(lv0.Z(), lv.Z(), "z");
ok+= compare(lv0.E(), lv.E(), "t");
if (ok == 0) std::cout << "\t OK " << std::endl;
else std::cout << std::endl;
std::cout << "Test Boost : ";
ok = 0;
ok+= compare(lvb.
X(), lvb2.
X(),
"x");
ok+= compare(lvb.
Y(), lvb2.
Y(),
"y");
ok+= compare(lvb.
Z(), lvb2.
Z(),
"z");
ok+= compare(lvb.
E(), lvb2.
E(),
"t");
ok+= compare(lvb.
M(), lv.M(),
"m",50);
lv0 = bst.Inverse() * lvb;
ok+= compare(lv0.X(), lv.X(), "x",5);
ok+= compare(lv0.Y(), lv.Y(), "y",5);
ok+= compare(lv0.Z(), lv.Z(), "z",3);
ok+= compare(lv0.E(), lv.E(), "t",3);
bst.SetComponents( brest.X(), brest.Y(), brest.Z() );
ok+= compare(lvr.
X(), 0.0,
"x",10);
ok+= compare(lvr.
Y(), 0.0,
"y",10);
ok+= compare(lvr.
Z(), 0.0,
"z",10);
ok+= compare(lvr.
M(), lv.M(),
"m",10);
if (ok == 0) std::cout << "\t OK " << std::endl;
else std::cout << std::endl;
return ok;
}
void mathcoreGenVector() {
testVector3D();
testPoint3D();
testLorentzVector();
testVectorUtil();
testRotation();
std::cout << "\n\nNumber of tests " << ntest << " failed = " << nfail << std::endl;
}