#include #include #include "TRandom3.h" #include "JTools/JQuantile.hh" #include "JGeometry3D/JPosition3D.hh" #include "JGeometry3D/JAngle3D.hh" #include "JGeometry3D/JRotation3D.hh" #include "JGeometry3D/JQuaternion3D.hh" #include "JMath/JMathToolkit.hh" #include "JMath/JConstants.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" #include "Jeep/JMessage.hh" /** * \file * * Example program to test rotations in 3D. * \author mdejong */ int main(int argc, char**argv) { using namespace std; int numberOfEvents; double precision; UInt_t seed; int debug; try { JParser<> zap("Example program to test rotations in 3D."); zap['n'] = make_field(numberOfEvents); zap['e'] = make_field(precision) = 1.0e-12; zap['S'] = make_field(seed) = 0; zap['d'] = make_field(debug) = 1; zap(argc, argv); } catch(const exception &error) { FATAL(error.what() << endl); } gRandom->SetSeed(seed); using namespace JPP; JQuantile Q("Distance"); const JVector3D v[] = { JVector3D(1,0,0), // unit vector along x-axis JVector3D(0,1,0), // unit vector along y-axis JVector3D(0,0,1) }; // unit vector along z-axis for (int i = 0; i != numberOfEvents; ++i) { STATUS("event: " << setw(10) << i << '\r'); DEBUG(endl); const double theta = gRandom->Uniform(0.0, PI); const double phi = gRandom->Uniform(0.0, 2*PI); const JAngle3D dir(theta, phi); // principal direction const JRotation3D R(dir); //const JRotation3Y Ry(theta); // rotation around y-axis //const JRotation3Z Rz(phi); // rotation around z-axis //const JQuaternion3D Ry(theta, JVersor3D(0,1,0)); // rotation around y-axis //const JQuaternion3D Rz(phi, JVersor3D(0,0,1)); // rotation around z-axis const JQuaternion3X Rx(0.0); // rotation around x-axis const JQuaternion3Y Ry(theta); // rotation around y-axis const JQuaternion3Z Rz(phi); // rotation around z-axis //const JQuaternion3D Rs = Rz * Ry; JQuaternion3D Rs(Rx, Ry, Rz); /* Rs.setIdentity(); Rs *= Rz; Rs *= Ry; */ for (int j = 0; j != sizeof(v)/sizeof(v[0]); ++j) { JPosition3D A(v[j]); JPosition3D B(A); //B.rotate(Ry); // rotates around y-axis //B.rotate(Rz); // rotates around z-axis B.rotate(Rs); B.rotate(R); // rotates principal direction to z-axis Q.put((B - A).getLength()); if ((B - A).getLengthSquared() > precision*precision) { DEBUG("theta " << theta << endl); DEBUG("phi " << phi << endl); //DEBUG("Ry" << endl << Ry); //DEBUG("Rz" << endl << Rz); DEBUG("R " << endl << R); ERROR("A " << FIXED(7,3) << A.getX() << ' ' << A.getY() << ' ' << A.getZ() << endl); ERROR("B " << FIXED(7,3) << B.getX() << ' ' << B.getY() << ' ' << B.getZ() << endl); } } } STATUS(endl); Q.print(cout); }