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