lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.7 -r1.8
--- HelixUtils.java 13 Oct 2008 01:05:58 -0000 1.7
+++ HelixUtils.java 30 Jun 2009 00:41:12 -0000 1.8
@@ -31,6 +31,8 @@
*/
public class HelixUtils {
private static double _minslope = 1.e-6;
+ private static double _maxpath = 2400; // todo, add a more "dynamic" maxpath
+ private static double _epsilon = 0.01; // value below witch a z vector Or an angle is considered as null
/**
* Creates a new instance of HelixUtils.
@@ -218,7 +220,202 @@
// Return the derivative matrix
return deriv;
}
+ /**
+ * return true if the helix is intercepting the bounded cylinder
+ *
+ * Check if at least one of the first ten intersection of the helix with the cylinder
+ * is between zmin and zmax...there might be a better way to do this
+ * @param helix
+ * @param r rarius of the cylinder
+ * @param zmin lower bound of the cylinder
+ * @param zmax higher bound of the cylinder
+ * @return
+ */
+ public static boolean isInterceptingBoundedCylinder(HelicalTrackFit helix,double r,double zmin,double zmax){
+ double minpath = PathToZPlane(helix, zmin);//not sure it's very efficient to calculate the maximum path
+ double maxpath = PathToZPlane(helix, zmax);
+ double path = Math.max(minpath, maxpath);
+ List<Double> pathlist = PathToCylinder(helix,r,path,10);
+ for(double s:pathlist){
+ double z = PointOnHelix(helix,s).z();
+ if(z<zmax && z> zmin)
+ return true;
+ }
+ return false;
+ }
+ /**
+ * return true if the helix is intercepting the given disk
+ * @param helix
+ * @param rmin
+ * @param rmax
+ * @param z
+ * @return
+ */
+ public static boolean isInterceptingZDisk(HelicalTrackFit helix,double rmin,double rmax, double z ){
+ double s = PathToZPlane(helix,z);
+ Hep3Vector point = PointOnHelix(helix,s);
+ double x = point.x();
+ double y = point.y();
+ double r = Math.sqrt(x*x+y*y);
+ if(r < rmax && r > rmin ){return true;}
+ return false;
+ }
+ /**
+ * return true if the point on the helix at coord z is intercepting a given ZPolygon, the
+ * method don't check if the polygone is parallel to a z plane
+ *
+ * this method check if the private methode insidePolygone return a value a less than .5 deg from 2PI
+ *
+ * @param helix
+ * @param z
+ * @param vertices
+ * @return true OR false wether the helix intercept or not the polygone
+ */
+ public static boolean isInterceptingZpolygone(HelicalTrackFit helix,double z,List<Hep3Vector> vertices){
+ //todo: check if all vertices are on same z
+ double epsilon = Math.PI/360.;
+ double s = PathToZPlane(helix,z);
+ if(s<0){return false;}
+ if(s>_maxpath){return false;}
+ Hep3Vector point = PointOnHelix(helix,s);
+ double angle = insidePolygone(point,vertices);
+ if(Math.abs(angle-2*Math.PI)<epsilon)
+ return true;
+ return false;
+ }
+
+
+
+ /**
+ * Check if the given helix is intercepting XY plan
+ * note, the z coordinate of the XY plane should be 0, but due to numerical
+ * approximation, it accept all z value < _epsilon
+ * @param helix the Helix
+ * @param normal a unitary normal vector to the plan
+ * @param orig one point of the plane
+ * @return true if the helix intersect at least once with the plane
+ */
+ public static boolean isInterceptingXYPlane(HelicalTrackFit helix, Hep3Vector normal,Hep3Vector orig){
+ if(normal.z()>_epsilon)
+ throw new UnsupportedOperationException("Not a XY plan !"+normal);
+ double x = -helix.x0()+orig.x();
+ double y = -helix.y0()+orig.y();
+ double z = -helix.z0()+orig.z();
+ double xn= normal.x();
+ double yn= normal.y();
+ double zn= normal.z();
+ double Phip = Math.atan2(yn, xn);
+ double dist = (xn*x+yn*y+zn*z);
+ double verif = Math.sin(Phip-helix.phi0())+helix.curvature()*dist;
+ if(normal.magnitude() < 0.99)
+ throw new UnsupportedOperationException("normal vector not unitary :"+normal.magnitude());
+ if(Math.abs(verif)>1)
+ return false;
+ else
+ return true ;
+ }
+
+ /**
+ * Check if the given helix is intercepting an XY plane bouded by a
+ * list of nodes
+ * @param helix
+ * @param normal normal vector to the plane
+ * @param nodes
+ * @return true if th e helix intercept the bounded plane
+ **/
+ public static boolean isInterceptingBoundedXYPlane(HelicalTrackFit helix,Hep3Vector normal,List<Hep3Vector> nodes){
+ Hep3Vector orig = nodes.get(0);
+ if(!isInterceptingXYPlane(helix,normal,orig)){return false;}
+ double s = PathToXYPlan(helix, normal, orig);
+ if(s<0) {return false;}
+ if(s>_maxpath){return false;}
+ Hep3Vector point = PointOnHelix(helix, s);
+ if(Math.abs(insidePolygone(point, nodes)-2*Math.PI)<_epsilon)
+ return true;
+ return false;
+ }
+ /**
+ * return one path length (projection on the XY plane) to an XY plane
+ * the methode is for now only used to check IF an helix is intercepting an
+ * XY plane
+ * note, the z coordinate of the XY plane should be 0, but due to numerical
+ * approximation, it accept all z value < _epsilon
+ * @param helix
+ * @param normal a UNITARY vector NORMAL to the plan
+ * @param origin one point of the plane
+ * @return path length
+ */
+ public static double PathToXYPlan(HelicalTrackFit helix,Hep3Vector normal,Hep3Vector origin){
+ if(normal.z()>_epsilon)
+ throw new UnsupportedOperationException("Not a XY plan ! normal vector is "+normal);
+ double x = -helix.x0()+origin.x();
+ double y = -helix.y0()+origin.y();
+ double z = -helix.z0()+origin.z();
+ double xn= normal.x();
+ double yn= normal.y();
+ double zn= normal.z();
+ double phip = Math.atan2(yn, xn);
+ double dist = (xn*x+yn*y+zn*z);
+ double sinPlan = Math.sin(phip-helix.phi0())+helix.curvature()*dist;
+ double Cs = (Math.asin(sinPlan)-phip+helix.phi0());
+ double s=dist/(Math.sqrt(xn*xn+yn*yn)*(2./Cs)*Math.sin(Cs/2.)*Math.cos(Cs/2+phip-helix.phi0())+zn*helix.slope());
+ return s;
+ }
+
+
+ /**
+ * adapted from http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
+ * find if a point is inside a given polygon in 3Dspace
+ * if the point is in the polygone the returned value is 2*PI
+ * if the point is outside of the polygone the returned value is 0
+ * if the point is not of the plan of the polygone the value is between 0 and 2PI
+ * The polygone HAVE to Be convex for the algorythme to work
+ *
+ * This fonction sum all the angle made by q and 2 consecutives point of the polygone
+ *
+ * as this method is generic, it might be usefull to make it public,and mabe in another
+ * file because it can be used for anything else than helix
+ *
+ * @param q the point to determine if it is or not in the polygon
+ * @param p List of edges of the polygone
+ * @return 2PI if inside, 0 if outside,intermediate values if not on the plane
+ */
+ private static double insidePolygone(Hep3Vector q,List<Hep3Vector> p)
+ {
+ int i;
+ double epsilon = _epsilon;
+ double m1,m2;
+ double anglesum=0,costheta;
+ double x1,x2,y1,y2,z1,z2,m1m2;
+ //might be usefull to check if all the p vector are on the same plane
+ Hep3Vector p1;
+ Hep3Vector p2;
+ int n = p.size();
+ for (i=0;i<n;i++) {
+ p1 = p.get(i);
+ p2 = p.get((i+1)%n); // the "last+1" vector is the first
+ x1 = (p1.x()-q.x()); // not optimal, operation can be made only once
+ y1 = (p1.y()-q.y()); //
+ z1 = (p1.z()-q.z());
+ x2 = (p2.x()-q.x());
+ y2 = (p2.y()-q.y());
+ z2 = (p2.z()-q.z());
+ m1m2 = Math.sqrt(x1*x1+y1*y1+z1*z1)*Math.sqrt(x2*x2+y2*y2+z2*z2);
+
+
+ if (m1m2 <= epsilon){
+ return 2*Math.PI; /* We are on a node, consider this inside */
+ }
+ else{
+ costheta = (x1*x2 + y1*y2 + z1*z2)/(m1m2);
+ anglesum += Math.acos(costheta);
+ }
+ }
+ return anglesum;
+ }
+
+
private static double RC(CircleFit cfit) {
return 1.0 / cfit.curvature();
}
@@ -238,4 +435,5 @@
private static double y0(CircleFit cfit) {
return cfit.yref() + cfit.dca() * Math.cos(cfit.phi());
}
+
}
\ No newline at end of file