Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
RunMPAlignment.java | +3 | -2 | 1.3 -> 1.4 |
AlignmentUtils.java | +32 | -7 | 1.2 -> 1.3 |
MPAlignmentParameters.java | +215 | -26 | 1.5 -> 1.6 |
+250 | -35 |
Updated track based alignment: new derivatives, new debugging plots and output. Fixed bug in phi calculation when calculating residuals, small effect.
diff -u -r1.3 -r1.4 --- RunMPAlignment.java 21 Aug 2012 19:34:18 -0000 1.3 +++ RunMPAlignment.java 29 Aug 2012 02:06:22 -0000 1.4 @@ -65,7 +65,7 @@
_minLayers = 8; //// if (_config.contains("Test")) //// flipSign = -1;
- ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug);
+
}
@@ -76,13 +76,14 @@
_config = config; if (_config.contains("Test")) flipSign = -1;
- ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug);
+
} public void detectorChanged(Detector detector) {
+ ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug);
//Initialize the res limits for(int i=1;i<=10;++i) {
diff -u -r1.2 -r1.3 --- AlignmentUtils.java 27 Aug 2012 18:40:10 -0000 1.2 +++ AlignmentUtils.java 29 Aug 2012 02:06:22 -0000 1.3 @@ -93,7 +93,7 @@
} public double dx_dR(double x, double xr, double yr, double d0, double phi0, double R) {
- return Math.sin(phi0) + deltayc(x,xr,yr,d0,phi0,R)*sign(R) + sign(R)*ddeltayc_dR(x,xr,d0,phi0,R);
+ return Math.sin(phi0) + deltayc(x,xr,yr,d0,phi0,R)*dsign_dR(R) + sign(R)*ddeltayc_dR(x,xr,d0,phi0,R);
} //-----------------------------------
@@ -125,7 +125,7 @@
//Simplify with xc double num = sign(R) * ( d0 + R + (-d0+R)*Math.cos(2*phi0) + 2*(x-xr)*Math.sin(phi0) ); double den = 2*Math.sqrt( R*R - Math.pow( x - xc(xr,d0,phi0,R), 2) );
- double C = Math.sqrt( R*R + Math.pow( x - xc(xr,d0,phi0,R), 2) )*sign(R);
+ double C = Math.sqrt( R*R + Math.pow( x - xc(xr,d0,phi0,R), 2) )*dsign_dR(R);
return -Math.cos(phi0) + num/den + C; }
@@ -174,8 +174,9 @@
return dphi; }
- -
+ public double dsign_dR(double R) { + return 0.0; + }
//Generic point y on circle/helix
@@ -244,7 +245,7 @@
// derivate of deltayc public double ddeltayc_dR(double x, double xr, double d0, double phi0, double R) { //simplifying using xc
- double num = ( 2*R - sign(R)*sign(R)*( 2*R + 2*Math.sin(phi0)*( x - xc(xr,d0,phi0,R)) ) - 2*sign(R)*( R*R - Math.pow(x-xc(xr,d0,phi0,R), 2))*sign(R) );
+ double num = ( 2*R - sign(R)*sign(R)*( 2*R + 2*Math.sin(phi0)*( x - xc(xr,d0,phi0,R)) ) - 2*sign(R)*( R*R - Math.pow(x-xc(xr,d0,phi0,R), 2))*dsign_dR(R) );
double den = 2*Math.sqrt( R*R - sign(R)*sign(R)*( R*R - Math.pow(x-xc(xr,d0,phi0,R), 2) ) ); return num/den; }
@@ -266,7 +267,7 @@
} // derivative of delta phi between point on helix (x,y) and point of closest approach (x0,y0)
- public double dtmpdphi_dR(double x, double xr, double yr, double d0, double phi0, double R) {
+ public double baddtmpdphi_dR(double x, double xr, double yr, double d0, double phi0, double R) {
// this calculation has been large "simplified" from the original by using x, x0, xc, y, y0, yc double term1 = ((yc(yr,d0,phi0,R) - y0(yr,d0,phi0))*Math.sin(phi0) - (x0(xr,d0,phi0) - xc(xr,d0,phi0,R))*Math.cos(phi0) ) / (R*R);
@@ -278,7 +279,31 @@
return term1 + term2; }
-
+ + // derivative of delta phi between point on helix (x,y) and point of closest approach (x0,y0) + public double dtmpdphi_dR(double x, double xr, double yr, double d0, double phi0, double R) { + // this calculation has been large "simplified" from the original by using x, x0, xc, y, y0, yc + + double num = R*sign(R)*(x-x0(xr,d0,phi0)); + + double den = (yc(yr,d0,phi0,R) - y0(yr,d0,phi0))*(2*(x-xc(xr,d0,phi0,R))-R*R); + + return num/den; + + +// double term1 = ((yc(yr,d0,phi0,R) - y0(yr,d0,phi0))*Math.sin(phi0) +R*Math.cos(phi0) ) / (R*R); +// +// double num1 = sign(R)*Math.sin(phi0)*(y(x, xr, yr, d0, phi0, R) - y0(yr,d0,phi0)); +// +// double num2 = sign(R)*(x-xc(xr,d0,phi0,R))*(R + Math.sin(phi0)*(x-xc(xr,d0,phi0,R)))/(y(x, xr, yr, d0, phi0, R)-yc(yr,d0,phi0,R)); +// +// double den = Math.pow(x-xc(xr,d0,phi0,R),2) + sign(R)*sign(R)*Math.pow(y(x, xr, yr, d0, phi0, R) - yc(yr,d0,phi0,R),2); +// +// double term2 = (num1 + num2)/den; +// +// return term1 + term2; + } +
// derivative of delta phi between point on helix (x,y) and point of closest approach (x0,y0) public double dtmpdphi_dd0(double x, double xr, double yr, double d0, double phi0, double R) {
diff -u -r1.5 -r1.6 --- MPAlignmentParameters.java 27 Aug 2012 18:40:10 -0000 1.5 +++ MPAlignmentParameters.java 29 Aug 2012 02:06:22 -0000 1.6 @@ -1,6 +1,7 @@
package org.lcsim.hps.users.phansson; import hep.aida.*;
+import hep.aida.ref.plotter.PlotterRegion;
import org.lcsim.hps.users.mgraham.alignment.*; import hep.physics.matrix.BasicMatrix; import hep.physics.matrix.MatrixOp;
@@ -23,13 +24,7 @@
import org.lcsim.event.RawTrackerHit; import org.lcsim.event.Track; import org.lcsim.event.TrackerHit;
-import org.lcsim.fit.helicaltrack.HelicalTrackCross; -import org.lcsim.fit.helicaltrack.HelicalTrackFit; -import org.lcsim.fit.helicaltrack.HelicalTrackHit; -import org.lcsim.fit.helicaltrack.HelicalTrackStrip; -import org.lcsim.fit.helicaltrack.HelixUtils; -import org.lcsim.fit.helicaltrack.MultipleScatter; -import org.lcsim.fit.helicaltrack.TrackDirection;
+import org.lcsim.fit.helicaltrack.*;
import org.lcsim.hps.event.HPSTransformations; import org.lcsim.hps.monitoring.AIDAFrame; import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
@@ -153,7 +148,7 @@
// IHistogram h = aida.histogram1D("res_"+direction[d]+"_layer" + (iLayer) + "_" + side[iSide] , 50, -10, 10); IHistogram h = aida.histogram1D("res_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, xmin, xmax); IHistogram h1 = aida.histogram1D("reserr_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, 0, 1);
- IHistogram h2 = aida.histogram1D("respull_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, -10, 10);
+ IHistogram h2 = aida.histogram1D("respull_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, -20, 20);
int region = (iLayer-1);//*2+iSide;
@@ -161,6 +156,14 @@
plotter_reserr.region(region).plot(h1); plotter_respull.region(region).plot(h2);
+ ((PlotterRegion) plotter_res.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_res.region(region)).getPlot().setAllowPopupMenus(true); + ((PlotterRegion) plotter_reserr.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_reserr.region(region)).getPlot().setAllowPopupMenus(true); + ((PlotterRegion) plotter_respull.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_respull.region(region)).getPlot().setAllowPopupMenus(true); + +
} //plotter_res.show();
@@ -170,6 +173,103 @@
} }
+ + for (int s=0;s<2;++s) { + + IPlotter plotter_hthresy = af.createPlotterFactory().create(); + plotter_hthresy.createRegions(5,1,0); + plotter_hthresy.setTitle("hthres_y " + side[s]); + IPlotter plotter_hthresypull = af.createPlotterFactory().create(); + plotter_hthresypull.createRegions(5,1,0); + plotter_hthresypull.setTitle("hthrespull_y " + side[s]); + IPlotter plotter_hthresyerr = af.createPlotterFactory().create(); + plotter_hthresyerr.createRegions(5,1,0); + plotter_hthresyerr.setTitle("hthreserr_y " + side[s]); + IPlotter plotter_hthresz = af.createPlotterFactory().create(); + plotter_hthresz.createRegions(5,1,0); + plotter_hthresz.setTitle("hthres_z " + side[s]); + IPlotter plotter_hthreszerr = af.createPlotterFactory().create(); + plotter_hthreszerr.createRegions(5,1,0); + plotter_hthreszerr.setTitle("hthreserr_z " + side[s]); + IPlotter plotter_hthreszpull = af.createPlotterFactory().create(); + plotter_hthreszpull.createRegions(5,1,0); + plotter_hthreszpull.setTitle("hthrespull_z " + side[s]); + IPlotter plotter_hthreszerrms = af.createPlotterFactory().create(); + plotter_hthreszerrms.createRegions(5,1,0); + plotter_hthreszerrms.setTitle("hthreserrms_z " + side[s]); + IPlotter plotter_hthreszerrres = af.createPlotterFactory().create(); + plotter_hthreszerrres.createRegions(5,1,0); + plotter_hthreszerrres.setTitle("hthreserrres_z " + side[s]); + + + + + + for (int iLayer=1;iLayer<6;++iLayer) { + + + IHistogram h = aida.histogram1D("hthres_y_layer" + (iLayer) + "_" + side[s] , 50, -2, 2); + IHistogram h1 = aida.histogram1D("hthreserr_y_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); + IHistogram h2 = aida.histogram1D("hthrespull_y_layer" + (iLayer) + "_" + side[s] , 50, -5, 5); + + IHistogram hz = aida.histogram1D("hthres_z_layer" + (iLayer) + "_" + side[s] , 50, -2, 2); + IHistogram hz1 = aida.histogram1D("hthreserr_z_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); + IHistogram hz2 = aida.histogram1D("hthrespull_z_layer" + (iLayer) + "_" + side[s] , 50, -5, 5); + IHistogram hz11 = aida.histogram1D("hthreserrms_z_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); + IHistogram hz111 = aida.histogram1D("hthreserrres_z_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); + + + + + + int region = (iLayer-1);//*2+iSide; + + plotter_hthresy.region(region).plot(h); + plotter_hthresyerr.region(region).plot(h1); + plotter_hthresypull.region(region).plot(h2); + + plotter_hthresz.region(region).plot(hz); + plotter_hthreszerr.region(region).plot(hz1); + plotter_hthreszpull.region(region).plot(hz2); + + plotter_hthreszerrms.region(region).plot(hz11); + plotter_hthreszerrres.region(region).plot(hz111); + + + ((PlotterRegion) plotter_hthresy.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_hthresy.region(region)).getPlot().setAllowPopupMenus(true); + ((PlotterRegion) plotter_hthresyerr.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_hthresyerr.region(region)).getPlot().setAllowPopupMenus(true); + ((PlotterRegion) plotter_hthresypull.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_hthresypull.region(region)).getPlot().setAllowPopupMenus(true); + + ((PlotterRegion) plotter_hthresz.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_hthresz.region(region)).getPlot().setAllowPopupMenus(true); + ((PlotterRegion) plotter_hthreszerr.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_hthreszerr.region(region)).getPlot().setAllowPopupMenus(true); + ((PlotterRegion) plotter_hthreszpull.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_hthreszpull.region(region)).getPlot().setAllowPopupMenus(true); + + ((PlotterRegion) plotter_hthreszerrms.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_hthreszerrms.region(region)).getPlot().setAllowPopupMenus(true); + ((PlotterRegion) plotter_hthreszerrres.region(region)).getPlot().setAllowUserInteraction(true); + ((PlotterRegion) plotter_hthreszerrres.region(region)).getPlot().setAllowPopupMenus(true); + + + + } + + plotterFrame.addPlotter(plotter_hthresy); + plotterFrame.addPlotter(plotter_hthresyerr); + plotterFrame.addPlotter(plotter_hthresypull); + plotterFrame.addPlotter(plotter_hthresz); + plotterFrame.addPlotter(plotter_hthreszerr); + plotterFrame.addPlotter(plotter_hthreszpull); + plotterFrame.addPlotter(plotter_hthreszerrms); + plotterFrame.addPlotter(plotter_hthreszerrres); + } + +
plotterFrame.pack(); plotterFrame.setVisible(true);
@@ -203,6 +303,7 @@
double msdrphi = msmap.get(htc).drphi(); double msdz = msmap.get(htc).dz(); double sHit = _trk.PathMap().get(htc);
+ calculateTrackHitResidual(htc,msmap);
HelicalTrackCross cross = (HelicalTrackCross) htc; List<HelicalTrackStrip> clusterlist = cross.getStrips(); TrackDirection trkdir = HelixUtils.CalculateTrackDirection(_trk, sHit);
@@ -249,7 +350,8 @@
BasicMatrix dfdqGlobal = FillMatrix(dfdq, 3, 5); BasicMatrix dfdqGlobalNew = _alignUtils.calculateLocalHelixDerivatives(_trk, strip); Hep3Matrix trkToStrip = getTrackToStripRotation(strip);
- _dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal);
+ //_dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal);
+ _dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobalNew);
if (_DEBUG) { double[] trackpars = {d0, z0, slope, phi0, R, s, xint};
@@ -258,19 +360,19 @@
System.out.println("trkToStrip Rotation:"); System.out.println(trkToStrip.toString()); printDerivatives(trackpars, dfdq);
+ //} + System.out.println("-----------------------------------------------------"); + System.out.println("Compare local derivatives for strip at "); + System.out.println(strip.origin()); + System.out.println("s xint"); + System.out.printf("%5.5f %5.5f\n", s, xint); + System.out.println(" d0 z0 slope phi0 R"); + System.out.printf("Values %5.5f %5.5f %5.5f %5.5f %5.5f\n", d0, z0, slope, phi0, R); + System.out.println("dfdqGlobal:"); + System.out.println(dfdqGlobal.toString()); + System.out.println("dfdqGlobalNew:"); + System.out.println(dfdqGlobalNew.toString());
}
- double[] trackpars = {d0, z0, slope, phi0, R, s, xint}; - System.out.println("-----------------------------------------------------"); - System.out.println("Compare local derivatives for strip at "); - System.out.println(strip.origin()); - System.out.println("s xint"); - System.out.printf("%5.5f %5.5f\n", s, xint); - System.out.println(" d0 z0 slope phi0 R"); - System.out.printf("Values %5.5f %5.5f %5.5f %5.5f %5.5f\n", d0, z0, slope, phi0, R); - System.out.println("dfdqGlobal:"); - System.out.println(dfdqGlobal.toString()); - System.out.println("dfdqGlobalNew:"); - System.out.println(dfdqGlobalNew.toString());
}
@@ -431,16 +533,19 @@
double R = _trk.R(); double xint = strip.origin().x(); double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
- double phi = s / R - phi0;
+ //double phi = s / R - phi0;
+ //corrected sign PHA
+ double phi = phi0 - s / R ;
Hep3Vector trkpos = HelixUtils.PointOnHelix(_trk, s); //System.out.println("trkpos = "+trkpos.toString()); //System.out.println("origin = "+corigin.toString());
-
+
Hep3Vector mserr = new BasicHep3Vector(msdrdphi * Math.sin(phi), msdrdphi * Math.sin(phi), msdz); Hep3Vector vdiffTrk = VecOp.sub(trkpos, corigin); Hep3Matrix trkToStrip = getTrackToStripRotation(strip); Hep3Vector vdiff = VecOp.mult(trkToStrip, vdiffTrk);
+ Hep3Vector mserrrot = VecOp.mult(trkToStrip, mserr);
double umc = vdiff.x(); double vmc = vdiff.y(); double wmc = vdiff.z();
@@ -468,12 +573,96 @@
System.out.println("u :"); System.out.println(u.toString()); System.out.println("umeas = " + umeas + "; umc = " + umc);
- System.out.println("udiff = " + _resid[0] + " +/- " + _error[0]); -
+ System.out.println("udiff = " + _resid[0] + " +/- " + _error[0] + " ( uError=" + uError + ", msuError=" + msuError + ")"); + System.out.println("MS: drdphi=" + msdrdphi + ", dz=" + msdz); + System.out.println("MS: phi=" + phi + " => msvec=" + mserr.toString()); + System.out.println("MS: msuError = " + msuError + " (msvec*u = " + mserr.toString() + " * " + u.toString() + ")" ); + System.out.println("MS: msvec rotated to strip? = " + mserrrot.toString()); +
} }
+ public void calculateTrackHitResidual(HelicalTrackHit hth,Map<HelicalTrackHit, MultipleScatter> msmap) { + + double msdrphi = msmap.get(hth).drphi(); + double msdz = msmap.get(hth).dz(); + + //Calculate the residuals that are being used in the track fit + + //Start with the bendplane y + double drphi_res = hth.drphi(); + double wrphi = Math.sqrt(drphi_res*drphi_res + msdrphi*msdrphi); + double s = _trk.PathMap().get(hth); + Hep3Vector posOnHelix = HelixUtils.PointOnHelix(_trk, s); + double resy = hth.y() - posOnHelix.y(); + double resypull = resy/wrphi; + + String side = "top"; + if (hth.getPosition()[2]<0) side = "bottom"; + int layer = hth.Layer(); + if(layer % 2 == 0) { + layer = layer/2; + } else { + layer = (layer+1)/2; + } + + aida.histogram1D("hthres_y_layer" + layer + "_" + side).fill(resy); + aida.histogram1D("hthrespull_y_layer" + layer + "_" + side).fill(resypull); + aida.histogram1D("hthreserr_y_layer" + layer + "_" + side).fill(wrphi); + + + //Now the residual for the "measurement" direction z + double resz = hth.z() - posOnHelix.z(); + //double dz = HitUtils.zres(hth, msmap, _trk); + + double slope = hth.z() / hth.r(); + // Don't use the helix slope if the magnitude of the slope is smaller than its uncertainty + if (Math.abs(_trk.slope()) > _trk.getSlopeError()) slope = _trk.slope(); + // Take the resolution uncertainty to be dr * |slope| + double dzres = hth.dr() * Math.abs(slope); + // Combine resolution and multiple scattering uncertainties in quadrature + double dz = Math.sqrt(dzres*dzres + msdz*msdz); + + aida.histogram1D("hthres_z_layer" + layer + "_" + side).fill(resz); + aida.histogram1D("hthrespull_z_layer" + layer + "_" + side).fill(resz/dz); + aida.histogram1D("hthreserr_z_layer" + layer + "_" + side).fill(dz); + aida.histogram1D("hthreserrms_z_layer" + layer + "_" + side).fill(msdz); + aida.histogram1D("hthreserrres_z_layer" + layer + "_" + side).fill(dzres); + + + // +// // Store the coordinates and errors for the line fit +// for (int i = 0; i < npix; i++) { +// HelicalTrackHit hit = pixel_hits.get(i); +// z[i] = hit.z(); +// dz[i] = HitUtils.zres(hit, msmap, oldhelix); +// s[i] = smap.get(hit); +// } +// +// + + + //From HelicalTrackFitter.java + // // Store the hit coordinates and weights in arrays for the circle fitter +// for (int i = 0; i < nc; i++) { +// HelicalTrackHit hit = circle_hits.get(i); +// // Store the hit position +// x[i] = hit.x(); +// y[i] = hit.y(); +// // Find the weight (= 1/uncertainty^2) for this hit +// // First get the multiple scattering uncertainty +// double drphi_ms = 0.; +// if (msmap.containsKey(hit)) { +// drphi_ms = msmap.get(hit).drphi(); +// } +// // Get the hit resolution and combine uncertainties in quadrature +// double drphi_res = hit.drphi(); +// wrphi[i] = 1. / (drphi_res * drphi_res + drphi_ms * drphi_ms); +// } + } + +
public double[] getResidual(Track track, int layer) { double[] res = new double[7]; SeedTrack st = (SeedTrack) track;
@@ -637,7 +826,7 @@
int iside = (df.intValue() % 10) - 1; //System.out.println("track is on " + s + " gl param is " + lbl + "(" + df + ","+iside+")" ); if(iside!=s) {
- System.out.println("ERROR track is on " + s + " while gl param is " + lbl + "(" + df + ")" );
+ System.out.println("WARNING track is on " + s + " while gl param is " + lbl + "(" + df + ")" );
System.exit(1); }
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