hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -N MPAlignmentParameters.java
--- MPAlignmentParameters.java 14 Nov 2012 22:08:23 -0000 1.18
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,1232 +0,0 @@
-package org.lcsim.hps.users.phansson;
-
-import hep.aida.*;
-import hep.aida.ref.plotter.PlotterRegion;
-import hep.physics.matrix.BasicMatrix;
-import hep.physics.matrix.MatrixOp;
-import hep.physics.vec.*;
-import java.io.FileWriter;
-import java.io.IOException;
-import java.io.PrintWriter;
-import java.util.ArrayList;
-import java.util.List;
-import java.util.Map;
-import java.util.logging.Level;
-import java.util.logging.Logger;
-import org.lcsim.detector.tracker.silicon.SiSensor;
-import org.lcsim.event.RawTrackerHit;
-import org.lcsim.event.Track;
-import org.lcsim.event.TrackerHit;
-import org.lcsim.fit.helicaltrack.*;
-import org.lcsim.hps.event.HPSTransformations;
-import org.lcsim.hps.monitoring.AIDAFrame;
-import org.lcsim.hps.recon.tracking.SvtUtils;
-import org.lcsim.hps.recon.tracking.TrackUtils;
-import org.lcsim.hps.recon.tracking.TrackerHitUtils;
-import org.lcsim.hps.users.mgraham.alignment.RunAlignment;
-import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
-import org.lcsim.recon.tracking.seedtracker.SeedTrack;
-import org.lcsim.util.aida.AIDA;
-
-/**
- * Class to calculate and print the residuals and derivatives
- * of the alignment parameters...used as input for MillePede
- * Notation follows the MillePede manual:
- * http://www.desy.de/~blobel/Mptwo.pdf
- *
- * the track is measured in the HelicalTrackFit frame
- * and residuals are in the sensor frame (u,v,w)
- *
- * ordering of track parameters is
- * double d0 = _trk.dca();
- * double z0 = _trk.z0();
- * double slope = _trk.slope();
- * double phi0 = _trk.phi0();
- * double R = _trk.R();
- *
- * @author mgraham
- */
-public class MPAlignmentParameters {
-
- 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];
- private double[] _error = new double[3];
- FileWriter fWriter;
- PrintWriter pWriter;
- private String _type = "LOCAL";
- boolean _DEBUG=false;
- private Hep3Vector _bfield;
- private ResLimit _resLimits;
- private HPSTransformations _detToTrk;
- private AlignmentUtils _alignUtils;
- private AlignmentUtils.OldAlignmentUtils _oldAlignUtils;
- private AlignmentUtils.NumDerivatives _numDerivatives;
- private TrackUtils trackUtil;
- private TrackerHitUtils trackerHitUtil;
-
- private boolean _includeMS = true;
-
- private AIDAFrame plotterFrame;
- private AIDAFrame plotterFrameSummary;
- private boolean hideFrame = false;
- private AIDA aida = AIDA.defaultInstance();
- private IAnalysisFactory af = aida.analysisFactory();
- IDataPointSet dps_t;
- IDataPointSet dps_b;
- IDataPointSet dps_pull_t;
- IDataPointSet dps_pull_b;
- IPlotter plotter_resuydiff_t;
- IPlotter plotter_resuydiff_b;
-
-
-
-
-
-
-
- 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.);
- try {
- fWriter = new FileWriter(outfile);
- pWriter = new PrintWriter(fWriter);
- } catch (IOException ex) {
- Logger.getLogger(RunAlignment.class.getName()).log(Level.SEVERE, null, ex);
- }
-
- makeAlignmentPlots();
-
- }
-
- public void setResLimits(int l,int d, double low,double high) {
- _resLimits.add(0,l,d, low,high); //top
- _resLimits.add(1,l,d, low,high); //bottom
- }
-
- public void setResLimits(int s, int l,int d, double low,double high) {
- _resLimits.add(s,l,d, low,high);
- }
-
- public void setHideFrame(boolean hide) {
- this.hideFrame = hide;
- }
-
- public void setDebug(boolean debug) {
- this._DEBUG = debug;
- _alignUtils.setDebug(debug);
- _numDerivatives.setDebug(debug);
- }
-
- public void setUniformZFieldStrength(double bfield) {
- _bfield = new BasicHep3Vector(0., 0., bfield);
- }
-
- public void setType(String t) {
- this._type = t;
- }
-
- public void setIncludeMS(boolean include) {
- this._includeMS = include;
- }
-
-
-
- public void PrintResidualsAndDerivatives(Track track, int itrack) {
-
- SeedTrack st = (SeedTrack) track;
- SeedCandidate seed = st.getSeedCandidate();
- Map<HelicalTrackHit, MultipleScatter> msmap = seed.getMSMap();
- _trk = seed.getHelix();
- List<TrackerHit> hitsOnTrack = track.getTrackerHits();
- String half = hitsOnTrack.get(0).getPosition()[2]>0 ? "top" : "bottom";
- pWriter.printf("TRACK %s (%d)\n",half,itrack);
- aida.cloud1D("Track Chi2 "+ half).fill(track.getChi2());
- aida.cloud1D("Track Chi2ndf "+ half).fill(track.getChi2()/track.getNDF());
- if(_DEBUG) System.out.printf("%s: track %d (chi2=%.2f ndf=%d) has %d hits\n",this.getClass().getSimpleName(),itrack,track.getChi2(),track.getNDF(),hitsOnTrack.size());
- for (TrackerHit hit : hitsOnTrack) {
-
- HelicalTrackHit htc = (HelicalTrackHit) hit;
-
- if(!(htc instanceof HelicalTrackCross)) {
- 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();
-
-
- // 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);
- }
-
- for (HelicalTrackStrip cl : clusterlist) {
- if(!"GLOBAL".equals(_type)) {
- CalculateResidual(cl, msdrphi, msdz);
- CalculateLocalDerivatives(cl);
- this.CalculateGlobalDerivativesSimple(cl);
- //CalculateGlobalDerivatives(cl);
- PrintStripResiduals(cl);
- } else {
- //CalculateLocalDerivativesGLOBAL(cl);
- CalculateGlobalDerivatives(cl);
- CalculateResidualGLOBAL(cl, msdrphi, msdz);
- }
- }
- }
-
- if(itrack%50==0) this.updatePlots();
- }
-
-
- private void CalculateLocalDerivatives(HelicalTrackStrip strip) {
- double xint = strip.origin().x();
- 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);
-
- if (_DEBUG) {
-
- //get track parameters.
- double d0 = _trk.dca();
- double z0 = _trk.z0();
- double slope = _trk.slope();
- double phi0 = _trk.phi0();
- double R = _trk.R();
- 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());
- System.out.printf("%s: Strip Origin %s \n",this.getClass().getSimpleName(),strip.origin().toString());
- 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",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",dfdqGlobalOld.toString());
- }
-
- }
-
-
-
-
- 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) {
-
- //Find interception with plane that the strips belongs to
- Hep3Vector x_vec = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()));
-
- if(Double.isNaN(x_vec.x())) {
- System.out.printf("%s: error this trkpos is wrong %s\n",this.getClass().getSimpleName(),x_vec.toString());
- System.out.printf("%s: origin %s trk \n%s\n",this.getClass().getSimpleName(),strip.origin(),_trk.toString());
- x_vec = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()),true);
- System.exit(1);
- }
- double slope = _trk.slope();
- double xr = 0.0;
- double yr = 0.0;
- double d0 = _trk.dca();
- double phi0 = _trk.phi0();
- double R = _trk.R();
- double z0 = _trk.z0();
- double s = HelixUtils.PathToXPlane(_trk, x_vec.x(), 0, 0).get(0);
- double phi = -s/R + phi0;
- double umeas = strip.umeas();
- Hep3Vector corigin = strip.origin();
- double vmeas = corigin.y(); //THis is wrong
- int layer = strip.layer();
- int side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? 10000 : 20000;
-
- if(_DEBUG) {
- System.out.printf("%s: --- Calculate global derivatives ---\n",this.getClass().getSimpleName());
- System.out.printf("%s: Side %d, layer %d, strip origin %s\n",this.getClass().getSimpleName(),side,layer,corigin.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(), d0, z0, slope, phi0, R,x_vec.x(),phi,x_vec.x(),s);
- }
-
- //1st index = alignment parameter (only u so far)
- //2nd index = residual coordinate (on du so far)
- //Naming scheme:
- //[Top]:
- // 10000 = top
- // 20000 = bottom
- //[Type]:
- // 1000 - translation
- // 2000 - rotation
- //[Direction] (tracker coord. frame)
- // 100 - x (beamline direction)
- // 200 - y (bend plane)
- // 300 - z (non-bend direction)
- // [Layer]
- // 1-10
-
-
- //Rotation matrix from the tracking fram to the sensor/strip frame
- Hep3Matrix T = trackerHitUtil.getTrackToStripRotation(strip);
-
-
-
- /*
- * Calculate jacobian da/db
- */
-
- BasicMatrix da_db = this._alignUtils.calculateJacobian(x_vec,T);
-
-
- if(_DEBUG) {
- this._alignUtils.printJacobianInfo(da_db);
- }
-
-
-
- /*
- * Invert the Jacobian
- */
-
-
- BasicMatrix db_da = (BasicMatrix) MatrixOp.inverse(da_db);
-
-
-
- if(_DEBUG) {
- System.out.printf("%s: Invert the Jacobian. We get da_db^-1=db_da: \n%s \n",this.getClass().getSimpleName(),db_da.toString());
- BasicMatrix prod = (BasicMatrix) MatrixOp.mult(db_da,da_db);
- System.out.printf("%s: Check the inversion i.e. db_da*da_db: \n%s \n",this.getClass().getSimpleName(),prod.toString());
- }
-
-
- /*
- * Calculate global derivative of the residual dz/db = dq_a^gl/db-dq_p^gl/db
- * dq_a^gl is the alignment corrected position in the global tracking frame
- * dq_p^gl_db is the predicted hit position (from the track model) in the global/tracking frame
- *
- */
-
-
-
- //****************************************************************************
- // First term in dz/db
-
- //3x6 matrix
- BasicMatrix dq_agl_db = _alignUtils.calculateGlobalHitPositionDers(x_vec);
-
-
- if(_DEBUG) {
- _alignUtils.printGlobalHitPositionDers(dq_agl_db);
- }
-
-
-
- /*
- * Second term in dz/db
- * dq_p^gl_db is the predicted hit position (from the track model) in the global/tracking frame
- *
- */
-
- //3x6 matrix
- BasicMatrix dq_pgl_db = _alignUtils.calculateGlobalPredictedHitPositionDers(_trk,x_vec);
-
- if(_DEBUG) {
- _alignUtils.printGlobalPredictedHitPositionDers(dq_pgl_db);
-
- System.out.printf("%s: Cross-check using numerical derivatives for dq_pgl/db (numder)\n",this.getClass().getSimpleName());
-
- BasicMatrix dq_pgl_db_numder = this._numDerivatives.calculateGlobalPredictedHitPositionDers(_trk,x_vec);
- _alignUtils.printGlobalPredictedHitPositionDers(dq_pgl_db_numder);
-
- }
-
- /*
- * Note that a shift in the position of the hit (q_agl) in in the
- * y/z plane is equivalent of the shift of the position of the prediction hit q_pgl
- * in the same plane. Thus we set these terms of the derivative to zero
- *
- */
-
- dq_pgl_db.setElement(0, 0, 0);
- dq_pgl_db.setElement(1, 1, 0);
- dq_pgl_db.setElement(2, 2, 0);
-
-
- /*
- * For the rotations in this global frame the rotation leads to a shift of the position
- * of the predicted hit. Since angles are small this rotation is the same as the
- * rotation of the hit position
- */
-
- dq_pgl_db.setElement(0, 3, 0);
- dq_pgl_db.setElement(1, 3, 0);
- dq_pgl_db.setElement(2, 3, 0);
- dq_pgl_db.setElement(0, 4, 0);
- dq_pgl_db.setElement(1, 4, 0);
- dq_pgl_db.setElement(2, 4, 0);
- dq_pgl_db.setElement(0, 5, 0);
- dq_pgl_db.setElement(1, 5, 0);
- dq_pgl_db.setElement(2, 5, 0);
-
-
- if(_DEBUG) {
- System.out.printf("%s: Fixing the predicted derivatives that are equivalent to the hit position derivatives. The result is below.\n",this.getClass().getSimpleName());
- _alignUtils.printGlobalPredictedHitPositionDers(dq_pgl_db);
- }
-
- /*
- * Putting them together i.e. subtract the predicted from the hit position derivative
- */
-
- //3x6 matrix
- BasicMatrix dz_db = (BasicMatrix) MatrixOp.sub(dq_agl_db, dq_pgl_db);
-
-
- if(_DEBUG) {
- _alignUtils.printGlobalResidualDers(dz_db);
- }
-
- /*
- * Convert the to the local frame using the Jacobian
- */
-
- //3x6 = 3x6*6x6 matrix
- BasicMatrix dz_da = (BasicMatrix) MatrixOp.mult(dz_db, db_da);
-
- if(_DEBUG) {
- _alignUtils.printLocalResidualDers(dz_da);
- }
-
-
-
- if(_DEBUG) {
- System.out.printf("%s: PELLE check manual one entry\n",this.getClass().getSimpleName());
- double dx_dx = dz_db.e(0, 0);
- double dx_dy = dz_db.e(0, 1);
- double dx_dz = dz_db.e(0, 2);
- double dx_da = dz_db.e(0, 3);
- double dx_db = dz_db.e(0, 4);
- double dx_dc = dz_db.e(0, 5);
- double dx_du = db_da.e(0, 0);
- double dy_du = db_da.e(1, 0);
- double dz_du = db_da.e(2, 0);
- double da_du = db_da.e(3, 0);
- double db_du = db_da.e(4, 0);
- double dc_du = db_da.e(5, 0);
- System.out.printf("%s: dx_dx*dx_du = %.3f ",this.getClass().getSimpleName(),dx_dx*dx_du);
- System.out.printf(" dx_dy*dy_du = %.3f ",dx_dy*dy_du);
- System.out.printf(" dx_dz*dz_du = %.3f (=%.3f*%.3f)\n",dx_dz*dz_du,dx_dz,dz_du);
- System.out.printf("%s: dx_da*da_du = %.3f ",this.getClass().getSimpleName(),dx_da*da_du);
- System.out.printf(" dx_db*db_du = %.3f ",dx_db*db_du);
- System.out.printf(" dx_dc*dc_du = %.3f \n",dx_dc*dc_du);
-
- double du_du = dx_dx*dx_du + dx_dy*dy_du + dx_dz*dz_du + dx_da*da_du + dx_db*db_du + dx_dc*dc_du;
- System.out.printf("%s: du_du = %.3f comapred to %.3f \n",this.getClass().getSimpleName(),du_du,dz_da.e(0, 0));
- }
-
-
-
-
- if(1==1) return;
-
-
- //Flag to tell if this hit is affected by the given global parameter
- boolean useGL = false;
-
- //Clear the old parameter list
- _glp.clear();
-
-
-
-
-
-
-
- BasicMatrix dqp_da_TRACK = new BasicMatrix(3,1);
-
-
- //Put it into a matrix to be able to transform easily
- //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1);
- //Get transformation matrix from tracking frame to sensor frame where the residuals are calculated
- Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip);
- if(_DEBUG) {
- System.out.println("Final transformation from tracking frame to strip frame:\n" + trkToStrip.toString());
- }
-
- //Transform to sensor frame!
- BasicMatrix dqp_da = (BasicMatrix) MatrixOp.mult(trkToStrip, dqp_da_TRACK);
- //Add it to the global parameter object
- GlobalParameter gp_tx = new GlobalParameter("Translation in x",side,layer,1000,100,true);
- gp_tx.setDfDp(dqp_da);
- _glp.add(gp_tx);
- if (_DEBUG) {
- gp_tx.print();
- System.out.println("Track frame dqp_da: " + dqp_da_TRACK);
- //System.out.printf("dqp_da = %5.5f %5.5f %5.5f GL%d name: %s\n", gp.dqp_da(0), gp.dqp_da(1), gp.dqp_da(2), gp.getLabel(),gp.getName());
- }
-
-
- //****************************************************************************
- //Derivatives of the predicted hit position qp for a translation of in y
- double y = -(R-d0)*Math.cos(phi0) + _alignUtils.sign(R) *Math.sqrt(Math.pow(R, 2)-Math.pow(x_vec.x()-(R-d0)*Math.sin(phi0), 2));
- dqp_da_TRACK.setElement(0, 0, _alignUtils.dx_dy(y, d0, phi0, R, phi));
- dqp_da_TRACK.setElement(1, 0, _alignUtils.dy_dy());
- dqp_da_TRACK.setElement(2, 0, _alignUtils.dz_dy(y, d0, phi0, slope, R, phi));
-
-
- //Put it into a matrix to be able to transform easily
- //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1);
- //Transform derivatives to sensor frame!
- dqp_da = (BasicMatrix) MatrixOp.mult(trkToStrip, dqp_da_TRACK);
- //Add it to the global parameter object
- GlobalParameter gp_ty = new GlobalParameter("Translation in y",side,layer,1000,200,true);
- gp_ty.setDfDp(dqp_da);
- _glp.add(gp_ty);
- if (_DEBUG) {
- gp_ty.print();
- System.out.println("Track frame dqp_da: " + dqp_da_TRACK);
- //System.out.printf("dfdp = %5.5f %5.5f %5.5f GL%d name: %s\n", gp.dfdp(0), gp.dfdp(1), gp.dfdp(2), gp.getLabel(),gp.getName());
- }
-
- //****************************************************************************
- //Derivatives of the predicted hit position qp for a translation of in z
-
- dqp_da_TRACK.setElement(0, 0, _alignUtils.dx_dz(slope, R, phi));
- dqp_da_TRACK.setElement(1, 0, _alignUtils.dy_dz(slope, R, phi));
- dqp_da_TRACK.setElement(2, 0, _alignUtils.dz_dz());
-
-
- //Put it into a matrix to be able to transform easily
- //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1);
- //Transform derivatives to sensor frame!
- dqp_da = (BasicMatrix) MatrixOp.mult(trkToStrip, dqp_da_TRACK);
- //Add it to the global parameter object
- GlobalParameter gp_tz = new GlobalParameter("Translation in z",side,layer,1000,300,true);
- gp_tz.setDfDp(dqp_da);
- _glp.add(gp_tz);
- if (_DEBUG) {
- gp_tz.print();
- System.out.println("Track frame dqp_da: " + dqp_da_TRACK);
- //System.out.printf("dfdp = %5.5f %5.5f %5.5f GL%d name: %s\n", gp.dfdp(0), gp.dfdp(1), gp.dfdp(2), gp.getLabel(),gp.getName());
- }
-
-
-
-
-
-
-
-
-
- }
-
-
-
-
-
-
- private void CalculateResidual(HelicalTrackStrip strip, double msdrdphi, double msdz) {
- if(_DEBUG) System.out.printf("%s: --- CalculateResidual ---\n",this.getClass().getSimpleName());
-
- Hep3Vector u = strip.u();
- Hep3Vector v = strip.v();
- Hep3Vector w = strip.w();
- Hep3Vector eta = VecOp.cross(u,v);
- Hep3Vector corigin = strip.origin();
- String side = SvtUtils.getInstance().isTopLayer((SiSensor)((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom";
-
- double bfield = Math.abs(this._bfield.z());
- 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, 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());
- }
- //Find interception with plane that the strips belongs to
- Hep3Vector trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, bfield);
-
- if(_DEBUG) {
- System.out.printf("%s: found interception point at %s \n",this.getClass().getSimpleName(),trkpos.toString());
- }
-
-
- if(Double.isNaN(trkpos.x()) || Double.isNaN(trkpos.y()) || Double.isNaN(trkpos.z())) {
- System.out.printf("%s: failed to get interception point (%s) \n",this.getClass().getSimpleName(),trkpos.toString());
- System.out.printf("%s: track params\n%s\n",this.getClass().getSimpleName(),_trk.toString());
- System.out.printf("%s: track pT=%.3f chi2=[%.3f][%.3f] \n",this.getClass().getSimpleName(),_trk.pT(bfield),_trk.chisq()[0],_trk.chisq()[1]);
- trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, bfield,true);
- System.exit(1);
- }
-
-
-
- double xint = trkpos.x();
- double phi0 = _trk.phi0();
- double R = _trk.R();
- double s = HelixUtils.PathToXPlane(_trk, xint, 0, 0).get(0);
- double phi = -s/R + phi0;
-
-
-
- Hep3Vector mserr = new BasicHep3Vector(msdrdphi * Math.sin(phi), msdrdphi * Math.sin(phi), msdz);
- Hep3Vector vdiffTrk = VecOp.sub(trkpos, corigin);
- Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip);
- Hep3Vector vdiff = VecOp.mult(trkToStrip, vdiffTrk);
-
-
-
- double umc = vdiff.x();
- double vmc = vdiff.y();
- double wmc = vdiff.z();
- double umeas = strip.umeas();
- double uError = strip.du();
- double msuError = VecOp.dot(mserr, u);
- double vmeas = 0;
- double vError = (strip.vmax() - strip.vmin()) / Math.sqrt(12);
- double wmeas = 0;
- double wError = 10.0/Math.sqrt(12); //0.001;
- _resid[0] = umeas - umc;
- _error[0] = this._includeMS ? Math.sqrt(uError * uError + msuError * msuError) : uError;
- _resid[1] = vmeas - vmc;
- _error[1] = vError;
- _resid[2] = wmeas - wmc;
- _error[2] = wError;
-
-
-
- //Fill residuals vs distrance from center of plane in the v directions
- aida.profile1D("res_u_vs_ydiff_layer_" + strip.layer() + "_" + side).fill(vdiffTrk.y(),_resid[0]);
-
- if (_DEBUG) {
- System.out.printf("%s: CalculateResidual Result ----\n",this.getClass().getSimpleName());
- 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());
- }
-
- }
-
-
-
-
-
-
-
-
-
- private void CalculateResidualGLOBAL(HelicalTrackStrip strip, double msdrdphi, double msdz) {
-
- /*
- * Calculate the residual in the GLOBAL tracking frame
- *
- */
-
- if(_DEBUG) System.out.println("--- CalculateResidualGLOBAL ---");
- if(_DEBUG) System.out.println("Strip origin = " + strip.origin().toString());
- Hep3Vector q_global = this.trackerHitUtil.getClusterPosition(strip,false);
- Hep3Vector q_global_error = trackerHitUtil.CalculateStripUncertaintyInGlobalFrame(strip,_trk, msdrdphi, msdz);
- double phi0 = _trk.phi0();
- double R = _trk.R();
- //double xint = strip.origin().x();
- double xint = trackUtil.calculateHelixInterceptXPlane(_trk, strip);
- 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);
- if(_DEBUG) System.out.println("posOnHelix = " + posOnHelix.toString());
- _resid[0] = q_global.x() - posOnHelix.x();
- _error[0] = q_global_error.x();
- _resid[1] = q_global.y() - posOnHelix.y();
- _error[1] = q_global_error.y();
- _resid[2] = q_global.z() - posOnHelix.z();
- _error[2] = q_global_error.z();
- if(_DEBUG) System.out.println("res x " + _resid[0] + " +- " + q_global_error.x());
- if(_DEBUG) System.out.println("res y " + _resid[1] + " +- " + q_global_error.y());
- if(_DEBUG) System.out.println("res z " + _resid[2] + " +- " + q_global_error.z());
-
- }
-
-
-
-
- private void PrintStripResiduals(HelicalTrackStrip strip) {
- String side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom";
- if (_DEBUG) {
- System.out.printf("%s: --- Alignment Results for this Strip ---\n",this.getClass().getSimpleName());
-
- String sensor_type = SvtUtils.getInstance().isAxial((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "axial" : "stereo";
- System.out.printf("%s: Strip layer %4d is %s in %s at origin %s\n",this.getClass().getSimpleName(), strip.layer(),sensor_type, side, strip.origin().toString());
- System.out.printf("%s: u=%s v=%s w=%s\n",this.getClass().getSimpleName(), strip.u().toString(),strip.v().toString(),strip.w().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]);
- System.out.printf("%s: Errors (u,v,w) : %5.5e %5.5e %5.5e\n",this.getClass().getSimpleName(), _error[0], _error[1], _error[2]);
- String[] q = {"d0", "z0", "slope", "phi0", "R"};
- System.out.printf("%s: track parameter derivatives:\n",this.getClass().getSimpleName());
- for (int i = 0; i < _nTrackParameters; i++) {
- System.out.printf("%s: %s %5.5e %5.5e %5.5e\n",this.getClass().getSimpleName(), q[i], _dfdq.e(0, i), _dfdq.e(1, i), _dfdq.e(2, i));
- }
- //String[] p = {"u-displacement"};
- System.out.printf("%s: Global parameter derivatives\n",this.getClass().getSimpleName());
- for (GlobalParameter gl : _glp) {
- 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: --- 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
- //String[] d = {"u","v","w"};
- String[] d = {"u"}; //phansson use only u direction!
-
- int s;
- for(int j=0;j<1;++j){
- side = "bottom";
- s = 1;
- if( strip.origin().z()>0) {
- side = "top";
- s = 0;
- }
-
-
- //if(!isAllowedResidual(s,strip.layer(),j,_resid[j])) {
- // if(_DEBUG) System.out.println("Layer " + strip.layer() + " in " + d[j] + " res " + _resid[j] + " +- " + _error[j] + " was outside allowed range");
- // continue;
- //}
-//// if(Math.abs(_resid[j]/_error[j])>3) {
-//// if(_DEBUG) System.out.println("Layer " + strip.layer() + " in " + d[j] + " res " + _resid[j] + " +- " + _error[j] + " had too high pull");
-//// continue;
-//// }
- pWriter.printf("res_%s %5.5e %5.5e \n", d[j], _resid[j], _error[j]);
- for (int i = 0; i < _nTrackParameters; i++) {
- pWriter.printf("lc_%s %5.5e \n", d[j], _dfdq.e(j, i));
- }
- for (GlobalParameter gl: _glp) {
- if(gl.active()){
- pWriter.printf("gl_%s %5.5e %5d\n", d[j], gl.dfdp(j), gl.getLabel());
- //x-check that side is correct
- int lbl = gl.getLabel();
- Double df = Math.floor(lbl/10000);
- 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("WARNING track is on " + s + " while gl param is " + lbl + "(" + df + ")" );
- System.exit(1);
- }
-
- }
- }
- if( _resid[j] < 9999999 ) {
- if (_DEBUG) System.out.println(this.getClass().getSimpleName() + ": filling ures with " + _resid[j]);
- aida.histogram1D("res_"+d[j]+"_layer" + strip.layer() + "_" + side).fill(_resid[j]);
- aida.histogram1D("reserr_"+d[j]+"_layer" + strip.layer() + "_" + side).fill(_error[j]);
- aida.histogram1D("respull_"+d[j]+"_layer" + strip.layer() + "_" + side).fill(_resid[j]/_error[j]);
-
-
- aida.histogram2D("respull_slope_"+d[j]+"_layer" + strip.layer() + "_" + side).fill(_trk.slope(),_resid[j]/_error[j]);
- }
-
- }
- }
-
-
-
-
-// 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 boolean isAllowedResidual(int side,int layer,int d,double res) {
-
- if(res <_resLimits.getMin(side,layer,d)) return false;
- if(res >_resLimits.getMax(side,layer,d)) return false;
- return true;
-
- }
-
-
- private void makeAlignmentPlots() {
-
-
-
- aida.tree().cd("/");
- plotterFrame = new AIDAFrame();
- plotterFrame.setTitle("Residuals");
- plotterFrameSummary = new AIDAFrame();
- plotterFrameSummary.setTitle("Summary");
-
- String[] side = {"top","bottom"};
-
- IPlotter plotter_chi2 = af.createPlotterFactory().create();
- plotter_chi2.createRegions(2,2);
- plotter_chi2.setTitle("Track Chi2");
- plotterFrame.addPlotter(plotter_chi2);
[truncated at 1000 lines; 236 more skipped]