Logo ROOT  

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