Geo.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file geo.h
3 /// \brief Supply basic geometry functions
4 /// \author messier@indiana.edu
5 /// \version
6 ////////////////////////////////////////////////////////////////////////
7 #include "GeometryObjects/Geo.h"
8 
9 #include <cmath>
10 #include <cassert>
11 #include <cfloat>
12 #include <iostream>
13 #include "TVector3.h"
14 
16 #include "Utilities/func/MathUtil.h"
17 #include "cetlib_except/exception.h"
18 
19 namespace geo {
20 
21  //......................................................................
22  /// \brief Project along a direction from a particular starting point to the
23  /// edge of a box
24  ///
25  /// \param xyz - The starting x,y,z location. Must be inside box.
26  /// \param dxyz - Direction vector
27  /// \param xlo - Low edge of box in x
28  /// \param xhi - Low edge of box in x
29  /// \param ylo - Low edge of box in y
30  /// \param yhi - Low edge of box in y
31  /// \param zlo - Low edge of box in z
32  /// \param zhi - Low edge of box in z
33  /// \param xyzout - On output, the position at the box edge
34  ///
35  /// Note: It should be safe to use the same array for input and
36  /// output.
37  ///
38  void ProjectToBoxEdge(const double xyz[],
39  const double dxyz[],
40  double xlo, double xhi,
41  double ylo, double yhi,
42  double zlo, double zhi,
43  double xyzout[])
44  {
45  // Make sure we're inside the box!
46 
47  if(xyz[0] < xlo || xyz[0] > xhi ||
48  xyz[1] < ylo || xyz[1] > yhi ||
49  xyz[2] < zlo || xyz[2] > zhi)
50  throw cet::exception("Geo") << "provided point is outside of the specified box";
51 
52  // Compute the distances to the x/y/z walls
53  double dx = 99.E99;
54  double dy = 99.E99;
55  double dz = 99.E99;
56  if (dxyz[0]>0.0) { dx = (xhi-xyz[0])/dxyz[0]; }
57  else if (dxyz[0]<0.0) { dx = (xlo-xyz[0])/dxyz[0]; }
58  if (dxyz[1]>0.0) { dy = (yhi-xyz[1])/dxyz[1]; }
59  else if (dxyz[1]<0.0) { dy = (ylo-xyz[1])/dxyz[1]; }
60  if (dxyz[2]>0.0) { dz = (zhi-xyz[2])/dxyz[2]; }
61  else if (dxyz[2]<0.0) { dz = (zlo-xyz[2])/dxyz[2]; }
62 
63  // Choose the shortest distance
64  double d = 0.0;
65  if (dx<dy && dx<dz) d = dx;
66  else if (dy<dz && dy<dx) d = dy;
67  else if (dz<dx && dz<dy) d = dz;
68 
69  // Make the step
70  for (int i=0; i<3; ++i) {
71  xyzout[i] = xyz[i] + dxyz[i]*d;
72  }
73  }
74 
75  int WhichWallofBox(const double xyz[],
76  double xlo, double xhi,
77  double ylo, double yhi,
78  double zlo, double zhi)
79  {
80  int wall = 0;
81  if(xyz[0] >= xlo-0.01 && xyz[0] <= xlo+0.01) wall = 1;
82  if(xyz[0] >= xhi-0.01 && xyz[0] <= xhi+0.01) wall = 2;
83  if(xyz[1] >= ylo-0.01 && xyz[1] <= ylo+0.01) wall = 3;
84  if(xyz[1] >= yhi-0.01 && xyz[1] <= yhi+0.01) wall = 4;
85  if(xyz[2] >= zlo-0.01 && xyz[2] <= zlo+0.01) wall = 5;
86  if(xyz[2] >= zhi-0.01 && xyz[2] <= zhi+0.01) wall = 6;
87 
88 // if(wall == 0) std::cout << xyz[0] << " " << xyz[1] << " " << xyz[2] << std::endl;
89 
90  return wall;
91 
92  }
93 
94  //......................................................................
95  /// \brief Determine if a particle starting at xyz with direction dxyz will
96  /// intersect a box defined by xlo, xhi, ylo, yhi, zlo, zhi
97  bool IntersectsBox(const double xyz[],
98  const double dxyz[],
99  double xlo, double xhi,
100  double ylo, double yhi,
101  double zlo, double zhi)
102  {
103  double xyzedge[3] = {0.};
104 
105  // there are 6 edges to a box, check them all. If a particle intersects
106  // any surface, then return true
107 
108  // top surface, y = yhi
109  geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 1, yhi, xyzedge);
110  if(xyzedge[0] > xlo && xyzedge[0] < xhi
111  && xyzedge[2] > zlo && xyzedge[2] < zhi) return true;
112 
113  // bottom surface, y = ylo
114  geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 1, ylo, xyzedge);
115  if(xyzedge[0] > xlo && xyzedge[0] < xhi
116  && xyzedge[2] > zlo && xyzedge[2] < zhi) return true;
117 
118  // x = xhi
119  geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 0, xhi, xyzedge);
120  if(xyzedge[1] > ylo && xyzedge[1] < yhi
121  && xyzedge[2] > zlo && xyzedge[2] < zhi) return true;
122 
123  // x = xlo
124  geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 0, xlo, xyzedge);
125  if(xyzedge[1] > ylo && xyzedge[1] < yhi
126  && xyzedge[2] > zlo && xyzedge[2] < zhi) return true;
127 
128  // z = zhi
129  geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 2, zhi, xyzedge);
130  if(xyzedge[1] > ylo && xyzedge[1] < yhi
131  && xyzedge[0] > xlo && xyzedge[0] < xhi) return true;
132 
133  // z = zlo
134  geo::ProjectToBoxEdgeFromOutside(xyz, dxyz, 2, zlo, xyzedge);
135  if(xyzedge[1] > ylo && xyzedge[1] < yhi
136  && xyzedge[0] > xlo && xyzedge[0] < xhi) return true;
137 
138  return false;
139 
140  }
141 
142  //......................................................................
143  /// \brief Project from a position outside of a box to an edge of the box
144  /// with coordinate value \a edge for the axis \a axis.
145  ///
146  /// axis = 0 --> x, axis = 1 --> y, axis = 2 --> z
147  /// \a xyzout is the position on the desired edge
148  void ProjectToBoxEdgeFromOutside(const double xyz[],
149  const double dxyz[],
150  int axis,
151  double edge,
152  double xyzout[])
153  {
154  double momDir = dxyz[axis];
155  double posDir = xyz[axis];
156  double momTot = sqrt(dxyz[0]*dxyz[0] + dxyz[1]*dxyz[1] + dxyz[2]*dxyz[2]);
157 
158  if(momTot < 1.e-9){
159  std::cerr << "no direction information supplied, do nothing" << std::endl;
160  return;
161  }
162 
163  double ddS = momDir/momTot; // direction cosine in the desired direction
164  double length1Dim = (edge - posDir); // length in desired direction to that edge;
165 
166  if(TMath::Abs(ddS) > 0.){
167  length1Dim /= ddS;
168  xyzout[0] = xyz[0] + length1Dim*dxyz[0]/momTot;
169  xyzout[1] = xyz[1] + length1Dim*dxyz[1]/momTot;
170  xyzout[2] = xyz[2] + length1Dim*dxyz[2]/momTot;
171  }
172 
173  return;
174  }
175 
176  //......................................................................
177  /// \brief Find the intersection between two line segments
178  ///
179  /// Given the two line segments (x0,y0)-(x1,y1), and (X0,Y0)-(X1,Y1), fills
180  /// (x,y) with the point that lies on both of them. The return value
181  /// indicates if the intersection of the (infinite) lines falls inside the
182  /// line segments. If it does not, the returned (x,y) is still meaningful,
183  /// unless the two input lines are parallel.
184  bool LineIntersection(double x0, double y0, double x1, double y1,
185  double X0, double Y0, double X1, double Y1,
186  double& x, double& y)
187  {
188  const double dx = x1-x0;
189  const double dy = y1-y0;
190  const double dX = X1-X0;
191  const double dY = Y1-Y0;
192 
193  const double denom = dX*dy-dY*dx;
194 
195  if(denom == 0){
196  // std::cerr << "Trying to find intercept of parallel lines" << std::endl;
197  x = x1; y = y1;
198  return false;
199  }
200 
201  // Distance of intercept point along each of the two lines
202  const double l = -(X0*dY-Y0*dX+dX*y0-dY*x0)/denom;
203  const double L = +(x0*dy-y0*dx+dx*Y0-dy*X0)/denom;
204 
205  x = X0+L*(X1-X0);
206  y = Y0+L*(Y1-Y0);
207 
208  return l >= 0 && l <= 1 && L >= 0 && L <= 1;
209  }
210 
211  //......................................................................
212  ///
213  /// \brief Find the distance of closest approach between point and line
214  ///
215  /// \param point - xyz coordinates of point
216  /// \param intercept - xyz coordinates of point on line
217  /// \param slopes - direction vector (need not be normalized)
218  /// \param closest - on output, point on line that is closest
219  ///
220  /// \returns distance from point to line
221  ///
222  double ClosestApproach(const double point[],
223  const double intercept[],
224  const double slopes[],
225  double closest[])
226  {
227  TVector3 closest_vec;
228  const double dist = ClosestApproach(TVector3(point),
229  TVector3(intercept), TVector3(slopes),
230  closest_vec);
231  closest_vec.GetXYZ(closest);
232  return dist;
233  }
234 
235  //......................................................................
236  ///
237  /// \brief Find the distance of closest approach between point and line
238  ///
239  /// \param point - xyz coordinates of point
240  /// \param intercept - xyz coordinates of point on line
241  /// \param slopes - direction vector (need not be normalized)
242  /// \param closest - on output, point on line that is closest
243  ///
244  /// \returns distance from point to line
245  ///
246  double ClosestApproach(TVector3 point, TVector3 intercept,
247  TVector3 slopes, TVector3& closest)
248  {
249  const double s = slopes.Dot(point-intercept);
250  const double sd = slopes.Mag2();
251 
252  if(sd > 0){
253  closest = intercept + (s/sd)*slopes;
254  }
255  else {
256  // How to handle this zero gracefully? Assume that the intercept
257  // is a particle vertex and "slopes" are momenta. In that case,
258  // the particle goes nowhere and the closest approach is the
259  // distance from the intercept to point
260  closest = intercept;
261  }
262 
263  return (point-closest).Mag();
264  }
265 
266  //....................................................................
267  /// \brief Find the distance of closest approach between two lines
268  /// which pass through points P0 and P1 and Q0 and Q1
269  ///
270  /// @param P0 - A point on line P
271  /// @param P1 - A second point on line P
272  /// @param Q0 - A point on line Q
273  /// @param Q1 - A second point on line Q
274  /// @param sc - distance from P0 to point on P which is closest
275  /// @param tc - distance from Q0 to point on Q which is closest
276  /// @param PC - point on line P which is closest
277  /// @param QC - point on line Q which is closest
278  ///
279  /// @returns The distance^2 between the points of closest approach.
280  ///
281  /// If any of the pointers sc, tc, PC, QC are zero, that part of the
282  /// calculation is skipped.
283  ///
284  double ClosestApproach(const TVector3& P0,
285  const TVector3& P1,
286  const TVector3& Q0,
287  const TVector3& Q1,
288  double* sc,
289  double* tc,
290  TVector3* PC,
291  TVector3* QC)
292  {
293  TVector3 w0(P0); w0 -= Q0;
294 
295  TVector3 u(P1); u -= P0; u = u.Unit();
296  TVector3 v(Q1); v -= Q0; v = v.Unit();
297 
298  double b = u.Dot(v);
299  double d = u.Dot(w0);
300  double e = v.Dot(w0);
301 
302  double denom = 1-b*b;
303  double sC, tC;
304  if (denom!=0.0) {
305  sC = (b*e-d)/denom;
306  tC = (e-b*d)/denom;
307  }
308  else {
309  sC = 0.0;
310  tC = d/b;
311  }
312  if (sc!=0) *sc = sC;
313  if (tc!=0) *tc = tC;
314 
315  TVector3 pC(P0+sC*u);
316  TVector3 qC(Q0+tC*v);
317  if (PC!=0) *PC = pC;
318  if (QC!=0) *QC = qC;
319 
320  TVector3 D(pC-qC);
321 
322  return D.Mag2();
323  }
324 
325  //......................................................................
326  ///
327  /// \brief In two dimensions, return the perpendicular distance from a
328  /// point (x0,y0) to a line defined by end points (x1,y1) and (x2,y2)
329  ///
330  /// @param x0 : x-value of point
331  /// @param y0 : y-value of point
332  /// @param x1 : x-value of point at one end of line
333  /// @param y1 : y-value of point at one end of line
334  /// @param x2 : x-value of point at other end of line
335  /// @param y2 : y-value of point at other end of line
336  /// @returns Squared distance from the point to the line
337  ///
338  double DsqrToLine(double x0, double y0,
339  double x1, double y1,
340  double x2, double y2)
341  {
342  using util::sqr;
343 
344  double d = 0.;
345 
346  // If we've been handed a point instead of a line, return distance
347  // to the point
348  if (x1==x2 && y1==y2) d = sqr(x0-x1)+sqr(y0-y1);
349  else d = sqr((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1)) / (sqr(x2-x1)+sqr(y2-y1));
350 
351  return d;
352 
353  }
354 
355  //......................................................................
356  ///
357  /// \brief Find the best-fit line to a collection of points in 2-D by
358  /// minimizing the squared vertical distance from the points to the line.
359  ///
360  /// To give the most robust fits and avoid cases where the slopes approach
361  /// infinity, the points are arranged to maximize the horizontal extent. So,
362  /// the fitter may fit x vs. y or y vs. x.
363  ///
364  /// @param x - input vector of x coordinates
365  /// @param y - input vector of y coordinates
366  /// @param w - input vector of weights for the points
367  /// @param x1 - output x coordinate at one end of line
368  /// @param y1 - output y coordinate at one end of line
369  /// @param x2 - output x coordinate at other end of line
370  /// @param y2 - output y coordinate at other end of line
371  /// @returns The chi^2 value defined by chi^2 = sum_i[w_i d^2_i]
372  ///
373  double LinFit(const std::vector<double>& x,
374  const std::vector<double>& y,
375  const std::vector<double>& w,
376  double* x1, double* y1,
377  double* x2, double* y2)
378  {
379  unsigned int i;
380 
381  // Before going ahead, make sure we have sensible arrays
382  if( x.size() != y.size() ||
383  x.size() != w.size() ||
384  x.size() < 2)
385  throw cet::exception("Geo") << "Input vectors to LinFit are not the "
386  << "correct length. x size: " << x.size()
387  << " y size: " << y.size()
388  << " w size: " << w.size();
389 
390  // Find the ranges of the input variables
391  double xlo = x[0];
392  double xhi = x[0];
393  double ylo = y[0];
394  double yhi = y[0];
395  for (i=1; i<w.size(); ++i) {
396  if (x[i]<xlo) xlo = x[i];
397  if (x[i]>xhi) xhi = x[i];
398  if (y[i]<ylo) ylo = y[i];
399  if (y[i]>yhi) yhi = y[i];
400  }
401  // Make sure the fit has some span
402  if (yhi-ylo==0.0 && xhi-xlo==0.0) {
403  std::cerr << __FILE__ << ":" << __LINE__
404  << " All points identical in line fit!" << std::endl;
405  *x1 = xlo;
406  *y1 = ylo;
407  *x2 = xhi;
408  *y2 = yhi;
409  return 0.0;
410  }
411 
412  //
413  // If the y extent is larger than the x extent, flip the variables
414  // and fit
415  //
416  if (yhi-ylo > xhi-xlo) {
417  return geo::LinFit(y,x,w,y1,x1,y2,x2);
418  }
419 
420  double m, b;
421  const double chi2 = util::LinFit(x, y, w, m, b);
422 
423  *x1 = xlo;
424  *x2 = xhi;
425  *y1 = m*(*x1)+b;
426  *y2 = m*(*x2)+b;
427 
428  return chi2;
429  }
430 
431  //......................................................................
432  ///
433  /// \brief Find the best-fit line to a collection of points in 2-D by
434  /// minimizing the squared perpendicular distance from the points to the
435  /// line.
436  ///
437  /// This is not the usual "linfit" which minimizes the squared
438  /// vertical distance
439  ///
440  /// @param x - input vector of x coordinates
441  /// @param y - input vector of y coordinates
442  /// @param w - input vector of weights for the points
443  /// @param x1 - output x coordinate at one end of line
444  /// @param y1 - output y coordinate at one end of line
445  /// @param x2 - output x coordinate at other end of line
446  /// @param y2 - output y coordinate at other end of line
447  /// @returns The chi^2 value defined by chi^2 = sum_i[w_i d^2_i]
448  ///
449  double LinFitMinDperp(const std::vector<double>& x,
450  const std::vector<double>& y,
451  const std::vector<double>& w,
452  double* x1, double* y1,
453  double* x2, double* y2)
454  {
455  // Before going ahead, make sure we have sensible arrays
456  if( x.size() != y.size() ||
457  x.size() != w.size() ||
458  x.size() < 2)
459  throw cet::exception("Geo") << "Input vectors to LinFit are not the "
460  << "correct length. x size: " << x.size()
461  << " y size: " << y.size()
462  << " w size: " << w.size();
463 
464  //--
465  // The following code implements the output of the following MAPLE
466  // script:
467  //--
468  // assume(x0,real);
469  // assume(y0,real);
470  // assume({M,B},real);
471  // assume({Dsqr},real);
472  // assume(Sw,real);
473  // assume(Swx,real);
474  // assume(Swx2,real);
475  // assume(Swy,real);
476  // assume(Swy2,real);
477  // assume(Swxy,real);
478  // y1 := M*x1 + B;
479  // y2 := M*x2 + B;
480  // Dsqr := (((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1))^2/((x2-x1)^2+(y2-y1)^2));
481  // eqn1 := simplify(expand(diff(Dsqr,M) = 0));
482  // eqn2 := simplify(expand(diff(Dsqr,B) = 0));
483  // Eqn1 := -2*B*Swy*M +B^2*M*Sw +Swy2*M -Swxy*M^2 -B*Swx +Swxy
484  // +B*Swx*M^2 -Swx2*M = 0;
485  // Eqn2 := B*Sw - Swy + M*Swx = 0;
486  // solve({Eqn1,Eqn2},{M,B});
487  // allvalues(%);
488  //--
489 
490  // register is deprecated in c++17 standard
491  //unsigned register int i;
492  unsigned int i;
493 
494  // Find the extremes of the inputs
495  double xlo = DBL_MAX;
496  double xhi = -DBL_MAX;
497  double ylo = DBL_MAX;
498  double yhi = -DBL_MAX;
499  for (i=0; i<w.size(); ++i) {
500  if(w[i] == 0) continue;
501  if (x[i]<xlo) xlo = x[i];
502  if (y[i]<ylo) ylo = y[i];
503  if (x[i]>xhi) xhi = x[i];
504  if (y[i]>yhi) yhi = y[i];
505  }
506 
507  // This algorithm really fails if the y values are identical
508  // return the low and hi values if it's too flat
509  if(fabs(yhi-ylo) < 0.1){
510  *x1 = xlo;
511  *y1 = ylo;
512  *x2 = xhi;
513  *y2 = yhi;
514  return 0.; //its a flat line, chi^2 is 0 by definition
515  }
516 
517  // For stability, fit in a way that avoids infinite slopes,
518  // exchanging the x and y variables as needed.
519  if (fabs(xlo - xhi) < 1.e-3) {
520  double chi2;
521  double xx1, yy1, xx2, yy2;
522  // Notice: (x,y) -> (y,x)
523  chi2 = LinFitMinDperp(y,x,w,&xx1, &yy1,&xx2, &yy2);
524  *x1 = xlo; // Infinite slope means that your initial x positions should be the input
525  *y1 = xx1;
526  *x2 = xhi;
527  *y2 = xx2;
528  return chi2;
529  }
530 
531  // Accumulate the sums needed to compute the fit line
532  double Sw = 0.0;
533  double Swx = 0.0;
534  double Swy = 0.0;
535  double Swxy= 0.0;
536  double Swy2= 0.0;
537  double Swx2= 0.0;
538  for (i=0; i<w.size(); ++i) {
539  Sw += w[i];
540  Swx += w[i]*x[i];
541  Swy += w[i]*y[i];
542  Swxy += w[i]*x[i]*y[i];
543  Swy2 += w[i]*y[i]*y[i];
544  Swx2 += w[i]*x[i]*x[i];
545  }
546 
547  if(Sw == 0.0)
548  throw cet::exception("Geo") << "sum of weights is 0 in LinFitMinDperp";
549 
550  // Precalculate reused squares
551  using util::sqr;
552  const double SwSq = sqr(Sw);
553  const double SwxSq = sqr(Swx);
554  const double SwySq = sqr(Swy);
555 
556  // C and D are two handy temporary variables
557  double C =
558  +sqr(SwxSq) // Swx^4
559  -2.0*SwxSq*Swx2*Sw
560  +2.0*SwxSq*SwySq
561  +2.0*SwxSq*Swy2*Sw
562  +sqr(Swx2)*SwSq
563  +2.0*Swx2*Sw*SwySq
564  -2.0*Swx2*SwSq*Swy2
565  +sqr(sqr(Swy)) // Swy^4
566  -2.0*SwySq*Swy2*Sw
567  +sqr(Swy2)*SwSq
568  -8.0*Swx*Swy*Swxy*Sw
569  +4.0*sqr(Swxy)*SwSq;
570  if (C>0.0) C = sqrt(C);
571  else C = 0.0;
572  double D = -Swx*Swy+Swxy*Sw;
573 
574  // The first of the two solutions. This one is the best fit through
575  // the points
576  double M1;
577  double B1;
578  if (D==0.0) {
579  // D could be zero, this might happen in a segmented detector if
580  // the track fit passed through only a single readout plane and
581  // all hits share the same x coordinates. I've tried to avoid
582  // this by checking if it makes more sense to fit with x and y
583  // exchanged above. This is an extra precaution.
584  D = 1.0E-9;
585  }
586  M1 = -(-SwxSq + Swx2*Sw + SwySq - Swy2*Sw - C)/(2.0*D);
587  B1 = -(-Swy+M1*Swx)/Sw;
588 
589  // This second solution is perpendicular to the first and represents
590  // the worst fit
591  // M2 = -(-SwxSq + Swx2*Sw + SwySq - Swy2*Sw + C)/(2.0*D);
592  // B2 = -(-Swy+M2*Swx)/Sw;
593 
594  *x1 = xlo;
595  *y1 = M1*xlo + B1;
596  *x2 = xhi;
597  *y2 = M1*xhi + B1;
598 
599  // Compute the chi^2 function for return
600  double chi2 = 0.0;
601  for (i=0; i<w.size(); ++i) {
602  chi2 += w[i]*geo::DsqrToLine(x[i],y[i],*x1,*y1,*x2, *y2);
603  }
604 
605  // test if linfit is better
606 
607  double chi2lin, linx1, linx2, liny1, liny2;
608 
609  chi2lin = LinFit(x,y,w,&linx1,&liny1,&linx2,&liny2);
610 
611  if(chi2lin < chi2) { //something is wrong with the perp fit, use the lin fit
612  *x1 = linx1;
613  *y1 = liny1;
614  *x2 = linx2;
615  *y2 = liny2;
616 
617  chi2 = chi2lin;
618  }
619 
620 
621  return chi2;
622  }
623 
624  //......................................................................
625  /// Find the distance from the given point to the edge of the detector
626  double DistToEdge(double *point,
627  double detHalfWidth,
628  double detHalfHeight,
629  double detLength)
630  {
631  // Distance to an edge in each direction
632  const double dx = fabs(detHalfWidth - fabs(point[0]));
633  const double dy = fabs(detHalfHeight - fabs(point[1]));
634  const double dz = std::min(fabs(point[2]), fabs(detLength-point[2]));
635 
636  return std::min(std::min(dx, dy), dz);
637  }
638 
639  //......................................................................
640  bool ClampRayToDetectorVolume(TVector3* p0, TVector3* p1,
641  const GeometryBase* geom)
642  {
643  const double Lx = geom->DetHalfWidth();
644  const double Ly = geom->DetHalfHeight();
645  const double Lz = geom->DetLength()/2;
646 
647  // Collect points in an array so we can loop
648  TVector3* ps[2] = {p0, p1};
649 
650  for(int n = 0; n < 2; ++n){
651  // Unit vector towards the other point
652  const TVector3 dp = (*ps[1-n]-*ps[n]).Unit();
653 
654  // Distance to detector edge along dp for each view. Zero if inside
655  const double distx = dp.X() ? (fabs(ps[n]->X() )-Lx)/fabs(dp.X()) : 0;
656  const double disty = dp.Y() ? (fabs(ps[n]->Y() )-Ly)/fabs(dp.Y()) : 0;
657  const double distz = dp.Z() ? (fabs(ps[n]->Z()-Lz)-Lz)/fabs(dp.Z()) : 0;
658  // Distance required to be inside the detector. Zero if already contained
659  const double dist = std::max(std::max(distx, disty), distz);
660 
661  if(dist > 0){
662  // Here we're assuming that the reconstructed track already goes
663  // through the detector, so that by moving towards the other point
664  // we'll hit it.
665  *ps[n] += dist*dp;
666  }
667 
668  // If that assumption was false, then we won't have ended up on the edge
669  // of the detector as we intended.
670  if(fabs(ps[n]->X()) > Lx+.1 ||
671  fabs(ps[n]->Y()) > Ly+.1 ||
672  ps[n]->Z() < -.1 ||
673  ps[n]->Z() > 2*Lz+.1){
674  return false;
675  }
676  } // end for n
677  return true;
678  }
679 
680  // Table for GetAvgPath2D. Filled by InitializePathTable().
681  // Floating point key is OK, because we only ever interpolate from this
682  // table, never try to look things up directly.
683  static std::map<double, double> gPathTable;
684 
685  //......................................................................
686  /// \brief Function to initialize \ref gPathTable for \ref GetAvgPath2D
687  ///
688  /// Trace paths through a 2D cell at all angles and all positions. Average
689  /// over the positions.
691  {
692  // Don't call me twice
693  if(!gPathTable.empty())
694  throw cet::exception("Geo") << "trying to non-empty initialize path table.";
695 
696  // Numbers off the ends don't mean anything, but they keep the
697  // interpolation happy for borderline cases
698  for(double dx = -.01; dx < 1.02; dx += .02){
699  const double dy = sqrt(std::max(1-dx*dx, 0.));
700 
701  // Dimensions of a cell, in cm. Convention is that x is the long axis
702  // and y the short axis.
703  // TODO: harcoding numbers is bad
704  const double Lx = 5.94;
705  const double Ly = 3.84;
706 
707  // Accumulate into S_tot. Count how many samples we took.
708  double S_tot = 0;
709  int samples = 0;
710 
711  // Step size for the integration, in cm
712  const double step = .005;
713 
714  // This is always a big enough line to integrate over that we'll cover
715  // the whole cell. Starting at zero and only doing half is justified
716  // by the symmetry.
717  const double U = util::pythag(Lx, Ly)/2;
718  for(double u = 0; u < U; u += step){
719  // There really should be exactly zero or two intersection times.
720  // We'll assert that in a minute. But avoid potential memory stomping.
721  double tInts[4];
722  // The place where we'll write the first time
723  double* tInt = tInts;
724 
725  // Integrate over points along a vector perpendicular to the track
726  const double x0 = -dy*u;
727  const double y0 = +dx*u;
728 
729  double x, y, t;
730 
731  // When do we intersect the +x wall? Where are we in y at the time?
732  x = +Lx/2;
733  t = (x-x0)/dx;
734  y = y0+t*dy;
735  // If this is within the cell: record it, and move onto the next slot
736  if(y > -Ly/2 && y < +Ly/2) *tInt++ = t;
737 
738  y = +Ly/2;
739  t = (y-y0)/dy;
740  x = x0+t*dx;
741  if(x > -Lx/2 && x < +Lx/2) *tInt++ = t;
742 
743  x = -Lx/2;
744  t = (x-x0)/dx;
745  y = y0+t*dy;
746  if(y > -Ly/2 && y < +Ly/2) *tInt++ = t;
747 
748  y = -Ly/2;
749  t = (y-y0)/dy;
750  x = x0+t*dx;
751  if(x > -Lx/2 && x < +Lx/2) *tInt++ = t;
752 
753  // Either the ray enters and exits, or it never hits the cell at all
754  if(tInt - tInts != 2 && tInt - tInts != 0)
755  throw cet::exception("Geo") << "ray enters the cell but does not exit, that is a problem";
756 
757  // In this case we missed
758  if(tInt == tInts) continue;
759 
760  // The direction is a unit vector, so the distance is just the times
761  S_tot += fabs(tInts[1]-tInts[0]);
762  ++samples;
763  } // end for u
764 
765  gPathTable[dx] = S_tot/samples;
766  } // end for dx
767  }
768 
769  // Need this to look up the position of a key in a list of (key,value) pairs
770  bool ltKey(std::pair<double, double> a, double b){return a.first < b;}
771 
772  //......................................................................
773  /// \brief Helper for \ref AverageCellPathLength
774  ///
775  /// Look up average length for a 2D ray in the table. dx is the component in
776  /// the long direction of the cell (detector dz).
777  double GetAvgPath2D(double dx)
778  {
779  if(gPathTable.empty()) InitializePathTable();
780 
781  dx = std::abs(dx);
782 
783  typedef std::map<double, double>::const_iterator it_t;
784 
785  // Find the point after dx
786  it_t itNext = std::lower_bound(gPathTable.begin(), gPathTable.end(), dx,
787  ltKey);
788 
789  if(itNext == gPathTable.end() || itNext == gPathTable.begin())
790  throw cet::exception("Geo") << "cannot find the point after dx: " << dx;
791 
792  // And the point before
793  it_t itPrev = itNext; --itPrev;
794 
795  // Interpolate
796  const double x0 = itPrev->first; const double y0 = itPrev->second;
797  const double x1 = itNext->first; const double y1 = itNext->second;
798  return ((x1-dx)*y0+(dx-x0)*y1)/(x1-x0);
799  }
800 
801  //......................................................................
802  /// \brief Mean path length of a ray with (unit) direction vector dx, dy, dz
803  /// through a cell in \a view, averaged over all transverse positions
805  double dx, double dy, double dz)
806  {
807  // Must be a unit vector
808  if(std::abs(util::pythag(dx, dy, dz)-1) > 1e-4)
809  throw cet::exception("Geo") << "provided values are not a unit vector - "
810  << "(dx, dy, dz) = (" << dx << " ," << dy << " ,"
811  << dz << ")";
812 
813  if(view == geo::kX){
814  // Straight down the cell, hardcode NDOS length...
815  if(std::abs(dy) > .99) return 393;
816  // Project down into xz 2D
817  const double scale = sqrt(1-dy*dy);
818  const double dist2D = GetAvgPath2D(dz/scale);
819  // And back up to 3D
820  return dist2D/scale;
821  }
822  else{
823  if(std::abs(dx) > .99) return 262;
824  const double scale = sqrt(1-dx*dx);
825  const double dist2D = GetAvgPath2D(dz/scale);
826  return dist2D/scale;
827  }
828  }
829 
830 } // namespace geo
831 ////////////////////////////////////////////////////////////////////////
bool ltKey(std::pair< double, double > a, double b)
Definition: Geo.cxx:770
T max(const caf::Proxy< T > &a, T b)
bool LineIntersection(double x0, double y0, double x1, double y1, double X0, double Y0, double X1, double Y1, double &x, double &y)
Find the intersection between two line segments.
Definition: Geo.cxx:184
static constexpr Double_t pC
Definition: Munits.h:237
bool IntersectsBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
Determine if a particle starting at xyz with direction dxyz will intersect a box defined by xlo...
Definition: Geo.cxx:97
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: Geo.cxx:373
Float_t y1[n_points_granero]
Definition: compare.C:5
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
double DistToEdge(double *point, double detHalfWidth, double detHalfHeight, double detLength)
Find the distance from the given point to the edge of the detector.
Definition: Geo.cxx:626
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double DetLength() const
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
Definition: Geo.cxx:640
double chi2()
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
float abs(float number)
Definition: d0nt_math.hpp:39
Float_t Y
Definition: plot.C:38
double dist
Definition: runWimpSim.h:113
Float_t Z
Definition: plot.C:38
void InitializePathTable()
Function to initialize gPathTable for GetAvgPath2D.
Definition: Geo.cxx:690
const double C
double sd(Eigen::VectorXd x)
void ProjectToBoxEdgeFromOutside(const double xyz[], const double dxyz[], int axis, double edge, double xyzout[])
Project from a position outside of a box to an edge of the box with coordinate value edge for the axi...
Definition: Geo.cxx:148
const XML_Char * s
Definition: expat.h:262
Double_t scale
Definition: plot.C:25
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
static std::map< double, double > gPathTable
Definition: Geo.cxx:683
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)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
Definition: Geo.cxx:449
static constexpr double L
Double_t X1
Definition: plot.C:264
Double_t Y1
Definition: plot.C:264
const double a
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
Definition: Geo.cxx:338
Collect Geo headers and supply basic geometry functions.
Float_t d
Definition: plot.C:236
double DetHalfHeight() const
double AverageCellPathLength(geo::View_t view, double dx, double dy, double dz)
Mean path length of a ray with (unit) direction vector dx, dy, dz through a cell in view...
Definition: Geo.cxx:804
TVector3 Unit() const
T sqr(T x)
Function to perform, and table to cache, timing fits to ADC values.
Definition: ADCShapeFit.cxx:23
#define PC
Definition: eph_manager.h:45
double DetHalfWidth() const
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void geom(int which=0)
Definition: geom.C:163
const hit & b
Definition: hits.cxx:21
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: MathUtil.cxx:36
double GetAvgPath2D(double dx)
Helper for AverageCellPathLength.
Definition: Geo.cxx:777
void ProjectToBoxEdge(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[])
Project along a direction from a particular starting point to the edge of a box.
Definition: Geo.cxx:38
The geometry of one entire detector (near, far, ipnd)
Definition: GeometryBase.h:49
T min(const caf::Proxy< T > &a, T b)
int WhichWallofBox(const double xyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
Definition: Geo.cxx:75
Float_t e
Definition: plot.C:35
Float_t X
Definition: plot.C:38
Helper for AttenCurve.
Definition: Path.h:10
Float_t w
Definition: plot.C:20