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/// 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
50using namespace ROOT::Math;
51
52int ntest = 0;
53int nfail = 0;
54int ok = 0;
55
56int 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
98int 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
232int 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 CLING
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
300int 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
424int 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
489int 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 CLING 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 transform 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
878void 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}
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
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
Dot
#define Dot(u, v)
Definition normal.c:49
v2
@ v2
Definition rootcling_impl.cxx:3690
v
@ v
Definition rootcling_impl.cxx:3687
v4
@ v4
Definition rootcling_impl.cxx:3692
v3
@ v3
Definition rootcling_impl.cxx:3691
v1
@ v1
Definition rootcling_impl.cxx:3689
lv
TLine lv
Definition textalign.C:5
m
TMarker m
Definition textangle.C:8