Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
mathcoreGenVector.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_math
3 /// \notebook -nodraw
4 /// Example macro testing available methods and operation of the GenVector classes.
5 /// The results are compared and check at the
6 /// numerical precision levels.
7 /// Some small discrepancy can appear when the macro
8 /// is executed on different architectures where it has been calibrated (Power PC G5)
9 /// The macro is divided in 4 parts:
10 /// - testVector3D : tests of the 3D Vector classes
11 /// - testPoint3D : tests of the 3D Point classes
12 /// - testLorentzVector : tests of the 4D LorentzVector classes
13 /// - testVectorUtil : tests of the utility functions of all the vector classes
14 ///
15 /// To execute the macro type in:
16 ///
17 /// ~~~{.cpp}
18 /// root[0] .x mathcoreGenVector.C
19 /// ~~~
20 ///
21 /// \macro_output
22 /// \macro_code
23 ///
24 /// \author Lorenzo Moneta
25 
26 #include "TMatrixD.h"
27 #include "TVectorD.h"
28 #include "TMath.h"
29 
30 #include "Math/Point3D.h"
31 #include "Math/Vector3D.h"
32 #include "Math/Vector4D.h"
33 #include "Math/GenVector/Rotation3D.h"
34 #include "Math/GenVector/EulerAngles.h"
35 #include "Math/GenVector/AxisAngle.h"
36 #include "Math/GenVector/Quaternion.h"
37 #include "Math/GenVector/RotationX.h"
38 #include "Math/GenVector/RotationY.h"
39 #include "Math/GenVector/RotationZ.h"
40 #include "Math/GenVector/RotationZYX.h"
41 #include "Math/GenVector/LorentzRotation.h"
42 #include "Math/GenVector/Boost.h"
43 #include "Math/GenVector/BoostX.h"
44 #include "Math/GenVector/BoostY.h"
45 #include "Math/GenVector/BoostZ.h"
46 #include "Math/GenVector/Transform3D.h"
47 #include "Math/GenVector/Plane3D.h"
48 #include "Math/GenVector/VectorUtil.h"
49 
50 using namespace ROOT::Math;
51 
52 int ntest = 0;
53 int nfail = 0;
54 int ok = 0;
55 
56 int compare( double v1, double v2, const char* name, double Scale = 1.0) {
57  ntest = ntest + 1;
58 
59  // numerical double limit for epsilon
60  double eps = Scale* 2.22044604925031308e-16;
61  int iret = 0;
62  double delta = v2 - v1;
63  double d = 0;
64  if (delta < 0 ) delta = - delta;
65  if (v1 == 0 || v2 == 0) {
66  if (delta > eps ) {
67  iret = 1;
68  }
69  }
70  // skip case v1 or v2 is infinity
71  else {
72  d = v1;
73 
74  if ( v1 < 0) d = -d;
75  // add also case when delta is small by default
76  if ( delta/d > eps && delta > eps )
77  iret = 1;
78  }
79 
80  if (iret == 0)
81  std::cout << ".";
82  else {
83  int pr = std::cout.precision (18);
84  int discr;
85  if (d != 0)
86  discr = int(delta/d/eps);
87  else
88  discr = int(delta/eps);
89 
90  std::cout << "\nDiscrepancy in " << name << "() : " << v1 << " != " << v2 << " discr = " << discr
91  << " (Allowed discrepancy is " << eps << ")\n";
92  std::cout.precision (pr);
93  nfail = nfail + 1;
94  }
95  return iret;
96 }
97 
98 int testVector3D() {
99  std::cout << "\n************************************************************************\n "
100  << " Vector 3D Test"
101  << "\n************************************************************************\n";
102 
103  XYZVector v1(0.01, 0.02, 16);
104 
105  std::cout << "Test Cartesian-Polar : " ;
106 
107  Polar3DVector v2(v1.R(), v1.Theta(), v1.Phi() );
108 
109  ok = 0;
110  ok+= compare(v1.X(), v2.X(), "x");
111  ok+= compare(v1.Y(), v2.Y(), "y");
112  ok+= compare(v1.Z(), v2.Z(), "z");
113  ok+= compare(v1.Phi(), v2.Phi(), "phi");
114  ok+= compare(v1.Theta(), v2.Theta(), "theta");
115  ok+= compare(v1.R(), v2.R(), "r");
116  ok+= compare(v1.Eta(), v2.Eta(), "eta");
117  ok+= compare(v1.Rho(), v2.Rho(), "rho");
118 
119  if (ok == 0) std::cout << "\t OK " << std::endl;
120 
121  std::cout << "Test Cartesian-CylindricalEta : ";
122 
123  RhoEtaPhiVector v3( v1.Rho(), v1.Eta(), v1.Phi() );
124 
125  ok = 0;
126  ok+= compare(v1.X(), v3.X(), "x");
127  ok+= compare(v1.Y(), v3.Y(), "y");
128  ok+= compare(v1.Z(), v3.Z(), "z");
129  ok+= compare(v1.Phi(), v3.Phi(), "phi");
130  ok+= compare(v1.Theta(), v3.Theta(), "theta");
131  ok+= compare(v1.R(), v3.R(), "r");
132  ok+= compare(v1.Eta(), v3.Eta(), "eta");
133  ok+= compare(v1.Rho(), v3.Rho(), "rho");
134 
135  if (ok == 0) std::cout << "\t OK " << std::endl;
136 
137  std::cout << "Test Cartesian-Cylindrical : ";
138 
139  RhoZPhiVector v4( v1.Rho(), v1.Z(), v1.Phi() );
140 
141  ok = 0;
142  ok+= compare(v1.X(), v4.X(), "x");
143  ok+= compare(v1.Y(), v4.Y(), "y");
144  ok+= compare(v1.Z(), v4.Z(), "z");
145  ok+= compare(v1.Phi(), v4.Phi(), "phi");
146  ok+= compare(v1.Theta(), v4.Theta(), "theta");
147  ok+= compare(v1.R(), v4.R(), "r");
148  ok+= compare(v1.Eta(), v4.Eta(), "eta");
149  ok+= compare(v1.Rho(), v4.Rho(), "rho");
150 
151  if (ok == 0) std::cout << "\t OK " << std::endl;
152 
153  std::cout << "Test Operations : " ;
154 
155  ok = 0;
156  double Dot = v1.Dot(v2);
157  ok+= compare( Dot, v1.Mag2(),"dot" );
158  XYZVector vcross = v1.Cross(v2);
159  ok+= compare( vcross.R(), 0,"cross" );
160 
161  XYZVector vscale1 = v1*10;
162  XYZVector vscale2 = vscale1/10;
163  ok+= compare( v1.R(), vscale2.R(), "scale");
164 
165  XYZVector vu = v1.Unit();
166  ok+= compare(v2.Phi(),vu.Phi(),"unit Phi");
167  ok+= compare(v2.Theta(),vu.Theta(),"unit Theta");
168  ok+= compare(1.0,vu.R(),"unit ");
169 
170  XYZVector q1 = v1;
171  RhoEtaPhiVector q2(1.0,1.0,1.0);
172 
173  XYZVector q3 = q1 + q2;
174  XYZVector q4 = q3 - q2;
175 
176  ok+= compare( q4.X(), q1.X(), "op X" );
177  ok+= compare( q4.Y(), q1.Y(), "op Y" );
178  ok+= compare( q4.Z(), q1.Z(), "op Z" );
179 
180  // test operator ==
181  XYZVector w1 = v1;
182  Polar3DVector w2 = v2;
183  RhoEtaPhiVector w3 = v3;
184  RhoZPhiVector w4 = v4;
185  ok+= compare( w1 == v1, static_cast<double>(true), "== XYZ");
186  ok+= compare( w2 == v2, static_cast<double>(true), "== Polar");
187  ok+= compare( w3 == v3, static_cast<double>(true), "== RhoEtaPhi");
188  ok+= compare( w4 == v4, static_cast<double>(true), "== RhoZPhi");
189 
190  if (ok == 0) std::cout << "\t OK " << std::endl;
191 
192  //test setters
193  std::cout << "Test Setters : " ;
194 
195  q2.SetXYZ(q1.X(), q1.Y(), q1.Z() );
196 
197  ok+= compare( q2.X(), q1.X(), "setXYZ X" );
198  ok+= compare( q2.Y(), q1.Y(), "setXYZ Y" );
199  ok+= compare( q2.Z(), q1.Z(), "setXYZ Z" );
200 
201  q2.SetCoordinates( 2.0*q1.Rho(), q1.Eta(), q1.Phi() );
202  XYZVector q1s = 2.0*q1;
203  ok+= compare( q2.X(), q1s.X(), "set X" );
204  ok+= compare( q2.Y(), q1s.Y(), "set Y" );
205  ok+= compare( q2.Z(), q1s.Z(), "set Z" );
206 
207  if (ok == 0) std::cout << "\t\t OK " << std::endl;
208 
209  std::cout << "Test Linear Algebra conversion: " ;
210 
211  XYZVector vxyz1(1.,2.,3.);
212 
213  TVectorD vla1(3);
214  vxyz1.Coordinates().GetCoordinates(vla1.GetMatrixArray() );
215 
216  TVectorD vla2(3);
217  vla2[0] = 1.; vla2[1] = -2.; vla2[2] = 1.;
218 
219  XYZVector vxyz2;
220  vxyz2.SetCoordinates(&vla2[0]);
221 
222  ok = 0;
223  double prod1 = vxyz1.Dot(vxyz2);
224  double prod2 = vla1*vla2;
225  ok+= compare( prod1, prod2, "la test" );
226 
227  if (ok == 0) std::cout << "\t\t OK " << std::endl;
228 
229  return ok;
230 }
231 
232 int testPoint3D() {
233  std::cout << "\n************************************************************************\n "
234  << " Point 3D Tests"
235  << "\n************************************************************************\n";
236 
237  XYZPoint p1(1.0, 2.0, 3.0);
238 
239  std::cout << "Test Cartesian-Polar : ";
240 
241  Polar3DPoint p2(p1.R(), p1.Theta(), p1.Phi() );
242 
243  ok = 0;
244  ok+= compare(p1.x(), p2.X(), "x");
245  ok+= compare(p1.y(), p2.Y(), "y");
246  ok+= compare(p1.z(), p2.Z(), "z");
247  ok+= compare(p1.phi(), p2.Phi(), "phi");
248  ok+= compare(p1.theta(), p2.Theta(), "theta");
249  ok+= compare(p1.r(), p2.R(), "r");
250  ok+= compare(p1.eta(), p2.Eta(), "eta");
251  ok+= compare(p1.rho(), p2.Rho(), "rho");
252 
253  if (ok == 0) std::cout << "\t OK " << std::endl;
254 
255  std::cout << "Test Polar-CylindricalEta : ";
256 
257  RhoEtaPhiPoint p3( p2.Rho(), p2.Eta(), p2.Phi() );
258 
259  ok = 0;
260  ok+= compare(p2.X(), p3.X(), "x");
261  ok+= compare(p2.Y(), p3.Y(), "y");
262  ok+= compare(p2.Z(), p3.Z(), "z",3);
263  ok+= compare(p2.Phi(), p3.Phi(), "phi");
264  ok+= compare(p2.Theta(), p3.Theta(), "theta");
265  ok+= compare(p2.R(), p3.R(), "r");
266  ok+= compare(p2.Eta(), p3.Eta(), "eta");
267  ok+= compare(p2.Rho(), p3.Rho(), "rho");
268 
269  if (ok == 0) std::cout << "\t OK " << std::endl;
270 
271  std::cout << "Test operations : ";
272 
273  //std::cout << "\nTest Dot and Cross products with Vectors : ";
274  Polar3DVector vperp(1.,p1.Theta() + TMath::PiOver2(),p1.Phi() );
275  double Dot = p1.Dot(vperp);
276  ok+= compare( Dot, 0.0,"dot", 10 );
277 
278  XYZPoint vcross = p1.Cross(vperp);
279  ok+= compare( vcross.R(), p1.R(),"cross mag" );
280  ok+= compare( vcross.Dot(vperp), 0.0,"cross dir" );
281 
282  XYZPoint pscale1 = 10*p1;
283  XYZPoint pscale2 = pscale1/10;
284  ok+= compare( p1.R(), pscale2.R(), "scale");
285 
286  // test operator ==
287  ok+= compare( p1 == pscale2, static_cast<double>(true), "== Point");
288 
289  //RhoEtaPhiPoint q1 = p1; ! constructor yet not working in CINT
290  RhoEtaPhiPoint q1; q1 = p1;
291  q1.SetCoordinates(p1.Rho(),2.0, p1.Phi() );
292 
293  Polar3DVector v2(p1.R(), p1.Theta(),p1.Phi());
294 
295  if (ok == 0) std::cout << "\t OK " << std::endl;
296 
297  return ok;
298 }
299 
300 int testLorentzVector() {
301  std::cout << "\n************************************************************************\n "
302  << " Lorentz Vector Tests"
303  << "\n************************************************************************\n";
304 
305  XYZTVector v1(1.0, 2.0, 3.0, 4.0);
306 
307  std::cout << "Test XYZT - PtEtaPhiE Vectors: ";
308 
309  PtEtaPhiEVector v2( v1.Rho(), v1.Eta(), v1.Phi(), v1.E() );
310 
311  ok = 0;
312  ok+= compare(v1.Px(), v2.X(), "x");
313  ok+= compare(v1.Py(), v2.Y(), "y");
314  ok+= compare(v1.Pz(), v2.Z(), "z", 2);
315  ok+= compare(v1.E(), v2.T(), "e");
316  ok+= compare(v1.Phi(), v2.Phi(), "phi");
317  ok+= compare(v1.Theta(), v2.Theta(), "theta");
318  ok+= compare(v1.Pt(), v2.Pt(), "pt");
319  ok+= compare(v1.M(), v2.M(), "mass", 5);
320  ok+= compare(v1.Et(), v2.Et(), "et");
321  ok+= compare(v1.Mt(), v2.Mt(), "mt", 3);
322 
323  if (ok == 0) std::cout << "\t OK " << std::endl;
324 
325  std::cout << "Test XYZT - PtEtaPhiM Vectors: ";
326 
327  PtEtaPhiMVector v3( v1.Rho(), v1.Eta(), v1.Phi(), v1.M() );
328 
329  ok = 0;
330  ok+= compare(v1.Px(), v3.X(), "x");
331  ok+= compare(v1.Py(), v3.Y(), "y");
332  ok+= compare(v1.Pz(), v3.Z(), "z", 2);
333  ok+= compare(v1.E(), v3.T(), "e");
334  ok+= compare(v1.Phi(), v3.Phi(), "phi");
335  ok+= compare(v1.Theta(), v3.Theta(), "theta");
336  ok+= compare(v1.Pt(), v3.Pt(), "pt");
337  ok+= compare(v1.M(), v3.M(), "mass", 5);
338  ok+= compare(v1.Et(), v3.Et(), "et");
339  ok+= compare(v1.Mt(), v3.Mt(), "mt", 3);
340 
341  if (ok == 0) std::cout << "\t OK " << std::endl;
342 
343  std::cout << "Test PtEtaPhiE - PxPyPzM Vect.: ";
344 
345  PxPyPzMVector v4( v3.X(), v3.Y(), v3.Z(), v3.M() );
346 
347  ok = 0;
348  ok+= compare(v4.Px(), v3.X(), "x");
349  ok+= compare(v4.Py(), v3.Y(), "y");
350  ok+= compare(v4.Pz(), v3.Z(), "z",2);
351  ok+= compare(v4.E(), v3.T(), "e");
352  ok+= compare(v4.Phi(), v3.Phi(), "phi");
353  ok+= compare(v4.Theta(), v3.Theta(), "theta");
354  ok+= compare(v4.Pt(), v3.Pt(), "pt");
355  ok+= compare(v4.M(), v3.M(), "mass",5);
356  ok+= compare(v4.Et(), v3.Et(), "et");
357  ok+= compare(v4.Mt(), v3.Mt(), "mt",3);
358 
359  if (ok == 0) std::cout << "\t OK " << std::endl;
360 
361  std::cout << "Test operations : ";
362 
363  ok = 0;
364  double Dot = v1.Dot(v2);
365  ok+= compare( Dot, v1.M2(),"dot" , 10 );
366 
367  XYZTVector vscale1 = v1*10;
368  XYZTVector vscale2 = vscale1/10;
369  ok+= compare( v1.M(), vscale2.M(), "scale");
370 
371  XYZTVector q1 = v1;
372  PtEtaPhiEVector q2(1.0,1.0,1.0,5.0);
373 
374  XYZTVector q3 = q1 + q2;
375  XYZTVector q4 = q3 - q2;
376 
377  ok+= compare( q4.x(), q1.X(), "op X" );
378  ok+= compare( q4.y(), q1.Y(), "op Y" );
379  ok+= compare( q4.z(), q1.Z(), "op Z" );
380  ok+= compare( q4.t(), q1.E(), "op E" );
381 
382  // test operator ==
383  XYZTVector w1 = v1;
384  PtEtaPhiEVector w2 = v2;
385  PtEtaPhiMVector w3 = v3;
386  PxPyPzMVector w4 = v4;
387  ok+= compare( w1 == v1, static_cast<double>(true), "== PxPyPzE");
388  ok+= compare( w2 == v2, static_cast<double>(true), "== PtEtaPhiE");
389  ok+= compare( w3 == v3, static_cast<double>(true), "== PtEtaPhiM");
390  ok+= compare( w4 == v4, static_cast<double>(true), "== PxPyPzM");
391 
392  // test gamma beta and boost
393  XYZVector b = q1.BoostToCM();
394  double beta = q1.Beta();
395  double gamma = q1.Gamma();
396 
397  ok += compare( b.R(), beta, "beta" );
398  ok += compare( gamma, 1./sqrt( 1 - beta*beta ), "gamma");
399 
400  if (ok == 0) std::cout << "\t OK " << std::endl;
401 
402  //test setters
403  std::cout << "Test Setters : " ;
404 
405  q2.SetXYZT(q1.Px(), q1.Py(), q1.Pz(), q1.E() );
406 
407  ok+= compare( q2.X(), q1.X(), "setXYZT X" );
408  ok+= compare( q2.Y(), q1.Y(), "setXYZT Y" );
409  ok+= compare( q2.Z(), q1.Z(), "setXYZT Z" ,2);
410  ok+= compare( q2.T(), q1.E(), "setXYZT E" );
411 
412  q2.SetCoordinates( 2.0*q1.Rho(), q1.Eta(), q1.Phi(), 2.0*q1.E() );
413  XYZTVector q1s = q1*2.0;
414  ok+= compare( q2.X(), q1s.X(), "set X" );
415  ok+= compare( q2.Y(), q1s.Y(), "set Y" );
416  ok+= compare( q2.Z(), q1s.Z(), "set Z" ,2);
417  ok+= compare( q2.T(), q1s.T(), "set E" );
418 
419  if (ok == 0) std::cout << "\t OK " << std::endl;
420 
421  return ok;
422 }
423 
424 int testVectorUtil() {
425  std::cout << "\n************************************************************************\n "
426  << " Utility Function Tests"
427  << "\n************************************************************************\n";
428 
429  std::cout << "Test Vector utility functions : ";
430 
431  XYZVector v1(1.0, 2.0, 3.0);
432  Polar3DVector v2pol(v1.R(), v1.Theta()+TMath::PiOver2(), v1.Phi() + 1.0);
433  // mixedmethods not yet impl.
434  XYZVector v2; v2 = v2pol;
435 
436  ok = 0;
437  ok += compare( VectorUtil::DeltaPhi(v1,v2), 1.0, "deltaPhi Vec");
438 
439  RhoEtaPhiVector v2cyl(v1.Rho(), v1.Eta() + 1.0, v1.Phi() + 1.0);
440  v2 = v2cyl;
441 
442  ok += compare( VectorUtil::DeltaR(v1,v2), sqrt(2.0), "DeltaR Vec");
443 
444  XYZVector vperp = v1.Cross(v2);
445  ok += compare( VectorUtil::CosTheta(v1,vperp), 0.0, "costheta Vec");
446  ok += compare( VectorUtil::Angle(v1,vperp), TMath::PiOver2(), "angle Vec");
447 
448  if (ok == 0) std::cout << "\t\t OK " << std::endl;
449 
450 
451  std::cout << "Test Point utility functions : ";
452 
453  XYZPoint p1(1.0, 2.0, 3.0);
454  Polar3DPoint p2pol(p1.R(), p1.Theta()+TMath::PiOver2(), p1.Phi() + 1.0);
455  // mixedmethods not yet impl.
456  XYZPoint p2; p2 = p2pol;
457 
458  ok = 0;
459  ok += compare( VectorUtil::DeltaPhi(p1,p2), 1.0, "deltaPhi Point");
460 
461  RhoEtaPhiPoint p2cyl(p1.Rho(), p1.Eta() + 1.0, p1.Phi() + 1.0);
462  p2 = p2cyl;
463  ok += compare( VectorUtil::DeltaR(p1,p2), sqrt(2.0), "DeltaR Point");
464 
465  XYZPoint pperp(vperp.X(), vperp.Y(), vperp.Z());
466  ok += compare( VectorUtil::CosTheta(p1,pperp), 0.0, "costheta Point");
467  ok += compare( VectorUtil::Angle(p1,pperp), TMath::PiOver2(), "angle Point");
468 
469  if (ok == 0) std::cout << "\t\t OK " << std::endl;
470 
471  std::cout << "LorentzVector utility funct.: ";
472 
473  XYZTVector q1(1.0, 2.0, 3.0,4.0);
474  PtEtaPhiEVector q2cyl(q1.Pt(), q1.Eta()+1.0, q1.Phi() + 1.0, q1.E() );
475  XYZTVector q2; q2 = q2cyl;
476 
477  ok = 0;
478  ok += compare( VectorUtil::DeltaPhi(q1,q2), 1.0, "deltaPhi LVec");
479  ok += compare( VectorUtil::DeltaR(q1,q2), sqrt(2.0), "DeltaR LVec");
480 
481  XYZTVector qsum = q1+q2;
482  ok += compare( VectorUtil::InvariantMass(q1,q2), qsum.M(), "InvMass");
483 
484  if (ok == 0) std::cout << "\t\t OK " << std::endl;
485 
486  return ok;
487 }
488 
489 int testRotation() {
490  std::cout << "\n************************************************************************\n "
491  << " Rotation and Transformation Tests"
492  << "\n************************************************************************\n";
493 
494  std::cout << "Test Vector Rotations : ";
495  ok = 0;
496 
497  XYZPoint v(1.,2,3.);
498 
499  double pi = TMath::Pi();
500  // initiate rotation with some non -trivial angles to test all matrix
501  EulerAngles r1( pi/2.,pi/4., pi/3 );
502  Rotation3D r2(r1);
503  // only operator= is in CINT for the other rotations
504  Quaternion r3; r3 = r2;
505  AxisAngle r4; r4 = r3;
506  RotationZYX r5; r5 = r2;
507 
508  XYZPoint v1 = r1 * v;
509  XYZPoint v2 = r2 * v;
510  XYZPoint v3 = r3 * v;
511  XYZPoint v4 = r4 * v;
512  XYZPoint v5 = r5 * v;
513 
514  ok+= compare(v1.X(), v2.X(), "x",2);
515  ok+= compare(v1.Y(), v2.Y(), "y",2);
516  ok+= compare(v1.Z(), v2.Z(), "z",2);
517 
518  ok+= compare(v1.X(), v3.X(), "x",2);
519  ok+= compare(v1.Y(), v3.Y(), "y",2);
520  ok+= compare(v1.Z(), v3.Z(), "z",2);
521 
522  ok+= compare(v1.X(), v4.X(), "x",5);
523  ok+= compare(v1.Y(), v4.Y(), "y",5);
524  ok+= compare(v1.Z(), v4.Z(), "z",5);
525 
526  ok+= compare(v1.X(), v5.X(), "x",2);
527  ok+= compare(v1.Y(), v5.Y(), "y",2);
528  ok+= compare(v1.Z(), v5.Z(), "z",2);
529 
530  // test with matrix
531  double rdata[9];
532  r2.GetComponents(rdata, rdata+9);
533  TMatrixD m(3,3,rdata);
534  double vdata[3];
535  v.GetCoordinates(vdata);
536  TVectorD q(3,vdata);
537  TVectorD q2 = m*q;
538 
539  XYZPoint v6;
540  v6.SetCoordinates( q2.GetMatrixArray() );
541 
542  ok+= compare(v1.X(), v6.X(), "x");
543  ok+= compare(v1.Y(), v6.Y(), "y");
544  ok+= compare(v1.Z(), v6.Z(), "z");
545 
546  if (ok == 0) std::cout << "\t OK " << std::endl;
547  else std::cout << std::endl;
548 
549  std::cout << "Test Axial Rotations : ";
550  ok = 0;
551 
552  RotationX rx( pi/3);
553  RotationY ry( pi/4);
554  RotationZ rz( 4*pi/5);
555 
556  Rotation3D r3x(rx);
557  Rotation3D r3y(ry);
558  Rotation3D r3z(rz);
559 
560  Quaternion qx; qx = rx;
561  Quaternion qy; qy = ry;
562  Quaternion qz; qz = rz;
563 
564  RotationZYX rzyx( rz.Angle(), ry.Angle(), rx.Angle() );
565 
566  XYZPoint vrot1 = rx * ry * rz * v;
567  XYZPoint vrot2 = r3x * r3y * r3z * v;
568 
569  ok+= compare(vrot1.X(), vrot2.X(), "x");
570  ok+= compare(vrot1.Y(), vrot2.Y(), "y");
571  ok+= compare(vrot1.Z(), vrot2.Z(), "z");
572 
573  vrot2 = qx * qy * qz * v;
574 
575  ok+= compare(vrot1.X(), vrot2.X(), "x",2);
576  ok+= compare(vrot1.Y(), vrot2.Y(), "y",2);
577  ok+= compare(vrot1.Z(), vrot2.Z(), "z",2);
578 
579  vrot2 = rzyx * v;
580 
581  ok+= compare(vrot1.X(), vrot2.X(), "x");
582  ok+= compare(vrot1.Y(), vrot2.Y(), "y");
583  ok+= compare(vrot1.Z(), vrot2.Z(), "z");
584 
585  // now inverse (first x then y then z)
586  vrot1 = rz * ry * rx * v;
587  vrot2 = r3z * r3y * r3x * v;
588 
589  ok+= compare(vrot1.X(), vrot2.X(), "x");
590  ok+= compare(vrot1.Y(), vrot2.Y(), "y");
591  ok+= compare(vrot1.Z(), vrot2.Z(), "z");
592 
593  XYZPoint vinv1 = rx.Inverse()*ry.Inverse()*rz.Inverse()*vrot1;
594 
595  ok+= compare(vinv1.X(), v.X(), "x",2);
596  ok+= compare(vinv1.Y(), v.Y(), "y");
597  ok+= compare(vinv1.Z(), v.Z(), "z");
598 
599  if (ok == 0) std::cout << "\t OK " << std::endl;
600  else std::cout << std::endl;
601 
602 
603  std::cout << "Test Rotations by a PI angle : ";
604  ok = 0;
605 
606  double b[4] = { 6,8,10,3.14159265358979323 };
607  AxisAngle arPi(b,b+4 );
608  Rotation3D rPi(arPi);
609  AxisAngle a1; a1 = rPi;
610  ok+= compare(arPi.Axis().X(), a1.Axis().X(),"x");
611  ok+= compare(arPi.Axis().Y(), a1.Axis().Y(),"y");
612  ok+= compare(arPi.Axis().Z(), a1.Axis().Z(),"z");
613  ok+= compare(arPi.Angle(), a1.Angle(),"angle");
614 
615  EulerAngles ePi; ePi=rPi;
616  EulerAngles e1; e1=Rotation3D(a1);
617  ok+= compare(ePi.Phi(), e1.Phi(),"phi");
618  ok+= compare(ePi.Theta(), e1.Theta(),"theta");
619  ok+= compare(ePi.Psi(), e1.Psi(),"ps1");
620 
621  if (ok == 0) std::cout << "\t\t OK " << std::endl;
622  else std::cout << std::endl;
623 
624  std::cout << "Test Inversions : ";
625  ok = 0;
626 
627  EulerAngles s1 = r1.Inverse();
628  Rotation3D s2 = r2.Inverse();
629  Quaternion s3 = r3.Inverse();
630  AxisAngle s4 = r4.Inverse();
631  RotationZYX s5 = r5.Inverse();
632 
633  // Euler angles not yet impl.
634  XYZPoint p = s2 * r2 * v;
635 
636  ok+= compare(p.X(), v.X(), "x",10);
637  ok+= compare(p.Y(), v.Y(), "y",10);
638  ok+= compare(p.Z(), v.Z(), "z",10);
639 
640 
641  p = s3 * r3 * v;
642 
643  ok+= compare(p.X(), v.X(), "x",10);
644  ok+= compare(p.Y(), v.Y(), "y",10);
645  ok+= compare(p.Z(), v.Z(), "z",10);
646 
647  p = s4 * r4 * v;
648  // axis angle inversion not very precise
649  ok+= compare(p.X(), v.X(), "x",1E9);
650  ok+= compare(p.Y(), v.Y(), "y",1E9);
651  ok+= compare(p.Z(), v.Z(), "z",1E9);
652 
653  p = s5 * r5 * v;
654 
655  ok+= compare(p.X(), v.X(), "x",10);
656  ok+= compare(p.Y(), v.Y(), "y",10);
657  ok+= compare(p.Z(), v.Z(), "z",10);
658 
659 
660  Rotation3D r6(r5);
661  Rotation3D s6 = r6.Inverse();
662 
663  p = s6 * r6 * v;
664 
665  ok+= compare(p.X(), v.X(), "x",10);
666  ok+= compare(p.Y(), v.Y(), "y",10);
667  ok+= compare(p.Z(), v.Z(), "z",10);
668 
669  if (ok == 0) std::cout << "\t OK " << std::endl;
670  else std::cout << std::endl;
671 
672  // test Rectify
673  std::cout << "Test rectify : ";
674  ok = 0;
675 
676  XYZVector u1(0.999498,-0.00118212,-0.0316611);
677  XYZVector u2(0,0.999304,-0.0373108);
678  XYZVector u3(0.0316832,0.0372921,0.998802);
679  Rotation3D rr(u1,u2,u3);
680  // check orto-normality
681  XYZPoint vrr = rr* v;
682  ok+= compare(v.R(), vrr.R(), "R",1.E9);
683 
684  if (ok == 0) std::cout << "\t\t OK " << std::endl;
685  else std::cout << std::endl;
686 
687  std::cout << "Test Transform3D : ";
688  ok = 0;
689 
690  XYZVector d(1.,-2.,3.);
691  Transform3D t(r2,d);
692 
693  XYZPoint pd = t * v;
694  // apply directly rotation
695  XYZPoint vd = r2 * v + d;
696 
697  ok+= compare(pd.X(), vd.X(), "x");
698  ok+= compare(pd.Y(), vd.Y(), "y");
699  ok+= compare(pd.Z(), vd.Z(), "z");
700 
701  // test with matrix
702  double tdata[12];
703  t.GetComponents(tdata);
704  TMatrixD mt(3,4,tdata);
705  double vData[4]; // needs a vector of dim 4
706  v.GetCoordinates(vData);
707  vData[3] = 1;
708  TVectorD q0(4,vData);
709 
710  TVectorD qt = mt*q0;
711 
712  ok+= compare(pd.X(), qt(0), "x");
713  ok+= compare(pd.Y(), qt(1), "y");
714  ok+= compare(pd.Z(), qt(2), "z");
715 
716  // test inverse
717 
718  Transform3D tinv = t.Inverse();
719 
720  p = tinv * t * v;
721 
722  ok+= compare(p.X(), v.X(), "x",10);
723  ok+= compare(p.Y(), v.Y(), "y",10);
724  ok+= compare(p.Z(), v.Z(), "z",10);
725 
726  // test construct inverse from translation first
727 
728  Transform3D tinv2 ( r2.Inverse(), r2.Inverse() *( -d) ) ;
729  p = tinv2 * t * v;
730 
731  ok+= compare(p.X(), v.X(), "x",10);
732  ok+= compare(p.Y(), v.Y(), "y",10);
733  ok+= compare(p.Z(), v.Z(), "z",10);
734 
735  // test from only rotation and only translation
736  Transform3D ta( EulerAngles(1.,2.,3.) );
737  Transform3D tb( XYZVector(1,2,3) );
738  Transform3D tc( Rotation3D(EulerAngles(1.,2.,3.)) , XYZVector(1,2,3) );
739  Transform3D td( ta.Rotation(), ta.Rotation() * XYZVector(1,2,3) ) ;
740 
741  ok+= compare( tc == tb*ta, static_cast<double>(true), "== Rot*Tra");
742  ok+= compare( td == ta*tb, static_cast<double>(true), "== Rot*Tra");
743 
744 
745  if (ok == 0) std::cout << "\t OK " << std::endl;
746  else std::cout << std::endl;
747 
748  std::cout << "Test Plane3D : ";
749  ok = 0;
750 
751  // test transfrom a 3D plane
752 
753  XYZPoint p1(1,2,3);
754  XYZPoint p2(-2,-1,4);
755  XYZPoint p3(-1,3,2);
756  Plane3D plane(p1,p2,p3);
757 
758  XYZVector n = plane.Normal();
759  // normal is perpendicular to vectors on the planes obtained from subtracting the points
760  ok+= compare(n.Dot(p2-p1), 0.0, "n.v12",10);
761  ok+= compare(n.Dot(p3-p1), 0.0, "n.v13",10);
762  ok+= compare(n.Dot(p3-p2), 0.0, "n.v23",10);
763 
764  Plane3D plane1 = t(plane);
765 
766  // transform the points
767  XYZPoint pt1 = t(p1);
768  XYZPoint pt2 = t(p2);
769  XYZPoint pt3 = t(p3);
770  Plane3D plane2(pt1,pt2,pt3);
771 
772  XYZVector n1 = plane1.Normal();
773  XYZVector n2 = plane2.Normal();
774 
775  ok+= compare(n1.X(), n2.X(), "a",10);
776  ok+= compare(n1.Y(), n2.Y(), "b",10);
777  ok+= compare(n1.Z(), n2.Z(), "c",10);
778  ok+= compare(plane1.HesseDistance(), plane2.HesseDistance(), "d",10);
779 
780  // check distances
781  ok += compare(plane1.Distance(pt1), 0.0, "distance",10);
782 
783  if (ok == 0) std::cout << "\t OK " << std::endl;
784  else std::cout << std::endl;
785 
786  std::cout << "Test LorentzRotation : ";
787  ok = 0;
788 
789  XYZTVector lv(1.,2.,3.,4.);
790 
791  // test from rotx (using boosts and 3D rotations not yet impl.)
792  // rx,ry and rz already defined
793  Rotation3D r3d = rx*ry*rz;
794 
795  LorentzRotation rlx(rx);
796  LorentzRotation rly(ry);
797  LorentzRotation rlz(rz);
798 
799  LorentzRotation rl0 = rlx*rly*rlz;
800  LorentzRotation rl1( r3d);
801 
802  XYZTVector lv0 = rl0 * lv;
803 
804  XYZTVector lv1 = rl1 * lv;
805 
806  XYZTVector lv2 = r3d * lv;
807 
808  ok+= compare(lv1== lv2,true,"V0==V2");
809  ok+= compare(lv1== lv2,true,"V1==V2");
810 
811  double rlData[16];
812  rl0.GetComponents(rlData);
813  TMatrixD ml(4,4,rlData);
814  double lvData[4];
815  lv.GetCoordinates(lvData);
816  TVectorD ql(4,lvData);
817 
818  TVectorD qlr = ml*ql;
819 
820  ok+= compare(lv1.X(), qlr(0), "x");
821  ok+= compare(lv1.Y(), qlr(1), "y");
822  ok+= compare(lv1.Z(), qlr(2), "z");
823  ok+= compare(lv1.E(), qlr(3), "t");
824 
825  // test inverse
826  lv0 = rl0 * rl0.Inverse() * lv;
827 
828  ok+= compare(lv0.X(), lv.X(), "x");
829  ok+= compare(lv0.Y(), lv.Y(), "y");
830  ok+= compare(lv0.Z(), lv.Z(), "z");
831  ok+= compare(lv0.E(), lv.E(), "t");
832 
833  if (ok == 0) std::cout << "\t OK " << std::endl;
834  else std::cout << std::endl;
835 
836  // test Boosts
837  std::cout << "Test Boost : ";
838  ok = 0;
839 
840  Boost bst( 0.3,0.4,0.5); // boost (must be <= 1)
841 
842  XYZTVector lvb = bst ( lv );
843 
844  LorentzRotation rl2 (bst);
845 
846  XYZTVector lvb2 = rl2 (lv);
847 
848  // test with lorentz rotation
849  ok+= compare(lvb.X(), lvb2.X(), "x");
850  ok+= compare(lvb.Y(), lvb2.Y(), "y");
851  ok+= compare(lvb.Z(), lvb2.Z(), "z");
852  ok+= compare(lvb.E(), lvb2.E(), "t");
853  ok+= compare(lvb.M(), lv.M(), "m",50); // m must stay constant
854 
855  // test inverse
856  lv0 = bst.Inverse() * lvb;
857 
858  ok+= compare(lv0.X(), lv.X(), "x",5);
859  ok+= compare(lv0.Y(), lv.Y(), "y",5);
860  ok+= compare(lv0.Z(), lv.Z(), "z",3);
861  ok+= compare(lv0.E(), lv.E(), "t",3);
862 
863  XYZVector brest = lv.BoostToCM();
864  bst.SetComponents( brest.X(), brest.Y(), brest.Z() );
865 
866  XYZTVector lvr = bst * lv;
867 
868  ok+= compare(lvr.X(), 0.0, "x",10);
869  ok+= compare(lvr.Y(), 0.0, "y",10);
870  ok+= compare(lvr.Z(), 0.0, "z",10);
871  ok+= compare(lvr.M(), lv.M(), "m",10);
872 
873  if (ok == 0) std::cout << "\t OK " << std::endl;
874  else std::cout << std::endl;
875  return ok;
876 }
877 
878 void mathcoreGenVector() {
879 
880 
881  testVector3D();
882  testPoint3D();
883  testLorentzVector();
884  testVectorUtil();
885  testRotation();
886 
887  std::cout << "\n\nNumber of tests " << ntest << " failed = " << nfail << std::endl;
888 }