Print

Print


Commit in lcsim/src/org/lcsim/fit/helicaltrack on MAIN
HelixUtils.java+1981.7 -> 1.8
Add some methods to org.lcsim.fit.helicaltrack.HelixUtils to know if and where an Helix is intersecting a surface in 3D Space.
The current methods are not changed

lcsim/src/org/lcsim/fit/helicaltrack
HelixUtils.java 1.7 -> 1.8
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
CVSspam 0.2.8