16 #include "Utilities/func/MathUtil.h" 17 #include "cetlib_except/exception.h" 42 double zlo,
double zhi,
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";
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]; }
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;
70 for (
int i=0;
i<3; ++
i) {
71 xyzout[
i] = xyz[
i] + dxyz[
i]*
d;
78 double zlo,
double zhi)
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;
101 double zlo,
double zhi)
103 double xyzedge[3] = {0.};
110 if(xyzedge[0] > xlo && xyzedge[0] < xhi
111 && xyzedge[2] > zlo && xyzedge[2] < zhi)
return true;
115 if(xyzedge[0] > xlo && xyzedge[0] < xhi
116 && xyzedge[2] > zlo && xyzedge[2] < zhi)
return true;
120 if(xyzedge[1] > ylo && xyzedge[1] < yhi
121 && xyzedge[2] > zlo && xyzedge[2] < zhi)
return true;
125 if(xyzedge[1] > ylo && xyzedge[1] < yhi
126 && xyzedge[2] > zlo && xyzedge[2] < zhi)
return true;
130 if(xyzedge[1] > ylo && xyzedge[1] < yhi
131 && xyzedge[0] > xlo && xyzedge[0] < xhi)
return true;
135 if(xyzedge[1] > ylo && xyzedge[1] < yhi
136 && xyzedge[0] > xlo && xyzedge[0] < xhi)
return true;
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]);
163 double ddS = momDir/momTot;
164 double length1Dim = (edge - posDir);
166 if(TMath::Abs(ddS) > 0.){
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;
185 double X0,
double Y0,
double X1,
double Y1,
186 double&
x,
double&
y)
188 const double dx = x1-x0;
189 const double dy = y1-y0;
190 const double dX = X1-X0;
191 const double dY = Y1-Y0;
193 const double denom = dX*dy-dY*
dx;
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;
208 return l >= 0 && l <= 1 && L >= 0 && L <= 1;
223 const double intercept[],
224 const double slopes[],
227 TVector3 closest_vec;
229 TVector3(intercept), TVector3(slopes),
231 closest_vec.GetXYZ(closest);
247 TVector3 slopes, TVector3& closest)
249 const double s = slopes.Dot(point-intercept);
250 const double sd = slopes.Mag2();
253 closest = intercept + (s/
sd)*slopes;
263 return (point-closest).Mag();
293 TVector3 w0(P0); w0 -= Q0;
295 TVector3
u(P1); u -= P0; u = u.Unit();
296 TVector3
v(Q1); v -= Q0; v = v.Unit();
299 double d = u.Dot(w0);
300 double e = v.Dot(w0);
302 double denom = 1-b*
b;
315 TVector3
pC(P0+sC*u);
316 TVector3 qC(Q0+tC*v);
339 double x1,
double y1,
340 double x2,
double y2)
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));
374 const std::vector<double>&
y,
375 const std::vector<double>&
w,
376 double*
x1,
double*
y1,
377 double*
x2,
double*
y2)
382 if( x.size() != y.size() ||
383 x.size() != w.size() ||
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();
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];
402 if (yhi-ylo==0.0 && xhi-xlo==0.0) {
404 <<
" All points identical in line fit!" <<
std::endl;
416 if (yhi-ylo > xhi-xlo) {
450 const std::vector<double>&
y,
451 const std::vector<double>&
w,
452 double*
x1,
double*
y1,
453 double*
x2,
double*
y2)
456 if( x.size() != y.size() ||
457 x.size() != w.size() ||
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();
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];
509 if(
fabs(yhi-ylo) < 0.1){
519 if (
fabs(xlo - xhi) < 1.e-3) {
521 double xx1, yy1, xx2, yy2;
538 for (i=0; i<w.size(); ++
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];
548 throw cet::exception(
"Geo") <<
"sum of weights is 0 in LinFitMinDperp";
552 const double SwSq =
sqr(Sw);
553 const double SwxSq =
sqr(Swx);
554 const double SwySq =
sqr(Swy);
570 if (C>0.0) C =
sqrt(C);
572 double D = -Swx*Swy+Swxy*Sw;
586 M1 = -(-SwxSq + Swx2*Sw + SwySq - Swy2*Sw -
C)/(2.0*D);
587 B1 = -(-Swy+M1*Swx)/Sw;
601 for (i=0; i<w.size(); ++
i) {
607 double chi2lin, linx1, linx2, liny1, liny2;
609 chi2lin =
LinFit(x,y,w,&linx1,&liny1,&linx2,&liny2);
628 double detHalfHeight,
632 const double dx =
fabs(detHalfWidth -
fabs(point[0]));
633 const double dy =
fabs(detHalfHeight -
fabs(point[1]));
648 TVector3*
ps[2] = {
p0, p1};
650 for(
int n = 0;
n < 2; ++
n){
652 const TVector3 dp = (*ps[1-
n]-*ps[
n]).
Unit();
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;
670 if(
fabs(ps[
n]->
X()) > Lx+.1 ||
671 fabs(ps[
n]->
Y()) > Ly+.1 ||
673 ps[
n]->
Z() > 2*Lz+.1){
693 if(!gPathTable.empty())
694 throw cet::exception(
"Geo") <<
"trying to non-empty initialize path table.";
698 for(
double dx = -.01;
dx < 1.02;
dx += .02){
704 const double Lx = 5.94;
705 const double Ly = 3.84;
712 const double step = .005;
718 for(
double u = 0;
u < U;
u +=
step){
723 double* tInt = tInts;
726 const double x0 = -dy*
u;
727 const double y0 = +dx*
u;
736 if(y > -Ly/2 && y < +Ly/2) *tInt++ =
t;
741 if(x > -Lx/2 && x < +Lx/2) *tInt++ =
t;
746 if(y > -Ly/2 && y < +Ly/2) *tInt++ =
t;
751 if(x > -Lx/2 && x < +Lx/2) *tInt++ =
t;
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";
758 if(tInt == tInts)
continue;
761 S_tot +=
fabs(tInts[1]-tInts[0]);
765 gPathTable[
dx] = S_tot/samples;
770 bool ltKey(std::pair<double, double>
a,
double b){
return a.first <
b;}
783 typedef std::map<double, double>::const_iterator it_t;
786 it_t itNext = std::lower_bound(gPathTable.begin(), gPathTable.end(),
dx,
789 if(itNext == gPathTable.end() || itNext == gPathTable.begin())
793 it_t itPrev = itNext; --itPrev;
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);
805 double dx,
double dy,
double dz)
809 throw cet::exception(
"Geo") <<
"provided values are not a unit vector - " 810 <<
"(dx, dy, dz) = (" << dx <<
" ," << dy <<
" ," bool ltKey(std::pair< double, double > a, double b)
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.
static constexpr Double_t pC
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...
fvar< T > fabs(const fvar< T > &x)
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...
Float_t y1[n_points_granero]
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.
Float_t x1[n_points_granero]
Vertical planes which measure X.
::xsd::cxx::tree::exception< char > exception
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
T sqr(T x)
More efficient square function than pow(x,2)
void InitializePathTable()
Function to initialize gPathTable for GetAvgPath2D.
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...
static std::map< double, double > gPathTable
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...
static constexpr double L
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...
Collect Geo headers and supply basic geometry functions.
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...
T sqr(T x)
Function to perform, and table to cache, timing fits to ADC values.
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.
double pythag(double x, double y)
2D Euclidean distance
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...
double GetAvgPath2D(double dx)
Helper for AverageCellPathLength.
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.
The geometry of one entire detector (near, far, ipnd)
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)