lcsim/sandbox/Partridge
diff -N ScatterLayer_1.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ScatterLayer_1.java 4 Feb 2008 20:53:19 -0000 1.1
@@ -0,0 +1,58 @@
+/*
+ * ScatterLayer.java
+ *
+ * Created on August 30, 2007, 1:41 PM
+ *
+ * To change this template, choose Tools | Template Manager
+ * and open the template in the editor.
+ */
+
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
+
+/**
+ *
+ * @author partridge
+ */
+public class ScatterLayer_1 {
+ private MaterialLayer _layer;
+ private double _s;
+ private Hep3Vector _pos;
+ private double _radlen;
+ private double _angle;
+
+ /** Creates a new instance of ScatterLayer */
+ public ScatterLayer_1(MaterialLayer layer, double p, double slope, double s, Hep3Vector pos) {
+ _layer = layer;
+ _s = s;
+ _pos = pos;
+ if (_layer.beflag() == BarrelEndcapFlag.BARREL) {
+ _radlen = _layer.thickness() * Math.sqrt(1+Math.pow(slope,2)) / _layer.X0();
+ } else if (_layer.beflag() == BarrelEndcapFlag.ENDCAP_NORTH || _layer.beflag() == BarrelEndcapFlag.ENDCAP_NORTH) {
+ _radlen = _layer.thickness() * Math.sqrt(1+Math.pow(slope,-2)) / _layer.X0();
+ }
+ // Theta = (0.0136 GeV)*(z/p)*(sqrt(radLength))*(1+0.038*log(radLength))
+ _angle = 0.0136 * (1.0 / p) * Math.sqrt(_radlen) * (1.0 + 0.038 * Math.log(_radlen));
+ }
+
+ public MaterialLayer layer() {
+ return _layer;
+ }
+
+ public double s() {
+ return _s;
+ }
+
+ public Hep3Vector position() {
+ return _pos;
+ }
+
+ public double radlen() {
+ return _radlen;
+ }
+
+ public double angle() {
+ return _angle;
+ }
+
+}
\ No newline at end of file
lcsim/sandbox/Partridge
diff -u -r1.3 -r1.4
--- HelicalTrackFitter.java 29 Nov 2007 19:42:40 -0000 1.3
+++ HelicalTrackFitter.java 4 Feb 2008 20:53:19 -0000 1.4
@@ -4,7 +4,7 @@
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFitter.java,v 1.3 2007/11/29 19:42:40 partridge Exp $
+ * $Id: HelicalTrackFitter.java,v 1.4 2008/02/04 20:53:19 partridge Exp $
*/
import hep.physics.matrix.SymmetricMatrix;
@@ -26,7 +26,7 @@
/**
* Fit a helix to a set of space points. First, a circle is fit to the x-y coordinates.
- * A straight-line fit is then performed on the s-z coordinates. If there are two few
+ * A straight-line fit is then performed on the s-z coordinates. If there are t0o few
* 3D hits to perform the s-z fit, the ZSegmentFitter is used to find the s-z parameters.
*
* The r-phi and z coordinate measurements are assumed to be uncorrelated. A block
@@ -129,7 +129,9 @@
boolean success = _cfitter.fit(x, y, wrphi, nc);
if(!success) return false;
- // Store the circle fit output
+ // Store the circle fit output, taking into that the circle fitter uses
+ // the "wrong" sign convention for the parameter d0. The off-diagonal
+ // covariance matrix elements involving d0 also have the wrong sign.
CircleFit cfit = _cfitter.getfit();
chisq[0] = cfit.chisq();
ndof[0] = nc - 3;
@@ -137,9 +139,9 @@
par[1] = cfit.phi();
par[2] = cfit.curvature();
cov.setElement(0, 0, cfit.cov()[0]);
- cov.setElement(0, 1, cfit.cov()[1]);
+ cov.setElement(0, 1, -cfit.cov()[1]);
cov.setElement(1, 1, cfit.cov()[2]);
- cov.setElement(0, 2, cfit.cov()[3]);
+ cov.setElement(0, 2, -cfit.cov()[3]);
cov.setElement(1, 2, cfit.cov()[4]);
cov.setElement(2, 2, cfit.cov()[5]);
System.out.println(cfit.toString());
lcsim/sandbox/Partridge/SegmentTest
diff -u -r1.1 -r1.2
--- SegmentTest.java 29 Nov 2007 19:43:36 -0000 1.1
+++ SegmentTest.java 4 Feb 2008 20:53:19 -0000 1.2
@@ -68,6 +68,8 @@
aida.histogram1D("Pull for z0", 200, -3., 3.).fill((z0f-z0)/z0err);
aida.histogram1D("Pull for slope", 200, -3., 3.).fill((slopef-slope)/slopeerr);
aida.histogram1D("Polygon sides", 15, 0, 15).fill(fit.getPolygon().size());
+ aida.cloud1D("Correlation").fill((slopef-slope)*(z0f-z0));
+ aida.cloud1D("Correlation pull").fill((slopef-slope)*(z0f-z0)/fit.getCovariance().e(0,1));
} else {System.out.println("Fit failed");
}
}
lcsim/src/org/lcsim/contrib/seedtracker
diff -N Scrap.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ Scrap.java 4 Feb 2008 20:53:19 -0000 1.1
@@ -0,0 +1,65 @@
+ public double MSError(double p, double d0, double z0, double r, double z) {
+ double slope = (r-d0) / (z-z0);
+ double mserrorsq = 0.;
+ for (MaterialLayer lyr : _matlyrs) {
+ if (inside(d0,z0,lyr) != inside(r,z,lyr)) {
+ if (lyr.beflag() == BarrelEndcapFlag.BARREL) {
+ double zint = z0 + (lyr.rmax()-d0) * slope;
+ if ((zint > lyr.zmin()) && (zint < lyr.zmax())) {
+ double pathsq = Math.pow(r-lyr.rmax(),2)+Math.pow(z-zint,2);
+ double radlen = lyr.thickness() * Math.sqrt(1+Math.pow(slope,2)) / lyr.X0();
+ mserrorsq += msanglesq(p,radlen) * pathsq;
+ }
+ } else if (lyr.beflag() == BarrelEndcapFlag.ENDCAP_NORTH) {
+ double rint = d0 + (lyr.zmax()-z0) / slope;
+ if ((rint > lyr.rmin()) && (rint < lyr.rmax())) {
+ double pathsq = Math.pow(r-rint,2)+Math.pow(z-lyr.zmax(),2);
+ double radlen = lyr.thickness() * Math.sqrt(1+Math.pow(slope,-2)) / lyr.X0();
+ mserrorsq += msanglesq(p,radlen) * pathsq;
+ }
+ } else if (lyr.beflag() == BarrelEndcapFlag.ENDCAP_SOUTH) {
+ double rint = d0 + (lyr.zmin()-z0) / slope;
+ if ((rint > lyr.rmin()) && (rint < lyr.rmax())) {
+ double pathsq = Math.pow(r-rint,2)+Math.pow(z-lyr.zmin(),2);
+ double radlen = lyr.thickness() * Math.sqrt(1+Math.pow(slope,-2)) / lyr.X0();
+ mserrorsq += msanglesq(p,radlen) * pathsq;
+ }
+ }
+ }
+ }
+ return Math.sqrt(mserrorsq);
+ }
+
+ public double MSError(SeedCandidate seed, double bfield, TrackerHit hit) {
+ double[] pars = seed.getHelix().parameters();
+ double p = 0.3 * bfield * Math.abs(pars[2]) / Math.sqrt(1+Math.pow(pars[4],-2));
+ double[] pos = hit.getPosition();
+ double r = Math.sqrt(Math.pow(pos[0],2)+Math.pow(pos[1],2));
+ return MSError(p,pars[0],pars[3],r,pos[2]);
+ }
+
+
+ private boolean inside(double r, double z, MaterialLayer lyr) {
+ if (lyr.beflag() == BarrelEndcapFlag.BARREL) {
+ return r < lyr.rmin();
+ } else if (lyr.beflag() == BarrelEndcapFlag.ENDCAP_NORTH) {
+ return z < lyr.zmin();
+ } else {
+ return z > lyr.zmax();
+ }
+ }
+
+ private boolean outside(double r, double z, MaterialLayer lyr) {
+ if (lyr.beflag() == BarrelEndcapFlag.BARREL) {
+ return r > lyr.rmax();
+ } else if (lyr.beflag() == BarrelEndcapFlag.ENDCAP_NORTH) {
+ return z > lyr.zmax();
+ } else {
+ return z < lyr.zmin();
+ }
+ }
+
+ private double msanglesq(double p, double radlength) {
+ // Theta = (0.0136 GeV)*(z/p)*(sqrt(radLength))*(1+0.038*log(radLength))
+ return radlength * Math.pow(0.0136*(1.0 + 0.038 * Math.log(radlength))/p,2);
+ }
