NOvA Offline: geo Namespace Reference

00025                            {
00026     kX,   
00027     kY,   
00028     kXorY 
00029   } View_t;

00018                         {
00019     kTop,    
00020     kWest,   
00021     kEast,   
00022     kAnySide 
00023   } Readout_t;

00026                      {
00027     kPLANE_NOT_FOUND = 999000
00028   };

00073   {
00074     double s = 
00075       (slopes[0]*(point[0]-intercept[0]) + 
00076        slopes[1]*(point[1]-intercept[1]) +
00077        slopes[2]*(point[2]-intercept[2]));
00078     double sd = 
00079       (slopes[0]*slopes[0] + 
00080        slopes[1]*slopes[1] +
00081        slopes[2]*slopes[2]);
00082     if (sd>0.0) {
00083       s /= sd;
00084       closest[0] = intercept[0] + s*slopes[0];
00085       closest[1] = intercept[1] + s*slopes[1];
00086       closest[2] = intercept[2] + s*slopes[2];
00087     }
00088     else {
00089       // How to handle this zero gracefully? Assume that the intercept
00090       // is a particle vertex and "slopes" are momenta. In that case,
00091       // the particle goes nowhere and the closest approach is the
00092       // distance from the intercept to point
00093       closest[0] = intercept[0];
00094       closest[1] = intercept[1];
00095       closest[2] = intercept[2];
00096     }
00097     return sqrt(pow((point[0]-closest[0]),2)+
00098                 pow((point[1]-closest[1]),2)+
00099                 pow((point[2]-closest[2]),2));
00100   }

00119   {
00120     // If we've been handed a point instead of a line, return distance
00121     // to the point
00122     if (x1==x2 && y1==y2) return pow(x0-x1,2)+pow(y0-y1,2);
00123     
00124     // Otherwise do the geometry for the line
00125     assert(x1!=x2 || y1!=y2);
00126     return
00127       pow((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1),2) / 
00128       (pow(x2-x1,2)+pow(y2-y1,2));
00129   }

00066   {
00067     char buff[8];
00068     const char* c;
00069     const char* start = nm;
00070     const char* end;
00071     int indx;
00072     while (1) {
00073       start = strchr(start, '_');
00074       if (start == NULL) break;
00075       ++start;
00076       end = strchr(start, '/');
00077       if (end==NULL) end = strchr(start, '\0');
00078       if (end==NULL) abort();
00079       
00080       for (indx=0, c=start; c<end; ++indx, ++c) buff[indx] = *c;
00081       buff[indx] = '\0';
00082       id.push_back(atoi(buff));
00083     }
00084   }

00041   {
00042     // Experience has shown that only depths below this "depth0"
00043     // contribute to a cell's 'uniqueness'.
00044     static const unsigned int depth0 = 3;
00045     
00046     // Now construct the id
00047     CellUniqueId id = 0ULL;
00048     for (unsigned int i=depth0; i<ids.size(); ++i) {
00049       // The following check is useful, although not really
00050       // required. The code in Geometry.cxx checks that the ID's are
00051       // unique which is all that really matters. Even if we happen to
00052       // overflow the storage limits of the CellUniqueId, the result
00053       // will always be the same for the same cell. So, as long as the
00054       // overflow is unique, all is OK. BTW, I've checked that as of
00055       // at least v1.6 2008/08/17 the CellUniqueID's do not overflow
00056       //
00057       // if (id>0 && ULLONG_MAX/id < 1000) abort();
00058       id = 1000*id + ids[i];
00059     }
00060     return id;
00061   }

