56 int compare(
double v1,
double v2,
const char*
name,
double Scale = 1.0) {
60 double eps = Scale* 2.22044604925031308e-16;
62 double delta = v2 - v1;
64 if (delta < 0 ) delta = - delta;
65 if (v1 == 0 || v2 == 0) {
76 if ( delta/d > eps && delta > eps )
83 int pr = std::cout.precision (18);
86 discr = int(delta/d/eps);
88 discr = int(delta/eps);
90 std::cout <<
"\nDiscrepancy in " << name <<
"() : " << v1 <<
" != " << v2 <<
" discr = " << discr
91 <<
" (Allowed discrepancy is " << eps <<
")\n";
92 std::cout.precision (pr);
99 std::cout <<
"\n************************************************************************\n " 101 <<
"\n************************************************************************\n";
107 std::cout <<
"Test Cartesian-Polar : " ;
112 ok+=
compare(v1.X(), v2.X(),
"x");
113 ok+=
compare(v1.Y(), v2.Y(),
"y");
114 ok+=
compare(v1.Z(), v2.Z(),
"z");
115 ok+=
compare(v1.Phi(), v2.Phi(),
"phi");
116 ok+=
compare(v1.Theta(), v2.Theta(),
"theta");
117 ok+=
compare(v1.R(), v2.R(),
"r");
118 ok+=
compare(v1.Eta(), v2.Eta(),
"eta");
119 ok+=
compare(v1.Rho(), v2.Rho(),
"rho");
121 if (ok == 0) std::cout <<
"\t OK " << std::endl;
123 std::cout <<
"Test Cartesian-CylindricalEta : ";
128 ok+=
compare(v1.X(), v3.X(),
"x");
129 ok+=
compare(v1.Y(), v3.Y(),
"y");
130 ok+=
compare(v1.Z(), v3.Z(),
"z");
131 ok+=
compare(v1.Phi(), v3.Phi(),
"phi");
132 ok+=
compare(v1.Theta(), v3.Theta(),
"theta");
133 ok+=
compare(v1.R(), v3.R(),
"r");
134 ok+=
compare(v1.Eta(), v3.Eta(),
"eta");
135 ok+=
compare(v1.Rho(), v3.Rho(),
"rho");
137 if (ok == 0) std::cout <<
"\t OK " << std::endl;
139 std::cout <<
"Test Cartesian-Cylindrical : ";
144 ok+=
compare(v1.X(), v4.X(),
"x");
145 ok+=
compare(v1.Y(), v4.Y(),
"y");
146 ok+=
compare(v1.Z(), v4.Z(),
"z");
147 ok+=
compare(v1.Phi(), v4.Phi(),
"phi");
148 ok+=
compare(v1.Theta(), v4.Theta(),
"theta");
149 ok+=
compare(v1.R(), v4.R(),
"r");
150 ok+=
compare(v1.Eta(), v4.Eta(),
"eta");
151 ok+=
compare(v1.Rho(), v4.Rho(),
"rho");
153 if (ok == 0) std::cout <<
"\t OK " << std::endl;
155 std::cout <<
"Test Operations : " ;
158 double Dot = v1.Dot(v2);
159 ok+=
compare( Dot, v1.Mag2(),
"dot" );
161 ok+=
compare( vcross.R(), 0,
"cross" );
165 ok+=
compare( v1.R(), vscale2.
R(),
"scale");
168 ok+=
compare(v2.Phi(),vu.Phi(),
"unit Phi");
169 ok+=
compare(v2.Theta(),vu.Theta(),
"unit Theta");
170 ok+=
compare(1.0,vu.R(),
"unit ");
178 ok+=
compare( q4.
X(), q1.X(),
"op X" );
179 ok+=
compare( q4.
Y(), q1.Y(),
"op Y" );
180 ok+=
compare( q4.
Z(), q1.Z(),
"op Z" );
187 ok+=
compare( w1 == v1, static_cast<double>(
true),
"== XYZ");
188 ok+=
compare( w2 == v2, static_cast<double>(
true),
"== Polar");
189 ok+=
compare( w3 == v3, static_cast<double>(
true),
"== RhoEtaPhi");
190 ok+=
compare( w4 == v4, static_cast<double>(
true),
"== RhoZPhi");
192 if (ok == 0) std::cout <<
"\t OK " << std::endl;
195 std::cout <<
"Test Setters : " ;
197 q2.SetXYZ(q1.X(), q1.Y(), q1.Z() );
199 ok+=
compare( q2.X(), q1.X(),
"setXYZ X" );
200 ok+=
compare( q2.Y(), q1.Y(),
"setXYZ Y" );
201 ok+=
compare( q2.Z(), q1.Z(),
"setXYZ Z" );
203 q2.SetCoordinates( 2.0*q1.Rho(), q1.Eta(), q1.Phi() );
205 ok+=
compare( q2.X(), q1s.
X(),
"set X" );
206 ok+=
compare( q2.Y(), q1s.
Y(),
"set Y" );
207 ok+=
compare( q2.Z(), q1s.
Z(),
"set Z" );
209 if (ok == 0) std::cout <<
"\t\t OK " << std::endl;
211 std::cout <<
"Test Linear Algebra conversion: " ;
216 vxyz1.Coordinates().GetCoordinates(vla1.GetMatrixArray() );
219 vla2[0] = 1.; vla2[1] = -2.; vla2[2] = 1.;
225 double prod1 = vxyz1.
Dot(vxyz2);
226 double prod2 = vla1*vla2;
227 ok+=
compare( prod1, prod2,
"la test" );
229 if (ok == 0) std::cout <<
"\t\t OK " << std::endl;
235 std::cout <<
"\n************************************************************************\n " 237 <<
"\n************************************************************************\n";
241 std::cout <<
"Test Cartesian-Polar : ";
255 if (ok == 0) std::cout <<
"\t OK " << std::endl;
257 std::cout <<
"Test Polar-CylindricalEta : ";
271 if (ok == 0) std::cout <<
"\t OK " << std::endl;
273 std::cout <<
"Test operations : ";
277 double Dot =
p1.Dot(vperp);
278 ok+=
compare( Dot, 0.0,
"dot", 10 );
281 ok+=
compare( vcross.
R(),
p1.R(),
"cross mag" );
282 ok+=
compare( vcross.
Dot(vperp), 0.0,
"cross dir" );
289 ok+=
compare(
p1 == pscale2, static_cast<double>(
true),
"== Point");
297 if (ok == 0) std::cout <<
"\t OK " << std::endl;
302 int testLorentzVector() {
303 std::cout <<
"\n************************************************************************\n " 304 <<
" Lorentz Vector Tests" 305 <<
"\n************************************************************************\n";
309 std::cout <<
"Test XYZT - PtEtaPhiE Vectors: ";
314 ok+=
compare(v1.Px(), v2.X(),
"x");
315 ok+=
compare(v1.Py(), v2.Y(),
"y");
316 ok+=
compare(v1.Pz(), v2.Z(),
"z", 2);
317 ok+=
compare(v1.E(), v2.T(),
"e");
318 ok+=
compare(v1.Phi(), v2.Phi(),
"phi");
319 ok+=
compare(v1.Theta(), v2.Theta(),
"theta");
320 ok+=
compare(v1.Pt(), v2.Pt(),
"pt");
321 ok+=
compare(v1.M(), v2.M(),
"mass", 5);
322 ok+=
compare(v1.Et(), v2.Et(),
"et");
323 ok+=
compare(v1.Mt(), v2.Mt(),
"mt", 3);
325 if (ok == 0) std::cout <<
"\t OK " << std::endl;
327 std::cout <<
"Test XYZT - PtEtaPhiM Vectors: ";
332 ok+=
compare(v1.Px(), v3.X(),
"x");
333 ok+=
compare(v1.Py(), v3.Y(),
"y");
334 ok+=
compare(v1.Pz(), v3.Z(),
"z", 2);
335 ok+=
compare(v1.E(), v3.T(),
"e");
336 ok+=
compare(v1.Phi(), v3.Phi(),
"phi");
337 ok+=
compare(v1.Theta(), v3.Theta(),
"theta");
338 ok+=
compare(v1.Pt(), v3.Pt(),
"pt");
339 ok+=
compare(v1.M(), v3.M(),
"mass", 5);
340 ok+=
compare(v1.Et(), v3.Et(),
"et");
341 ok+=
compare(v1.Mt(), v3.Mt(),
"mt", 3);
343 if (ok == 0) std::cout <<
"\t OK " << std::endl;
345 std::cout <<
"Test PtEtaPhiE - PxPyPzM Vect.: ";
350 ok+=
compare(v4.Px(), v3.X(),
"x");
351 ok+=
compare(v4.Py(), v3.Y(),
"y");
352 ok+=
compare(v4.Pz(), v3.Z(),
"z",2);
353 ok+=
compare(v4.E(), v3.T(),
"e");
354 ok+=
compare(v4.Phi(), v3.Phi(),
"phi");
355 ok+=
compare(v4.Theta(), v3.Theta(),
"theta");
356 ok+=
compare(v4.Pt(), v3.Pt(),
"pt");
357 ok+=
compare(v4.M(), v3.M(),
"mass",5);
358 ok+=
compare(v4.Et(), v3.Et(),
"et");
359 ok+=
compare(v4.Mt(), v3.Mt(),
"mt",3);
361 if (ok == 0) std::cout <<
"\t OK " << std::endl;
363 std::cout <<
"Test operations : ";
366 double Dot = v1.Dot(v2);
367 ok+=
compare( Dot, v1.M2(),
"dot" , 10 );
371 ok+=
compare( v1.M(), vscale2.
M(),
"scale");
379 ok+=
compare( q4.
x(), q1.X(),
"op X" );
380 ok+=
compare( q4.
y(), q1.Y(),
"op Y" );
381 ok+=
compare( q4.
z(), q1.Z(),
"op Z" );
382 ok+=
compare( q4.
t(), q1.E(),
"op E" );
389 ok+=
compare( w1 == v1, static_cast<double>(
true),
"== PxPyPzE");
390 ok+=
compare( w2 == v2, static_cast<double>(
true),
"== PtEtaPhiE");
391 ok+=
compare( w3 == v3, static_cast<double>(
true),
"== PtEtaPhiM");
392 ok+=
compare( w4 == v4, static_cast<double>(
true),
"== PxPyPzM");
396 double beta = q1.Beta();
397 double gamma = q1.Gamma();
400 ok +=
compare( gamma, 1./
sqrt( 1 - beta*beta ),
"gamma");
402 if (ok == 0) std::cout <<
"\t OK " << std::endl;
405 std::cout <<
"Test Setters : " ;
407 q2.SetXYZT(q1.Px(), q1.Py(), q1.Pz(), q1.E() );
409 ok+=
compare( q2.X(), q1.X(),
"setXYZT X" );
410 ok+=
compare( q2.Y(), q1.Y(),
"setXYZT Y" );
411 ok+=
compare( q2.Z(), q1.Z(),
"setXYZT Z" ,2);
412 ok+=
compare( q2.T(), q1.E(),
"setXYZT E" );
414 q2.SetCoordinates( 2.0*q1.Rho(), q1.Eta(), q1.Phi(), 2.0*q1.E() );
416 ok+=
compare( q2.X(), q1s.
X(),
"set X" );
417 ok+=
compare( q2.Y(), q1s.
Y(),
"set Y" );
418 ok+=
compare( q2.Z(), q1s.
Z(),
"set Z" ,2);
419 ok+=
compare( q2.T(), q1s.
T(),
"set E" );
421 if (ok == 0) std::cout <<
"\t OK " << std::endl;
427 std::cout <<
"\n************************************************************************\n " 428 <<
" Utility Function Tests" 429 <<
"\n************************************************************************\n";
431 std::cout <<
"Test Vector utility functions : ";
450 if (ok == 0) std::cout <<
"\t\t OK " << std::endl;
453 std::cout <<
"Test Point utility functions : ";
471 if (ok == 0) std::cout <<
"\t\t OK " << std::endl;
473 std::cout <<
"LorentzVector utility funct.: ";
486 if (ok == 0) std::cout <<
"\t\t OK " << std::endl;
492 std::cout <<
"\n************************************************************************\n " 493 <<
" Rotation and Transformation Tests" 494 <<
"\n************************************************************************\n";
496 std::cout <<
"Test Vector Rotations : ";
534 r2.GetComponents(rdata, rdata+9);
537 v.GetCoordinates(vdata);
548 if (ok == 0) std::cout <<
"\t OK " << std::endl;
549 else std::cout << std::endl;
551 std::cout <<
"Test Axial Rotations : ";
566 RotationZYX rzyx( rz.Angle(), ry.Angle(), rx.Angle() );
571 ok+=
compare(vrot1.X(), vrot2.
X(),
"x");
572 ok+=
compare(vrot1.Y(), vrot2.
Y(),
"y");
573 ok+=
compare(vrot1.Z(), vrot2.
Z(),
"z");
575 vrot2 = qx * qy * qz *
v;
577 ok+=
compare(vrot1.X(), vrot2.X(),
"x",2);
578 ok+=
compare(vrot1.Y(), vrot2.Y(),
"y",2);
579 ok+=
compare(vrot1.Z(), vrot2.Z(),
"z",2);
583 ok+=
compare(vrot1.X(), vrot2.X(),
"x");
584 ok+=
compare(vrot1.Y(), vrot2.Y(),
"y");
585 ok+=
compare(vrot1.Z(), vrot2.Z(),
"z");
588 vrot1 = rz * ry * rx *
v;
589 vrot2 = r3z * r3y * r3x *
v;
591 ok+=
compare(vrot1.X(), vrot2.X(),
"x");
592 ok+=
compare(vrot1.Y(), vrot2.Y(),
"y");
593 ok+=
compare(vrot1.Z(), vrot2.Z(),
"z");
595 XYZPoint vinv1 = rx.Inverse()*ry.Inverse()*rz.Inverse()*vrot1;
597 ok+=
compare(vinv1.X(), v.X(),
"x",2);
598 ok+=
compare(vinv1.Y(), v.Y(),
"y");
599 ok+=
compare(vinv1.Z(), v.Z(),
"z");
601 if (ok == 0) std::cout <<
"\t OK " << std::endl;
602 else std::cout << std::endl;
605 std::cout <<
"Test Rotations by a PI angle : ";
608 double b[4] = { 6,8,10,3.14159265358979323 };
623 if (ok == 0) std::cout <<
"\t\t OK " << std::endl;
624 else std::cout << std::endl;
626 std::cout <<
"Test Inversions : ";
645 ok+=
compare(p.X(), v.X(),
"x",10);
646 ok+=
compare(p.Y(), v.Y(),
"y",10);
647 ok+=
compare(p.Z(), v.Z(),
"z",10);
651 ok+=
compare(p.X(), v.X(),
"x",1E9);
652 ok+=
compare(p.Y(), v.Y(),
"y",1E9);
653 ok+=
compare(p.Z(), v.Z(),
"z",1E9);
657 ok+=
compare(p.X(), v.X(),
"x",10);
658 ok+=
compare(p.Y(), v.Y(),
"y",10);
659 ok+=
compare(p.Z(), v.Z(),
"z",10);
667 ok+=
compare(p.X(), v.X(),
"x",10);
668 ok+=
compare(p.Y(), v.Y(),
"y",10);
669 ok+=
compare(p.Z(), v.Z(),
"z",10);
671 if (ok == 0) std::cout <<
"\t OK " << std::endl;
672 else std::cout << std::endl;
675 std::cout <<
"Test rectify : ";
678 XYZVector u1(0.999498,-0.00118212,-0.0316611);
680 XYZVector u3(0.0316832,0.0372921,0.998802);
684 ok+=
compare(v.R(), vrr.
R(),
"R",1.E9);
686 if (ok == 0) std::cout <<
"\t\t OK " << std::endl;
687 else std::cout << std::endl;
689 std::cout <<
"Test Transform3D : ";
705 t.GetComponents(tdata);
708 v.GetCoordinates(vData);
724 ok+=
compare(p.X(), v.X(),
"x",10);
725 ok+=
compare(p.Y(), v.Y(),
"y",10);
726 ok+=
compare(p.Z(), v.Z(),
"z",10);
733 ok+=
compare(p.X(), v.X(),
"x",10);
734 ok+=
compare(p.Y(), v.Y(),
"y",10);
735 ok+=
compare(p.Z(), v.Z(),
"z",10);
743 ok+=
compare( tc == tb*ta, static_cast<double>(
true),
"== Rot*Tra");
744 ok+=
compare( td == ta*tb, static_cast<double>(
true),
"== Rot*Tra");
747 if (ok == 0) std::cout <<
"\t OK " << std::endl;
748 else std::cout << std::endl;
750 std::cout <<
"Test Plane3D : ";
780 ok+=
compare(plane1.HesseDistance(), plane2.HesseDistance(),
"d",10);
783 ok +=
compare(plane1.Distance(pt1), 0.0,
"distance",10);
785 if (ok == 0) std::cout <<
"\t OK " << std::endl;
786 else std::cout << std::endl;
788 std::cout <<
"Test LorentzRotation : ";
810 ok+=
compare(lv1== lv2,
true,
"V0==V2");
811 ok+=
compare(lv1== lv2,
true,
"V1==V2");
817 lv.GetCoordinates(lvData);
830 ok+=
compare(lv0.X(), lv.X(),
"x");
831 ok+=
compare(lv0.Y(), lv.Y(),
"y");
832 ok+=
compare(lv0.Z(), lv.Z(),
"z");
833 ok+=
compare(lv0.E(), lv.E(),
"t");
835 if (ok == 0) std::cout <<
"\t OK " << std::endl;
836 else std::cout << std::endl;
839 std::cout <<
"Test Boost : ";
842 Boost bst( 0.3,0.4,0.5);
855 ok+=
compare(lvb.
M(), lv.M(),
"m",50);
858 lv0 = bst.Inverse() * lvb;
860 ok+=
compare(lv0.X(), lv.X(),
"x",5);
861 ok+=
compare(lv0.Y(), lv.Y(),
"y",5);
862 ok+=
compare(lv0.Z(), lv.Z(),
"z",3);
863 ok+=
compare(lv0.E(), lv.E(),
"t",3);
866 bst.SetComponents( brest.X(), brest.Y(), brest.Z() );
873 ok+=
compare(lvr.
M(), lv.M(),
"m",10);
875 if (ok == 0) std::cout <<
"\t OK " << std::endl;
876 else std::cout << std::endl;
880 void mathcoreGenVector() {
889 std::cout <<
"\n\nNumber of tests " << ntest <<
" failed = " << nfail << std::endl;
T Dot(const SVector< T, D > &lhs, const SVector< T, D > &rhs)
Vector dot product.
Scalar R() const
Polar R, converting if necessary from internal coordinate system.
Class describing a geometrical plane in 3 dimensions.
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
DisplacementVector3D< CoordSystem, Tag > & SetCoordinates(const Scalar src[])
Set internal data based on a C-style array of 3 Scalar numbers.
Scalar R() const
Polar R, converting if necessary from internal coordinate system.
static double p3(double t, double a, double b, double c, double d)
Rotation class representing a 3D rotation about the Z axis by the angle of rotation.
int compare(double v1, double v2, const std::string &name="", double scale=1.0)
Class describing a generic position vector (point) in 3 dimensions.
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Rotation class with the (3D) rotation represented by angles describing first a rotation of an angle p...
DisplacementVector3D Unit() const
return unit vector parallel to this
Rotation3D Inverse() const
Return inverse of a rotation.
Rotation class with the (3D) rotation represented by a unit quaternion (u, i, j, k).
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
double beta(double x, double y)
Calculates the beta function.
Scalar Phi() const
Return Phi Euler angle.
AxisAngle class describing rotation represented with direction axis (3D Vector) and an angle of rotat...
Scalar X() const
Cartesian X, converting if necessary from internal coordinate system.
LorentzRotation Inverse() const
Return inverse of a rotation.
Scalar Psi() const
Return Psi Euler angle.
Quaternion Inverse() const
Return inverse of a rotation.
Rotation class representing a 3D rotation about the Y axis by the angle of rotation.
unsigned int r3[N_CITIES]
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...
Scalar Z() const
Cartesian Z, converting if necessary from internal coordinate system.
Scalar Y() const
Cartesian Y, converting if necessary from internal coordinate system.
Lorentz boost class with the (4D) transformation represented internally by a 4x4 orthosymplectic matr...
Scalar Theta() const
Return Theta Euler angle.
static double p2(double t, double a, double b, double c)
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...
PositionVector3D< CoordSystem, Tag > & SetCoordinates(const Scalar src[])
Set internal data based on a C-style array of 3 Scalar numbers.
Scalar E() const
return 4-th component (time, or energy for a 4-momentum vector)
AxisVector Axis() const
accesss to rotation axis
Element * GetMatrixArray()
Lorentz transformation class with the (4D) transformation represented by a 4x4 orthosymplectic matrix...
Scalar Z() const
Cartesian Z, converting if necessary from internal coordinate system.
Vector1::Scalar DeltaPhi(const Vector1 &v1, const Vector2 &v2)
Find aximutal Angle difference between two generic vectors ( v2.Phi() - v1.Phi() ) The only requireme...
Scalar Angle() const
access to rotation angle
R__EXTERN TSystem * gSystem
double Angle(const Vector1 &v1, const Vector2 &v2)
Find Angle between two vectors.
Scalar X() const
Cartesian X, converting if necessary from internal coordinate system.
unsigned int r1[N_CITIES]
Rotation class representing a 3D rotation about the X axis by the angle of rotation.
DisplacementVector3D Cross(const DisplacementVector3D< OtherCoords, Tag > &v) const
Return vector (cross) product of two displacement vectors, as a vector in the coordinate system of th...
Rotation class with the (3D) rotation represented by a 3x3 orthogonal matrix.
static double p1(double t, double a, double b)
void GetComponents(Foreign4Vector &v1, Foreign4Vector &v2, Foreign4Vector &v3, Foreign4Vector &v4) const
Get components into four 4-vectors which will be the (orthosymplectic) columns of the rotation matrix...
EulerAngles class describing rotation as three angles (Euler Angles).
Scalar M() const
return magnitude (mass) using the (-,-,-,+) metric.
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZVector
3D Vector based on the cartesian coordinates x,y,z in double precision
AxisAngle Inverse() const
Return inverse of an AxisAngle rotation.
Scalar Y() const
Cartesian Y, converting if necessary from internal coordinate system.
RotationZYX Inverse() const
Return inverse of a rotation.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Scalar Dot(const DisplacementVector3D< OtherCoords, Tag > &v) const
Return the scalar (Dot) product of this with a displacement vector in any coordinate system...
double CosTheta(const Vector1 &v1, const Vector2 &v2)
Find CosTheta Angle between two generic 3D vectors pre-requisite: vectors implement the X()...
unsigned int r2[N_CITIES]
Scalar Dot(const DisplacementVector3D< OtherCoords, Tag > &v) const
Return the scalar (dot) product of two displacement vectors.