Logo ROOT  
Reference Guide
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"
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 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
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 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;
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
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}
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define s1(x)
Definition: RSha256.hxx:91
char name[80]
Definition: TGX11.cxx:109
float * q
Definition: THbookFile.cxx:87
double sqrt(double)
AxisAngle class describing rotation represented with direction axis (3D Vector) and an angle of rotat...
Definition: AxisAngle.h:41
AxisVector Axis() const
accesss to rotation axis
Definition: AxisAngle.h:177
Scalar Angle() const
access to rotation angle
Definition: AxisAngle.h:182
AxisAngle Inverse() const
Return inverse of an AxisAngle rotation.
Definition: AxisAngle.h:260
Lorentz boost class with the (4D) transformation represented internally by a 4x4 orthosymplectic matr...
Definition: Boost.h:46
Scalar Theta() const
Polar theta, converting if necessary from internal coordinate system.
Scalar R() const
Polar R, converting if necessary from internal coordinate system.
Scalar X() const
Cartesian X, converting if necessary from internal coordinate system.
Scalar Y() const
Cartesian Y, converting if necessary from internal coordinate system.
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...
Scalar Rho() const
Cylindrical transverse component rho.
DisplacementVector3D< CoordSystem, Tag > & SetCoordinates(const Scalar src[])
Set internal data based on a C-style array of 3 Scalar numbers.
Scalar Phi() const
Polar phi, converting if necessary from internal coordinate system.
Scalar Eta() const
Polar eta, converting if necessary from internal coordinate system.
Scalar Z() const
Cartesian Z, converting if necessary from internal coordinate system.
EulerAngles class describing rotation as three angles (Euler Angles).
Definition: EulerAngles.h:43
Scalar Psi() const
Return Psi Euler angle.
Definition: EulerAngles.h:226
Scalar Theta() const
Return Theta Euler angle.
Definition: EulerAngles.h:216
Scalar Phi() const
Return Phi Euler angle.
Definition: EulerAngles.h:206
Class describing a geometrical plane in 3 dimensions.
Definition: Plane3D.h:51
Vector Normal() const
Return normal vector to the plane as Cartesian DisplacementVector.
Definition: Plane3D.h:156
Scalar HesseDistance() const
Return the Hesse Distance (distance from the origin) of the plane or the d coefficient expressed in n...
Definition: Plane3D.h:162
Scalar Distance(const Point &p) const
Return the signed distance to a Point.
Definition: Plane3D.h:170
Basic 3D Transformation class describing a rotation and then a translation The internal data are a 3D...
Definition: Transform3D.h:78
Transform3D< T > Inverse() const
Return the inverse of the transformation.
Definition: Transform3D.h:875
Lorentz transformation class with the (4D) transformation represented by a 4x4 orthosymplectic matrix...
LorentzRotation Inverse() const
Return inverse of a rotation.
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...
Scalar E() const
return 4-th component (time, or energy for a 4-momentum vector)
BetaVector BoostToCM() const
The beta vector for the boost that would bring this vector into its center of mass frame (zero moment...
Scalar M() const
return magnitude (mass) using the (-,-,-,+) metric.
Scalar Pt() const
return the transverse spatial component sqrt ( X**2 + Y**2 )
Scalar Beta() const
Return beta scalar value.
Scalar Eta() const
pseudorapidity
Scalar Py() const
spatial Y component
Scalar Pz() const
spatial Z component
Scalar Phi() const
azimuthal Angle
Scalar Gamma() const
Return Gamma scalar value.
Scalar Px() const
spatial X component
Class describing a generic position vector (point) in 3 dimensions.
Scalar Y() const
Cartesian Y, converting if necessary from internal coordinate system.
Scalar Z() const
Cartesian Z, converting if necessary from internal coordinate system.
Scalar X() const
Cartesian X, converting if necessary from internal coordinate system.
PositionVector3D< CoordSystem, Tag > & SetCoordinates(const Scalar src[])
Set internal data based on a C-style array of 3 Scalar numbers.
Scalar Dot(const DisplacementVector3D< OtherCoords, Tag > &v) const
Return the scalar (Dot) product of this with a displacement vector in any coordinate system,...
Scalar R() const
Polar R, converting if necessary from internal coordinate system.
Rotation class with the (3D) rotation represented by a unit quaternion (u, i, j, k).
Definition: Quaternion.h:47
Quaternion Inverse() const
Return inverse of a rotation.
Definition: Quaternion.h:253
Rotation class with the (3D) rotation represented by a 3x3 orthogonal matrix.
Definition: Rotation3D.h:65
Rotation3D Inverse() const
Return inverse of a rotation.
Definition: Rotation3D.h:414
Rotation class representing a 3D rotation about the X axis by the angle of rotation.
Definition: RotationX.h:43
Rotation class representing a 3D rotation about the Y axis by the angle of rotation.
Definition: RotationY.h:43
Rotation class with the (3D) rotation represented by angles describing first a rotation of an angle p...
Definition: RotationZYX.h:61
RotationZYX Inverse() const
Return inverse of a rotation.
Definition: RotationZYX.h:267
Rotation class representing a 3D rotation about the Z axis by the angle of rotation.
Definition: RotationZ.h:43
Element * GetMatrixArray()
Definition: TVectorT.h:78
double beta(double x, double y)
Calculates the beta function.
const Int_t n
Definition: legend1.C:16
double gamma(double x)
Vector1::Scalar DeltaR(const Vector1 &v1, const Vector2 &v2)
Find difference in pseudorapidity (Eta) and Phi betwen two generic vectors The only requirements on t...
Definition: VectorUtil.h:95
Vector1::Scalar DeltaPhi(const Vector1 &v1, const Vector2 &v2)
Find aximutal Angle difference between two generic vectors ( v2.Phi() - v1.Phi() ) The only requireme...
Definition: VectorUtil.h:59
double Angle(const Vector1 &v1, const Vector2 &v2)
Find Angle between two vectors.
Definition: VectorUtil.h:138
double CosTheta(const Vector1 &v1, const Vector2 &v2)
Find CosTheta Angle between two generic 3D vectors pre-requisite: vectors implement the X(),...
Definition: VectorUtil.h:112
Vector1::Scalar InvariantMass(const Vector1 &v1, const Vector2 &v2)
return the invariant mass of two LorentzVector The only requirement on the LorentzVector is that they...
Definition: VectorUtil.h:215
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZVector
3D Vector based on the cartesian coordinates x,y,z in double precision
Definition: Vector3Dfwd.h:34
static constexpr double pi
constexpr Double_t PiOver2()
Definition: TMath.h:52
constexpr Double_t Pi()
Definition: TMath.h:38
#define Dot(u, v)
Definition: normal.c:49
auto * lv
Definition: textalign.C:5
auto * m
Definition: textangle.C:8