00154   {
00155     // Before going ahead, make sure we have sensible arrays
00156     assert(x.size()==y.size());
00157     assert(x.size()==w.size());
00158     assert(x.size()>=2);
00159     //--
00160     // The following code implements the output of the following MAPLE
00161     // script:
00162     //--
00163     // assume(x0,real);
00164     // assume(y0,real);
00165     // assume({M,B},real);
00166     // assume({Dsqr},real);
00167     // assume(Sw,real);
00168     // assume(Swx,real);
00169     // assume(Swx2,real);
00170     // assume(Swy,real);
00171     // assume(Swy2,real);
00172     // assume(Swxy,real);
00173     // y1 := M*x1 + B;
00174     // y2 := M*x2 + B;
00175     // Dsqr := (((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1))^2/((x2-x1)^2+(y2-y1)^2));
00176     // eqn1 := simplify(expand(diff(Dsqr,M) = 0));
00177     // eqn2 := simplify(expand(diff(Dsqr,B) = 0));
00178     // Eqn1 := -2*B*Swy*M +B^2*M*Sw +Swy2*M -Swxy*M^2 -B*Swx +Swxy
00179     // +B*Swx*M^2 -Swx2*M = 0;
00180     // Eqn2 := B*Sw - Swy + M*Swx = 0;
00181     // solve({Eqn1,Eqn2},{M,B});
00182     // allvalues(%);
00183     //--
00184     unsigned register int i;
00185     
00186     // Find the extremes of the inputs
00187     double xlo = x[0];
00188     double xhi = x[0];
00189     double ylo = y[0];
00190     double yhi = y[0];
00191     for (i=1; i<w.size(); ++i) {
00192       if (x[i]<xlo) xlo = x[i];
00193       if (y[i]<ylo) ylo = y[i];
00194       if (x[i]>xhi) xhi = x[i];
00195       if (y[i]>yhi) yhi = y[i];
00196     }
00197     
00198     // For stability, fit in a way that avoids infinite slopes,
00199     // exchanging the x and y variables as needed.
00200     if ((yhi-ylo)>(xhi-xlo)) {
00201       double chi2;
00202       double xx1, yy1, xx2, yy2;
00203       // Notice: (x,y) -> (y,x)
00204       chi2 = LinFitMinDperp(y,x,w,&xx1, &yy1,&xx2, &yy2);
00205       *x1 = yy1;
00206       *y1 = xx1;
00207       *x2 = yy2;
00208       *y2 = xx2;
00209       return chi2;
00210     }
00211     
00212     // Accumulate the sums needed to compute the fit line
00213     double Sw  = 0.0;
00214     double Swx = 0.0;
00215     double Swy = 0.0;
00216     double Swxy= 0.0;
00217     double Swy2= 0.0;
00218     double Swx2= 0.0;
00219     for (i=0; i<w.size(); ++i) {
00220       Sw   += w[i];
00221       Swx  += w[i]*x[i];
00222       Swy  += w[i]*y[i];
00223       Swxy += w[i]*x[i]*y[i];
00224       Swy2 += w[i]*y[i]*y[i];
00225       Swx2 += w[i]*x[i]*x[i];
00226     }
00227     assert(Sw>0.0);
00228     
00229     // C and D are two handy temporary variables
00230     double C =
00231       +pow(Swx,4)
00232       -2.0*pow(Swx,2)*Swx2*Sw
00233       +2.0*pow(Swx,2)*pow(Swy,2)
00234       +2.0*pow(Swx,2)*Swy2*Sw
00235       +pow(Swx2,2)*pow(Sw,2)
00236       +2.0*Swx2*Sw*pow(Swy,2)
00237       -2.0*Swx2*pow(Sw,2)*Swy2
00238       +pow(Swy,4)
00239       -2.0*pow(Swy,2)*Swy2*Sw
00240       +pow(Swy2,2)*pow(Sw,2)
00241       -8.0*Swx*Swy*Swxy*Sw
00242       +4.0*pow(Swxy,2)*pow(Sw,2);
00243     if (C>0.0) C = sqrt(C);
00244     else       C = 0.0;
00245     double D = -Swx*Swy+Swxy*Sw;
00246     
00247     // The first of the two solutions. This one is the best fit through
00248     // the points
00249     double M1;
00250     double B1;
00251     if (D==0.0) {
00252       // D could be zero, this might happen in a segmented detector if
00253       // the track fit passed through only a single readout plane and
00254       // all hits share the same x coordinates. I've tried to avoid
00255       // this by checking if it makes more sense to fit with x and y
00256       // exchanged above. This is an extra precaution.
00257       D = 1.0E-9;
00258     }
00259     M1 = -(-pow(Swx,2) + Swx2*Sw + pow(Swy,2) - Swy2*Sw - C)/(2.0*D);
00260     B1 = -(-Swy+M1*Swx)/Sw;
00261     
00262     // This second solution is perpendicular to the first and represents
00263     // the worst fit
00264     // M2 = -(-pow(Swx,2) + Swx2*Sw + pow(Swy,2) - Swy2*Sw + C)/(2.0*D);
00265     // B2 = -(-Swy+M2*Swx)/Sw;
00266     
00267     *x1 = xlo;
00268     *y1 = M1*xlo + B1;
00269     *x2 = xhi;
00270     *y2 = M1*xhi + B1;
00271     
00272     // Compute the chi^2 function for return
00273     double chi2 = 0.0;
00274     for (i=0; i<w.size(); ++i) {
00275       chi2 += w[i]*geo::DsqrToLine(x[i],y[i],*x1,*y1,*x2, *y2);
00276     }
00277     
00278     return chi2;
00279   }

