lcsim/src/org/lcsim/util/swim
diff -N HelixSwim.java
--- HelixSwim.java 24 May 2006 17:48:47 -0000 1.12
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,481 +0,0 @@
-package org.lcsim.util.swim;
-
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-
-import java.util.ArrayList;
-import java.util.List;
-import org.lcsim.event.Track;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.SpacePoint;
-
-/**
- * Given a particle of initial momentum, position and charge (p, r0 and iq),
- * calculates its trajectory in a uniform magnetic field in the z-direction. The
- * trajectory is parametrized as function of the path length alpha. <br>
- * This implementation works for charged and neutral particles alike.
- *
- * @author W.Walkowiak
- * @version $Id: HelixSwim.java,v 1.12 2006/05/24 17:48:47 jstrube Exp $
- * @deprecated Please use {@link org.lcsim.util.swim.HelixSwimmer#HelixSwimmer(double)} and {@link org.lcsim.util.swim.Helix} instead
- * @see org.lcsim.util.swim.HelixSwimmer
- */
-@Deprecated public class HelixSwim {
- /**
- * Create a helix swimmmer.
- *
- * @param B field strength; uniform, solenoidal, directed along z-axis
- */
- public HelixSwim(double B) {
- isValid = false;
- centerOK = false;
- bFieldZ = B;
- }
-
- /**
- * Create a helix swimmmer set B field and cylindrical geometry parallel to
- * z-axis.
- *
- * @param B field strength; uniform, solenoidal, directed along z-axis
- * @param rhoMin minimum radius
- * @param rhoMax maximum radius
- * @param zMax maximum extend in +/-z
- * @see #setCylinderGeometry
- */
- public HelixSwim(double B, double rhoMin, double rhoMax, double zMax) {
- isValid = false;
- centerOK = false;
- bFieldZ = B;
- setCylinderGeometry(rhoMin, rhoMax, zMax);
- }
-
- /**
- * Create and initialize a helix swimmmer.
- *
- * @param B field strength; uniform, solenoidal, directed along z-axis
- * @param p 3-momentum (px,py,pz)
- * @param r0 initial position(x0,y0,z0)
- * @param iq charge iq = q/|e| = +1/0/-1
- */
- public HelixSwim(double B, double[] p, double[] r0, int iq) {
- isValid = false;
- centerOK = false;
- bFieldZ = B;
- setTrack(p, r0, iq);
- }
-
- /**
- * Sets parameters for helix swimmmer.
- *
- * @param p 3-momentum (px,py,pz)
- * @param r0 initial position(x0,y0,z0)
- * @param iq charge iq = q/|e| = +1/0/-1
- */
- public void setTrack(double[] p, double[] r0, int iq) {
- isValid = false;
- centerOK = false;
-
- momentum = new double[3];
- referencePoint = new double[3];
- for (int i = 0; i < 3; i++) {
- momentum[i] = p[i];
- referencePoint[i] = r0[i];
- }
- m_iq = iq;
-
- m_pt = Math.sqrt(p[0] * p[0] + p[1] * p[1]);
- m_ptot = Math.sqrt(m_pt * m_pt + p[2] * p[2]);
- m_pz = p[2];
- m_phi = Math.atan2(p[1], p[0]);
-
- calcCenter();
- }
-
- /**
- * Sets the track to swim
- *
- * @param t the Track
- */
- public void setTrack(Track t) {
- isValid = false;
- centerOK = false;
-
- referencePoint = t.getReferencePoint();
- momentum = t.getMomentum();
- m_iq = t.getCharge();
- m_pt = Math.sqrt(momentum[0] * momentum[0] + momentum[1] * momentum[1]);
- m_ptot = Math.sqrt(m_pt * m_pt + momentum[2] * momentum[2]);
- m_pz = momentum[2];
- m_phi = Math.atan2(momentum[1], momentum[0]);
-
- calcCenter();
- }
-
- /**
- * Swims the track to a point at a distance alpha
- *
- * @param alpha The length parameter
- * @return The point at a length alpha from the origin
- */
- // FIXME the parameter alpha is not well defined
- public SpacePoint getPointAtLength(double alpha) {
- return new CartesianPoint(swimBy(alpha));
- }
-
- /**
- * swim along this track to point at distance alpha from it's origin (in
- * positive track direction)
- *
- * @deprecated in favor of
- * @link #getPointAtLength
- */
- @Deprecated
- public double[] swimBy(double alpha) {
- if (alpha != m_alpha || !centerOK)
- isValid = false;
-
- if (!centerOK)
- calcCenter();
-
- if (!isValid) {
-
- m_alpha = alpha;
-
- // now swim it
-
- m_x = new double[3];
-
- if (m_iq != 0 && bFieldZ != 0) { // charged tracks in field
-
- double alp = -m_iq * alpha;
-
- m_x[0] = m_xc + m_rc * Math.cos(m_alpha0 + alp);
- m_x[1] = m_yc + m_rc * Math.sin(m_alpha0 + alp);
- m_x[2] = m_z0 - m_iq * m_c * m_pz * alp;
-
- } else { // neutrals or no field
-
- m_x[0] = referencePoint[0] + alpha * momentum[0] / m_ptot;
- m_x[1] = referencePoint[1] + alpha * momentum[1] / m_ptot;
- m_x[2] = referencePoint[2] + alpha * momentum[2] / m_ptot;
- }
- }
- isValid = true;
-
- return m_x;
- }
-
- /**
- * report the circle's radius in the rho-phi plane
- */
- public double getRc() {
- return m_rc;
- }
-
- /**
- * report tan lambda
- *
- * @return tangent of the dip angle lambda
- */
- public double getTanL() {
- double tanL = m_pz / PTINY;
- if (m_pt > PTINY)
- tanL = m_pz / m_pt;
-
- return tanL;
- }
-
- /**
- * report sin lambda
- *
- * @return sine of the dip angle lambda
- */
- public double getSinL() {
- double sinL = m_pz / PTINY;
- if (m_ptot > PTINY)
- sinL = m_pz / m_ptot;
-
- return sinL;
- }
-
- /**
- * report cos lambda
- *
- * @return cosine of the dip angle lambda
- */
- public double getCosL() {
- double cosL = m_pt / PTINY;
- if (m_ptot > PTINY)
- cosL = m_pt / m_ptot;
-
- return cosL;
- }
-
- /**
- * Set minimum point number on a full circle.
- *
- * @param nPoints number of points on a full circle
- * @see #asPoints
- */
- public void setPointDensity(int nPoints) {
- m_nPoints = Math.abs(nPoints);
- }
-
- /**
- * Set point distance along trajectory: <br>
- *
- * @param ptDist distance between two points
- * @see #asPoints
- */
- public void setPointDistance(double ptDist) {
- m_ptDist = Math.abs(ptDist);
- }
-
- /**
- * Swims the helix along its trajectory to the
- * point of closest approach to the given SpacePoint.
- * The equation to solve for s is an O(2) Taylor-expansion.
- * The parameterization is as follows:<br />
- * p<sub>x</sub> = p<sub>0x</sub> cos(&rho s) - p<sub>0y</sub> sin(&rho s)<br />
- * p<sub>y</sub> = p<sub>0y</sub> cos(&rho s) + p<sub>0x</sub> sin(&rho s)<br />
- * p<sub>z</sub> = p<sub>0z</sub><br />
- * x = x<sub>0</sub> + p<sub>0x</sub>/a sin(&rho s) - p<sub>0y</sub>/a (1-cos(&rho s))<br />
- * y = y<sub>0</sub> + p<sub>0y</sub>/a sin(&rho s) + p<sub>0x</sub>/a (1-cos(&rho s))<br />
- * z = z<sub>0</sub> + p<sub>z</sub>/p s<br />
- * a = -0.299792458 B q, where [B] = T and [q] = e;<br />
- * &rh0 = a / p;<br />
- * @param point Point in Space to swim to.
- * @return the length Parameter s
- */
- public double swimToPoca(SpacePoint point) {
-
- return 0;
- }
-
- /**
- * Set cylindrical geometry parallel to z-axis.
- *
- * @param rhoMin minimum radius
- * @param rhoMax maximum radius
- * @param zMax maximum extend in +/-z
- * @see #asPoints
- */
- public void setCylinderGeometry(double rhoMin, double rhoMax, double zMax) {
- if (Math.abs(rhoMax) < rhoMin)
- throw new IllegalArgumentException("rhoMax < rhoMin is illegal!");
-
- isValid = false;
- m_rhoMin = rhoMin;
- m_rhoMax = Math.abs(rhoMax);
- m_zMax = Math.abs(zMax);
- }
-
-// private enum Status {PROBLEM, PLANE_PLUSZ, PLANE_MINUSZ, PLANE_R};
-
- // FIXME This function duplicates some already existing code
- public Hep3Vector swimToRadialIntersect(double r_cyl, double zMax) {
- double s, sin_lambda, cos_lambda, sin_phi0, cos_phi0, szhit, srhit, sin_lam_min;
- double Rc, phic, Rq, darg, diff, tdiff;
-// Status i_hit = Status.PROBLEM;
-
- Hep3Vector r0 = new BasicHep3Vector(referencePoint[0], referencePoint[1], referencePoint[2]);
- // calculate useful quantities: lambda is dip angle (w/rt x-y plane;
- // phi0 is initial azimuthal angle in x-y plane w/rt x-axis;
- // R is radius of circle in x-y plane;
- // Rc**2=xc**2 + yc**2 is distance to center of circle in x-y plane
- sin_lambda = m_pz / m_ptot;
- cos_lambda = m_pt / m_ptot;
- cos_phi0 = momentum[0] / m_pt;
- sin_phi0 = momentum[1] / m_pt;
- //
- // First, check intercept with planes at z=+/- z_cyl
- // Eqn of motion z(s)=z0 + s*sin(lambda)
- // Check +z_cyl plane:
- //
- final double sbig = 99999.;
- final double pt_tiny = 0.010;
-
- szhit = sbig;
- sin_lam_min = zMax / sbig;
- if (sin_lambda > sin_lam_min) {
-// i_hit = Status.PLANE_PLUSZ;
- szhit = (zMax - r0.z()) / sin_lambda;
- }
- // likewise check -z_cyl plane
- else if (sin_lambda < -sin_lam_min) {
-// i_hit = Status.PLANE_MINUSZ;
- szhit = (-zMax - r0.z()) / sin_lambda;
- }
- //
- // Now check intercept of x-y plane path with cylinder
- //
- // First, check that pt is finite
- if (m_pt < pt_tiny) {
- return new BasicHep3Vector(0, 0, r0.z() + szhit * sin_lambda);
- }
- //
- // calculate center of circle and circle radius=R
- //
- Rq = ONE_METER*m_pt / (0.3 * bFieldZ) / -m_iq;
- // FIXME calcCenter already does something like that
- double xc = r0.x() + (Rq) * sin_phi0;
- double yc = r0.y() - (Rq) * cos_phi0;
- Rc = Math.sqrt(xc * xc + yc * yc);
- phic = Math.atan2(yc, xc);
- //
- // now calculate path length at intersection point
- //
- darg = r_cyl * r_cyl / (2. * Rq * Rc) - Rc / (2. * Rq) - Rq / (2. * Rc);
- if (Math.abs(darg) > 1.0) {
- srhit = sbig;
- } else {
- diff = Math.asin(darg) + m_phi - phic;
- tdiff = Math.tan(diff);
- srhit = (Rq / cos_lambda) * Math.atan(tdiff);
- }
- //
- // Determine if cyl or z-plane hit first (smaller path length)
- //
- if (Math.abs(srhit) <= szhit) {
- s = srhit;
-// i_hit = Status.PLANE_R;
- } else {
- s = szhit;
- }
- //
- // Calculate position at this intercept point
- //
- darg = s * cos_lambda / Rq - m_phi;
- return new BasicHep3Vector(xc + Rq * Math.sin(darg)
- , yc + Rq * Math.cos(darg)
- , r0.x() + s * sin_lambda);
- }
-
- /**
- * Provide points along the particle's track within cylindrical volume.
- *
- * @param alphaMax : maximum path length along trajectory, alphaMax < 0 :
- * alphaMax is automatically set according to the cylinders
- * extensions.
- * @return List of double[3] containing points on trajectory.
- * @see #setCylinderGeometry
- * @see #setPointDistance
- * @see #setPointDensity
- */
- public List<double[]> asPoints(double alphaMax) {
- // FIXME Tony, please take a look at this
- double refRc = 100; // 1 m
-
- List<double[]> vp = new ArrayList<double[]>();
-
- if (alphaMax < 0) {
- double tanL = getTanL();
- double sinL = getSinL();
- if (Math.abs(momentum[2]) > HelixSwim.PTINY && tanL != 0. && m_iq != 0) {
- alphaMax = Math.abs((m_zMax + Math.abs(referencePoint[2])) / tanL);
- } else if (sinL != 0. && m_iq == 0) {
- alphaMax = Math.abs((m_zMax + Math.abs(referencePoint[2])) / sinL);
- } else {
- alphaMax = 2. * Math.PI;
- }
- // System.out.println("alphaMax: "+alphaMax+" tanL: "+tanL+
- // " SinL: "+sinL+" zMax: "+m_zMax+" iq: "+m_iq);
- }
- double cosL = getCosL();
- double rc = Math.abs(getRc());
- double dAlpha = m_ptDist;
-
- if (m_iq != 0) {
- double dAlphaMin1 = (m_pt > PTINY) ? dAlpha * Math.abs(cosL) / rc : dAlpha;
- double dAlphaMin2 = (rc < refRc) ? 2. * Math.PI / m_nPoints : refRc / rc * 2. * Math.PI
- / m_nPoints;
- dAlpha = Math.min(dAlphaMin1, dAlphaMin2);
- // System.out.println("rc: "+rc+" dAlphaMin1: "+dAlphaMin1+
- // " dAlphaMin2: "+dAlphaMin2);
- }
- // System.out.println("dAlpha: "+dAlpha);
- boolean hitBounds = false;
- boolean wasVisible = false;
- double alpha = 0;
-
- if (m_rhoMax > m_rhoMin) {
- while (alpha <= alphaMax && !hitBounds) {
- double[] x = swimBy(alpha);
- double rho = Math.sqrt(x[0] * x[0] + x[1] * x[1]);
- if ((m_zMax != 0 && Math.abs(x[2]) > m_zMax) || (m_rhoMax > 0 && rho > m_rhoMax)) {
- hitBounds = wasVisible;
- } else if (rho < m_rhoMin) {
- } else {
- vp.add(new double[] { x[0], x[1], x[2] });
- wasVisible = true;
- }
- alpha += dAlpha;
- }
- // if ( hitBounds )
- // System.out.println("BANG: alpha: "+alpha+" NVert: "+vp.size());
- } else {
- System.out.println("HelixSwim: asPoints: rhoMax < rhoMin !!");
- }
- return vp;
- }
-
- // calculates the center of the helix hypothesis
- private void calcCenter() {
- isValid = false;
- m_alpha0 = m_phi + 0.5 * m_iq * Math.PI;
- if (m_iq != 0 && bFieldZ != 0.) {
- // FIXME hardcoded number needs explanation, please
- m_c = ONE_METER / (0.3 * bFieldZ);
- m_rc = m_c * m_pt;
-
- m_xc = referencePoint[0] - m_rc * Math.cos(m_alpha0);
- m_yc = referencePoint[1] - m_rc * Math.sin(m_alpha0);
- } else {
- m_rc = -1.;
- }
- m_z0 = referencePoint[2];
- centerOK = true;
- }
-
- /**
- * public static constants
- */
- private final static double PTINY = 1.e-6;
- // one meter in org.lcsim standard units
- private final static int ONE_METER = 1000;
-
- // input: environment and track data
- private double bFieldZ;
- private double[] momentum;
- private double[] referencePoint;
- // the charge in multiples of the positron charge
- private int m_iq;
-
- // cylinder geomtry
- private double m_rhoMin = 0;
- private double m_rhoMax = 1.e6;
- private double m_zMax = 1.e6;
-
- // No. of points on circle with radius of 1 m or less.
- // Scaled if radius is bigger.
- private double m_nPoints = 200;
- private double m_ptDist = 5.; // 5 mm
- // swimming
- private double m_alpha;
- private double m_xc, m_yc;
- private double m_rc;
- private double m_z0;
- private double m_pt, m_pz, m_ptot;
- private double m_phi;
- private double m_alpha0;
- // unit conversion constant
- // FIXME this needs a real definition
- private double m_c;
-
- // resulting position at distance alpha
- private double[] m_x;
-
- // status
- private boolean isValid = false;
- private boolean centerOK = false;
-}
