hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -u -r1.16 -r1.17
--- MPAlignmentParameters.java 1 Nov 2012 15:42:34 -0000 1.16
+++ MPAlignmentParameters.java 7 Nov 2012 20:55:23 -0000 1.17
@@ -181,7 +181,8 @@
if(!"GLOBAL".equals(_type)) {
CalculateResidual(cl, msdrphi, msdz);
CalculateLocalDerivatives(cl);
- CalculateGlobalDerivatives(cl);
+ this.CalculateGlobalDerivativesSimple(cl);
+ //CalculateGlobalDerivatives(cl);
PrintStripResiduals(cl);
} else {
//CalculateLocalDerivativesGLOBAL(cl);
@@ -228,9 +229,203 @@
}
+
+
+
+ private void CalculateGlobalDerivativesSimple(HelicalTrackStrip strip) {
+
+ /*
+ * Residual in local sensor frame is defined as r = m_a-p
+ * with m_a as the alignment corrected position (u_a,v_a,w_a)
+ * and p is the predicted hit position in the local sensor frame
+ *
+ * Calcualte the derivative of dr/da where
+ * a=(delta_u,delta_v,delta_w,alpha,beta,gamma)
+ * are the alingment paramerers in the local sensor frame
+ *
+ * Factorize: dr/da=dr/dm_a*dm_a/da
+ *
+ * Start with dr/dma
+ */
+
+ //Find interception with the plane the hit belongs to i.e. the predicted hit position
+ Hep3Vector p = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()));
+ double pathLengthToInterception = HelixUtils.PathToXPlane(_trk, p.x(),0,0).get(0);
+ //Find the unit vector of the track direction
+ TrackDirection trkdir = HelixUtils.CalculateTrackDirection(_trk, pathLengthToInterception);
+ Hep3Vector t_TRACK = trkdir.Direction();
+ Hep3Vector n_TRACK = strip.w();
+ Hep3Matrix T = trackerHitUtil.getTrackToStripRotation(strip);
+ Hep3Vector t = VecOp.mult(T, t_TRACK);
+ Hep3Vector n = VecOp.mult(T, n_TRACK);
+
+
+ /*
+ * Measured position, either
+ * 1. use the measured hit position on the plane w.r.t. an origin that is centered on the sensor plane.
+ * 2. use the predicted hit position in order to get a better measure of the unmeasured directions
+ */
+ double umeas,vmeas,wmeas;
+ boolean useMeasuredHitPosition = false;
+ if(useMeasuredHitPosition) {
+ if(_DEBUG) System.out.printf("%s: using measured hit position as \"m\"\n",this.getClass().getSimpleName());
+ umeas = strip.umeas();
+ vmeas = 0.; //the un-measured (along the strip) is set to zero
+ wmeas = 0.; // the hit is on the surface of the plane
+ } else {
+ if(_DEBUG) System.out.printf("%s: using predicted hit position at %s as \"m\"\n",this.getClass().getSimpleName(),p.toString());
+ Hep3Vector p_local = VecOp.sub(p,strip.origin()); //subtract the center of the sensor in tracking frame
+ if(_DEBUG) System.out.printf("%s: vector between origin of sensor and predicted position in tracking frame: %s \n",this.getClass().getSimpleName(),p_local.toString());
+ p_local = VecOp.mult(T, p_local); //rotate to local frame
+ if(_DEBUG) System.out.printf("%s: predicted position in local frame: %s \n",this.getClass().getSimpleName(),p_local.toString());
+ umeas = p_local.x();
+ vmeas = p_local.y();
+ wmeas = p_local.z();
+ Hep3Vector res_tmp = new BasicHep3Vector(strip.umeas()-p_local.x(),0.-p_local.y(),0.-p_local.z());
+ if(_DEBUG) System.out.printf("%s: cross check residuals %s\n",this.getClass().getSimpleName(),res_tmp.toString());
+
+ }
+
+ /*
+ * Calculate the dma_da
+ */
+
+ BasicMatrix dma_da = new BasicMatrix(3,6);
+ //dma_ddeltau
+ dma_da.setElement(0, 0, 1);
+ dma_da.setElement(1, 0, 0);
+ dma_da.setElement(2, 0, 0);
+ //dma_ddeltav
+ dma_da.setElement(0, 1, 0);
+ dma_da.setElement(1, 1, 1);
+ dma_da.setElement(2, 1, 0);
+ //dma_ddeltau
+ dma_da.setElement(0, 2, 0);
+ dma_da.setElement(1, 2, 0);
+ dma_da.setElement(2, 2, 1);
+ //dma_dalpha
+ dma_da.setElement(0, 3, 0);
+ dma_da.setElement(1, 3, wmeas);
+ dma_da.setElement(2, 3, -vmeas);
+ //dma_dbeta
+ dma_da.setElement(0, 4, -wmeas);
+ dma_da.setElement(1, 4, 0);
+ dma_da.setElement(2, 4, umeas);
+ //dma_dgamma
+ dma_da.setElement(0, 5, vmeas);
+ dma_da.setElement(1, 5, -umeas);
+ dma_da.setElement(2, 5, 0);
+
+
+
+ /*
+ * Calculate the dr_dma
+ * e_ij = delta(i,j) - t_i*t_j/(t*n)
+ */
+ BasicMatrix dr_dma = new BasicMatrix(3,3);
+ double tn = VecOp.dot(t, n);
+ double[] t_arr = t.v();
+ double[] n_arr = n.v();
+ for(int i=0;i<3;++i) {
+ for(int j=0;j<3;++j) {
+ double delta_ij = i==j ? 1 : 0;
+ double elem = delta_ij - t_arr[i]*n_arr[j]/tn;
+ dr_dma.setElement(i, j, elem);
+ }
+ }
+
+
+
+ /*
+ * Calculate the dr/da=dr/dma*dma/da
+ */
+
+ BasicMatrix dr_da = (BasicMatrix) MatrixOp.mult(dr_dma, dma_da);
+
+ int layer = strip.layer();
+ int side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? 10000 : 20000;
+
+
+ if(_DEBUG) {
+ double phi = -pathLengthToInterception/_trk.R()+_trk.phi0();
+ System.out.printf("%s: --- Calculate SIMPLE global derivatives ---\n",this.getClass().getSimpleName());
+ System.out.printf("%s: Side %d, layer %d, strip origin %s\n",this.getClass().getSimpleName(),side,layer,strip.origin().toString());
+ System.out.printf("%s: %10s%10s%10s%10s%10s%10s%10s%10s%10s\n",this.getClass().getSimpleName(),"d0","z0","slope","phi0","R","xint","phi", "xint","s");
+ System.out.printf("%s: %10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n",this.getClass().getSimpleName(), _trk.dca(), _trk.z0(), _trk.slope(), _trk.phi0(), _trk.R(),p.x(),phi,p.x(),pathLengthToInterception);
+ System.out.printf("%s: Track direction t = %s\n",this.getClass().getSimpleName(),t.toString());
+ System.out.printf("%s: Plane normal n = %s\n",this.getClass().getSimpleName(),n.toString());
+ System.out.printf("%s: m=(umeas,vmeas,wmeas)=(%.3f,%.3f,%.3f)\n",this.getClass().getSimpleName(),umeas,vmeas,wmeas);
+ System.out.printf("dma_da=\n%s\ndr_dma=\n%s\ndr_da=\n%s\n",dma_da.toString(),dr_dma.toString(),dr_da.toString());
+
+ }
+
+
+
+
+ /*
+ * Prepare the derivatives for output
+ */
+
+
+
+ _glp.clear();
+
+ GlobalParameter gp_tu = new GlobalParameter("Translation in u",side,layer,1000,100,true);
+ Hep3Vector mat = new BasicHep3Vector(dr_da.e(0, 0),dr_da.e(1, 0),dr_da.e(2, 0));
+ gp_tu.setDfDp(mat);
+ _glp.add(gp_tu);
+ if (_DEBUG) {
+ gp_tu.print();
+ }
+
+ GlobalParameter gp_tv = new GlobalParameter("Translation in v",side,layer,1000,200,true);
+ mat = new BasicHep3Vector(dr_da.e(0, 1),dr_da.e(1, 1),dr_da.e(2, 1));
+ gp_tv.setDfDp(mat);
+ _glp.add(gp_tv);
+ if (_DEBUG) {
+ gp_tv.print();
+ }
+
+
+ GlobalParameter gp_tw = new GlobalParameter("Translation in w",side,layer,1000,300,true);
+ mat = new BasicHep3Vector(dr_da.e(0, 2),dr_da.e(1, 2),dr_da.e(2, 2));
+ gp_tw.setDfDp(mat);
+ _glp.add(gp_tw);
+ if (_DEBUG) {
+ gp_tw.print();
+ }
+
+ GlobalParameter gp_ta = new GlobalParameter("Rotation alpha",side,layer,2000,100,true);
+ mat = new BasicHep3Vector(dr_da.e(0, 3),dr_da.e(1, 3),dr_da.e(2, 3));
+ gp_ta.setDfDp(mat);
+ _glp.add(gp_ta);
+ if (_DEBUG) {
+ gp_ta.print();
+ }
+
+ GlobalParameter gp_tb = new GlobalParameter("Rotation beta",side,layer,2000,200,true);
+ mat = new BasicHep3Vector(dr_da.e(0, 4),dr_da.e(1, 4),dr_da.e(2, 4));
+ gp_tb.setDfDp(mat);
+ _glp.add(gp_tb);
+ if (_DEBUG) {
+ gp_tb.print();
+ }
+
+ GlobalParameter gp_tg = new GlobalParameter("Rotation gamma",side,layer,2000,300,true);
+ mat = new BasicHep3Vector(dr_da.e(0, 5),dr_da.e(1, 5),dr_da.e(2, 5));
+ gp_tg.setDfDp(mat);
+ _glp.add(gp_tg);
+ if (_DEBUG) {
+ gp_tg.print();
+ }
+
+
+
+
+
+ }
-
private void CalculateGlobalDerivatives(HelicalTrackStrip strip) {
@@ -753,7 +948,7 @@
private void PrintStripResiduals(HelicalTrackStrip strip) {
if (_DEBUG) {
- System.out.printf("%s: --- Alignment Results for this Strip ---",this.getClass().getSimpleName());
+ System.out.printf("%s: --- Alignment Results for this Strip ---\n",this.getClass().getSimpleName());
String side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom";
System.out.printf("%s: Strip layer %4d in %s at origin %s\n",this.getClass().getSimpleName(), strip.layer(), side, strip.origin().toString());
System.out.printf("%s: Residuals (u,v,w) : %5.5e %5.5e %5.5e\n",this.getClass().getSimpleName(), _resid[0], _resid[1], _resid[2]);
@@ -766,10 +961,10 @@
//String[] p = {"u-displacement"};
System.out.printf("%s: Global parameter derivatives\n",this.getClass().getSimpleName());
for (GlobalParameter gl : _glp) {
- System.out.printf("%s: %s %5.5e %5.5e %5.5e %5d %10s\n",this.getClass().getSimpleName(), "", gl.dfdp(0), gl.dfdp(1), gl.dfdp(2), gl.getLabel(),gl.getName());
+ System.out.printf("%s: %s %8.3e %5.8e %8.3e %5d \"%10s\"\n",this.getClass().getSimpleName(), "", gl.dfdp(0), gl.dfdp(1), gl.dfdp(2), gl.getLabel(),gl.getName());
//System.out.printf("%s %5.5e %5.5e %5.5e %5d\n", p[j], _dfdp.e(0, j), _dfdp.e(1, j), _dfdp.e(2, j), _globalLabel[j]);
}
- System.out.printf("%s: --- Alignment Results for this Strip ---\n",this.getClass().getSimpleName());
+ System.out.printf("%s: --- END Alignment Results for this Strip ---\n",this.getClass().getSimpleName());
}
pWriter.printf("%d\n", strip.layer());
// Loop over the three directions u,v,w and print residuals and derivatives