JGeometry3D.cc

Go to the documentation of this file.

 
#include <iostream>
#include <iomanip>
 
#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;
  ULong_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);
}