Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
MPAlignmentParameters.java | +200 | -5 | 1.16 -> 1.17 |
New gl der calculation.
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
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