lcsim-contrib/src/main/java/org/lcsim/contrib/onoprien/util/swim

diff -u -r1.1.1.1 -r1.2
--- Helix.java	10 Dec 2008 22:03:06 -0000	1.1.1.1
+++ Helix.java	22 Jan 2009 21:01:58 -0000	1.2
@@ -1,7 +1,5 @@

 package org.lcsim.contrib.onoprien.util.swim;

-import java.util.*;
-

 import hep.physics.matrix.SymmetricMatrix;
 import hep.physics.vec.Hep3Vector;

@@ -11,187 +9,237 @@

  * Helical {@link Trajectory} with its axis parallel to Z.
  *
  * @author D. Onoprienko

- * @version $Id: Helix.java,v 1.1.1.1 2008/12/10 22:03:06 jeremy Exp $

+ * @version $Id: Helix.java,v 1.2 2009/01/22 21:01:58 onoprien Exp $

  */
 public class Helix extends AbstractTrajectory {

-  /** 
-   * Emumeration of helix parameters in C-representation.
-   */

+  /**  Emumeration of helix parameters in C-representation. */

   public enum CRep {

-    /** D0 (signed impact parameter in XY plane). */

+
+    /**
+     * D0 (signed impact parameter in XY plane).<br>
+     * D0 = - (X0-Xr)*sin(PHI0) + (Y0-Yr)*cos(PHI0)
+     */

D0,

-    /** Phi0. */

+
+    /**
+     * Azimuthal angle of the trajectory direction at the point of closest approach
+     * in XY plane to the reference point.
+     */

     PHI0,

-    /**  Omega (signed curvature). */
-    OMEGA,
-    /** Z0. */

+
+    /**
+     * Signed curvature. Positive if rotation in XY plane is clock-wise as viewed from positive Z.
+     * If magnetic field is along positive Z direction, positively charged particles
+     * will leave tracks with positive curvature.
+     */
+    C,
+
+    /** Z of the point of closest approach in XY plane to the reference point. */

Z0,

-    /** Tangent of the dip angle. */

+
+    /** Tangent of the dip angle (dZ/sqrt(dX**2+dY**2)). */

     LAMBDA
   }

-  /** 
-   * Emumeration of helix parameters in W-representation.
-   */
-  public enum WRep {

+  /** Emumeration of helix parameters in V-representation. */
+  public enum VRep {
+

     /** X of the helix origin. */
     X,
     /** Y of the helix origin. */
     Y,
     /** Z of the helix origin. */
     Z,

     /** X component of the unit vector in the direction of the helix at its origin. */
     DX,
     /** Y component of the unit vector in the direction of the helix at its origin. */
     DY,
     /** Z component of the unit vector in the direction of the helix at its origin. */
     DZ,

-    /** Amplitude of the transverse component of the unit vector in the direction of 
-     * the helix at its origin, devided by the helix radius. */
-    RO

+
+    /**
+     * Signed curvature. Positive if rotation in XY plane is clock-wise as viewed from positive Z.
+     * If magnetic field is along positive Z direction, positively charged particles
+     * will leave tracks with positive curvature.
+     */
+    C

   }
 
 // -- Constructors :  ----------------------------------------------------------

+
+  /**
+   * Construct a Helix given origin, direction, and signed curvature.
+   * Reference point is not set.
+   */

   public Helix(Hep3Vector origin, Hep3Vector direction, double curvature) {

     if (origin instanceof ConstHep3Vector) {
       _orig = (ConstHep3Vector) origin;
     } else {
       _orig = new ConstHep3Vector(origin);
     }

     double mag = direction.magnitude();
     _dir = new ConstHep3Vector(direction.x()/mag, direction.y()/mag, direction.z()/mag);

-    
-    _ro = - curvature * Math.hypot(_dir.x(), _dir.y());

+    _c = - curvature;
+
+    _ro = Math.hypot(_dir.x(), _dir.y()) * _c;

+
+  /**
+   * Construct a helix from parameter vector in C-representation and reference point.
+   * The origin is set to the point of closest approach in XY plane.
+   */

   public Helix(ParVector<CRep> parameters, Hep3Vector referencePoint) {
     _ref = new ConstHep3Vector(referencePoint);

-    _parC = new ParVector<CRep>(parameters);
-    repCtoW();

+    _d0 = parameters.get(CRep.D0);
+    _phi0 = parameters.get(CRep.PHI0);
+    _c = parameters.get(CRep.C);
+    _z0 = parameters.get(CRep.Z0);
+    _tan = parameters.get(CRep.LAMBDA);
+    repCtoV();

-  public Helix(ParVector<WRep> parameters) {
-    _orig = new ConstHep3Vector(parameters.get(WRep.X), parameters.get(WRep.Y), parameters.get(WRep.Z));
-    _dir = new ConstHep3Vector(parameters.get(WRep.DX), parameters.get(WRep.DY), parameters.get(WRep.DZ));
-    _ro = parameters.get(WRep.RO);

+  /**
+   * Construct a helix from parameter vector in V-representation.
+   * Reference point is not set.
+   */
+  public Helix(ParVector<VRep> parameters) {
+    _orig = new ConstHep3Vector(parameters.get(VRep.X), parameters.get(VRep.Y), parameters.get(VRep.Z));
+    _dir = new ConstHep3Vector(parameters.get(VRep.DX), parameters.get(VRep.DY), parameters.get(VRep.DZ));
+    _c = parameters.get(VRep.C);
+    _ro = Math.hypot(_dir.x(), _dir.y()) * _c;

-  /** Copy constructor. */

+  /**
+   * Copy constructor.
+   */

   public Helix(Helix helix) {

     _orig = helix._orig;
     _dir = helix._dir;

-    _ro = helix._ro;

+    _c = helix._c;
+

     _ref = helix._ref;

+    _d0 = helix._d0;
+    _phi0 = helix._phi0;
+    _z0 = helix._z0;
+    _tan = helix._tan;
+
+    _ro = helix._ro;

   }
 
 // -- Get helix parameters :  --------------------------------------------------
   
   /** 
    * Returns helix parameters in the requested representation.

+   * An attempt to ask for parameters in C-representation before setting reference point
+   * results in <tt>IllegalArgumentException</tt>.

    *
    * @param representation  Enum class that defines helix parameters representation.

-   *                        Currently either {@link CRep} or {@link WRep}.

+   *                        Currently either {@link CRep} or {@link VRep}.

    */
   public <T extends Enum<T>> ParVector<T> getParameters(Class<T> representation) {
     if (CRep.class.isAssignableFrom(representation)) {

-      if (_parC == null) repWtoC();
-      return (ParVector) _parC;
-    } else if (WRep.class.isAssignableFrom(representation)) {
-      return new ParVector<T>(representation, _orig.x(),_orig.y(),_orig.z(),_dir.x(),_dir.y(),_dir.z(),_ro);

+      if (_ref == null) throw new IllegalArgumentException(ERR1);
+      return new ParVector<T>(representation, _d0, _phi0, _c, _z0, _tan);
+    } else if (VRep.class.isAssignableFrom(representation)) {
+      return new ParVector<T>(representation, _orig.x(),_orig.y(),_orig.z(),_dir.x(),_dir.y(),_dir.z(),_c);

     } else {
       throw new IllegalArgumentException();
     }
   }

   /** 
    * Returns covariance matrix for helix parameters in the requested representation.
    *
    * @param representation  Enum class that defines helix parameters representation.

-   *                        Currently either {@link CRep} or {@link WRep}.

+   *                        Currently either {@link CRep} or {@link VRep}.

    */
   public <T extends Enum<T>> ParCovMatrix<T> getCovMatrix(Class<T> representation) {
     throw new UnsupportedOperationException("Covariance matrix operations are not yet implemented");  // FIXME
   }

+
+  /** Returns helix period in Z. */

   public double getZPeriod() {
     return Math.abs(_dir.z()*TWOPI / _ro);
   }

+
+  /** Returns helix period in path length. */

   public double getLengthPeriod() {
     return Math.abs(TWOPI / _ro);
   }

+  /** Returns helix radius. */
+  public double getRadius() {
+    return 1./Math.abs(_c);
+  }
+

/**

-   * Returns reference point used to calculate helix parameters in C-representation.

+   * Returns reference point used to calculate helix parameters in VRep-representation.
+   * Throws <tt>IllegalStateException</tt> if the reference point has not been set.

    */
   public ConstHep3Vector getReferencePoint() {

-    return (_ref == null) ? _orig : _ref ;

+    if (_ref == null) throw new IllegalStateException(ERR1) ;
+    return _ref;

+
+
+// -- Setting ref point and calculating parameters in C-representation :  ------

/**

-   * Sets the reference point for C-representation of the helix parameters.
-   * Once set by this method, the reference point will not be affected by calls to {@link #swim}.
-   * Calling this method with <tt>null</tt> vector resets the reference point to
-   * the origin of this trajectory, and allows subsequent calls to {@link #swim} to
-   * modify it.

+   * Sets the reference point for VRep-representation of the helix parameters.

    */
   public void setReferencePoint(Hep3Vector referencePoint) {

-    if (referencePoint == null) {
-      _ref = null;
-    } else {
-      _ref = new ConstHep3Vector(referencePoint);
-    }
-    _parC = null;

+    _ref = (referencePoint instanceof ConstHep3Vector) ? (ConstHep3Vector)referencePoint : new ConstHep3Vector(referencePoint);
+    Intersection inter = intersect(new PocaXY(referencePoint));
+    if (!inter.hasNext()) inter.stepBack();
+    Hep3Vector pos = inter.getPosition();
+    Hep3Vector d = inter.getDirection();
+    _phi0 = Math.atan2(d.y(), d.x());
+    _d0 = (pos.y()-_ref.y())*Math.cos(_phi0) - (pos.x()-_ref.x())*Math.sin(_phi0);
+    _z0 = pos.z();
+    _tan = _dir.z()/Math.hypot(_dir.x(), _dir.y());

   }
   
 // -- Position at the given point :  -------------------------------------------

-  /**
-   * Returns the position of the origin of the trajectory.
-   */
-  public ConstHep3Vector getPoint() {

+  /** Returns the position of the origin of the trajectory. */
+  public ConstHep3Vector getPosition() {

     return _orig;
   }

-  /**
-   * Returns the position at given path length along the trajectory.
-   */
-  public ConstHep3Vector getPoint(double pathLength) {

+  /** Returns the position at given path length along the trajectory. */
+  public ConstHep3Vector getPosition(double pathLength) {
+    double angle = Math.abs(_ro)*pathLength;

     double nx = _dir.x()/_ro;
     double ny = _dir.y()/_ro;

-    double x, y, z;

+    double x, y;

     if (Double.isInfinite(nx) || Double.isInfinite(ny)) {
       x = _orig.x() + _dir.x() * pathLength;
       y = _orig.y() + _dir.y() * pathLength;
     } else {

-      double sin = Math.sin(_ro*pathLength);
-      double cos1 = 1. - Math.cos(_ro*pathLength);
-      x = _orig.x() + nx * sin - ny * cos1;
-      y = _orig.y() + ny * sin + nx * cos1;

+      double sin = Math.sin(angle);
+      double cos1 = 1. - Math.cos(angle);
+      x = _orig.x() + ny * cos1 - nx * sin;
+      y = _orig.y() - nx * cos1 - ny * sin;

-    z = _orig.z() + _dir.z() * pathLength;

+    double z = _orig.z() + _dir.z() * pathLength;

     return new ConstHep3Vector(x,y,z);
   }
   
 // -- Direction at the given point :  ------------------------------------------

-  /**
-   * Returns the direction unit vector at the origin.
-   */

+  /** Returns the direction unit vector at the origin. */

   public ConstHep3Vector getDirection() {
     return _dir;
   }

-  /**
-   * Returns the direction unit vector at given path length along the trajectory.
-   */

+  /** Returns the direction unit vector at given path length along the trajectory. */

   public ConstHep3Vector getDirection(double pathLength) {

-    double sin = Math.sin(_ro*pathLength);
-    double cos = Math.cos(_ro*pathLength);

+    double angle = Math.abs(_ro)*pathLength;
+    double sin = Math.sin(angle);
+    double cos = Math.cos(angle);

     double dx = _dir.x() * cos - _dir.y() * sin;
     double dy = _dir.y() * cos + _dir.x() * sin;
     return new ConstHep3Vector(dx, dy, _dir.z());

@@ -214,36 +262,51 @@

   }
   
   public Intersection intersect(Poca poca) {

-    Hep3Vector point = poca.getPoint();
-    double s = 0;   // FIXME
-    return new IntersectionSingle(s, this, poca);

+//    Hep3Vector point = poca.getPoint();
+//    double s = 0;   // FIXME
+//    return new IntersectionSingle(s, this, poca);
+    throw new UnsupportedOperationException();

+
+  /**
+   * Returns intersection with <tt>PocaXY</tt>.
+   * Z coordinate of the specified reference point is ignored - might want to change this behavior later.
+   */

   public Intersection intersect(PocaXY poca) {

+    double roa = Math.abs(_ro);

     Hep3Vector point = poca.getPoint();

-    double s = 0;   // FIXME

+    double deltaX = _orig.x() - point.x();
+    double deltaY = _orig.y() - point.y();
+    double a1 = 1. - _dir.z()*_dir.z() - roa*(_dir.y()*deltaX-_dir.x()*deltaY);
+    double a2 = _dir.x()*deltaX + _dir.y()*deltaY;
+    double a12 = Math.sqrt(a1*a1 + roa*roa*a2*a2);
+    double cos = a1/a12;
+    double sin = - (roa*a2)/a12;
+    double s = getLength(sin,cos);

     return new IntersectionSingle(s, this, poca);
   }

+
+  /** Returns intersection with <tt>ZPlane</tt>. */

   public Intersection intersect(ZPlane zPlane) {

-    double s = getLengthToZ(zPlane._z);

+    double s = findLengthToZ(zPlane._z);

     return (Double.isInfinite(s)) ? new IntersectionNone(this, zPlane) : new IntersectionSingle(s, this, zPlane) ;
   }

+
+  /** Returns intersection with <tt>ZDisk</tt>. */

   public Intersection intersect(ZDisk zDisk) {

-    double s = getLengthToZ(zDisk._z);

+    double s = findLengthToZ(zDisk._z);

     if (Double.isInfinite(s)) return new IntersectionNone(this, zDisk);

-    Hep3Vector pos = getPoint(s);

+    Hep3Vector pos = getPosition(s);

     double x = pos.x();
     double y = pos.y();
     double r2 = x*x + y*y;
     return (r2 > zDisk._rMin2 && r2 < zDisk._rMax2) ? new IntersectionSingle(s, this, zDisk) : new IntersectionNone(this, zDisk);
   }

-  public Intersection intersectZCylinder(ZCylinder zCylinder) { // FIXME

+  public Intersection intersect(ZCylinder zCylinder) { // FIXME

-    double ro2 = _ro*_ro;
-    double a1 = 1. - _dir.z()*_dir.z() - _ro*(_orig.x()*_dir.y() + _orig.y()*_dir.x());

+    double ro2 = _c*_c;
+    double a1 = 1. - _dir.z()*_dir.z() - _c*(_orig.x()*_dir.y() + _orig.y()*_dir.x());

     double a2 = _orig.x()*_dir.x() - _orig.y()*_dir.y();
     double deltaR = _orig.x()*_orig.x() + _orig.y()*_orig.y() - zCylinder._r*zCylinder._r;
     double a12 = a1*a1 + ro2*a2*a2;

@@ -254,8 +317,8 @@

     double cos1 = 1. - (ro2* (2.*a2*a2 - deltaR*a1 + a2*d) ) / (2.*a12);
     double cos2 = 1. - (ro2* (2.*a2*a2 - deltaR*a1 - a2*d) ) / (2.*a12);

-    double sin1 = (- _ro* (2.*a1*a2 + ro2*deltaR*a2 + a1*d)) / (2.*a12);
-    double sin2 = (- _ro* (2.*a1*a2 + ro2*deltaR*a2 - a1*d)) / (2.*a12);

+    double sin1 = (- _c* (2.*a1*a2 + ro2*deltaR*a2 + a1*d)) / (2.*a12);
+    double sin2 = (- _c* (2.*a1*a2 + ro2*deltaR*a2 - a1*d)) / (2.*a12);

     
     double s1 = getLength(sin1, cos1);
     double s2 = getLength(sin2, cos2);

@@ -274,32 +337,44 @@

         s = new double[]{s1,s2};
       }
     }

-    s1 = getLengthToZ(zCylinder.getMinZ());
-    s2 = getLengthToZ(zCylinder.getMaxZ());

+    s1 = findLengthToZ(zCylinder.getMinZ());
+    s2 = findLengthToZ(zCylinder.getMaxZ());

     return new IntersectionPeriodicBound(false, s1, s2, s, getLengthPeriod(), this, zCylinder);
   }

+
+
+// -- Helix-specific intersection methods :  -----------------------------------
+
+  /** Returns path length to Z plane specified by the argument. */
+  double findLengthToZ(double z) {
+    return (z - _orig.z()) / _dir.z();
+  }
+
+

 // -- Modifying the trajectory :  ----------------------------------------------

-  
-  /**

+
+  /**

    * Reverse the direction of the trajectory.
    */
   public void reverse() {

-    throw new UnsupportedOperationException();  // FIXME

+    _dir = new ConstHep3Vector(-_dir.x(), -_dir.y(), -_dir.z());
+    _c = - _c;
+    _ro = -_ro;
+    if (_ref != null) setReferencePoint(_ref);

   /**
    * Moves the origin by the specified path length along the trajectory.

-   * Returns <tt>null</tt> if an attemt to move the origin by the requested length would place
-   * it beyond the boundary of this trajectory.
-   * Unless the reference point that defines C-representation of helix parameters has
-   * been set explicitly in a constructor or through a call to {@link #setReferencePoint(Hep3Vector)}
-   * method, the new origin also becomes the new reference point.

+   * Returns <tt>true</tt>.

    */
   public boolean swim(double pathLength) {

-    throw new UnsupportedOperationException();  // FIXME

+    ConstHep3Vector orig = getPosition(pathLength);
+    _dir = getDirection(pathLength);
+    _orig = orig;
+    return true;

+
+

 // -- Helper methods :  --------------------------------------------------------
   
   double getLength(double sin, double cos) {

@@ -307,39 +382,27 @@

     if (sin < 0.) {
       ros = - ros;
     }

-    return ros/_ro;

+    return ros/Math.abs(_ro);

-  /** Returns path length to Z plane specified by the argument. */
-  double getLengthToZ(double z) {
-    return (z - _orig.z()) / _dir.z();
-  }
-  
-  private void repWtoC() {
-    throw new UnsupportedOperationException();  // FIXME
-  }
-  
-  private void repCtoW() {

+  private void repCtoV() {

     throw new UnsupportedOperationException();  // FIXME
   }
 
 // -- Private parts :  ---------------------------------------------------------

-  // W-representation :
-  
-  ConstHep3Vector _orig;
-  ConstHep3Vector _dir;
-  double _ro;
-  SymmetricMatrix _covW;
-  
-  // C-representation :
-  
-  ConstHep3Vector _ref;    // Ref.point 
-  ParVector<CRep> _parC;   // D0, Phi0, Omega, Z0, Tan Lambda at XY POCA to ref.point
-  SymmetricMatrix _covC;
-  
-  static private final double TWOPI = 2. * Math.PI;

+  ConstHep3Vector _orig; // trajectory origin
+  ConstHep3Vector _dir; // direction at origin
+  double _c; // signed curvature
+  
+  ConstHep3Vector _ref;    // Ref.point
+  double _d0, _phi0, _z0, _tan; // parameters in CRep representation
+ 
+  SymmetricMatrix _covW; // covariance matrix in VRep representation
+  SymmetricMatrix _covC; // covariance matrix in CRep representation

+  double _ro; // sqrt(dx**2+dy**2) * _c

+  static private final double TWOPI = 2. * Math.PI;
+  static private final String ERR1 = "Reference point not set";

* after the last point. The cursor can be moved by calling <tt>step(int)</tt>, * <tt>stepForward()</tt>, and <tt>stepBack()</tt> methods. <tt>hasPrevious()</tt> and * <tt>hasNext()</tt> methods check whether thare is an intersection point upstream and

* <tt>getDirection()</tt> methods return trajectory parameters at the point downstream * from the cursor. An attempt to move the cursor beyond the ends of the point list, or to * obtain trajectory parameters when the cursor is downstream from the last intersection point,

*/ public <T extends Enum<T>> ParCovMatrix<T> getCovMatrix(Class<T> representation) { throw new UnsupportedOperationException("Covariance matrix operations are not yet implemented"); // FIXME

return new ConstHep3Vector(_orig.x() + pathLength * _dir.x(), _orig.y() + pathLength * _dir.y(), _orig.z() + pathLength * _dir.z());

double x = pos.x(); double y = pos.y(); double r2 = x*x + y*y; return (r2 > zDisk._rMin2 && r2 < zDisk._rMax2) ? new IntersectionSingle(s, this, zDisk) : new IntersectionNone(this, zDisk); }

Commit in `lcsim-contrib/src/main/java/org/lcsim/contrib/onoprien/util/swim` on MAIN
`AbstractIntersection.java`	+3	-3	1.1.1.1 -> 1.2
`AbstractSurface.java`	+2	-2	1.1.1.1 -> 1.2
`AbstractTrajectory.java`	+4	-4	1.1.1.1 -> 1.2
`BField.java`	+3	-6	1.1.1.1 -> 1.2
`Helix.java`	+195	-132	1.1.1.1 -> 1.2
`Intersection.java`	+3	-3	1.1.1.1 -> 1.2
`IntersectionHelixZ.java`	+4	-4	1.1.1.1 -> 1.2
`IntersectionMultiple.java`	+2	-3	1.1.1.1 -> 1.2
`IntersectionNone.java`	+4	-4	1.1.1.1 -> 1.2
`IntersectionPeriodic.java`	+2	-2	1.1.1.1 -> 1.2
`IntersectionPeriodicBound.java`	+2	-2	1.1.1.1 -> 1.2
`IntersectionSingle.java`	+3	-3	1.1.1.1 -> 1.2
`Line.java`	+11	-12	1.1.1.1 -> 1.2
`ParCovMatrix.java`	+23	-4	1.1.1.1 -> 1.2
`ParVector.java`	+9	-2	1.1.1.1 -> 1.2
`Poca.java`	+6	-3	1.1.1.1 -> 1.2
`PocaXY.java`	+5	-1	1.1.1.1 -> 1.2
`Surface.java`	+1	-3	1.1.1.1 -> 1.2
`Trajectory.java`	+3	-3	1.1.1.1 -> 1.2
`ZCylinder.java`	+3	-5	1.1.1.1 -> 1.2
`ZDisk.java`	+3	-7	1.1.1.1 -> 1.2
`ZPlane.java`	+3	-5	1.1.1.1 -> 1.2
`package-info.java`	+4		1.1.1.1 -> 1.2
	+298	-213