JCylinder3D.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <vector>
 
#include "TRandom3.h"
 
#include "JGeometry3D/JAxis3D.hh"
#include "JGeometry3D/JCylinder3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to test intersection of straight line with cylinder.
 * \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 intersection of straight line with cylinder.");
 
    zap['n'] = make_field(numberOfEvents)      = 1000;
    zap['e'] = make_field(precision)           = 1.0e-12;
    zap['S'] = make_field(seed)                = 0;
    zap['d'] = make_field(debug)               = 3;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  gRandom->SetSeed(seed);
 
 
  using namespace JPP;
 
  const JCylinder3D cylinder(JCircle2D(JVector2D(0.0,0.0), 2.0), -1.0, +1.0);
  
 
  int number_of_errors = 0;
  int number_of_events = 0;
 
  for (int i = 0; i != numberOfEvents; ++i) {
 
    STATUS(setw(4) << i << '\r'); DEBUG(endl);
 
    const double x  = gRandom->Uniform(-1.0, +1.0);
    const double y  = gRandom->Uniform(-1.0, +1.0);
    const double z  = gRandom->Uniform(-1.0, +1.0);
 
    double dx;
    double dy;
    double dz;
 
    gRandom->Sphere(dx, dy, dz, 1.0);
 
    JAxis3D axis(JVector3D(x, y, z), JVersor3D(dx, dy, dz));
 
    JCylinder3D::intersection_type intersection = cylinder.getIntersection(axis);
 
    const pair<double, bool> path[] = { 
      make_pair(intersection.first  - 1*precision,                      false),
      make_pair(2*precision,                                            true),
      make_pair(intersection.second - intersection.first - 2*precision, true),
      make_pair(2*precision,                                            false)
    };
 
    ++number_of_events;
    
    bool ok = true;
 
    const int N = sizeof(path)/sizeof(path[0]);
    
    JVector3D pos[N];
 
    for (int i = 0; i != N; ++i) {
                    
      axis.move(path[i].first);
      
      if ((getDistance(cylinder, axis.getPosition()) < 1e-3*precision) != path[i].second) {
        ok = false;
      }
 
      pos[i] = axis.getPosition();
    }
 
    if (!ok) {
 
      for (int i = 0; i != N; ++i) {
        cout << endl << "point[" << noshowpos << setw(1) << i << "] = " 
             << showpos << FIXED(9,5) << pos[i].getX()  << ' '
             << showpos << FIXED(9,5) << pos[i].getY()  << ' '
             << showpos << FIXED(9,5) << pos[i].getZ()  << ' '
             << SCIENTIFIC(12,3) << cylinder.getDistance(pos[i]) << ' ' 
             << (cylinder.getDistance(pos[i]) < 1e-3*precision ? "in" : "out") << endl;
      }
    }
 
    if (!ok) {
      ++number_of_errors;
    }
  }
  STATUS(endl);
 
  NOTICE("Number of errors: " << number_of_errors << " out of " << number_of_events << " events." << endl);
 
  ASSERT(number_of_events >  0);
  ASSERT(number_of_errors == 0);
 
  return 0;
}