 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
| StripMPAlignmentInput.java | +1114 | added 1.1 | |
| MPAlignmentInputCalculator.java | +137 | added 1.1 | |
| ModuleMPAlignmentInput.java | +709 | added 1.1 | |
| RunMPAlignment.java | +9 | -10 | 1.13 -> 1.14 |
| +1969 | -10 | ||
Split up into super and derived classes.
diff -N StripMPAlignmentInput.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ StripMPAlignmentInput.java 19 Nov 2012 21:51:07 -0000 1.1 @@ -0,0 +1,1114 @@
+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.Arrays;
+import java.util.List;
+import java.util.Map;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.constants.Constants;
+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 StripMPAlignmentInput extends MPAlignmentInputCalculator {
+
+ private double[] _resid = new double[3];
+ private double[] _error = new double[3];
+ private final int _nTrackParameters = 5; //the five track parameters
+
+
+ private AIDAFrame plotterFrame;
+ private AIDAFrame plotterFrameSummary;
+ private IDataPointSet dps_t;
+ private IDataPointSet dps_b;
+ private IDataPointSet dps_pull_t;
+ private IDataPointSet dps_pull_b;
+ private IPlotter plotter_resuydiff_t;
+ private IPlotter plotter_resuydiff_b;
+
+
+
+
+
+
+
+ public StripMPAlignmentInput(String outfile,String type) {
+ super(outfile,type);
+
+ makeAlignmentPlots();
+
+ }
+
+ @Override
+ public void setResLimits(int l,int d, double low,double high) {
+
+ }
+
+ public void setResLimits(int s, int l,int d, double low,double high) {
+
+ }
+
+
+
+ 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";
+ this.addMilleInputLine(String.format("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) {
+ CalculateResidual(cl, msdrphi, msdz);
+ CalculateLocalDerivatives(cl);
+ CalculateGlobalDerivatives(cl);
+ PrintStripResiduals(cl);
+
+ }
+
+ }
+
+ 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 CalculateGlobalDerivatives(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 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] = _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 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));
+ 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: --- END Alignment Results for this Strip ---\n",this.getClass().getSimpleName());
+ }
+
+
+
+ this.addMilleInputLine(String.format("%d\n", strip.layer()));
+
+
+ List<String> direction = new ArrayList<String>(Arrays.asList("u")); //use only u direction!
+
+ int s;
+ for(int d=0;d<direction.size();++d){
+ side = "bottom";
+ s = 1;
+ if( strip.origin().z()>0) {
+ side = "top";
+ s = 0;
+ }
+
+
+ //if(!isAllowedResidual(s,strip.layer(),j,_residirection.get(d))) {
+ // if(_DEBUG) System.out.println("Layer " + strip.layer() + " in " + direction.get(d) + " res " + _residirection.get(d) + " +- " + _error[d] + " was outside allowed range");
+ // continue;
+ //}
+
+ this.addMilleInputLine(String.format("res_%s %5.5e %5.5e \n", direction.get(d), _resid[d], _error[d]));
+ for (int i = 0; i < _nTrackParameters; i++) {
+ this.addMilleInputLine(String.format("lc_%s %5.5e \n", direction.get(d), _dfdq.e(d, i)));
+ }
+ for (GlobalParameter gl: _glp) {
+ if(gl.active()){
+ this.addMilleInputLine(String.format("gl_%s %5.5e %5d\n", direction.get(d), gl.dfdp(d), 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 (_DEBUG) System.out.println(this.getClass().getSimpleName() + ": filling ures with " + _resid[d]);
+ aida.histogram1D("res_"+direction.get(d)+"_layer" + strip.layer() + "_" + side).fill(_resid[d]);
+ aida.histogram1D("reserr_"+direction.get(d)+"_layer" + strip.layer() + "_" + side).fill(_error[d]);
+ aida.histogram1D("respull_"+direction.get(d)+"_layer" + strip.layer() + "_" + side).fill(_resid[d]/_error[d]);
+ aida.histogram2D("respull_slope_"+direction.get(d)+"_layer" + strip.layer() + "_" + side).fill(_trk.slope(),_resid[d]/_error[d]);
+
+ }
+ }
+
+
+
+
+
+
+ @Override
+ protected void makeAlignmentPlots() {
+
+
+
+ aida.tree().cd("/");
+ plotterFrame = new AIDAFrame();
+ plotterFrame.setTitle("Residuals");
+ plotterFrameSummary = new AIDAFrame();
+ plotterFrameSummary.setTitle("Summary");
+
+ List<String> sides = new ArrayList<String>(Arrays.asList("top","bottom"));
+
+ IPlotter plotter_chi2 = af.createPlotterFactory().create();
+ plotter_chi2.createRegions(2,2);
+ plotter_chi2.setTitle("Track Chi2");
+ plotterFrame.addPlotter(plotter_chi2);
+ ICloud1D hchi2_top = aida.cloud1D("Track Chi2 top");
+ ICloud1D hchi2_bot = aida.cloud1D("Track Chi2 bottom");
+ ICloud1D hchi2ndf_top = aida.cloud1D("Track Chi2ndf top");
+ ICloud1D hchi2ndf_bot = aida.cloud1D("Track Chi2ndf bottom");
+ plotter_chi2.region(0).plot(hchi2_top);
+ plotter_chi2.region(1).plot(hchi2_bot);
+ plotter_chi2.region(2).plot(hchi2ndf_top);
+ plotter_chi2.region(3).plot(hchi2ndf_bot);
+
+
+ List<String> direction = "MODULE".equals(this._type) ? new ArrayList<String>(Arrays.asList("x","y","z")) : new ArrayList<String>(Arrays.asList("u"));
+
+
+
+ double xbins_u_res[][] = {
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01}
+ {-0.15,0.15},
+ {-0.4,0.4},
+ {-0.7,0.7},
+ {-0.5,0.5},
+ {-0.5,0.5},
+ {-0.5,0.5},
+ {-2,2},
+ {-2,2},
+ {-2,2},
+ {-2,2}
+ };
+ double xbins_u_reserr[][] = {
+ {0,0.04},
+ {0,0.04},
+ {0,0.5},
+ {0,0.5},
+ {0,0.7},
+ {0,0.7},
+ {0,1.2},
+ {0,1.2},
+ {0,1.5},
+ {0,1.5}
+ };
+ double xbins_u_respull[][] = {
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5}
+ };
+
+ for (int d=0;d<direction.size();++d) {
+
+ for (int s=0;s<2;++s) {
+ // if(iSide==1) continue;
+ IPlotter plotter_res = af.createPlotterFactory().create();
+ plotter_res.createRegions(5,2,0);
+ plotter_res.setTitle("res_" + direction.get(d) + " " + sides.get(s));
+ IPlotter plotter_reserr = af.createPlotterFactory().create();
+ plotter_reserr.createRegions(5,2,0);
+ plotter_reserr.setTitle("reserr_" + direction.get(d) + " " + sides.get(s));
+ IPlotter plotter_respull = af.createPlotterFactory().create();
+ plotter_respull.createRegions(5,2,0);
+ plotter_respull.setTitle("respull_" + direction.get(d) + " " + sides.get(s));
+ IPlotter plotter_respull_slope = af.createPlotterFactory().create();
+ plotter_respull_slope.createRegions(5,2,0);
+ plotter_respull_slope.setTitle("respull_slope_" + direction.get(d) + " " + sides.get(s));
+
+ for (int iLayer=1;iLayer<11;++iLayer) {
+
+ IHistogram h = aida.histogram1D("res_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) , 50, xbins_u_res[iLayer-1][0], xbins_u_res[iLayer-1][1]);
+ IHistogram h1 = aida.histogram1D("reserr_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) , 50, xbins_u_reserr[iLayer-1][0], xbins_u_reserr[iLayer-1][1]);
+ IHistogram h2 = aida.histogram1D("respull_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) , 50, xbins_u_respull[iLayer-1][0], xbins_u_respull[iLayer-1][1]);
+
+ IHistogram hslope;
+ if(sides.get(s)=="top") hslope = aida.histogram2D("respull_slope_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) ,50, 0 , 0.1, 50, xbins_u_respull[iLayer-1][0], xbins_u_respull[iLayer-1][1]);
+ else hslope = aida.histogram2D("respull_slope_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) ,50, -0.1 , 0, 50, xbins_u_respull[iLayer-1][0], xbins_u_respull[iLayer-1][1]);
+
+
+ int region = (iLayer-1);//*2+iSide;
+
+ plotter_res.region(region).plot(h);
+ plotter_reserr.region(region).plot(h1);
+ plotter_respull.region(region).plot(h2);
+ plotter_respull_slope.region(region).plot(hslope);
+
+ ((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);
+ ((PlotterRegion) plotter_respull_slope.region(region)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_respull_slope.region(region)).getPlot().setAllowPopupMenus(true);
+ plotter_respull_slope.style().dataStyle().fillStyle().setParameter("colorMapScheme","rainbow");
+ ((PlotterRegion) plotter_respull_slope.region(region)).getPlot().setAllowPopupMenus(true);
+
+
+ }
+
+ //plotter_res.show();
+ plotterFrame.addPlotter(plotter_res);
+ plotterFrame.addPlotter(plotter_reserr);
+ plotterFrame.addPlotter(plotter_respull);
+ plotterFrame.addPlotter(plotter_respull_slope);
+ }
+
+
+
+
+ }
+
+
+
+
+ IPlotter plotter_prf = af.createPlotterFactory().create();
+ plotter_prf.createRegions(1,2,0);
+ plotter_prf.setTitle("<Residuals>");
+ IDataPointSetFactory dpsf = aida.analysisFactory().createDataPointSetFactory(null);
+ dps_t = dpsf.create("dps_t", "Mean of u residuals top",2);
+ plotter_prf.region(0).plot(dps_t);
+ dps_b = dpsf.create("dps_b", "Mean of u residuals bottom",2);
+ plotter_prf.region(0).plot(dps_b);
+
+ dps_pull_t = dpsf.create("dps_pull_t", "Mean of u pulls top",2);
+ plotter_prf.region(1).plot(dps_pull_t);
+ dps_pull_b = dpsf.create("dps_pull_b", "Mean of u pulls bottom",2);
+ plotter_prf.region(1).plot(dps_pull_b);
+
+ for(int region=0;region<2;++region) {
+ ((PlotterRegion) plotter_prf.region(region)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_prf.region(region)).getPlot().setAllowPopupMenus(true);
+ }
+
+ plotterFrameSummary.addPlotter(plotter_prf);
+
+
+ plotter_resuydiff_t = af.createPlotterFactory().create();
+ plotter_resuydiff_t.createRegions(5,2,0);
+ plotter_resuydiff_t.setTitle("res u vs ydiff top");
+ plotter_resuydiff_b = af.createPlotterFactory().create();
+ plotter_resuydiff_b.createRegions(5,2,0);
+ plotter_resuydiff_b.setTitle("res u vs ydiff bot");
+ //plotter_resuydiff_t_b.createRegions(5,2,0);
+ //plotter_resuydiff_b.setTitle("distr: res u vs ydiff");
+ for(int iLayer=1;iLayer<11;++iLayer) {
+ IProfile1D prf_t = aida.profile1D("res_u_vs_ydiff_layer_"+iLayer+"_top",10,-50,50);
+ IProfile1D prf_b = aida.profile1D("res_u_vs_ydiff_layer_"+iLayer+"_bottom",10,-50,50);
+ plotter_resuydiff_t.region(iLayer-1).plot(prf_t);
+ plotter_resuydiff_b.region(iLayer-1).plot(prf_b);
+ ((PlotterRegion) plotter_resuydiff_t.region(iLayer-1)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_resuydiff_t.region(iLayer-1)).getPlot().setAllowPopupMenus(true);
+ ((PlotterRegion) plotter_resuydiff_b.region(iLayer-1)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_resuydiff_b.region(iLayer-1)).getPlot().setAllowPopupMenus(true);
+ }
+
+ plotterFrameSummary.addPlotter(plotter_resuydiff_t);
+ plotterFrameSummary.addPlotter(plotter_resuydiff_b);
+
+
+ plotterFrame.pack();
+ plotterFrame.setVisible(!hideFrame);
+ plotterFrameSummary.pack();
+ plotterFrameSummary.setVisible(!hideFrame);
+
+
+ }
+
+
+ @Override
+ public void updatePlots() {
+ dps_t.clear();
+ dps_b.clear();
+ dps_pull_t.clear();
+ dps_pull_b.clear();
+
+ final int nLayers = 10;
+ final List<String> direction = this._type.equals("MODULE") ? new ArrayList<String>(Arrays.asList("x","y","z")) : new ArrayList<String>(Arrays.asList("u"));
+ for(int i=1;i<nLayers+1;++i) {
+
+ for(String d : direction) {
+
+ double mean = aida.histogram1D("res_"+d+"_layer"+i+"_top").mean();
+ double stddev = aida.histogram1D("res_"+d+"_layer"+i+"_top").rms();
+ double N = aida.histogram1D("res_"+d+"_layer"+i+"_top").entries();
+ double error = N >0 ? stddev/Math.sqrt(N) : 0;
+ dps_t.addPoint();
+ dps_t.point(i-1).coordinate(1).setValue(mean);
+ dps_t.point(i-1).coordinate(1).setErrorPlus(error);
+ dps_t.point(i-1).coordinate(0).setValue(i);
+ dps_t.point(i-1).coordinate(0).setErrorPlus(0);
+
+ mean = aida.histogram1D("res_"+d+"_layer"+i+"_bottom").mean();
+ stddev = aida.histogram1D("res_"+d+"_layer"+i+"_bottom").rms();
+ N = aida.histogram1D("res_"+d+"_layer"+i+"_bottom").entries();
+ error = N >0 ? stddev/Math.sqrt(N) : 0;
+ dps_b.addPoint();
+ dps_b.point(i-1).coordinate(1).setValue(mean);
+ dps_b.point(i-1).coordinate(1).setErrorPlus(error);
+ dps_b.point(i-1).coordinate(0).setValue(i);
+ dps_b.point(i-1).coordinate(0).setErrorPlus(0);
+
+ mean = aida.histogram1D("respull_"+d+"_layer"+i+"_top").mean();
+ stddev = aida.histogram1D("respull_"+d+"_layer"+i+"_top").rms();
+ N = aida.histogram1D("respull_"+d+"_layer"+i+"_top").entries();
+ error = N >0 ? stddev/Math.sqrt(N) : 0;
+ dps_pull_t.addPoint();
+ dps_pull_t.point(i-1).coordinate(1).setValue(mean);
+ dps_pull_t.point(i-1).coordinate(1).setErrorPlus(error);
+ dps_pull_t.point(i-1).coordinate(0).setValue(i);
+ dps_pull_t.point(i-1).coordinate(0).setErrorPlus(0);
+
+ mean = aida.histogram1D("respull_"+d+"_layer"+i+"_bottom").mean();
+ stddev = aida.histogram1D("respull_"+d+"_layer"+i+"_bottom").rms();
+ N = aida.histogram1D("respull_"+d+"_layer"+i+"_bottom").entries();
+ error = N >0 ? stddev/Math.sqrt(N) : 0;
+ dps_pull_b.addPoint();
+ dps_pull_b.point(i-1).coordinate(1).setValue(mean);
+ dps_pull_b.point(i-1).coordinate(1).setErrorPlus(error);
+ dps_pull_b.point(i-1).coordinate(0).setValue(i);
+ dps_pull_b.point(i-1).coordinate(0).setErrorPlus(0);
+
+ }
+ }
+
+
+
+ }
+
+
+
+
+
+ private void CalculateGlobalDerivativesWRONG(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);
[truncated at 1000 lines; 118 more skipped]
diff -N MPAlignmentInputCalculator.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MPAlignmentInputCalculator.java 19 Nov 2012 21:51:07 -0000 1.1 @@ -0,0 +1,137 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.hps.users.phansson;
+
+import hep.aida.IAnalysisFactory;
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.event.Track;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.hps.event.HPSTransformations;
+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.util.aida.AIDA;
+
+/**
+ *
+ * @author phansson
+ */
+public abstract class MPAlignmentInputCalculator {
+ protected final int _nTrackParameters = 5; //the five track parameters
+ protected boolean _DEBUG = false;
+ private String _outfile;
+ protected List<GlobalParameter> _glp;
+ protected BasicMatrix _dfdq;
+ protected HelicalTrackFit _trk;
+ protected String _type = "LOCAL";
+ protected Hep3Vector _bfield;
+ protected ResLimit _resLimits;
+ protected HPSTransformations _detToTrk;
+ protected AlignmentUtils _alignUtils;
+ protected AlignmentUtils.OldAlignmentUtils _oldAlignUtils;
+ protected AlignmentUtils.NumDerivatives _numDerivatives;
+ protected TrackUtils trackUtil;
+ protected TrackerHitUtils trackerHitUtil;
+ protected boolean hideFrame = false;
+ protected AIDA aida = AIDA.defaultInstance();
+ protected IAnalysisFactory af = aida.analysisFactory();
+
+ protected boolean _includeMS = true;
+ public abstract void PrintResidualsAndDerivatives(Track track, int itrack);
+ protected abstract void makeAlignmentPlots();
+ public abstract void updatePlots();
+ private FileWriter fWriter;
+ private PrintWriter pWriter;
+
+ public MPAlignmentInputCalculator(String outfile,String type) {
+ _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.);
+ _type = type;
+ _outfile = outfile;
+ openFile();
+ }
+
+
+
+ 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 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;
+
+ }
+
+
+ public void closeFile() {
+ try {
+ pWriter.close();
+ fWriter.close();
+ } catch(IOException ex) {
+ Logger.getLogger(RunAlignment.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ private void openFile() {
+ try {
+ fWriter = new FileWriter(_outfile);
+ pWriter = new PrintWriter(fWriter);
+ } catch (IOException ex) {
+ Logger.getLogger(RunAlignment.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void addMilleInputLine(String line) {
+ this.pWriter.print(line);
+ }
+
+ public void setDebug(boolean debug) {
+ this._DEBUG = debug;
+ _alignUtils.setDebug(debug);
+ _numDerivatives.setDebug(debug);
+ }
+
+ public void setHideFrame(boolean hide) {
+ this.hideFrame = hide;
+ }
+
+ public void setIncludeMS(boolean include) {
+ this._includeMS = include;
+ }
+
+
+ public void setUniformZFieldStrength(double bfield) {
+ _bfield = new BasicHep3Vector(0.0, 0.0, bfield);
+ }
+
+
+
+}
diff -N ModuleMPAlignmentInput.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ ModuleMPAlignmentInput.java 19 Nov 2012 21:51:07 -0000 1.1 @@ -0,0 +1,709 @@
+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.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+import java.util.Map;
+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.monitoring.AIDAFrame;
+import org.lcsim.hps.recon.tracking.SvtUtils;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+
+/**
+ * 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 tracking frame (x,y,z)
+ *
+ * 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 final class ModuleMPAlignmentInput extends MPAlignmentInputCalculator {
+
+ private double[] _resid = new double[3];
+ private double[] _error = new double[3];
+
+
+ private AIDAFrame plotterFrame;
+ private AIDAFrame plotterFrameSummary;
+ private IDataPointSet dps_t;
+ private IDataPointSet dps_b;
+ private IDataPointSet dps_pull_t;
+ private IDataPointSet dps_pull_b;
+ private IPlotter plotter_resuydiff_t;
+ private IPlotter plotter_resuydiff_b;
+
+
+
+
+
+
+
+ public ModuleMPAlignmentInput(String outfile,String type) {
+ super(outfile,type);
+
+ makeAlignmentPlots();
+
+ }
+
+
+
+ @Override
+ 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";
+ this.addMilleInputLine(String.format("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);
+
+
+ CalculateResidual(htc);
+ CalculateLocalDerivatives(htc);
+ CalculateGlobalDerivatives(htc);
+ PrintHitResiduals(htc);
+
+
+ }
+
+ if(itrack%50==0) this.updatePlots();
+ }
+
+
+
+
+ private void CalculateLocalDerivatives(HelicalTrackHit hit) {
+ double xint = hit.x();
+ BasicMatrix dfdqGlobalOld = _oldAlignUtils.calculateLocalHelixDerivatives(_trk, xint);
+ BasicMatrix dfdqGlobal = _alignUtils.calculateLocalHelixDerivatives(_trk, xint);
+ BasicMatrix dfdqGlobalNum = this._numDerivatives.calculateLocalHelixDerivatives(_trk, xint);
+ _dfdq = 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 HelicalTrackHit Result --- \n",this.getClass().getSimpleName());
+ System.out.printf("%s: Hit position: %s \n",this.getClass().getSimpleName(),hit.getCorrectedPosition().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 CalculateGlobalDerivatives(HelicalTrackHit hit) {
+
+ /*
+ * Residual in tracking 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 tracking frame
+ *
+ * Calcualte the derivative of dr/da where
+ * a=(delta_u,delta_v,delta_w,alpha,beta,gamma)
+ * are the alingment paramerers in the tracking frame
+ *
+ * Factorize: dr/da=dr/dm_a*dm_a/da
+ *
+ * Start with dr/dma
+ */
+
+ //Find the unit vector of the track direction
+ TrackDirection trkdir = HelixUtils.CalculateTrackDirection(_trk, _trk.PathMap().get(hit));
+ Hep3Vector t = trkdir.Direction();
+ //Find the unit vector of the plane where the hit happened
+ if(!(hit instanceof HelicalTrackCross) ) {
+ throw new RuntimeException("Error: gl ders for HTH only works if they are crosses!");
+ }
+
+
+ HelicalTrackCross cross = (HelicalTrackCross) hit;
+ HelicalTrackStrip strip = cross.getStrips().get(0);
+ Hep3Vector n = strip.w();
+
+
+
+ /*
+ * use the measured hit position
+ */
+ double umeas,vmeas,wmeas;
+ umeas = hit.x();
+ vmeas = hit.y(); //the un-measured (along the strip) is set to zero
+ wmeas = hit.z(); // the hit is on the surface of the plane
+
+ /*
+ * 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 = hit.Layer();
+ int side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? 10000 : 20000;
+
+
+ if(_DEBUG) {
+ double phi = -_trk.PathMap().get(hit)/_trk.R()+_trk.phi0();
+ System.out.printf("%s: --- Calculate SIMPLE global derivatives for helical track hit ---\n",this.getClass().getSimpleName());
+ System.out.printf("%s: Side %d, layer %d, hit position %s\n",this.getClass().getSimpleName(),side,layer,hit.getCorrectedPosition().toString());
+ System.out.printf("%s: %10s%10s%10s%10s%10s%10s%10s%10s\n",this.getClass().getSimpleName(),"d0","z0","slope","phi0","R","phi","xint","s");
+ System.out.printf("%s: %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(),phi,hit.x(),_trk.PathMap().get(hit));
+ 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 x",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 y",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 z",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 CalculateResidual(HelicalTrackHit hit) {
+ if(_DEBUG) System.out.printf("%s: --- Calculate Residual for HelicalTrackhit ---\n",this.getClass().getSimpleName());
+
+ Map<String,Double> res_map = this.trackUtil.calculateTrackHitResidual(hit, _trk, _includeMS);
+ _resid[0] = 0.;
+ _error[0] = 0.320/Math.sqrt(12);
+ _resid[1] = res_map.get("resy");
+ _error[1] = res_map.get("erry");
+ _resid[2] = res_map.get("resz");
+ _error[2] = res_map.get("errz");
+
+
+
+ //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 HelicalTrackHit Result ----\n",this.getClass().getSimpleName());
+ System.out.printf("%s: Hit positon: %s\n",this.getClass().getSimpleName(),hit.getCorrectedPosition().toString());
+ Hep3Vector trkpos = HelixUtils.PointOnHelix(_trk, HelixUtils.PathToXPlane(_trk, hit.x(), 0, 0).get(0));
+ System.out.printf("%s: Predicted position: %s\n",this.getClass().getSimpleName(),trkpos.toString());
+ System.out.printf("%s: drphi=%,4f dz=%.4f MS: drdphi=%.4f, dz=%.4f\n",this.getClass().getSimpleName(),hit.drphi(),hit.dr()*Math.abs(_trk.slope()),_trk.ScatterMap().get(hit).drphi(),_trk.ScatterMap().get(hit).dz());
+ System.out.printf("%s: %12s %12s %12s\n","resx","resy","resz");
+ System.out.printf("%s: %5.5f+-%5.5f %5.5f+-%5.5f %5.5f+-%5.5f\n",_resid[0],_error[0],_resid[1],_error[1],_resid[2],_error[2]);
+ }
+
+ }
+
+
+
+
+
+
+
+
+ private void PrintHitResiduals(HelicalTrackHit hit) {
+ HelicalTrackStrip strip = ((HelicalTrackCross) hit).getStrips().get(0);
+ String side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom";
+ if (_DEBUG) {
+ System.out.printf("%s: --- Alignment Results for this helical track hit ---\n",this.getClass().getSimpleName());
+
+ String sensor_type = SvtUtils.getInstance().isAxial((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "axial" : "stereo";
+ System.out.printf("%s: Layer %d in %s at position %s\n",this.getClass().getSimpleName(), hit.Layer(),side, hit.getCorrectedPosition().toString());
+ System.out.printf("%s: Residuals (x,y,z) : %5.5e %5.5e %5.5e\n",this.getClass().getSimpleName(), _resid[0], _resid[1], _resid[2]);
+ System.out.printf("%s: Errors (x,y,z) : %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: --- END Alignment Results for this helical track hit ---\n",this.getClass().getSimpleName());
+ }
+ this.addMilleInputLine(String.format("%d\n", hit.Layer()));
+ // Loop over the three directions u,v,w and print residuals and derivatives
+
+ List<String> direction = new ArrayList<String>(Arrays.asList("x","y","z"));
+
+ int s = "bottom".equals(side) ? 1 : 0;
+ for(int j=0;j<direction.size();++j){
+
+
+
+ //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;
+//// }
+ //This should really be x,y,z but the downstream parser likes u,v,w -> FIX THIS!
+ String d = direction.get(j).equals("x") ? "u" : direction.get(j).equals("y") ? "v" : "w";
+ this.addMilleInputLine(String.format("res_%s %5.5e %5.5e \n", d, _resid[j], _error[j]));
+ for (int i = 0; i < _nTrackParameters; i++) {
+ this.addMilleInputLine(String.format("lc_%s %5.5e \n", d, _dfdq.e(j, i)));
+ }
+ for (GlobalParameter gl: _glp) {
+ if(gl.active()){
+ this.addMilleInputLine(String.format("gl_%s %5.5e %5d\n", d, gl.dfdp(j), gl.getLabel()));
+
+ //Cross 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_"+direction.get(j)+"_layer" + strip.layer() + "_" + side).fill(_resid[j]);
+ aida.histogram1D("reserr_"+direction.get(j)+"_layer" + strip.layer() + "_" + side).fill(_error[j]);
+ aida.histogram1D("respull_"+direction.get(j)+"_layer" + strip.layer() + "_" + side).fill(_resid[j]/_error[j]);
+
+
+ aida.histogram2D("respull_slope_"+direction.get(j)+"_layer" + strip.layer() + "_" + side).fill(_trk.slope(),_resid[j]/_error[j]);
+ }
+
+ }
+ }
+
+
+
+
+
+ @Override
+ public void makeAlignmentPlots() {
+
+
+
+ aida.tree().cd("/");
+ plotterFrame = new AIDAFrame();
+ plotterFrame.setTitle("Residuals");
+ plotterFrameSummary = new AIDAFrame();
+ plotterFrameSummary.setTitle("Summary");
+
+ List<String> sides = new ArrayList<String>(Arrays.asList("top","bottom"));
+
+ IPlotter plotter_chi2 = af.createPlotterFactory().create();
+ plotter_chi2.createRegions(2,2);
+ plotter_chi2.setTitle("Track Chi2");
+ plotterFrame.addPlotter(plotter_chi2);
+ ICloud1D hchi2_top = aida.cloud1D("Track Chi2 top");
+ ICloud1D hchi2_bot = aida.cloud1D("Track Chi2 bottom");
+ ICloud1D hchi2ndf_top = aida.cloud1D("Track Chi2ndf top");
+ ICloud1D hchi2ndf_bot = aida.cloud1D("Track Chi2ndf bottom");
+ plotter_chi2.region(0).plot(hchi2_top);
+ plotter_chi2.region(1).plot(hchi2_bot);
+ plotter_chi2.region(2).plot(hchi2ndf_top);
+ plotter_chi2.region(3).plot(hchi2ndf_bot);
+
+
+ List<String> direction = new ArrayList<String>(Arrays.asList("x","y","z"));
+
+
+
+ double xbins_u_res[][] = {
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01},
+// {-0.01,0.01}
+ {-0.15,0.15},
+ {-0.4,0.4},
+ {-0.7,0.7},
+ {-0.5,0.5},
+ {-0.5,0.5},
+ {-0.5,0.5},
+ {-2,2},
+ {-2,2},
+ {-2,2},
+ {-2,2}
+ };
+ double xbins_u_reserr[][] = {
+ {0,0.04},
+ {0,0.04},
+ {0,0.5},
+ {0,0.5},
+ {0,0.7},
+ {0,0.7},
+ {0,1.2},
+ {0,1.2},
+ {0,1.5},
+ {0,1.5}
+ };
+ double xbins_u_respull[][] = {
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5},
+ {-5,5}
+ };
+
+ for (int d=0;d<direction.size();++d) {
+
+ for (int s=0;s<2;++s) {
+ // if(iSide==1) continue;
+ IPlotter plotter_res = af.createPlotterFactory().create();
+ plotter_res.createRegions(5,2,0);
+ plotter_res.setTitle("res_" + direction.get(d) + " " + sides.get(s));
+ IPlotter plotter_reserr = af.createPlotterFactory().create();
+ plotter_reserr.createRegions(5,2,0);
+ plotter_reserr.setTitle("reserr_" + direction.get(d) + " " + sides.get(s));
+ IPlotter plotter_respull = af.createPlotterFactory().create();
+ plotter_respull.createRegions(5,2,0);
+ plotter_respull.setTitle("respull_" + direction.get(d) + " " + sides.get(s));
+ IPlotter plotter_respull_slope = af.createPlotterFactory().create();
+ plotter_respull_slope.createRegions(5,2,0);
+ plotter_respull_slope.setTitle("respull_slope_" + direction.get(d) + " " + sides.get(s));
+
+ for (int iLayer=1;iLayer<11;++iLayer) {
+
+ IHistogram h = aida.histogram1D("res_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) , 50, xbins_u_res[iLayer-1][0], xbins_u_res[iLayer-1][1]);
+ IHistogram h1 = aida.histogram1D("reserr_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) , 50, xbins_u_reserr[iLayer-1][0], xbins_u_reserr[iLayer-1][1]);
+ IHistogram h2 = aida.histogram1D("respull_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) , 50, xbins_u_respull[iLayer-1][0], xbins_u_respull[iLayer-1][1]);
+
+ IHistogram hslope;
+ if(sides.get(s)=="top") hslope = aida.histogram2D("respull_slope_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) ,50, 0 , 0.1, 50, xbins_u_respull[iLayer-1][0], xbins_u_respull[iLayer-1][1]);
+ else hslope = aida.histogram2D("respull_slope_"+direction.get(d)+"_layer" + (iLayer) + "_" + sides.get(s) ,50, -0.1 , 0, 50, xbins_u_respull[iLayer-1][0], xbins_u_respull[iLayer-1][1]);
+
+
+ int region = (iLayer-1);//*2+iSide;
+
+ plotter_res.region(region).plot(h);
+ plotter_reserr.region(region).plot(h1);
+ plotter_respull.region(region).plot(h2);
+ plotter_respull_slope.region(region).plot(hslope);
+
+ ((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);
+ ((PlotterRegion) plotter_respull_slope.region(region)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_respull_slope.region(region)).getPlot().setAllowPopupMenus(true);
+ plotter_respull_slope.style().dataStyle().fillStyle().setParameter("colorMapScheme","rainbow");
+ ((PlotterRegion) plotter_respull_slope.region(region)).getPlot().setAllowPopupMenus(true);
+
+
+ }
+
+ //plotter_res.show();
+ plotterFrame.addPlotter(plotter_res);
+ plotterFrame.addPlotter(plotter_reserr);
+ plotterFrame.addPlotter(plotter_respull);
+ plotterFrame.addPlotter(plotter_respull_slope);
+ }
+
+
+
+
+ }
+
+
+
+
+ IPlotter plotter_prf = af.createPlotterFactory().create();
+ plotter_prf.createRegions(1,2,0);
+ plotter_prf.setTitle("<Residuals>");
+ IDataPointSetFactory dpsf = aida.analysisFactory().createDataPointSetFactory(null);
+ dps_t = dpsf.create("dps_t", "Mean of u residuals top",2);
+ plotter_prf.region(0).plot(dps_t);
+ dps_b = dpsf.create("dps_b", "Mean of u residuals bottom",2);
+ plotter_prf.region(0).plot(dps_b);
+
+ dps_pull_t = dpsf.create("dps_pull_t", "Mean of u pulls top",2);
+ plotter_prf.region(1).plot(dps_pull_t);
+ dps_pull_b = dpsf.create("dps_pull_b", "Mean of u pulls bottom",2);
+ plotter_prf.region(1).plot(dps_pull_b);
+
+ for(int region=0;region<2;++region) {
+ ((PlotterRegion) plotter_prf.region(region)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_prf.region(region)).getPlot().setAllowPopupMenus(true);
+ }
+
+ plotterFrameSummary.addPlotter(plotter_prf);
+
+
+ plotter_resuydiff_t = af.createPlotterFactory().create();
+ plotter_resuydiff_t.createRegions(5,2,0);
+ plotter_resuydiff_t.setTitle("res u vs ydiff top");
+ plotter_resuydiff_b = af.createPlotterFactory().create();
+ plotter_resuydiff_b.createRegions(5,2,0);
+ plotter_resuydiff_b.setTitle("res u vs ydiff bot");
+ //plotter_resuydiff_t_b.createRegions(5,2,0);
+ //plotter_resuydiff_b.setTitle("distr: res u vs ydiff");
+ for(int iLayer=1;iLayer<11;++iLayer) {
+ IProfile1D prf_t = aida.profile1D("res_u_vs_ydiff_layer_"+iLayer+"_top",10,-50,50);
+ IProfile1D prf_b = aida.profile1D("res_u_vs_ydiff_layer_"+iLayer+"_bottom",10,-50,50);
+ plotter_resuydiff_t.region(iLayer-1).plot(prf_t);
+ plotter_resuydiff_b.region(iLayer-1).plot(prf_b);
+ ((PlotterRegion) plotter_resuydiff_t.region(iLayer-1)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_resuydiff_t.region(iLayer-1)).getPlot().setAllowPopupMenus(true);
+ ((PlotterRegion) plotter_resuydiff_b.region(iLayer-1)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_resuydiff_b.region(iLayer-1)).getPlot().setAllowPopupMenus(true);
+ }
+
+ plotterFrameSummary.addPlotter(plotter_resuydiff_t);
+ plotterFrameSummary.addPlotter(plotter_resuydiff_b);
+
+
+ plotterFrame.pack();
+ plotterFrame.setVisible(!hideFrame);
+ plotterFrameSummary.pack();
+ plotterFrameSummary.setVisible(!hideFrame);
+
+
+ }
+
+
+ @Override
+ public void updatePlots() {
+ dps_t.clear();
+ dps_b.clear();
+ dps_pull_t.clear();
+ dps_pull_b.clear();
+
+ final int nLayers = 10;
+ final List<String> direction = new ArrayList<String>(Arrays.asList("x","y","z"));
+ for(int i=1;i<nLayers+1;++i) {
+
+ for(String d : direction) {
+
+ double mean = aida.histogram1D("res_"+d+"_layer"+i+"_top").mean();
+ double stddev = aida.histogram1D("res_"+d+"_layer"+i+"_top").rms();
+ double N = aida.histogram1D("res_"+d+"_layer"+i+"_top").entries();
+ double error = N >0 ? stddev/Math.sqrt(N) : 0;
+ dps_t.addPoint();
+ dps_t.point(i-1).coordinate(1).setValue(mean);
+ dps_t.point(i-1).coordinate(1).setErrorPlus(error);
+ dps_t.point(i-1).coordinate(0).setValue(i);
+ dps_t.point(i-1).coordinate(0).setErrorPlus(0);
+
+ mean = aida.histogram1D("res_"+d+"_layer"+i+"_bottom").mean();
+ stddev = aida.histogram1D("res_"+d+"_layer"+i+"_bottom").rms();
+ N = aida.histogram1D("res_"+d+"_layer"+i+"_bottom").entries();
+ error = N >0 ? stddev/Math.sqrt(N) : 0;
+ dps_b.addPoint();
+ dps_b.point(i-1).coordinate(1).setValue(mean);
+ dps_b.point(i-1).coordinate(1).setErrorPlus(error);
+ dps_b.point(i-1).coordinate(0).setValue(i);
+ dps_b.point(i-1).coordinate(0).setErrorPlus(0);
+
+ mean = aida.histogram1D("respull_"+d+"_layer"+i+"_top").mean();
+ stddev = aida.histogram1D("respull_"+d+"_layer"+i+"_top").rms();
+ N = aida.histogram1D("respull_"+d+"_layer"+i+"_top").entries();
+ error = N >0 ? stddev/Math.sqrt(N) : 0;
+ dps_pull_t.addPoint();
+ dps_pull_t.point(i-1).coordinate(1).setValue(mean);
+ dps_pull_t.point(i-1).coordinate(1).setErrorPlus(error);
+ dps_pull_t.point(i-1).coordinate(0).setValue(i);
+ dps_pull_t.point(i-1).coordinate(0).setErrorPlus(0);
+
+ mean = aida.histogram1D("respull_"+d+"_layer"+i+"_bottom").mean();
+ stddev = aida.histogram1D("respull_"+d+"_layer"+i+"_bottom").rms();
+ N = aida.histogram1D("respull_"+d+"_layer"+i+"_bottom").entries();
+ error = N >0 ? stddev/Math.sqrt(N) : 0;
+ dps_pull_b.addPoint();
+ dps_pull_b.point(i-1).coordinate(1).setValue(mean);
+ dps_pull_b.point(i-1).coordinate(1).setErrorPlus(error);
+ dps_pull_b.point(i-1).coordinate(0).setValue(i);
+ dps_pull_b.point(i-1).coordinate(0).setErrorPlus(0);
+
+ }
+ }
+
+
+
+ }
+
+
+
+
+
+}
diff -u -r1.13 -r1.14 --- RunMPAlignment.java 14 Nov 2012 22:08:23 -0000 1.13 +++ RunMPAlignment.java 19 Nov 2012 21:51:07 -0000 1.14 @@ -28,7 +28,7 @@
int nevt = 0;
double[] beamsize = {0.001, 0.02, 0.02};
String _config = "";
- MPAlignmentParameters ap;
+ StripMPAlignmentInput ap;
int totalTracks=0;
int totalTracksProcessed=0;
private String _resLimitFileName="";
@@ -49,6 +49,9 @@
public void setDebug(boolean v) {
this._debug = v;
}
+ public void setType(String type) {
+ this._type = type;
+ }
public void setMilleFile(String filename) {
milleFile = filename;
}
@@ -72,9 +75,8 @@
@Override
public void detectorChanged(Detector detector) {
- ap = new MPAlignmentParameters(milleFile);
+ ap = new StripMPAlignmentInput(milleFile,_type);
ap.setDebug(_debug);
- ap.setType(this._type);
ap.setHideFrame(hideFrame);
double bfield = detector.getFieldMap().getField(new BasicHep3Vector(0., 0., 1.)).y();
System.out.printf("%s: B-field in z %.3f\n",this.getClass().getSimpleName(),bfield);
@@ -130,13 +132,7 @@
@Override
public void endOfData() {
ap.updatePlots();
- try {
- System.out.println(this.getClass().getSimpleName() + ": Total Number of Tracks Found = "+totalTracks);
- System.out.println(this.getClass().getSimpleName() + ": Total Number of Tracks Processed = "+totalTracksProcessed);
- ap.closeFile();
- } catch (IOException ex) {
- Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE, null, ex);
- }
+ ap.closeFile();
if (!"".equals(outputPlotFileName)) {
try {
aida.saveAs(outputPlotFileName);
@@ -144,6 +140,9 @@
Logger.getLogger(TrigRateDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);
}
}
+ System.out.println(this.getClass().getSimpleName() + ": Total Number of Tracks Found = "+totalTracks); + System.out.println(this.getClass().getSimpleName() + ": Total Number of Tracks Processed = "+totalTracksProcessed); +
}
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