Geo.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file    geo.h
/// \brief   Supply basic geometry functions
/// \author  messier@indiana.edu
/// \version
////////////////////////////////////////////////////////////////////////
#include "GeometryObjects/Geo.h"

#include <cmath>
#include <cassert>
#include <cfloat>
#include <iostream>
#include "TVector3.h"

#include "GeometryObjects/CellGeo.h"
#include "Utilities/func/MathUtil.h"
#include "cetlib_except/exception.h"

namespace geo {

  //......................................................................
  /// \brief Project along a direction from a particular starting point to the
  /// edge of a box
  ///
  /// \param xyz    - The starting x,y,z location. Must be inside box.
  /// \param dxyz   - Direction vector
  /// \param xlo    - Low edge of box in x
  /// \param xhi    - Low edge of box in x
  /// \param ylo    - Low edge of box in y
  /// \param yhi    - Low edge of box in y
  /// \param zlo    - Low edge of box in z
  /// \param zhi    - Low edge of box in z
  /// \param xyzout - On output, the position at the box edge
  ///
  /// Note: It should be safe to use the same array for input and
  /// output.
  ///
  void 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("Geo") << "provided point is outside of the specified box";
    
    // 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;
    }
  }

  int WhichWallofBox(const double xyz[],
                        double xlo, double xhi,
                        double ylo, double yhi,
                        double zlo, double zhi)
  {
    int wall = 0;
    if(xyz[0] >= xlo-0.01 && xyz[0] <= xlo+0.01) wall = 1;
    if(xyz[0] >= xhi-0.01 && xyz[0] <= xhi+0.01) wall = 2;
    if(xyz[1] >= ylo-0.01 && xyz[1] <= ylo+0.01) wall = 3;
    if(xyz[1] >= yhi-0.01 && xyz[1] <= yhi+0.01) wall = 4;
    if(xyz[2] >= zlo-0.01 && xyz[2] <= zlo+0.01) wall = 5;
    if(xyz[2] >= zhi-0.01 && xyz[2] <= zhi+0.01) wall = 6;

//    if(wall == 0) std::cout << xyz[0] << " " << xyz[1] << " " << xyz[2] << std::endl;

    return wall;
    
  }

  //......................................................................
  /// \brief Determine if a particle starting at xyz with direction dxyz will
  /// intersect a box defined by xlo, xhi, ylo, yhi, zlo, zhi
  bool IntersectsBox(const double xyz[],
                     const double dxyz[],
                     double xlo, double xhi,
                     double ylo, double yhi,
                     double zlo, double zhi)
  {
    double xyzedge[3] = {0.};

    // there are 6 edges to a box, check them all.  If a particle intersects
    // any surface, then return true

    // top surface, y = yhi
    geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 1, yhi, xyzedge);
    if(xyzedge[0] > xlo && xyzedge[0] < xhi
       && xyzedge[2] > zlo && xyzedge[2] < zhi) return true;

    // bottom surface, y = ylo
    geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 1, ylo, xyzedge);
    if(xyzedge[0] > xlo && xyzedge[0] < xhi
       && xyzedge[2] > zlo && xyzedge[2] < zhi) return true;

    // x = xhi
    geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 0, xhi, xyzedge);
    if(xyzedge[1] > ylo && xyzedge[1] < yhi
       && xyzedge[2] > zlo && xyzedge[2] < zhi) return true;

    // x = xlo
    geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 0, xlo, xyzedge);
    if(xyzedge[1] > ylo && xyzedge[1] < yhi
       && xyzedge[2] > zlo && xyzedge[2] < zhi) return true;

    // z = zhi
    geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 2, zhi, xyzedge);
    if(xyzedge[1] > ylo && xyzedge[1] < yhi
       && xyzedge[0] > xlo && xyzedge[0] < xhi) return true;

    // z = zlo
    geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 2, zlo, xyzedge);
    if(xyzedge[1] > ylo && xyzedge[1] < yhi
       && xyzedge[0] > xlo && xyzedge[0] < xhi) return true;

    return false;

  }

  //......................................................................
  /// \brief Project from a position outside of a box to an edge of the box
  /// with coordinate value \a edge for the axis \a axis.
  ///
  /// axis = 0 --> x, axis = 1 --> y, axis = 2 --> z
  /// \a xyzout is the position on the desired edge
  void ProjectToBoxEdgeFromOutside(const double xyz[],
                                   const double dxyz[],
                                   int    axis,
                                   double edge,
                                   double xyzout[])
  {
    double momDir = dxyz[axis];
    double posDir = xyz[axis];
    double momTot = sqrt(dxyz[0]*dxyz[0] + dxyz[1]*dxyz[1] + dxyz[2]*dxyz[2]);

    if(momTot < 1.e-9){
      std::cerr << "no direction information supplied, do nothing" << std::endl;
      return;
    }

    double ddS        = momDir/momTot;   // direction cosine in the desired direction
    double length1Dim = (edge - posDir); // length in desired direction to that edge;

    if(TMath::Abs(ddS) > 0.){
      length1Dim /= ddS;
      xyzout[0] = xyz[0] + length1Dim*dxyz[0]/momTot;
      xyzout[1] = xyz[1] + length1Dim*dxyz[1]/momTot;
      xyzout[2] = xyz[2] + length1Dim*dxyz[2]/momTot;
    }

    return;
  }

  //......................................................................
  /// \brief Find the intersection between two line segments
  ///
  /// Given the two line segments (x0,y0)-(x1,y1), and (X0,Y0)-(X1,Y1), fills
  /// (x,y) with the point that lies on both of them. The return value
  /// indicates if the intersection of the (infinite) lines falls inside the
  /// line segments. If it does not, the returned (x,y) is still meaningful,
  /// unless the two input lines are parallel.
  bool LineIntersection(double x0, double y0, double x1, double y1,
                        double X0, double Y0, double X1, double Y1,
                        double& x, double& y)
  {
    const double dx = x1-x0;
    const double dy = y1-y0;
    const double dX = X1-X0;
    const double dY = Y1-Y0;

    const double denom = dX*dy-dY*dx;

    if(denom == 0){
      // std::cerr << "Trying to find intercept of parallel lines" << std::endl;
      x = x1; y = y1;
      return false;
    }

    // Distance of intercept point along each of the two lines
    const double l = -(X0*dY-Y0*dX+dX*y0-dY*x0)/denom;
    const double L = +(x0*dy-y0*dx+dx*Y0-dy*X0)/denom;

    x = X0+L*(X1-X0);
    y = Y0+L*(Y1-Y0);

    return l >= 0 && l <= 1 && L >= 0 && L <= 1;
  }

  //......................................................................
  ///
  /// \brief Find the distance of closest approach between point and line
  ///
  /// \param point - xyz coordinates of point
  /// \param intercept - xyz coordinates of point on line
  /// \param slopes - direction vector (need not be normalized)
  /// \param closest - on output, point on line that is closest
  ///
  /// \returns distance from point to line
  ///
  double ClosestApproach(const double point[],
                         const double intercept[],
                         const double slopes[],
                         double closest[])
  {
    TVector3 closest_vec;
    const double dist = ClosestApproach(TVector3(point),
                                        TVector3(intercept), TVector3(slopes),
                                        closest_vec);
    closest_vec.GetXYZ(closest);
    return dist;
  }

  //......................................................................
  ///
  /// \brief Find the distance of closest approach between point and line
  ///
  /// \param point - xyz coordinates of point
  /// \param intercept - xyz coordinates of point on line
  /// \param slopes - direction vector (need not be normalized)
  /// \param closest - on output, point on line that is closest
  ///
  /// \returns distance from point to line
  ///
  double ClosestApproach(TVector3 point, TVector3 intercept,
                         TVector3 slopes, TVector3& closest)
  {
    const double s = slopes.Dot(point-intercept);
    const double sd = slopes.Mag2();

    if(sd > 0){
      closest = intercept + (s/sd)*slopes;
    }
    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 = intercept;
    }

    return (point-closest).Mag();
  }

  //....................................................................
  /// \brief Find the distance of closest approach between two lines
  /// which pass through points P0 and P1 and Q0 and Q1
  /// 
  /// @param P0 - A point on line P
  /// @param P1 - A second point on line P
  /// @param Q0 - A point on line Q
  /// @param Q1 - A second point on line Q
  /// @param sc - distance from P0 to point on P which is closest
  /// @param tc - distance from Q0 to point on Q which is closest
  /// @param PC - point on line P which is closest
  /// @param QC - point on line Q which is closest
  ///
  /// @returns The distance^2 between the points of closest approach.
  ///
  /// If any of the pointers sc, tc, PC, QC are zero, that part of the
  /// calculation is skipped.
  ///
  double ClosestApproach(const TVector3& P0,
                         const TVector3& P1,
                         const TVector3& Q0,
                         const TVector3& Q1,
                         double* sc,
                         double* tc,
                         TVector3* PC,
                         TVector3* QC) 
  {
    TVector3 w0(P0); w0 -= Q0;

    TVector3 u(P1); u -= P0; u = u.Unit();
    TVector3 v(Q1); v -= Q0; v = v.Unit();
    
    double b = u.Dot(v);
    double d = u.Dot(w0);
    double e = v.Dot(w0);
    
    double denom = 1-b*b;
    double sC, tC;
    if (denom!=0.0) {
      sC = (b*e-d)/denom;
      tC = (e-b*d)/denom;
    }
    else {
      sC = 0.0;
      tC = d/b;
    }
    if (sc!=0) *sc = sC;
    if (tc!=0) *tc = tC;
    
    TVector3 pC(P0+sC*u);
    TVector3 qC(Q0+tC*v);
    if (PC!=0) *PC = pC;
    if (QC!=0) *QC = qC;
    
    TVector3 D(pC-qC);
    
    return D.Mag2();
  }

  //......................................................................
  ///
  /// \brief In two dimensions, return the perpendicular distance from a
  /// point (x0,y0) to a line defined by end points (x1,y1) and (x2,y2)
  ///
  /// @param x0 : x-value of point
  /// @param y0 : y-value of point
  /// @param x1 : x-value of point at one end of line
  /// @param y1 : y-value of point at one end of line
  /// @param x2 : x-value of point at other end of line
  /// @param y2 : y-value of point at other end of line
  /// @returns Squared distance from the point to the line
  ///
  double DsqrToLine(double x0, double y0,
                    double x1, double y1,
                    double x2, double y2)
  {
    using util::sqr;

    double d = 0.;

    // If we've been handed a point instead of a line, return distance
    // to the point
    if (x1==x2 && y1==y2) d = sqr(x0-x1)+sqr(y0-y1);
    else d = sqr((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1)) / (sqr(x2-x1)+sqr(y2-y1));

    return d;
      
  }

  //......................................................................
  ///
  /// \brief Find the best-fit line to a collection of points in 2-D by
  /// minimizing the squared vertical distance from the points to the line.
  ///
  /// To give the most robust fits and avoid cases where the slopes approach
  /// infinity, the points are arranged to maximize the horizontal extent. So,
  /// the fitter may fit x vs. y or y vs. x.
  ///
  /// @param x  - input vector of x coordinates
  /// @param y  - input vector of y coordinates
  /// @param w  - input vector of weights for the points
  /// @param x1 - output x coordinate at one end of line
  /// @param y1 - output y coordinate at one end of line
  /// @param x2 - output x coordinate at other end of line
  /// @param y2 - output y coordinate at other end of line
  /// @returns The chi^2 value defined by chi^2 = sum_i[w_i d^2_i]
  ///
  double LinFit(const std::vector<double>& x,
                const std::vector<double>& y,
                const std::vector<double>& w,
                double* x1, double* y1,
                double* x2, double* y2)
  {
    unsigned int i;

    // Before going ahead, make sure we have sensible arrays
    if( x.size() != y.size() ||
        x.size() != w.size() ||
        x.size() < 2)
      throw cet::exception("Geo") << "Input vectors to LinFit are not the "
                                  << "correct length. x size: " << x.size()
                                  << " y size: " << y.size()
                                  << " w size: " << w.size();

    // Find the ranges of the input variables
    double xlo = x[0];
    double xhi = x[0];
    double ylo = y[0];
    double yhi = y[0];
    for (i=1; i<w.size(); ++i) {
      if (x[i]<xlo) xlo = x[i];
      if (x[i]>xhi) xhi = x[i];
      if (y[i]<ylo) ylo = y[i];
      if (y[i]>yhi) yhi = y[i];
    }
    // Make sure the fit has some span
    if (yhi-ylo==0.0 && xhi-xlo==0.0) {
      std::cerr << __FILE__ << ":" << __LINE__
                << " All points identical in line fit!" << std::endl;
      *x1 = xlo;
      *y1 = ylo;
      *x2 = xhi;
      *y2 = yhi;
      return 0.0;
    }

    //
    // If the y extent is larger than the x extent, flip the variables
    // and fit
    //
    if (yhi-ylo > xhi-xlo) {
      return geo::LinFit(y,x,w,y1,x1,y2,x2);
    }

    double m, b;
    const double chi2 = util::LinFit(x, y, w, m, b);

    *x1 = xlo;
    *x2 = xhi;
    *y1 = m*(*x1)+b;
    *y2 = m*(*x2)+b;

    return chi2;
  }

  //......................................................................
  ///
  /// \brief Find the best-fit line to a collection of points in 2-D by
  /// minimizing the squared perpendicular distance from the points to the
  /// line.
  ///
  /// This is not the usual "linfit" which minimizes the squared
  /// vertical distance
  ///
  /// @param x  - input vector of x coordinates
  /// @param y  - input vector of y coordinates
  /// @param w  - input vector of weights for the points
  /// @param x1 - output x coordinate at one end of line
  /// @param y1 - output y coordinate at one end of line
  /// @param x2 - output x coordinate at other end of line
  /// @param y2 - output y coordinate at other end of line
  /// @returns The chi^2 value defined by chi^2 = sum_i[w_i d^2_i]
  ///
  double LinFitMinDperp(const std::vector<double>& x,
                        const std::vector<double>& y,
                        const std::vector<double>& w,
                        double* x1, double* y1,
                        double* x2, double* y2)
  {
    // Before going ahead, make sure we have sensible arrays
    if( x.size() != y.size() ||
        x.size() != w.size() ||
        x.size() < 2)
      throw cet::exception("Geo") << "Input vectors to LinFit are not the "
                                  << "correct length. x size: " << x.size()
                                  << " y size: " << y.size()
                                  << " w size: " << w.size();

    //--
    // The following code implements the output of the following MAPLE
    // script:
    //--
    // assume(x0,real);
    // assume(y0,real);
    // assume({M,B},real);
    // assume({Dsqr},real);
    // assume(Sw,real);
    // assume(Swx,real);
    // assume(Swx2,real);
    // assume(Swy,real);
    // assume(Swy2,real);
    // assume(Swxy,real);
    // y1 := M*x1 + B;
    // y2 := M*x2 + B;
    // Dsqr := (((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1))^2/((x2-x1)^2+(y2-y1)^2));
    // eqn1 := simplify(expand(diff(Dsqr,M) = 0));
    // eqn2 := simplify(expand(diff(Dsqr,B) = 0));
    // Eqn1 := -2*B*Swy*M +B^2*M*Sw +Swy2*M -Swxy*M^2 -B*Swx +Swxy
    // +B*Swx*M^2 -Swx2*M = 0;
    // Eqn2 := B*Sw - Swy + M*Swx = 0;
    // solve({Eqn1,Eqn2},{M,B});
    // allvalues(%);
    //--

    // register is deprecated in c++17 standard
    //unsigned register int i;
    unsigned int i;

    // Find the extremes of the inputs
    double xlo = DBL_MAX;
    double xhi = -DBL_MAX;
    double ylo = DBL_MAX;
    double yhi = -DBL_MAX;
    for (i=0; i<w.size(); ++i) {
      if(w[i] == 0) continue;
      if (x[i]<xlo) xlo = x[i];
      if (y[i]<ylo) ylo = y[i];
      if (x[i]>xhi) xhi = x[i];
      if (y[i]>yhi) yhi = y[i];
    }

    // This algorithm really fails if the y values are identical
    // return the low and hi values if it's too flat
    if(fabs(yhi-ylo) < 0.1){
      *x1 = xlo;
      *y1 = ylo;
      *x2 = xhi;
      *y2 = yhi;
      return 0.; //its a flat line, chi^2 is 0 by definition
    }

    // For stability, fit in a way that avoids infinite slopes,
    // exchanging the x and y variables as needed.
    if (fabs(xlo - xhi) < 1.e-3) {
      double chi2;
      double xx1, yy1, xx2, yy2;
      // Notice: (x,y) -> (y,x)
      chi2 = LinFitMinDperp(y,x,w,&xx1, &yy1,&xx2, &yy2);
      *x1 = xlo; // Infinite slope means that your initial x positions should be the input
      *y1 = xx1;
      *x2 = xhi;
      *y2 = xx2;
      return chi2;
    }

    // Accumulate the sums needed to compute the fit line
    double Sw  = 0.0;
    double Swx = 0.0;
    double Swy = 0.0;
    double Swxy= 0.0;
    double Swy2= 0.0;
    double Swx2= 0.0;
    for (i=0; i<w.size(); ++i) {
      Sw   += w[i];
      Swx  += w[i]*x[i];
      Swy  += w[i]*y[i];
      Swxy += w[i]*x[i]*y[i];
      Swy2 += w[i]*y[i]*y[i];
      Swx2 += w[i]*x[i]*x[i];
    }

    if(Sw == 0.0)
      throw cet::exception("Geo") << "sum of weights is 0 in LinFitMinDperp";

    // Precalculate reused squares
    using util::sqr;
    const double SwSq  = sqr(Sw);
    const double SwxSq = sqr(Swx);
    const double SwySq = sqr(Swy);

    // C and D are two handy temporary variables
    double C =
      +sqr(SwxSq) // Swx^4
      -2.0*SwxSq*Swx2*Sw
      +2.0*SwxSq*SwySq
      +2.0*SwxSq*Swy2*Sw
      +sqr(Swx2)*SwSq
      +2.0*Swx2*Sw*SwySq
      -2.0*Swx2*SwSq*Swy2
      +sqr(sqr(Swy)) // Swy^4
      -2.0*SwySq*Swy2*Sw
      +sqr(Swy2)*SwSq
      -8.0*Swx*Swy*Swxy*Sw
      +4.0*sqr(Swxy)*SwSq;
    if (C>0.0) C = sqrt(C);
    else       C = 0.0;
    double D = -Swx*Swy+Swxy*Sw;

    // The first of the two solutions. This one is the best fit through
    // the points
    double M1;
    double B1;
    if (D==0.0) {
      // D could be zero, this might happen in a segmented detector if
      // the track fit passed through only a single readout plane and
      // all hits share the same x coordinates. I've tried to avoid
      // this by checking if it makes more sense to fit with x and y
      // exchanged above. This is an extra precaution.
      D = 1.0E-9;
    }
    M1 = -(-SwxSq + Swx2*Sw + SwySq - Swy2*Sw - C)/(2.0*D);
    B1 = -(-Swy+M1*Swx)/Sw;

    // This second solution is perpendicular to the first and represents
    // the worst fit
    // M2 = -(-SwxSq + Swx2*Sw + SwySq - Swy2*Sw + C)/(2.0*D);
    // B2 = -(-Swy+M2*Swx)/Sw;

    *x1 = xlo;
    *y1 = M1*xlo + B1;
    *x2 = xhi;
    *y2 = M1*xhi + B1;

    // Compute the chi^2 function for return
    double chi2 = 0.0;
    for (i=0; i<w.size(); ++i) {
      chi2 += w[i]*geo::DsqrToLine(x[i],y[i],*x1,*y1,*x2, *y2);
    }

    // test if linfit is better

    double chi2lin, linx1, linx2, liny1, liny2;

    chi2lin = LinFit(x,y,w,&linx1,&liny1,&linx2,&liny2);

    if(chi2lin < chi2) { //something is wrong with the perp fit, use the lin fit
      *x1 = linx1;
      *y1 = liny1;
      *x2 = linx2;
      *y2 = liny2;

      chi2 = chi2lin;
    }


    return chi2;
  }

  //......................................................................
  /// Find the distance from the given point to the edge of the detector
  double DistToEdge(double *point,
                    double detHalfWidth,
                    double detHalfHeight,
                    double detLength)
  {
    // Distance to an edge in each direction
    const double dx = fabs(detHalfWidth  - fabs(point[0]));
    const double dy = fabs(detHalfHeight - fabs(point[1]));
    const double dz = std::min(fabs(point[2]), fabs(detLength-point[2]));

    return std::min(std::min(dx, dy), dz);
  }

  //......................................................................
  bool ClampRayToDetectorVolume(TVector3* p0, TVector3* p1,
                                const GeometryBase* geom)
  {
    const double Lx = geom->DetHalfWidth();
    const double Ly = geom->DetHalfHeight();
    const double Lz = geom->DetLength()/2;

    // Collect points in an array so we can loop
    TVector3* ps[2] = {p0, p1};

    for(int n = 0; n < 2; ++n){
      // Unit vector towards the other point
      const TVector3 dp = (*ps[1-n]-*ps[n]).Unit();

      // Distance to detector edge along dp for each view. Zero if inside
      const double distx = dp.X() ? (fabs(ps[n]->X()   )-Lx)/fabs(dp.X()) : 0;
      const double disty = dp.Y() ? (fabs(ps[n]->Y()   )-Ly)/fabs(dp.Y()) : 0;
      const double distz = dp.Z() ? (fabs(ps[n]->Z()-Lz)-Lz)/fabs(dp.Z()) : 0;
      // Distance required to be inside the detector. Zero if already contained
      const double dist = std::max(std::max(distx, disty), distz);

      if(dist > 0){
        // Here we're assuming that the reconstructed track already goes
        // through the detector, so that by moving towards the other point
        // we'll hit it.
        *ps[n] += dist*dp;
      }

      // If that assumption was false, then we won't have ended up on the edge
      // of the detector as we intended.
      if(fabs(ps[n]->X()) > Lx+.1 ||
         fabs(ps[n]->Y()) > Ly+.1 ||
         ps[n]->Z() < -.1 ||
         ps[n]->Z() > 2*Lz+.1){
        return false;
      }
    } // end for n
    return true;
  }

  // Table for GetAvgPath2D. Filled by InitializePathTable().
  // Floating point key is OK, because we only ever interpolate from this
  // table, never try to look things up directly.
  static std::map<double, double> gPathTable;

  //......................................................................
  /// \brief Function to initialize \ref gPathTable for \ref GetAvgPath2D
  ///
  /// Trace paths through a 2D cell at all angles and all positions. Average
  /// over the positions.
  void InitializePathTable()
  {
    // Don't call me twice
    if(!gPathTable.empty()) 
      throw cet::exception("Geo") << "trying to non-empty initialize path table.";

    // Numbers off the ends don't mean anything, but they keep the
    // interpolation happy for borderline cases
    for(double dx = -.01; dx < 1.02; dx += .02){
      const double dy = sqrt(std::max(1-dx*dx, 0.));

      // Dimensions of a cell, in cm. Convention is that x is the long axis
      // and y the short axis.
      // TODO: harcoding numbers is bad
      const double Lx = 5.94;
      const double Ly = 3.84;

      // Accumulate into S_tot. Count how many samples we took.
      double S_tot = 0;
      int samples = 0;

      // Step size for the integration, in cm
      const double step = .005;

      // This is always a big enough line to integrate over that we'll cover
      // the whole cell. Starting at zero and only doing half is justified
      // by the symmetry.
      const double U = util::pythag(Lx, Ly)/2;
      for(double u = 0; u < U; u += step){
        // There really should be exactly zero or two intersection times.
        // We'll assert that in a minute. But avoid potential memory stomping.
        double tInts[4];
        // The place where we'll write the first time
        double* tInt = tInts;

        // Integrate over points along a vector perpendicular to the track
        const double x0 = -dy*u;
        const double y0 = +dx*u;

        double x, y, t;

        // When do we intersect the +x wall? Where are we in y at the time?
        x = +Lx/2;
        t  = (x-x0)/dx;
        y = y0+t*dy;
        // If this is within the cell: record it, and move onto the next slot
        if(y > -Ly/2 && y < +Ly/2) *tInt++ = t;

        y = +Ly/2;
        t = (y-y0)/dy;
        x = x0+t*dx;
        if(x > -Lx/2 && x < +Lx/2) *tInt++ = t;

        x = -Lx/2;
        t  = (x-x0)/dx;
        y = y0+t*dy;
        if(y > -Ly/2 && y < +Ly/2) *tInt++ = t;

        y = -Ly/2;
        t  = (y-y0)/dy;
        x = x0+t*dx;
        if(x > -Lx/2 && x < +Lx/2) *tInt++ = t;

        // Either the ray enters and exits, or it never hits the cell at all
        if(tInt - tInts != 2 && tInt - tInts != 0)
          throw cet::exception("Geo") << "ray enters the cell but does not exit, that is a problem";

        // In this case we missed
        if(tInt == tInts) continue;

        // The direction is a unit vector, so the distance is just the times
        S_tot += fabs(tInts[1]-tInts[0]);
        ++samples;
      } // end for u

      gPathTable[dx] = S_tot/samples;
    } // end for dx
  }

  // Need this to look up the position of a key in a list of (key,value) pairs
  bool ltKey(std::pair<double, double> a, double b){return a.first < b;}

  //......................................................................
  /// \brief Helper for \ref AverageCellPathLength
  ///
  /// Look up average length for a 2D ray in the table. dx is the component in
  /// the long direction of the cell (detector dz).
  double GetAvgPath2D(double dx)
  {
    if(gPathTable.empty()) InitializePathTable();

    dx = std::abs(dx);

    typedef std::map<double, double>::const_iterator it_t;

    // Find the point after dx
    it_t itNext = std::lower_bound(gPathTable.begin(), gPathTable.end(), dx,
                                   ltKey);

    if(itNext == gPathTable.end() || itNext == gPathTable.begin()) 
      throw cet::exception("Geo") << "cannot find the point after dx: " << dx;

    // And the point before
    it_t itPrev = itNext; --itPrev;

    // Interpolate
    const double x0 = itPrev->first; const double y0 = itPrev->second;
    const double x1 = itNext->first; const double y1 = itNext->second;
    return ((x1-dx)*y0+(dx-x0)*y1)/(x1-x0);
  }

  //......................................................................
  /// \brief Mean path length of a ray with (unit) direction vector dx, dy, dz
  /// through a cell in \a view, averaged over all transverse positions
  double AverageCellPathLength(geo::View_t view,
                               double dx, double dy, double dz)
  {
    // Must be a unit vector
    if(std::abs(util::pythag(dx, dy, dz)-1) > 1e-4)
      throw cet::exception("Geo") << "provided values are not a unit vector - "
                                  << "(dx, dy, dz) = (" << dx << " ," << dy << " ,"
                                  << dz << ")";

    if(view == geo::kX){
      // Straight down the cell, hardcode NDOS length...
      if(std::abs(dy) > .99) return 393;
      // Project down into xz 2D
      const double scale = sqrt(1-dy*dy);
      const double dist2D = GetAvgPath2D(dz/scale);
      // And back up to 3D
      return dist2D/scale;
    }
    else{
      if(std::abs(dx) > .99) return 262;
      const double scale = sqrt(1-dx*dx);
      const double dist2D = GetAvgPath2D(dz/scale);
      return dist2D/scale;
    }
  }

} // namespace geo
////////////////////////////////////////////////////////////////////////