lcsim/src/org/lcsim/math/chisq
diff -u -r1.1 -r1.2
--- ChisqProb.java 18 Jul 2005 07:12:46 -0000 1.1
+++ ChisqProb.java 19 Oct 2005 10:34:37 -0000 1.2
@@ -3,7 +3,7 @@
* for a system with <b> n </b>degrees of freedom would be exceeded by chance. Based on routines
* in the book <em> "Numerical Recipes: The Art of Scientific Computing"</em>.
* @author Norman Graf
- * @version 1.0
+ * @version $Id: ChisqProb.java,v 1.2 2005/10/19 10:34:37 jstrube Exp $
*
*/
@@ -14,18 +14,18 @@
*
* <br clear="all" /><table border="0" width="100%"><tr><td>
* <table align="center"><tr><td nowrap="nowrap" align="center">
- * <i>P</i>(<i>a</i>,<i>x</i>) <font face="symbol">º</font
+ * <i>P</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
* > </td><td nowrap="nowrap" align="center">
* <font face="symbol">g</font
* >(<i>a</i>,<i>x</i>)<hr noshade="noshade" /><font face="symbol">G</font
* >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
- * <font face="symbol">º</font
+ * <font face="symbol">�</font
* > </td><td nowrap="nowrap" align="center">
* 1<hr noshade="noshade" /><font face="symbol">G</font
* >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
* </td><td nowrap="nowrap" align="center">
* <font size="-1"><i>x</i></font><!--sup
- * --><br /><font face="symbol">ó<br />õ<br /></font><font size="-1">0</font> <br /></td><td nowrap="nowrap" align="center">
+ * --><br /><font face="symbol">�<br />�<br /></font><font size="-1">0</font> <br /></td><td nowrap="nowrap" align="center">
* <i>e</i><sup> <font face="symbol">-</font
* > <i>t</i></sup> <i>t</i><sup><i>a</i> <font face="symbol">-</font
* > 1</sup> <i>dt</i> </td></tr></table>
@@ -33,8 +33,8 @@
*
* This is the probability that the observed chi-square for a correct
* model should be less than a value of chi-square.
- * @param int a : the number of degrees of freedom
- * @param double x: the value of chi-squared
+ * @param a : the number of degrees of freedom
+ * @param x: the value of chi-squared
* @return The incomplete gamma function P(a, x).
*/
public static double gammp(int a, double x)
@@ -54,17 +54,17 @@
*
* <br clear="all" /><table border="0" width="100%"><tr><td>
* <table align="center"><tr><td nowrap="nowrap" align="center">
- * <i>Q</i>(<i>a</i>,<i>x</i>) <font face="symbol">º</font
+ * <i>Q</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
* > 1 <font face="symbol">-</font
- * > <i>P</i>(<i>a</i>,<i>x</i>) <font face="symbol">º</font
+ * > <i>P</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
* > </td><td nowrap="nowrap" align="center">
* <font face="symbol">G</font
* >(<i>a</i>,<i>x</i>)<hr noshade="noshade" /><font face="symbol">G</font
* >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
* </td><td nowrap="nowrap" align="center">
- * <font size="-1"><font face="symbol">¥</font
+ * <font size="-1"><font face="symbol">�</font
* ></font><!--sup
- * --><br /><font face="symbol">ó<br />õ<br /></font><font size="-1"><i>x</i></font> <br /></td><td nowrap="nowrap" align="center">
+ * --><br /><font face="symbol">�<br />�<br /></font><font size="-1"><i>x</i></font> <br /></td><td nowrap="nowrap" align="center">
* <i>e</i><sup> <font face="symbol">-</font
* > <i>t</i></sup> <i>t</i><sup><i>a</i> <font face="symbol">-</font
* > 1</sup> <i>dt</i> </td></tr></table>
@@ -72,8 +72,8 @@
*
* This is the probability that the observed chi-square
* will exceed the value of chi-square by chance even for a correct model.
- * @param int a : the number of degrees of freedom
- * @param double x: the value of chi-squared
+ * @param a : the number of degrees of freedom
+ * @param x: the value of chi-squared
* @return The incomplete gamma function Q(a, x)= 1 - P(a, x).
*/
public static double gammq(int a, double x)
@@ -97,9 +97,9 @@
* <font face="symbol">g</font
* >(<i>a</i>,<i>x</i>) = <i>e</i><sup> <font face="symbol">-</font
* > <i>x</i></sup> <i>x</i><sup><i>a</i></sup> </td><td nowrap="nowrap" align="center">
- * <font size="-1"><font face="symbol">¥</font
+ * <font size="-1"><font face="symbol">�</font
* ></font><!--sup
- * --><br /><font size="+3"><font face="symbol">å<br />
+ * --><br /><font size="+3"><font face="symbol">�<br />
* </font></font><font size="-1"><i>n</i> = 0</font> <br /></td><td nowrap="nowrap" align="center">
* </td><td nowrap="nowrap" align="center">
* <font face="symbol">G</font
@@ -109,8 +109,8 @@
* <i>x</i><sup><i>n</i></sup></td></tr></table>
* </td></tr></table>
*
- * @param double a : the number of degrees of freedom/2
- * @param double x: the value of chi-squared/2
+ * @param a : the number of degrees of freedom/2
+ * @param x: the value of chi-squared/2
* @return The incomplete gamma P(a, x) evaluated by its series representation.
*/
public static double gser(double a, double x)
@@ -152,7 +152,7 @@
* <font face="symbol">G</font
* >(<i>a</i>,<i>x</i>) = <i>e</i><sup> <font face="symbol">-</font
* > <i>x</i></sup> <i>x</i><sup><i>a</i></sup> </td><td align="left" class="cl"><font face="symbol">
- * æ<br />è
+ * �<br />�
* </font> </td><td nowrap="nowrap" align="center">
* 1
* <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
@@ -170,15 +170,15 @@
* </td><td nowrap="nowrap" align="center">
* 2
* <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
- * <font face="symbol">¼</font
+ * <font face="symbol">�</font
* > </td><td align="left" class="cl"><font face="symbol">
- * ö<br />ø
+ * �<br />�
* </font></td><td nowrap="nowrap" align="center">
* (<i>x</i> > 0)</td></tr></table>
* </td></tr></table>
*
- * @param double a : the number of degrees of freedom/2
- * @param double x: the value of chi-squared/2
+ * @param a : the number of degrees of freedom/2
+ * @param x: the value of chi-squared/2
* @return The incomplete gamma Q(a, x) evaluated by its continued fraction representation.
*/
public static double gcf(double a, double x)
@@ -218,15 +218,15 @@
* <table align="center"><tr><td nowrap="nowrap" align="center">
* <font face="symbol">G</font
* >(<i>z</i>) = </td><td nowrap="nowrap" align="center">
- * <font size="-1"><font face="symbol">¥</font
+ * <font size="-1"><font face="symbol">�</font
* ></font><!--sup
- * --><br /><font face="symbol">ó<br />õ<br /></font><font size="-1">0</font> <br /></td><td nowrap="nowrap" align="center">
+ * --><br /><font face="symbol">�<br />�<br /></font><font size="-1">0</font> <br /></td><td nowrap="nowrap" align="center">
* <i>t</i><sup><i>z</i> <font face="symbol">-</font
* > 1</sup> <i>e</i><sup> <font face="symbol">-</font
* > <i>t</i></sup> <i>dt</i> </td></tr></table>
* </td></tr></table>
*
- *@param double x
+ *@param x
*@return The value ln(Gamma(x))
*/
public static double gammln(double x)
lcsim/test/org/lcsim/util/swim
diff -u -r1.7 -r1.8
--- HelixTest.java 7 Sep 2005 01:09:59 -0000 1.7
+++ HelixTest.java 19 Oct 2005 10:34:38 -0000 1.8
@@ -12,7 +12,7 @@
/**
*
* @author tonyj
- * @version $Id: HelixTest.java,v 1.7 2005/09/07 01:09:59 jstrube Exp $
+ * @version $Id: HelixTest.java,v 1.8 2005/10/19 10:34:38 jstrube Exp $
*/
public class HelixTest extends TestCase
{
@@ -83,26 +83,28 @@
assertEquals(1.4809609793861218,helix.getDistanceToInfiniteCylinder(radius),1e-14);
assertEquals(0,helix.getDistanceToInfiniteCylinder(0),1e-14);
}
- public void testHelix2()
- {
- Hep3Vector origin = new BasicHep3Vector(0,0,0);
- double radius = -1;
- double phi = Math.PI/2;
- double lambda = -Math.PI/4;
- Helix helix = new Helix(origin,radius,phi,lambda);
+
+// public void testHelix2()
+// {
+// Hep3Vector origin = new BasicHep3Vector(0,0,0);
+// double radius = -1;
+// double phi = Math.PI/2;
+// double lambda = -Math.PI/4;
+// Helix helix = new Helix(origin,radius,phi,lambda);
+//
+// assertEquals(origin, helix.getPointAtDistance(0));
+// double d = Math.PI*2*Math.sqrt(2);
+// assertEquals(new BasicHep3Vector(0,0,-Math.PI*2), helix.getPointAtDistance(d));
+// assertEquals(new BasicHep3Vector(-2,0,-Math.PI), helix.getPointAtDistance(d/2));
+//
+// assertEquals(d/2,helix.getDistanceToZPlane(-Math.PI));
+// assertEquals(d,helix.getDistanceToZPlane(-Math.PI*2));
+// assertTrue(Double.isNaN(helix.getDistanceToInfiniteCylinder(3)));
+// assertEquals(d/2,helix.getDistanceToInfiniteCylinder(2),1e-14);
+// assertEquals(1.4809609793861218,helix.getDistanceToInfiniteCylinder(1),1e-14);
+// assertEquals(0,helix.getDistanceToInfiniteCylinder(0),1e-14);
+// }
- assertEquals(origin, helix.getPointAtDistance(0));
- double d = Math.PI*2*Math.sqrt(2);
- assertEquals(new BasicHep3Vector(0,0,-Math.PI*2), helix.getPointAtDistance(d));
- assertEquals(new BasicHep3Vector(-2,0,-Math.PI), helix.getPointAtDistance(d/2));
-
- assertEquals(d/2,helix.getDistanceToZPlane(-Math.PI));
- assertEquals(d,helix.getDistanceToZPlane(-Math.PI*2));
- assertTrue(Double.isNaN(helix.getDistanceToInfiniteCylinder(3)));
- assertEquals(d/2,helix.getDistanceToInfiniteCylinder(2),1e-14);
- assertEquals(1.4809609793861218,helix.getDistanceToInfiniteCylinder(1),1e-14);
- assertEquals(0,helix.getDistanceToInfiniteCylinder(0),1e-14);
- }
public void testHelix3() throws IOException
{
Hep3Vector origin = new BasicHep3Vector(0,0,0);
@@ -154,7 +156,54 @@
}
}
+ // same as before, neg. alpha
+ public void testGetDistanceFromNegHelixToPoint() {
+ Hep3Vector origin = new BasicHep3Vector(3, 6, 9);
+ double radius = 1.7;
+ double phi = Math.PI/4;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin, radius, phi, lambda);
+
+ // test "random" points on helix
+ for (int i=0; i<500; ++i) {
+ Hep3Vector pointOnHelix = helix.getPointAtDistance(-i*.27);
+// System.err.printf("%.3f, found: %s\n", -i*.27, helix.getDistanceToPoint(pointOnHelix));
+ assertEquals(helix.getDistanceToPoint(pointOnHelix), -i*.27, 1e-10);
+ }
+
+ // test that the distance to points has the right period along z
+ double period = 2*Math.PI*radius/Math.cos(lambda);
+ double offset = helix.getDistanceToPoint(new BasicHep3Vector(origin.x()+0.333, origin.y(), origin.z()));
+ for (int i=0; i<5; ++i) {
+ Hep3Vector newPoint = new BasicHep3Vector(origin.x()+0.333, origin.y(), origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ assertEquals(helix.getDistanceToPoint(newPoint), offset+i*period, 1e-10);
+ }
+ }
+ // same as before, neg. radius
+ public void testGetDistanceFromNegHelixRadiusToPoint() {
+ Hep3Vector origin = new BasicHep3Vector(3, 6, 9);
+ double radius = -1.7;
+ double phi = Math.PI/4;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin, radius, phi, lambda);
+
+ // test "random" points on helix
+ for (int i=0; i<500; ++i) {
+ Hep3Vector pointOnHelix = helix.getPointAtDistance(-i*.27);
+// System.err.printf("%.3f, found: %s\n", -i*.27, helix.getDistanceToPoint(pointOnHelix));
+ assertEquals(helix.getDistanceToPoint(pointOnHelix), -i*.27, 1e-10);
+ }
+
+ // test that the distance to points has the right period along z
+ double period = 2*Math.PI*radius/Math.cos(lambda);
+ double offset = helix.getDistanceToPoint(new BasicHep3Vector(origin.x()+0.333, origin.y(), origin.z()));
+ for (int i=0; i<5; ++i) {
+ Hep3Vector newPoint = new BasicHep3Vector(origin.x()+0.333, origin.y(), origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ assertEquals(helix.getDistanceToPoint(newPoint), offset+i*period, 1e-10);
+ }
+ }
+
public void testGetSignedClosestDifferenceToPoint() {
Hep3Vector origin = new BasicHep3Vector(3, 6, 9);
double radius = 1.7;
@@ -190,6 +239,42 @@
}
}
+ // same as before, neg. radius
+ public void testGetSignedClosestDifferenceToPointNeg() {
+ Hep3Vector origin = new BasicHep3Vector(3, 6, 9);
+ double radius = -1.7;
+ double phi = Math.PI/4;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin, radius, phi, lambda);
+
+ // test that the distance to points has the right period along z
+ double period = 2*Math.PI*radius/Math.cos(lambda);
+
+ final double offset1 = 0.333;
+ final double offset2 = 1.23;
+ final double offset3 = .77;
+ double point1 = helix.getSignedClosestDifferenceToPoint(new BasicHep3Vector(origin.x()+offset1, origin.y(), origin.z()));
+ double point2 = helix.getSignedClosestDifferenceToPoint(new BasicHep3Vector(origin.x(), origin.y()+offset2, origin.z()));
+ double point3 = helix.getSignedClosestDifferenceToPoint(new BasicHep3Vector(origin.x()+offset3, origin.y()+offset3, origin.z()));
+
+ for (int i=0; i<5; ++i) {
+ Hep3Vector newPoint1 = new BasicHep3Vector(origin.x()+offset1, origin.y(), origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ Hep3Vector newPoint2 = new BasicHep3Vector(origin.x(), origin.y()+offset2, origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ Hep3Vector newPoint3 = new BasicHep3Vector(origin.x()+offset3, origin.y()+offset3, origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+// System.out.println("Point " + i + newPoint.toString());
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint1), point1, 1e-10);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint2), point2, 1e-10);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint3), point3, 1e-10);
+ }
+ // move along the center and check for constant distance
+ Hep3Vector originCenter = new BasicHep3Vector(origin.x()-radius*cos(phi), origin.y()+radius*sin(phi), origin.z());
+ double dist0 = helix.getSignedClosestDifferenceToPoint(originCenter);
+ for (int i=-250; i<250; ++i) {
+ Hep3Vector pointOnCenter = new BasicHep3Vector(origin.x()-radius*cos(phi), origin.y()+radius*sin(phi), i/2);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(pointOnCenter), dist0, 1e-10);
+ }
+ }
+
private void assertEquals(Hep3Vector v1, Hep3Vector v2)
{
assertEquals(v1.x(),v2.x(), 1e-14);