geo.h

Go to the documentation of this file.

#ifndef GEO_GEO_H
#define GEO_GEO_H

#include "cetlib_except/exception.h"

#include "stdexcept"


namespace geo {
  void ProjectToBoxEdge(const double xyz[],
                        const double dxyz[],
                        double xlo, double xhi,
                        double ylo, double yhi,
                        double zlo, double zhi,
                        double xyzout[]);
  double ClosestApproach(const double point[],
                         const double intercept[],
                         const double slopes[],
                         double closest[]);
  bool CrossesBoundary ( double x0[],
                         double gradient[],
                         double x_lo, double x_hi,
                         double y_lo, double y_hi,
                         double z_lo, double z_hi,
                         double point[] );
}

inline void geo::ProjectToBoxEdge(const double xyz[],
                                  const double dxyz[],
                                  double xlo, double xhi,
                                  double ylo, double yhi,
                                  double zlo, double zhi,
                                  double xyzout[]) 
{
  // Make sure we're inside the box!
  if( !(xyz[0]>=xlo && xyz[0]<=xhi) ||
      !(xyz[1]>=ylo && xyz[1]<=yhi) ||
      !(xyz[2]>=zlo && xyz[2]<=zhi)  )
    throw cet::exception("ProjectToBoxEdge") << "desired point is not"
                                             << " in the specififed box\n"; 
  
  // Compute the distances to the x/y/z walls
  double dx = 99.E99;
  double dy = 99.E99;
  double dz = 99.E99;
  if      (dxyz[0]>0.0) { dx = (xhi-xyz[0])/dxyz[0]; }
  else if (dxyz[0]<0.0) { dx = (xlo-xyz[0])/dxyz[0]; }
  if      (dxyz[1]>0.0) { dy = (yhi-xyz[1])/dxyz[1]; }
  else if (dxyz[1]<0.0) { dy = (ylo-xyz[1])/dxyz[1]; }
  if      (dxyz[2]>0.0) { dz = (zhi-xyz[2])/dxyz[2]; }
  else if (dxyz[2]<0.0) { dz = (zlo-xyz[2])/dxyz[2]; }
  
  // Choose the shortest distance
  double d = 0.0;
  if      (dx<dy && dx<dz) d = dx;
  else if (dy<dz && dy<dx) d = dy;
  else if (dz<dx && dz<dy) d = dz;
  
  // Make the step
  for (int i=0; i<3; ++i) {
    xyzout[i] = xyz[i] + dxyz[i]*d;
  }
}

inline double geo::ClosestApproach(const double point[],
                                   const double intercept[],
                                   const double slopes[],
                                   double closest[]) 
{
  double s = 
    (slopes[0]*(point[0]-intercept[0]) + 
     slopes[1]*(point[1]-intercept[1]) +
     slopes[2]*(point[2]-intercept[2]));
  double sd = 
    (slopes[0]*slopes[0] + 
     slopes[1]*slopes[1] +
     slopes[2]*slopes[2]);
  if (sd>0.0) {
    s /= sd;
    closest[0] = intercept[0] + s*slopes[0];
    closest[1] = intercept[1] + s*slopes[1];
    closest[2] = intercept[2] + s*slopes[2];
  }
  else {
    // How to handle this zero gracefully? Assume that the intercept
    // is a particle vertex and "slopes" are momenta. In that case,
    // the particle goes nowhere and the closest approach is the
    // distance from the intercept to point
    closest[0] = intercept[0];
    closest[1] = intercept[1];
    closest[2] = intercept[2];
  }
  return std::sqrt(pow((point[0]-closest[0]),2)+
              pow((point[1]-closest[1]),2)+
              pow((point[2]-closest[2]),2));
}


inline bool geo::CrossesBoundary ( double x0[],          // initial particle position
                                   double gradient[],    // initial particle gradient
                                   double x_lo,          // -
                                   double x_hi,          //  |
                                   double y_lo,          //  |- box coordinates
                                   double y_hi,          //  |  (make into vectors?)
                                   double z_lo,          //  |
                                   double z_hi,          // -
                                   double point[] )
{
  
  double distance[3]; // distance to plane
  
  // puts box coordinates into more useful vectors (for loop later)
  double lo[3] = { x_lo , y_lo , z_lo };
  double hi[3] = { x_hi , y_hi , z_hi };
  
  // puts box coordinates into more useful vector (for loop later)
  double facecoord[6] = { lo[0] , hi[0] ,
                          lo[1] , hi[1] ,
                          lo[2] , hi[2] };
  
  int intersect[6]={0,0,0,0,0,0}; // initialize intersection tally vector
  int count=0;
  // iterates through spatial axes (0,1,2) = (x,y,z)
  for(int i=0; i<3; i++) {
    // try both planes with normal parallel to axis
    for(int p=0; p<2; p++ ) {
      
      point[i] = facecoord[count];           // point on face
      distance[i] = point[i] - x0[i];        // calculate x-coordinate of track distance
      
      double C_dg = distance[i] / gradient[i];  // calculate correlation b/w gradient and distance
      
      for(int m=0; m<3; m++){ distance[m] = C_dg * gradient[m];     }
      for(int n=0; n<3; n++){    point[n] = x0[n] + distance[n]; }
      
      int j, k;
      switch (i) {
        case 0: j=1; k=2; break;
        case 1: j=2; k=0; break;
        case 2: j=0; k=1; break;
        default:
          throw std::logic_error("Big trouble");
      } // switch
      
      // now want to check to see if the point is in the right plane
      if ( lo[j] < point[j] && point[j] < hi[j]
           && lo[k] < point[k] && point[k] < hi[k] ) {
        
        //           double length = std::sqrt( distance[0]*distance[0]
        //                               + distance[1]*distance[1]
        //                               + distance[2]*distance[2] );
        
        // direction of motion w.r.t. start point
        int direction = distance[i]*gradient[i]
          / std::sqrt( (distance[i]*distance[i]) * (gradient[i]*gradient[i]) );
        bool directed = ( direction + 1 ) / 2;
        
        // checks if particle passes through face
        // and also checks to see whether it passes inward or outward
        //if ( track_length > length && directed ) {
        if ( directed ) {
          int normal = pow( -1 , count + 1 );
          int thru = normal * gradient[i] / std::sqrt(gradient[i]*gradient[i]) ;
          intersect[count]=thru;
        }
      }
      count++;
    }
  }
  
  // count faces it passes through, 
  // ... not necessary now, but maybe useful in the future
  int passes=0;
  for ( int face=0; face<6; ++face ) { 
    passes+=(intersect[face]*intersect[face]); 
  }
  
  if ( passes==0 ) {
    return 0;
  } else {
    return 1;
  }
}

#endif