lcsim/src/org/lcsim/util/swim
diff -u -r1.5 -r1.6
--- HelixSwim.java 3 Aug 2005 22:53:37 -0000 1.5
+++ HelixSwim.java 15 Aug 2005 05:23:52 -0000 1.6
@@ -13,7 +13,7 @@
* This implementation works for charged and neutral particles alike.
*
* @author W.Walkowiak
- * @version $Id: HelixSwim.java,v 1.5 2005/08/03 22:53:37 jstrube Exp $
+ * @version $Id: HelixSwim.java,v 1.6 2005/08/15 05:23:52 jstrube Exp $
*/
public class HelixSwim
{
@@ -25,7 +25,7 @@
{
isValid = false;
centerOK = false;
- m_B = B;
+ bFieldZ = B;
}
/**
@@ -41,7 +41,7 @@
{
isValid = false;
centerOK = false;
- m_B = B;
+ bFieldZ = B;
setCylinderGeometry(rhoMin,rhoMax,zMax);
}
@@ -50,13 +50,13 @@
* @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
+ * @param iq charge iq = q/|e| = +1/0/-1
*/
public HelixSwim(double B, double[] p, double [] r0, int iq)
{
isValid = false;
centerOK = false;
- m_B = B;
+ bFieldZ = B;
setTrack(p, r0, iq);
}
@@ -65,18 +65,18 @@
*
* @param p 3-momentum (px,py,pz)
* @param r0 initial position(x0,y0,z0)
- * @param iq charge iq = q/|e| = +/-1
+ * @param iq charge iq = q/|e| = +1/0/-1
*/
public void setTrack(double[] p, double[] r0, int iq)
{
isValid = false;
centerOK = false;
- m_p = new double[3];
- m_r0 = new double[3];
+ momentum = new double[3];
+ referencePoint = new double[3];
for (int i=0; i<3; i++ ) {
- m_p[i] = p[i];
- m_r0[i] = r0[i];
+ momentum[i] = p[i];
+ referencePoint[i] = r0[i];
}
m_iq = iq;
@@ -97,13 +97,13 @@
isValid = false;
centerOK = false;
- m_r0 = t.getReferencePoint();
- m_p = t.getMomentum();
+ referencePoint = t.getReferencePoint();
+ momentum = t.getMomentum();
m_iq = t.getCharge();
- m_pt = Math.sqrt(m_p[0]*m_p[0]+m_p[1]*m_p[1]);
- m_ptot = Math.sqrt(m_pt*m_pt+m_p[2]*m_p[2]);
- m_pz = m_p[2];
- m_phi = Math.atan2(m_p[1],m_p[0]);
+ 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();
}
@@ -122,7 +122,7 @@
/**
* swim along this track to point at distance alpha from it's origin
* (in positive track direction)
- * @deprecated in favor of @see getPointAtLength
+ * @deprecated in favor of @link #getPointAtLength
*/
@Deprecated
public double[] swimBy(double alpha)
@@ -139,7 +139,7 @@
m_x = new double[3];
- if ( m_iq != 0 && m_B != 0 ) { // charged tracks in field
+ if ( m_iq != 0 && bFieldZ != 0 ) { // charged tracks in field
double alp = - m_iq * alpha;
@@ -149,9 +149,9 @@
} else { // neutrals or no field
- m_x[0] = m_r0[0] + alpha*m_p[0]/m_ptot;
- m_x[1] = m_r0[1] + alpha*m_p[1]/m_ptot;
- m_x[2] = m_r0[2] + alpha*m_p[2]/m_ptot;
+ 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;
@@ -212,7 +212,7 @@
/**
* Set minimum point number on a full circle.
* @param nPoints number of points on a full circle
- * @see #asPoints
+ * @see #asPoints
*/
public void setPointDensity(int nPoints)
{
@@ -261,112 +261,93 @@
*/
public List<double[]> asPoints(double alphaMax)
{
- double refRc = 100; // 1 m
-
- List<double[]> vp = new ArrayList<double[]>();
-
- if ( alphaMax < 0 ) {
- double tanL = getTanL();
- double sinL = getSinL();
- if ( Math.abs(m_p[2]) > HelixSwim.PTINY
- && tanL != 0. && m_iq != 0 ) {
- alphaMax = Math.abs((m_zMax+Math.abs(m_r0[2]))/tanL);
- } else if ( sinL != 0. && m_iq == 0 ) {
- alphaMax = Math.abs((m_zMax+Math.abs(m_r0[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;
- }
-
- /*
- * private methods
- */
- private void calcCenter()
- {
- isValid = false;
-
- m_alpha0 = m_phi + 0.5*m_iq*Math.PI;
-
- if ( m_iq != 0 && m_B != 0. ) {
-
- m_c = mm/(0.3*m_B);
- m_rc = m_c*m_pt;
+ // 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 = m_r0[0] - m_rc*Math.cos(m_alpha0);
- m_yc = m_r0[1] - m_rc*Math.sin(m_alpha0);
-
- } else {
- m_rc = -1.;
- }
- m_z0 = m_r0[2];
-
- centerOK = true;
+ 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;
- private final static int mm = 1000;
+ // one meter in org.lcsim standard units
+ private final static int ONE_METER = 1000;
- /*
- * private members
- */
// input: environment and track data
- private double m_B;
- private double[] m_p;
- private double[] m_r0;
+ private double bFieldZ;
+ private double[] momentum;
+ private double[] referencePoint;
+ // the charge in multiples of the positron charge
private int m_iq;
// cylinder geomtry
@@ -374,7 +355,8 @@
private double m_rhoMax = 1.e6;
private double m_zMax = 1.e6;
- // no points on circle with radius of 1 m
+ // 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
@@ -385,6 +367,8 @@
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
