hps-java/src/main/java/org/lcsim/hps/users/phansson
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);
- }
-
- }
- *
- */
-
}