Logo ROOT  

Reference Guide
mathcoreGenVector.C File Reference
Tutorials » Math tutorials

Detailed Description

View in nbviewer Open in SWAN Example macro testing available methods and operation of the GenVector classes.

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:

  • testVector3D : tests of the 3D Vector classes
  • testPoint3D : tests of the 3D Point classes
  • testLorentzVector : tests of the 4D LorentzVector classes
  • testVectorUtil : tests of the utility functions of all the vector classes

To execute the macro type in:

root[0] .x mathcoreGenVector.C
************************************************************************
Vector 3D Test
************************************************************************
Test Cartesian-Polar : ........ OK
Test Cartesian-CylindricalEta : ........ OK
Test Cartesian-Cylindrical : ........ OK
Test Operations : ............. OK
Test Setters : ...... OK
Test Linear Algebra conversion: . OK
************************************************************************
Point 3D Tests
************************************************************************
Test Cartesian-Polar : ........ OK
Test Polar-CylindricalEta : ........ OK
Test operations : ..... OK
************************************************************************
Lorentz Vector Tests
************************************************************************
Test XYZT - PtEtaPhiE Vectors: .......... OK
Test XYZT - PtEtaPhiM Vectors: .......... OK
Test PtEtaPhiE - PxPyPzM Vect.: .......... OK
Test operations : ............ OK
Test Setters : ........ OK
************************************************************************
Utility Function Tests
************************************************************************
Test Vector utility functions : .... OK
Test Point utility functions : .... OK
LorentzVector utility funct.: ... OK
************************************************************************
Rotation and Transformation Tests
************************************************************************
Test Vector Rotations : ............... OK
Test Axial Rotations : ............... OK
Test Rotations by a PI angle : ....... OK
Test Inversions : ............... OK
Test rectify : . OK
Test Transform3D : .............. OK
Test Plane3D : ........ OK
Test LorentzRotation : .......... OK
Test Boost : ............. OK
Number of tests 224 failed = 0
#include "TMatrixD.h"
#include "TVectorD.h"
#include "TMath.h"
#include "Math/Point3D.h"
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/GenVector/Rotation3D.h"
#include "Math/GenVector/EulerAngles.h"
#include "Math/GenVector/AxisAngle.h"
#include "Math/GenVector/Quaternion.h"
#include "Math/GenVector/RotationX.h"
#include "Math/GenVector/RotationY.h"
#include "Math/GenVector/RotationZ.h"
#include "Math/GenVector/RotationZYX.h"
#include "Math/GenVector/LorentzRotation.h"
#include "Math/GenVector/Boost.h"
#include "Math/GenVector/BoostX.h"
#include "Math/GenVector/BoostY.h"
#include "Math/GenVector/BoostZ.h"
#include "Math/GenVector/Transform3D.h"
#include "Math/GenVector/Plane3D.h"
#include "Math/GenVector/VectorUtil.h"
using namespace ROOT::Math;
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;
// numerical double limit for epsilon
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;
}
}
// skip case v1 or v2 is infinity
else {
d = v1;
if ( v1 < 0) d = -d;
// add also case when delta is small by default
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";
XYZVector v1(0.01, 0.02, 16);
std::cout << "Test Cartesian-Polar : " ;
Polar3DVector v2(v1.R(), v1.Theta(), v1.Phi() );
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 : ";
RhoEtaPhiVector v3( v1.Rho(), v1.Eta(), v1.Phi() );
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 : ";
RhoZPhiVector v4( v1.Rho(), v1.Z(), v1.Phi() );
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;
double Dot = v1.Dot(v2);
ok+= compare( Dot, v1.Mag2(),"dot" );
XYZVector vcross = v1.Cross(v2);
ok+= compare( vcross.R(), 0,"cross" );
XYZVector vscale1 = v1*10;
XYZVector vscale2 = vscale1/10;
ok+= compare( v1.R(), vscale2.R(), "scale");
XYZVector vu = v1.Unit();
ok+= compare(v2.Phi(),vu.Phi(),"unit Phi");
ok+= compare(v2.Theta(),vu.Theta(),"unit Theta");
ok+= compare(1.0,vu.R(),"unit ");
XYZVector q1 = v1;
RhoEtaPhiVector q2(1.0,1.0,1.0);
XYZVector q3 = q1 + q2;
XYZVector q4 = q3 - q2;
ok+= compare( q4.X(), q1.X(), "op X" );
ok+= compare( q4.Y(), q1.Y(), "op Y" );
ok+= compare( q4.Z(), q1.Z(), "op Z" );
// test operator ==
XYZVector w1 = v1;
Polar3DVector w2 = v2;
RhoEtaPhiVector w3 = v3;
RhoZPhiVector w4 = v4;
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;
//test setters
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() );
XYZVector q1s = 2.0*q1;
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: " ;
XYZVector vxyz1(1.,2.,3.);
TVectorD vla1(3);
vxyz1.Coordinates().GetCoordinates(vla1.GetMatrixArray() );
TVectorD vla2(3);
vla2[0] = 1.; vla2[1] = -2.; vla2[2] = 1.;
XYZVector vxyz2;
vxyz2.SetCoordinates(&vla2[0]);
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";
XYZPoint p1(1.0, 2.0, 3.0);
std::cout << "Test Cartesian-Polar : ";
Polar3DPoint p2(p1.R(), p1.Theta(), p1.Phi() );
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 : ";
RhoEtaPhiPoint p3( p2.Rho(), p2.Eta(), p2.Phi() );
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 : ";
//std::cout << "\nTest Dot and Cross products with Vectors : ";
Polar3DVector vperp(1.,p1.Theta() + TMath::PiOver2(),p1.Phi() );
double Dot = p1.Dot(vperp);
ok+= compare( Dot, 0.0,"dot", 10 );
XYZPoint vcross = p1.Cross(vperp);
ok+= compare( vcross.R(), p1.R(),"cross mag" );
ok+= compare( vcross.Dot(vperp), 0.0,"cross dir" );
XYZPoint pscale1 = 10*p1;
XYZPoint pscale2 = pscale1/10;
ok+= compare( p1.R(), pscale2.R(), "scale");
// test operator ==
ok+= compare( p1 == pscale2, static_cast<double>(true), "== Point");
//RhoEtaPhiPoint q1 = p1; ! constructor yet not working in CINT
RhoEtaPhiPoint q1; q1 = p1;
q1.SetCoordinates(p1.Rho(),2.0, p1.Phi() );
Polar3DVector v2(p1.R(), p1.Theta(),p1.Phi());
if (ok == 0) std::cout << "\t OK " << std::endl;
return ok;
}
int testLorentzVector() {
std::cout << "\n************************************************************************\n "
<< " Lorentz Vector Tests"
<< "\n************************************************************************\n";
XYZTVector v1(1.0, 2.0, 3.0, 4.0);
std::cout << "Test XYZT - PtEtaPhiE Vectors: ";
PtEtaPhiEVector v2( v1.Rho(), v1.Eta(), v1.Phi(), v1.E() );
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: ";
PtEtaPhiMVector v3( v1.Rho(), v1.Eta(), v1.Phi(), v1.M() );
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.: ";
PxPyPzMVector v4( v3.X(), v3.Y(), v3.Z(), v3.M() );
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;
double Dot = v1.Dot(v2);
ok+= compare( Dot, v1.M2(),"dot" , 10 );
XYZTVector vscale1 = v1*10;
XYZTVector vscale2 = vscale1/10;
ok+= compare( v1.M(), vscale2.M(), "scale");
XYZTVector q1 = v1;
PtEtaPhiEVector q2(1.0,1.0,1.0,5.0);
XYZTVector q3 = q1 + q2;
XYZTVector q4 = q3 - q2;
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" );
// test operator ==
XYZTVector w1 = v1;
PtEtaPhiEVector w2 = v2;
PtEtaPhiMVector w3 = v3;
PxPyPzMVector w4 = v4;
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");
// test gamma beta and boost
XYZVector b = q1.BoostToCM();
double beta = q1.Beta();
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;
//test setters
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() );
XYZTVector q1s = q1*2.0;
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 : ";
XYZVector v1(1.0, 2.0, 3.0);
Polar3DVector v2pol(v1.R(), v1.Theta()+TMath::PiOver2(), v1.Phi() + 1.0);
// mixedmethods not yet impl.
XYZVector v2; v2 = v2pol;
ok = 0;
ok += compare( VectorUtil::DeltaPhi(v1,v2), 1.0, "deltaPhi Vec");
RhoEtaPhiVector v2cyl(v1.Rho(), v1.Eta() + 1.0, v1.Phi() + 1.0);
v2 = v2cyl;
ok += compare( VectorUtil::DeltaR(v1,v2), sqrt(2.0), "DeltaR Vec");
XYZVector vperp = v1.Cross(v2);
ok += compare( VectorUtil::CosTheta(v1,vperp), 0.0, "costheta Vec");
ok += compare( VectorUtil::Angle(v1,vperp), TMath::PiOver2(), "angle Vec");
if (ok == 0) std::cout << "\t\t OK " << std::endl;
std::cout << "Test Point utility functions : ";
XYZPoint p1(1.0, 2.0, 3.0);
Polar3DPoint p2pol(p1.R(), p1.Theta()+TMath::PiOver2(), p1.Phi() + 1.0);
// mixedmethods not yet impl.
XYZPoint p2; p2 = p2pol;
ok = 0;
ok += compare( VectorUtil::DeltaPhi(p1,p2), 1.0, "deltaPhi Point");
RhoEtaPhiPoint p2cyl(p1.Rho(), p1.Eta() + 1.0, p1.Phi() + 1.0);
p2 = p2cyl;
ok += compare( VectorUtil::DeltaR(p1,p2), sqrt(2.0), "DeltaR Point");
XYZPoint pperp(vperp.X(), vperp.Y(), vperp.Z());
ok += compare( VectorUtil::CosTheta(p1,pperp), 0.0, "costheta Point");
ok += compare( VectorUtil::Angle(p1,pperp), TMath::PiOver2(), "angle Point");
if (ok == 0) std::cout << "\t\t OK " << std::endl;
std::cout << "LorentzVector utility funct.: ";
XYZTVector q1(1.0, 2.0, 3.0,4.0);
PtEtaPhiEVector q2cyl(q1.Pt(), q1.Eta()+1.0, q1.Phi() + 1.0, q1.E() );
XYZTVector q2; q2 = q2cyl;
ok = 0;
ok += compare( VectorUtil::DeltaPhi(q1,q2), 1.0, "deltaPhi LVec");
ok += compare( VectorUtil::DeltaR(q1,q2), sqrt(2.0), "DeltaR LVec");
XYZTVector qsum = q1+q2;
ok += compare( VectorUtil::InvariantMass(q1,q2), qsum.M(), "InvMass");
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;
XYZPoint v(1.,2,3.);
double pi = TMath::Pi();
// initiate rotation with some non -trivial angles to test all matrix
EulerAngles r1( pi/2.,pi/4., pi/3 );
Rotation3D r2(r1);
// only operator= is in CINT for the other rotations
Quaternion r3; r3 = r2;
AxisAngle r4; r4 = r3;
RotationZYX r5; r5 = r2;
XYZPoint v1 = r1 * v;
XYZPoint v2 = r2 * v;
XYZPoint v3 = r3 * v;
XYZPoint v4 = r4 * v;
XYZPoint v5 = r5 * v;
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);
// test with matrix
double rdata[9];
r2.GetComponents(rdata, rdata+9);
TMatrixD m(3,3,rdata);
double vdata[3];
v.GetCoordinates(vdata);
TVectorD q(3,vdata);
TVectorD q2 = m*q;
XYZPoint v6;
v6.SetCoordinates( q2.GetMatrixArray() );
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;
RotationX rx( pi/3);
RotationY ry( pi/4);
RotationZ rz( 4*pi/5);
Rotation3D r3x(rx);
Rotation3D r3y(ry);
Rotation3D r3z(rz);
Quaternion qx; qx = rx;
Quaternion qy; qy = ry;
Quaternion qz; qz = rz;
RotationZYX rzyx( rz.Angle(), ry.Angle(), rx.Angle() );
XYZPoint vrot1 = rx * ry * rz * v;
XYZPoint vrot2 = r3x * r3y * r3z * v;
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);
vrot2 = rzyx * v;
ok+= compare(vrot1.X(), vrot2.X(), "x");
ok+= compare(vrot1.Y(), vrot2.Y(), "y");
ok+= compare(vrot1.Z(), vrot2.Z(), "z");
// now inverse (first x then y then 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 };
AxisAngle arPi(b,b+4 );
Rotation3D rPi(arPi);
AxisAngle a1; a1 = rPi;
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");
EulerAngles ePi; ePi=rPi;
EulerAngles e1; e1=Rotation3D(a1);
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;
EulerAngles s1 = r1.Inverse();
Rotation3D s2 = r2.Inverse();
Quaternion s3 = r3.Inverse();
AxisAngle s4 = r4.Inverse();
RotationZYX s5 = r5.Inverse();
// Euler angles not yet impl.
XYZPoint p = s2 * r2 * v;
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
p = s3 * r3 * v;
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
p = s4 * r4 * v;
// axis angle inversion not very precise
ok+= compare(p.X(), v.X(), "x",1E9);
ok+= compare(p.Y(), v.Y(), "y",1E9);
ok+= compare(p.Z(), v.Z(), "z",1E9);
p = s5 * r5 * v;
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
Rotation3D r6(r5);
Rotation3D s6 = r6.Inverse();
p = s6 * r6 * v;
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;
// test Rectify
std::cout << "Test rectify : ";
ok = 0;
XYZVector u1(0.999498,-0.00118212,-0.0316611);
XYZVector u2(0,0.999304,-0.0373108);
XYZVector u3(0.0316832,0.0372921,0.998802);
Rotation3D rr(u1,u2,u3);
// check orto-normality
XYZPoint vrr = rr* v;
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;
XYZVector d(1.,-2.,3.);
Transform3D t(r2,d);
XYZPoint pd = t * v;
// apply directly rotation
XYZPoint vd = r2 * v + d;
ok+= compare(pd.X(), vd.X(), "x");
ok+= compare(pd.Y(), vd.Y(), "y");
ok+= compare(pd.Z(), vd.Z(), "z");
// test with matrix
double tdata[12];
t.GetComponents(tdata);
TMatrixD mt(3,4,tdata);
double vData[4]; // needs a vector of dim 4
v.GetCoordinates(vData);
vData[3] = 1;
TVectorD q0(4,vData);
TVectorD qt = mt*q0;
ok+= compare(pd.X(), qt(0), "x");
ok+= compare(pd.Y(), qt(1), "y");
ok+= compare(pd.Z(), qt(2), "z");
// test inverse
Transform3D tinv = t.Inverse();
p = tinv * t * v;
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
// test construct inverse from translation first
Transform3D tinv2 ( r2.Inverse(), r2.Inverse() *( -d) ) ;
p = tinv2 * t * v;
ok+= compare(p.X(), v.X(), "x",10);
ok+= compare(p.Y(), v.Y(), "y",10);
ok+= compare(p.Z(), v.Z(), "z",10);
// test from only rotation and only translation
Transform3D ta( EulerAngles(1.,2.,3.) );
Transform3D tb( XYZVector(1,2,3) );
Transform3D tc( Rotation3D(EulerAngles(1.,2.,3.)) , XYZVector(1,2,3) );
Transform3D td( ta.Rotation(), ta.Rotation() * XYZVector(1,2,3) ) ;
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;
// test transfrom a 3D plane
XYZPoint p1(1,2,3);
XYZPoint p2(-2,-1,4);
XYZPoint p3(-1,3,2);
Plane3D plane(p1,p2,p3);
XYZVector n = plane.Normal();
// normal is perpendicular to vectors on the planes obtained from subtracting the points
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);
Plane3D plane1 = t(plane);
// transform the points
XYZPoint pt1 = t(p1);
XYZPoint pt2 = t(p2);
XYZPoint pt3 = t(p3);
Plane3D plane2(pt1,pt2,pt3);
XYZVector n1 = plane1.Normal();
XYZVector n2 = plane2.Normal();
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);
// check distances
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;
XYZTVector lv(1.,2.,3.,4.);
// test from rotx (using boosts and 3D rotations not yet impl.)
// rx,ry and rz already defined
Rotation3D r3d = rx*ry*rz;
LorentzRotation rlx(rx);
LorentzRotation rly(ry);
LorentzRotation rlz(rz);
LorentzRotation rl0 = rlx*rly*rlz;
LorentzRotation rl1( r3d);
XYZTVector lv0 = rl0 * lv;
XYZTVector lv1 = rl1 * lv;
XYZTVector lv2 = r3d * lv;
ok+= compare(lv1== lv2,true,"V0==V2");
ok+= compare(lv1== lv2,true,"V1==V2");
double rlData[16];
rl0.GetComponents(rlData);
TMatrixD ml(4,4,rlData);
double lvData[4];
lv.GetCoordinates(lvData);
TVectorD ql(4,lvData);
TVectorD qlr = ml*ql;
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");
// test inverse
lv0 = rl0 * rl0.Inverse() * lv;
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;
// test Boosts
std::cout << "Test Boost : ";
ok = 0;
Boost bst( 0.3,0.4,0.5); // boost (must be <= 1)
XYZTVector lvb = bst ( lv );
LorentzRotation rl2 (bst);
XYZTVector lvb2 = rl2 (lv);
// test with lorentz rotation
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); // m must stay constant
// test inverse
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);
XYZVector brest = lv.BoostToCM();
bst.SetComponents( brest.X(), brest.Y(), brest.Z() );
XYZTVector lvr = bst * lv;
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;
}
d
#define d(i)
Definition: RSha256.hxx:102
b
#define b(i)
Definition: RSha256.hxx:100
s1
#define s1(x)
Definition: RSha256.hxx:91
name
char name[80]
Definition: TGX11.cxx:109
q
float * q
Definition: THbookFile.cxx:87
sqrt
double sqrt(double)
ROOT::Math::DisplacementVector3D::SetCoordinates
DisplacementVector3D< CoordSystem, Tag > & SetCoordinates(const Scalar src[])
Set internal data based on a C-style array of 3 Scalar numbers.
Definition: DisplacementVector3D.h:198
TVectorT::GetMatrixArray
Element * GetMatrixArray()
Definition: TVectorT.h:78
ROOT::Math::beta
double beta(double x, double y)
Calculates the beta function.
Definition: SpecFuncMathCore.cxx:111
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::Cephes::gamma
double gamma(double x)
Definition: SpecFuncCephes.cxx:339
ROOT::Math::VectorUtil::DeltaR
Vector1::Scalar DeltaR(const Vector1 &v1, const Vector2 &v2)
Find difference in pseudorapidity (Eta) and Phi betwen two generic vectors The only requirements on t...
Definition: VectorUtil.h:95
ROOT::Math::VectorUtil::DeltaPhi
Vector1::Scalar DeltaPhi(const Vector1 &v1, const Vector2 &v2)
Find aximutal Angle difference between two generic vectors ( v2.Phi() - v1.Phi() ) The only requireme...
Definition: VectorUtil.h:59
ROOT::Math::VectorUtil::Angle
double Angle(const Vector1 &v1, const Vector2 &v2)
Find Angle between two vectors.
Definition: VectorUtil.h:138
ROOT::Math::VectorUtil::CosTheta
double CosTheta(const Vector1 &v1, const Vector2 &v2)
Find CosTheta Angle between two generic 3D vectors pre-requisite: vectors implement the X(),...
Definition: VectorUtil.h:112
ROOT::Math::VectorUtil::InvariantMass
Vector1::Scalar InvariantMass(const Vector1 &v1, const Vector2 &v2)
return the invariant mass of two LorentzVector The only requirement on the LorentzVector is that they...
Definition: VectorUtil.h:215
ROOT::Math::RhoZPhiVector
DisplacementVector3D< Cylindrical3D< double >, DefaultCoordinateSystemTag > RhoZPhiVector
3D Vector based on the cylindrical coordinates rho, z, phi in double precision.
Definition: Vector3Dfwd.h:70
ROOT::Math::PtEtaPhiEVector
LorentzVector< PtEtaPhiE4D< double > > PtEtaPhiEVector
LorentzVector based on the cylindrical coordinates Pt, eta, phi and E (rho, eta, phi,...
Definition: Vector4Dfwd.h:61
ROOT::Math::PtEtaPhiMVector
LorentzVector< PtEtaPhiM4D< double > > PtEtaPhiMVector
LorentzVector based on the cylindrical coordinates pt, eta, phi and Mass in double precision.
Definition: Vector4Dfwd.h:66
ROOT::Math::Plane3D
Impl::Plane3D< double > Plane3D
Definition: Plane3D.h:302
ROOT::Math::RhoEtaPhiPoint
PositionVector3D< CylindricalEta3D< double >, DefaultCoordinateSystemTag > RhoEtaPhiPoint
3D Point based on the eta based cylindrical coordinates rho, eta, phi in double precision.
Definition: Point3Dfwd.h:49
ROOT::Math::XYZPoint
PositionVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZPoint
3D Point based on the cartesian coordinates x,y,z in double precision
Definition: Point3Dfwd.h:33
ROOT::Math::Transform3D
Impl::Transform3D< double > Transform3D
Definition: Transform3D.h:1312
ROOT::Math::XYZVector
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZVector
3D Vector based on the cartesian coordinates x,y,z in double precision
Definition: Vector3Dfwd.h:34
ROOT::Math::Polar3DVector
DisplacementVector3D< Polar3D< double >, DefaultCoordinateSystemTag > Polar3DVector
3D Vector based on the polar coordinates rho, theta, phi in double precision.
Definition: Vector3Dfwd.h:60
ROOT::Math::RhoEtaPhiVector
DisplacementVector3D< CylindricalEta3D< double >, DefaultCoordinateSystemTag > RhoEtaPhiVector
3D Vector based on the eta based cylindrical coordinates rho, eta, phi in double precision.
Definition: Vector3Dfwd.h:50
ROOT::Math::PxPyPzMVector
LorentzVector< PxPyPzM4D< double > > PxPyPzMVector
LorentzVector based on the x, y, z, and Mass in double precision.
Definition: Vector4Dfwd.h:56
ROOT::Math::XYZTVector
LorentzVector< PxPyPzE4D< double > > XYZTVector
LorentzVector based on x,y,x,t (or px,py,pz,E) coordinates in double precision with metric (-,...
Definition: Vector4Dfwd.h:33
ROOT::Math::Polar3DPoint
PositionVector3D< Polar3D< double >, DefaultCoordinateSystemTag > Polar3DPoint
3D Point based on the polar coordinates rho, theta, phi in double precision.
Definition: Point3Dfwd.h:59
TGeant4Unit::pi
static constexpr double pi
Definition: TGeant4SystemOfUnits.h:67
TMath::PiOver2
constexpr Double_t PiOver2()
Definition: TMath.h:52
TMath::Pi
constexpr Double_t Pi()
Definition: TMath.h:38
Dot
#define Dot(u, v)
Definition: normal.c:49
v2
@ v2
Definition: rootcling_impl.cxx:3627
v
@ v
Definition: rootcling_impl.cxx:3624
v4
@ v4
Definition: rootcling_impl.cxx:3629
v3
@ v3
Definition: rootcling_impl.cxx:3628
v1
@ v1
Definition: rootcling_impl.cxx:3626
lv
auto * lv
Definition: textalign.C:5
m
auto * m
Definition: textangle.C:8
Author
Lorenzo Moneta

Definition in file mathcoreGenVector.C.