Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
RunMPAlignment.java | +1 | -1 | 1.9 -> 1.10 |
MPAlignmentParameters.java | +54 | -232 | 1.14 -> 1.15 |
+55 | -233 |
Added numerical der cross checks. Added cleaner debug code. Cleaned up old code.
diff -u -r1.9 -r1.10 --- RunMPAlignment.java 23 Oct 2012 01:31:17 -0000 1.9 +++ RunMPAlignment.java 24 Oct 2012 06:47:54 -0000 1.10 @@ -78,7 +78,7 @@
ap.setHideFrame(hideFrame); double bfield = detector.getFieldMap().getField(new BasicHep3Vector(0., 0., 1.)).y(); System.out.printf("%s: B-field in z %.3f\n",this.getClass().getSimpleName(),bfield);
- ap.setBfield(bfield);
+ ap.setUniformZFieldStrength(bfield);
loadResidualLimits(); }
diff -u -r1.14 -r1.15 --- MPAlignmentParameters.java 23 Oct 2012 01:31:17 -0000 1.14 +++ MPAlignmentParameters.java 24 Oct 2012 06:47:54 -0000 1.15 @@ -15,7 +15,6 @@
import java.util.logging.Logger; import org.lcsim.detector.tracker.silicon.SiSensor; import org.lcsim.event.RawTrackerHit;
-import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.Track; import org.lcsim.event.TrackerHit; import org.lcsim.fit.helicaltrack.*;
@@ -49,9 +48,8 @@
*/ public class MPAlignmentParameters {
- private int _nlc = 5; //the five track parameters - private int _ngl = 1; //delta(u) and delta(gamma) for each plane - private List<GlobalParameter> _glp = new ArrayList<GlobalParameter>();
+ private final int _nTrackParameters = 5; //the five track parameters + private List<GlobalParameter> _glp;
private BasicMatrix _dfdq; private HelicalTrackFit _trk; private double[] _resid = new double[3];
@@ -60,14 +58,14 @@
PrintWriter pWriter; private String _type = "LOCAL"; boolean _DEBUG=false;
- double smax = 1e3; - Hep3Vector _bfield = null; - ResLimit _resLimits = new ResLimit();
+ private Hep3Vector _bfield; + private ResLimit _resLimits;
private HPSTransformations _detToTrk;
- AlignmentUtils _alignUtils; - TrackUtils trackUtil; - TrackerHitUtils trackerHitUtil; - WTrackUtils wutils = new WTrackUtils();
+ private AlignmentUtils _alignUtils; + private AlignmentUtils.OldAlignmentUtils _oldAlignUtils; + private AlignmentUtils.NumDerivatives _numDerivatives; + private TrackUtils trackUtil; + private TrackerHitUtils trackerHitUtil;
private boolean _includeMS = true;
@@ -82,8 +80,6 @@
IDataPointSet dps_pull_b; IPlotter plotter_resuydiff_t; IPlotter plotter_resuydiff_b;
- IDataPointSet[] dps_resuydiff_t = new IDataPointSet[10]; - IDataPointSet[] dps_resuydiff_b = new IDataPointSet[10];
@@ -92,13 +88,15 @@
public MPAlignmentParameters(String outfile) {
-
+ _glp = new ArrayList<GlobalParameter>(); + _resLimits = new ResLimit();
_detToTrk = new HPSTransformations(); _alignUtils = new AlignmentUtils(_DEBUG);
+ _oldAlignUtils = new AlignmentUtils(_DEBUG).new OldAlignmentUtils(); + _numDerivatives = new AlignmentUtils(_DEBUG).new NumDerivatives();
trackUtil = new TrackUtils(); trackerHitUtil = new TrackerHitUtils(_DEBUG); _bfield = new BasicHep3Vector(0., 0., 1.);
- this.wutils.setDebug(this._DEBUG);
try { fWriter = new FileWriter(outfile); pWriter = new PrintWriter(fWriter);
@@ -108,9 +106,6 @@
makeAlignmentPlots();
- - -
} public void setResLimits(int l,int d, double low,double high) {
@@ -130,7 +125,7 @@
this._DEBUG = debug; }
- public void setBfield(double bfield) {
+ public void setUniformZFieldStrength(double bfield) {
_bfield = new BasicHep3Vector(0., 0., bfield); }
@@ -163,16 +158,19 @@
if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": this hit is not a cross"); continue; }
+ + // Update the hit position to be sure it's the latest
HelicalTrackCross cross = (HelicalTrackCross) htc; cross.setTrackDirection(_trk);
+ // Get the MS errors
double msdrphi = msmap.get(htc).drphi(); double msdz = msmap.get(htc).dz();
- double sHit = _trk.PathMap().get(htc); -
+
+ // Find the strip clusters for this 3D hit
List<HelicalTrackStrip> clusterlist = cross.getStrips();
-
+
if(_DEBUG) { System.out.printf("%s: This hit has %d clusters msdrphi=%.4f msdz=%.4f\n",this.getClass().getSimpleName(),clusterlist.size(),msdrphi,msdz); }
@@ -194,39 +192,12 @@
if(itrack%50==0) this.updatePlots(); }
- private BasicMatrix CalculateLocalDerivativesOld(double xint) { - //get track parameters. - double d0 = _trk.dca(); - double z0 = _trk.z0(); - double slope = _trk.slope(); - double phi0 = _trk.phi0(); - double R = _trk.R(); -//strip origin is defined in the tracking coordinate system (x=beamline) - double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0); - double phi = -s/R + phi0; - double[][] dfdq = new double[3][5]; - //dx/dq - //these are wrong for X, but for now it doesn't matter - dfdq[0][0] = Math.sin(phi0); - dfdq[0][1] = 0; - dfdq[0][2] = 0; - dfdq[0][3] = d0 * Math.cos(phi0) + R * Math.sin(phi0) - s * Math.cos(phi0); - dfdq[0][4] = (phi - phi0) * Math.cos(phi0); - double[] mydydq = dydq(R, d0, phi0, xint, s); - double[] mydzdq = dzdq(R, d0, phi0, xint, slope, s); - for (int i = 0; i < 5; i++) { - dfdq[1][i] = mydydq[i]; - dfdq[2][i] = mydzdq[i]; - } - - BasicMatrix dfdqGlobal = FillMatrix(dfdq, 3, 5); - return dfdqGlobal; - }
private void CalculateLocalDerivatives(HelicalTrackStrip strip) { double xint = strip.origin().x();
- BasicMatrix dfdqGlobalOld = CalculateLocalDerivativesOld(xint);
+ BasicMatrix dfdqGlobalOld = _oldAlignUtils.calculateLocalHelixDerivatives(_trk, xint);
BasicMatrix dfdqGlobal = _alignUtils.calculateLocalHelixDerivatives(_trk, xint);
+ BasicMatrix dfdqGlobalNum = this._numDerivatives.calculateLocalHelixDerivatives(_trk, xint);
Hep3Matrix trkToStrip = trackerHitUtil.getTrackToStripRotation(strip); _dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal);
@@ -238,7 +209,7 @@
double slope = _trk.slope(); double phi0 = _trk.phi0(); double R = _trk.R();
- double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
+ double s = HelixUtils.PathToXPlane(_trk, xint, 0, 0).get(0);
double phi = -s/R + phi0; double[] trackpars = {d0, z0, slope, phi0, R, s, xint}; System.out.printf("%s: --- CalculateLocalDerivatives Result --- \n",this.getClass().getSimpleName());
@@ -246,9 +217,11 @@
System.out.printf("%s: %10s%10s%10s%10s%10s%10s%10s\n",this.getClass().getSimpleName(),"d0","z0","slope","phi0","R", "xint", "s"); System.out.printf("%s: %10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n",this.getClass().getSimpleName(), d0, z0, slope, phi0, R,xint,s); System.out.printf("%s: Local derivatives:\n",this.getClass().getSimpleName());
- System.out.printf("%s: \n%s\n",this.getClass().getSimpleName(),dfdqGlobal.toString());
+ System.out.printf("%s\n",dfdqGlobal.toString()); + System.out.printf("%s: Numerical Local derivatives:\n",this.getClass().getSimpleName()); + System.out.printf("%s\n",dfdqGlobalNum.toString());
System.out.printf("%s: OLD Local derivatives:\n",this.getClass().getSimpleName());
- System.out.printf("%s: \n%s\n",this.getClass().getSimpleName(),dfdqGlobalOld.toString());
+ System.out.printf("%s\n",dfdqGlobalOld.toString());
} }
@@ -771,7 +744,7 @@
if(_DEBUG) { System.out.printf("%s: Finding interception point for residual calculation (B=%s)\n",this.getClass().getSimpleName(),this._bfield.toString()); double xint_wrong = trackUtil.calculateHelixInterceptXPlane(_trk, strip);
- double s_wrong = HelixUtils.PathToXPlane(_trk, xint_wrong, smax, _nlc).get(0);
+ double s_wrong = HelixUtils.PathToXPlane(_trk, xint_wrong, 0, 0).get(0);
Hep3Vector trkpos_wrong= HelixUtils.PointOnHelix(_trk, s_wrong); System.out.printf("%s: using wrong method xint_wrong=%.3f s_wrong=%.3f -> trkpos_wrong=%s\n",this.getClass().getSimpleName(),xint_wrong,s_wrong,trkpos_wrong.toString()); }
@@ -796,7 +769,7 @@
double xint = trkpos.x(); double phi0 = _trk.phi0(); double R = _trk.R();
- double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
+ double s = HelixUtils.PathToXPlane(_trk, xint, 0, 0).get(0);
double phi = -s/R + phi0;
@@ -832,16 +805,16 @@
if (_DEBUG) { System.out.printf("%s: ---- CalculateResidual ----\n",this.getClass().getSimpleName());
- System.out.printf("%s: Strip Origin: %s\n",this.getClass(),corigin.toString()); - System.out.printf("%s: Position on the track (tracking coordinates) at the strip: %s\n",this.getClass(),trkpos.toString()); - System.out.printf("%s: vdiffTrk %s\n",this.getClass(),vdiffTrk.toString()); - System.out.printf("%s: vdiff %s\n",this.getClass(),vdiff.toString()); - System.out.printf("%s: u %s\n",this.getClass(),u.toString()); - System.out.printf("%s: umeas = %.4f umc = %.4f\n",this.getClass(),umeas,umc); - System.out.printf("%s: udiff = %.4f +/- %.4f ( uError=%.4f , msuError=%.4f\n",this.getClass(),_resid[0],_error[0],uError,msuError); - System.out.printf("%s: MS: drdphi=%.4f, dz=%.4f\n",this.getClass(),msdrdphi,msdz); - System.out.printf("%s: MS: phi=%.4f => msvec=%s\n",this.getClass(),phi,mserr.toString()); - System.out.printf("%s: MS: msuError = %.4f (msvec*u = %s * %s\n",this.getClass(),msuError,mserr.toString(),u.toString());
+ System.out.printf("%s: Strip Origin: %s\n",this.getClass().getSimpleName(),corigin.toString()); + System.out.printf("%s: Position on the track (tracking coordinates) at the strip: %s\n",this.getClass().getSimpleName(),trkpos.toString()); + System.out.printf("%s: vdiffTrk %s\n",this.getClass().getSimpleName(),vdiffTrk.toString()); + System.out.printf("%s: vdiff %s\n",this.getClass().getSimpleName(),vdiff.toString()); + System.out.printf("%s: u %s\n",this.getClass().getSimpleName(),u.toString()); + System.out.printf("%s: umeas = %.4f umc = %.4f\n",this.getClass().getSimpleName(),umeas,umc); + System.out.printf("%s: udiff = %.4f +/- %.4f ( uError=%.4f , msuError=%.4f\n",this.getClass().getSimpleName(),_resid[0],_error[0],uError,msuError); + System.out.printf("%s: MS: drdphi=%.4f, dz=%.4f\n",this.getClass().getSimpleName(),msdrdphi,msdz); + System.out.printf("%s: MS: phi=%.4f => msvec=%s\n",this.getClass().getSimpleName(),phi,mserr.toString()); + System.out.printf("%s: MS: msuError = %.4f (msvec*u = %s * %s\n",this.getClass().getSimpleName(),msuError,mserr.toString(),u.toString());
} }
@@ -869,7 +842,7 @@
double R = _trk.R(); //double xint = strip.origin().x(); double xint = trackUtil.calculateHelixInterceptXPlane(_trk, strip);
- double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
+ double s = HelixUtils.PathToXPlane(_trk, xint, 0, 0).get(0);
double phi = -s/R + phi0; Hep3Vector posOnHelix = HelixUtils.PointOnHelix(_trk, s); if(_DEBUG) System.out.println("phi0 " + phi0 + " R " + R + " xint " + xint + " s " + s + " phi " + phi);
@@ -898,7 +871,7 @@
System.out.printf("Errors (u,v,w) : %5.5e %5.5e %5.5e\n", _error[0], _error[1], _error[2]); String[] q = {"d0", "z0", "slope", "phi0", "R"}; System.out.println("track parameter derivatives");
- for (int i = 0; i < _nlc; i++) {
+ for (int i = 0; i < _nTrackParameters; i++) {
System.out.printf("%s %5.5e %5.5e %5.5e\n", q[i], _dfdq.e(0, i), _dfdq.e(1, i), _dfdq.e(2, i)); } //String[] p = {"u-displacement"};
@@ -933,7 +906,7 @@
//// continue; //// } pWriter.printf("res_%s %5.5e %5.5e \n", d[j], _resid[j], _error[j]);
- for (int i = 0; i < _nlc; i++) {
+ for (int i = 0; i < _nTrackParameters; i++) {
pWriter.printf("lc_%s %5.5e \n", d[j], _dfdq.e(j, i)); } for (GlobalParameter gl: _glp) {
@@ -961,7 +934,7 @@
// double R = _trk.R(); // double xint = trackUtil.calculateHelixInterceptXPlane(_trk, strip); // //double xint = strip.origin().x();
-// double pathLength = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
+// double pathLength = HelixUtils.PathToXPlane(_trk, xint, 0, 0).get(0);
// double phi = -s/R + phi0; aida.histogram2D("respull_slope_"+d[j]+"_layer" + strip.layer() + "_" + side).fill(_trk.slope(),_resid[j]/_error[j]);
@@ -973,114 +946,21 @@
- private BasicMatrix FillMatrix(double[][] array, int nrow, int ncol) { - BasicMatrix retMat = new BasicMatrix(nrow, ncol); - for (int i = 0; i < nrow; i++) { - for (int j = 0; j < ncol; j++) { - retMat.setElement(i, j, array[i][j]); - } - } - return retMat; - }
+// private BasicMatrix FillMatrix(double[][] array, int nrow, int ncol) { +// BasicMatrix retMat = new BasicMatrix(nrow, ncol); +// for (int i = 0; i < nrow; i++) { +// for (int j = 0; j < ncol; j++) { +// retMat.setElement(i, j, array[i][j]); +// } +// } +// return retMat; +// }
public void closeFile() throws IOException { pWriter.close(); fWriter.close(); }
- - private double dsdR(double R, double d0, double phi0, double xint) { - double sqrtTerm = Sqrt(R * R - Math.pow(((d0 - R) * Sin(phi0) + xint), 2)); - - double rsign = Math.signum(R); - double dsdr = (1 / sqrtTerm) * ((-rsign * xint) + (-rsign) * d0 * Sin(phi0) - + ArcTan(R * Cos(phi0), (-R) * Sin(phi0)) - * sqrtTerm - - ArcTan(rsign * sqrtTerm, xint + (d0 - R) * Sin(phi0)) - * sqrtTerm); - - -// if (_DEBUG) -// System.out.println("xint = " + xint + "; dsdr = " + dsdr); - return dsdr; - - } - - private double dsdphi(double R, double d0, double phi0, double xint) { - double sqrtTerm = Sqrt(R * R - Math.pow(((d0 - R) * Sin(phi0) + xint), 2)); - double rsign = Math.signum(R); - double dsdphi = R * (sqrtTerm + rsign * d0 * Cos(phi0) - rsign * R * Cos(phi0)) / sqrtTerm; -// if (_DEBUG) -// System.out.println("xint = " + xint + "; dsdphi = " + dsdphi); - return dsdphi; - } - - private double dsdd0(double R, double d0, double phi0, double xint) { - double sqrtTerm = Sqrt(R * R - Math.pow(((d0 - R) * Sin(phi0) + xint), 2)); - double rsign = Math.signum(R); - double dsdd0 = rsign * (R * Sin(phi0)) / sqrtTerm; -// if (_DEBUG) -// System.out.println("xint = " + xint + "; dsdd0 = " + dsdd0); - return dsdd0; - } - - private double[] dydq(double R, double d0, double phi0, double xint, double s) { - double[] dy = new double[5]; -// dy[0] = Cos(phi0) + Cot(phi0 - s / R) * Csc(phi0 - s / R) * dsdd0(R, d0, phi0, xint); - dy[0] = Cos(phi0) - Sec(phi0 - s / R) * Tan(phi0 - s / R) * dsdd0(R, d0, phi0, xint); - dy[1] = 0; - dy[2] = 0; -// dy[3] = (-(d0 - R)) * Sin(phi0) - R * Cot(phi0 - s / R) * Csc(phi0 - s / R) * (1 - dsdphi(R, d0, phi0, xint) / R); - dy[3] = (-(d0 - R)) * Sin(phi0) + Sec(phi0 - s / R) * Tan(phi0 - s / R) * (R - dsdphi(R, d0, phi0, xint)); - // dy[4] = -Cos(phi0) + Csc(phi0 - s / R) - R * Cot(phi0 - s / R) * Csc(phi0 - s / R) * (s / (R * R) - dsdR(R, d0, phi0, xint) / R); - dy[4] = -Cos(phi0) + Sec(phi0 - s / R) + (1 / R) * Sec(phi0 - s / R) * Tan(phi0 - s / R) * (s - R * dsdR(R, d0, phi0, xint)); - return dy; - } - - private double[] dzdq(double R, double d0, double phi0, double xint, double slope, double s) { - double[] dz = new double[5]; - dz[0] = slope * dsdd0(R, d0, phi0, xint); - dz[1] = 1; - dz[2] = s; - dz[3] = slope * dsdphi(R, d0, phi0, xint); - dz[4] = slope * dsdR(R, d0, phi0, xint); - return dz; - } - - private double Csc(double val) { - return 1 / Math.sin(val); - } - - private double Cot(double val) { - return 1 / Math.tan(val); - } - - private double Sec(double val) { - return 1 / Math.cos(val); - } - - private double Sin(double val) { - return Math.sin(val); - } - - private double Cos(double val) { - return Math.cos(val); - } - - private double Tan(double val) { - return Math.tan(val); - } - - private double ArcTan(double val1, double val2) { - return Math.atan2(val1, val2); - } - - private double Sign(double val) { - return Math.signum(val); - } - - private double Sqrt(double val) { - return Math.sqrt(val); - }
+
private boolean isAllowedResidual(int side,int layer,int d,double res) {
@@ -1220,11 +1100,6 @@
- - - - -
}
@@ -1263,48 +1138,25 @@
//plotter_resuydiff_b.setTitle("distr: res u vs ydiff"); for(int iLayer=1;iLayer<11;++iLayer) { IProfile1D prf_t = aida.profile1D("res_u_vs_ydiff_layer_"+iLayer+"_top",10,-50,50);
- IProfile1D prf_b = aida.profile1D("res_u_vs_ydiff_layer_"+iLayer+"_bottom",10,-50,50); -// for(int i=-2;i<3;++i) { -// IHistogram ht = aida.histogram1D("res_u_vs_ydiff_"+i+"_layer_"+iLayer+"_top",50,-0.5,0.5); -// IHistogram hb = aida.histogram1D("res_u_vs_ydiff_"+i+"_layer_"+iLayer+"_bottom",50,-0.5,0.5); -// if(i==0) plotter_resuydiff_b.region(iLayer-1).plot(ht); -// else if(i==99990) plotter_resuydiff_b.region(iLayer-1).plot(ht,"mode=overlay"); -// } -
+ IProfile1D prf_b = aida.profile1D("res_u_vs_ydiff_layer_"+iLayer+"_bottom",10,-50,50);
plotter_resuydiff_t.region(iLayer-1).plot(prf_t); plotter_resuydiff_b.region(iLayer-1).plot(prf_b);
-
((PlotterRegion) plotter_resuydiff_t.region(iLayer-1)).getPlot().setAllowUserInteraction(true); ((PlotterRegion) plotter_resuydiff_t.region(iLayer-1)).getPlot().setAllowPopupMenus(true); ((PlotterRegion) plotter_resuydiff_b.region(iLayer-1)).getPlot().setAllowUserInteraction(true); ((PlotterRegion) plotter_resuydiff_b.region(iLayer-1)).getPlot().setAllowPopupMenus(true);
- - //plotter_resuydiff_t.style().statisticsBoxStyle().setVisible(false); -// dps_resuydiff_t[iLayer-1] = dpsf.create("dps_resuydiff_layer_"+iLayer+"_t", "L"+iLayer+"u residual vs sensor v position (top)",2); -// plotter_resuydiff.region(iLayer-1).plot(dps_resuydiff_t[iLayer-1]); -// dps_resuydiff_b[iLayer-1] = dpsf.create("dps_resuydiff_layer_"+iLayer+"_b", "L"+iLayer+"u residual vs sensor v position (bot)",2); -// plotter_resuydiff.region(iLayer-1).plot(dps_resuydiff_b[iLayer-1]); -// plotter_resuydiff.style().statisticsBoxStyle().setVisible(false);
}
- - -
plotterFrameSummary.addPlotter(plotter_resuydiff_t); plotterFrameSummary.addPlotter(plotter_resuydiff_b);
- - -
plotterFrame.pack(); plotterFrame.setVisible(!hideFrame); plotterFrameSummary.pack(); plotterFrameSummary.setVisible(!hideFrame);
- -
}
@@ -1359,36 +1211,6 @@
}
- - /* - for(int iLayer=1;iLayer<11;++iLayer) { - dps_resuydiff_t[iLayer-1].clear(); - dps_resuydiff_b[iLayer-1].clear(); - for(int i=-2;i<3;++i) { - IHistogram1D h = aida.histogram1D("res_u_vs_ydiff_"+i+"_layer_"+iLayer+"_top"); - dps_resuydiff_t[iLayer-1].addPoint(); - dps_resuydiff_t[iLayer-1].point(i+2).coordinate(1).setValue(h.mean()); - double N = h.entries(); - double error = N>0 ? h.rms()/Math.sqrt(N) : 0; - double vdiff = -i*25.0; - dps_resuydiff_t[iLayer-1].point(i+2).coordinate(1).setValue(h.mean()); - dps_resuydiff_t[iLayer-1].point(i+2).coordinate(1).setErrorPlus(error); - dps_resuydiff_t[iLayer-1].point(i+2).coordinate(0).setValue(vdiff); - - IHistogram1D hb = aida.histogram1D("res_u_vs_ydiff_"+i+"_layer_"+iLayer+"_bottom"); - dps_resuydiff_b[iLayer-1].addPoint(); - dps_resuydiff_b[iLayer-1].point(i+2).coordinate(1).setValue(hb.mean()); - N = hb.entries(); - error = N>0 ? hb.rms()/Math.sqrt(N) : 0; - dps_resuydiff_b[iLayer-1].point(i+2).coordinate(1).setValue(hb.mean()); - dps_resuydiff_b[iLayer-1].point(i+2).coordinate(1).setErrorPlus(error); - dps_resuydiff_b[iLayer-1].point(i+2).coordinate(0).setValue(vdiff); - } - - } - * - */ -
}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1