lcsim/sandbox/RobKutschke
diff -N cyl_int.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ cyl_int.java 30 Aug 2007 14:13:09 -0000 1.1
@@ -0,0 +1,124 @@
+import java.lang.Math;
+
+/**
+ *
+ * Find the interesection of a track, specified by the CLEO paramters
+ * with the a cylinder of radius r, centered on the origin.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: cyl_int.java,v 1.1 2007/08/30 14:13:09 kutschke Exp $
+ *
+ * Date $Date: 2007/08/30 14:13:09 $
+ *
+ */
+
+public class cyl_int{
+
+ private double x;
+ private double y;
+ private double psi;
+
+ private double k;
+ private double d0;
+ private double phi;
+ private double z;
+ private double t;
+
+ private double q;
+ private double rho;
+
+ public cyl_int( double k,
+ double d0,
+ double phi,
+ double z,
+ double t,
+ double r
+ ){
+
+ this.k = k;
+ this.d0 = d0;
+ this.phi = phi;
+ this.z = z;
+ this.t = t;
+
+ // Some other quantities derived from the track parameters.
+
+ q = (k > 0 ) ? 1.0 : -1.0;
+ rho = Math.abs( 0.5/k );
+
+
+ double u0 = Math.cos(phi);
+ double v0 = Math.sin(phi);
+ double x0 = -d0 * v0;
+ double y0 = d0 * u0;
+
+ double x0t = x0-q*rho*v0;
+ double y0t = y0+q*rho*u0;
+
+ double dpsi = Math.atan2(-y0t,x0t);
+
+ double sang = (r*r-x0t*x0t-y0t*y0t-rho*rho)/
+ (2.*q*rho*Math.sqrt(x0t*x0t+y0t*y0t));
+ if ( sang > 1.0 ){
+ sang = 1.0;
+ }else if ( sang < -1.0 ){
+ sang = -1.0;
+ }
+
+ double psi1=(Math.asin(sang)-phi-dpsi);
+ double psi2=(Math.PI-Math.asin(sang)-phi-dpsi);
+
+ double twopi = 2.*Math.PI;
+
+ // Put angles in the correct quadrant.
+ if (psi1 > twopi){
+ psi1 -= twopi;
+ }
+ if (psi1 < 0.0) {
+ psi1 += twopi;
+ }
+ if (psi1 < 0.0) {
+ psi1 += twopi;
+ }
+
+ if (psi2 > twopi) {
+ psi2 -= twopi;
+ }
+ if (psi2 < 0.0) {
+ psi2 += twopi;
+ }
+ if (psi2 < 0.0) {
+ psi2 += twopi;
+ }
+
+ psi = ( psi2 < psi1) ? psi2 : psi1;
+
+ x = x0t + q*rho*Math.sin(phi+q*psi);
+ y = y0t - q*rho*Math.cos(phi+q*psi);
+ }
+
+ double getX(){
+ return x;
+ }
+
+ double getY(){
+ return y;
+ }
+
+ double getZ(){
+ double cz = t / Math.sqrt( 1. + t*t );
+ return z + rho*psi*t;
+ }
+
+
+ double getPsi(){
+ return psi;
+ }
+
+
+}
+
+
+
+
+