00019   {
00020     static std::vector<int> ids;
00021     ids.clear();
00022     for (unsigned int i=0; i<=depth; ++i) {
00023       GetNodeNumbers(n[i]->GetName(), ids);
00024     }
00025     return IdsToUniqueId(ids);
00026   }

00031   {
00032     static std::vector<int> ids;
00033     ids.clear();
00034     GetNodeNumbers(path, ids);
00035     return IdsToUniqueId(ids);
00036   }

00029   {
00030     // Make sure we're inside the box!
00031     assert(xyz[0]>=xlo && xyz[0]<=xhi);
00032     assert(xyz[1]>=ylo && xyz[1]<=yhi);
00033     assert(xyz[2]>=zlo && xyz[2]<=zhi);
00034     
00035     // Compute the distances to the x/y/z walls
00036     double dx = 99.E99;
00037     double dy = 99.E99;
00038     double dz = 99.E99;
00039     if      (dxyz[0]>0.0) { dx = (xhi-xyz[0])/dxyz[0]; }
00040     else if (dxyz[0]<0.0) { dx = (xlo-xyz[0])/dxyz[0]; }
00041     if      (dxyz[1]>0.0) { dy = (yhi-xyz[1])/dxyz[1]; }
00042     else if (dxyz[1]<0.0) { dy = (ylo-xyz[1])/dxyz[1]; }
00043     if      (dxyz[2]>0.0) { dz = (zhi-xyz[2])/dxyz[2]; }
00044     else if (dxyz[2]<0.0) { dz = (zlo-xyz[2])/dxyz[2]; }
00045     
00046     // Choose the shortest distance
00047     double d = 0.0;
00048     if      (dx<dy && dx<dz) d = dx;
00049     else if (dy<dz && dy<dx) d = dy;
00050     else if (dz<dx && dz<dy) d = dz;
00051     
00052     // Make the step
00053     for (int i=0; i<3; ++i) {
00054       xyzout[i] = xyz[i] + dxyz[i]*d;
00055     }
00056   }


Classes
class	geo::CellGeo
	Encapsulate the cell geometry. More...
class	geo::Exception
	Exceptions thrown by the geometry package. More...
class	geo::Geometry
	The geometry of one entire detector (near, far, ipnd). More...
class	geo::PlaneGeo
	Geometry information for a single readout plane. More...
Typedefs
typedef unsigned long long int	CellUniqueId
typedef enum geo::_readout	Readout_t
	Enumerate the possible locations of the cell readout.
typedef enum geo::_plane_proj	View_t
	Enumerate the possible plane projections.
Enumerations
enum	_return_codes { kPLANE_NOT_FOUND = 999000 }
	Return codes which flag errors and the like. More...
enum	_readout { kTop, kWest, kEast, kAnySide }
	Enumerate the possible locations of the cell readout. More...
enum	_plane_proj { kX, kY, kXorY }
	Enumerate the possible plane projections. More...
Functions
CellUniqueId	NodesToUniqueId (const std::vector< const TGeoNode * > &n, unsigned int depth)
CellUniqueId	PathToUniqueId (const char *path)
CellUniqueId	IdsToUniqueId (const std::vector< int > &ids)
void	GetNodeNumbers (const char *nm, std::vector< int > &id)
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[])
double	DsqrToLine (double x0, double y0, double x1, double y1, double x2, double y2)
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)

geo Namespace Reference

Classes

Typedefs

Enumerations

Functions

Detailed Description

Typedef Documentation

Enumeration Type Documentation

Function Documentation