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;
  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);
}