 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
| HelicalTrackHitResidualsDriver.java | -315 | 1.3 removed | |
| MPAlignmentParameters.java | -1232 | 1.18 removed | |
| WTrackUtils.java | -259 | 1.4 removed | |
| -1806 | |||
Remove old versions
diff -N HelicalTrackHitResidualsDriver.java --- HelicalTrackHitResidualsDriver.java 14 Nov 2012 00:51:25 -0000 1.3 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,315 +0,0 @@
-package org.lcsim.hps.users.phansson;
-
-import hep.aida.*;
-import hep.aida.ref.plotter.PlotterRegion;
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-import java.io.BufferedReader;
-import java.io.FileNotFoundException;
-import java.io.FileReader;
-import java.io.IOException;
-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.*;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelicalTrackHit;
-import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
-import org.lcsim.fit.helicaltrack.MultipleScatter;
-import org.lcsim.geometry.Detector;
-import org.lcsim.hps.monitoring.AIDAFrame;
-import org.lcsim.hps.recon.tracking.SvtUtils;
-import org.lcsim.hps.recon.tracking.TrackUtils;
-import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
-import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
-import org.lcsim.recon.tracking.seedtracker.SeedTrack;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-
-/**
- *
- * @author mgraham
- */
-public class HelicalTrackHitResidualsDriver extends Driver {
-
- private AIDA aida = AIDA.defaultInstance();
- int totalTracks=0;
- int totalTracksProcessed=0;
- private String outputPlotFileName="";
- private String trackCollectionName="MatchedTracks";
- private boolean hideFrame = false;
- private boolean _debug = false;
- private TrackUtils _trackUtils = new TrackUtils();
-
- private AIDAFrame plotterFrame;
-// ICloud1D[] _h_resz_track_top = new ICloud1D[5];
-// ICloud1D[] _h_resz_track_bottom = new ICloud1D[5];
-// ICloud1D[] _h_resy_track_top = new ICloud1D[5];
-// ICloud1D[] _h_resy_track_bottom = new ICloud1D[5];
- IHistogram1D[] _h_resz_track_top = new IHistogram1D[5];
- IHistogram1D[] _h_resz_track_bottom = new IHistogram1D[5];
- IHistogram1D[] _h_resy_track_top = new IHistogram1D[5];
- IHistogram1D[] _h_resy_track_bottom = new IHistogram1D[5];
- IDataPointSet dps_hth_y_b;
- IDataPointSet dps_hth_y_t;
- IDataPointSet dps_hth_z_b;
- IDataPointSet dps_hth_z_t;
- IPlotter _plotter_resz_top;
- IPlotter _plotter_resy_top;
- IPlotter _plotter_resz_bottom;
- IPlotter _plotter_resy_bottom;
- IPlotter _plotter_mean_res;
-
-
-
-
- public void setDebug(boolean v) {
- this._debug = v;
- }
-
- public void setTrackCollectionName(String trackCollectionName) {
- this.trackCollectionName = trackCollectionName;
- }
-
- public void setOutputPlotFileName(String filename) {
- outputPlotFileName = filename;
- }
-
- public void setHideFrame(boolean hide) {
- hideFrame = hide;
- }
-
- public HelicalTrackHitResidualsDriver() {
-
- }
-
-
-
- @Override
- public void detectorChanged(Detector detector) {
-
- Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
- Hep3Vector _bfield = new BasicHep3Vector(0,0,detector.getFieldMap().getField(IP).y());
-
- makePlots();
-
-
- }
-
-
-
- @Override
- public void process(EventHeader event) {
-
- List<Track> tracklist = new ArrayList<Track>();
- if(event.hasCollection(Track.class,this.trackCollectionName)) {
- tracklist = event.get(Track.class, this.trackCollectionName);
- if(_debug) {
- System.out.println(this.getClass().getSimpleName() + ": Number of Tracks = " + tracklist.size());
- }
- }
-
-
- for (Track track : tracklist) {
- SeedTrack st = (SeedTrack) track;
- SeedCandidate seed = st.getSeedCandidate();
- HelicalTrackFit trk = seed.getHelix();
- List<TrackerHit> hitsOnTrack = track.getTrackerHits();
- for(TrackerHit hit : hitsOnTrack) {
- HelicalTrackHit hth = (HelicalTrackHit) hit;
- Map<String,Double> res_track = _trackUtils.calculateTrackHitResidual(hth, trk, false);
- //Map<String,Double> res_trackgen = _trackUtils.calculateTrackHitResidual(hth, trkgen, false);
- boolean isTop = false;
- if(SvtUtils.getInstance().isTopLayer((SiSensor)((RawTrackerHit)hth.getRawHits().get(0)).getDetectorElement())) {
- isTop = true;
- }
- int layer = hth.Layer();
- if(_debug) System.out.println(this.getClass().getSimpleName() + ": residual for hit at " + hth.toString() + " and layer " + layer);
- if(layer%2==0) {
- System.out.println(this.getClass().getSimpleName() + ": HTH layer is not odd!" + layer);
- System.exit(1);
- }
- int layer_idx = (layer-1)/2;
- if(isTop) {
- this._h_resz_track_top[layer_idx].fill(res_track.get("resz"));
- this._h_resy_track_top[layer_idx].fill(res_track.get("resy"));
- } else {
- this._h_resz_track_bottom[layer_idx].fill(res_track.get("resz"));
- this._h_resy_track_bottom[layer_idx].fill(res_track.get("resy"));
- }
-
-// if(Math.abs(res_track.get("resy"))>0.02) {
-// System.out.println(this.getClass().getSimpleName() + ": this has large y res = " + res_track.get("resy"));
-// System.exit(1);
-// }
-
- }
-
- }
- totalTracks++;
- totalTracksProcessed++;
-
-
-
- if(totalTracks%50==0) this.updatePlots();
-
-
-
- }
-
- public void endOfData() {
- this.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);
-// } catch (IOException ex) {
-// Logger.getLogger(CmpGenToFittedTracksDriver.class.getName()).log(Level.SEVERE, null, ex);
-// }
- if (!"".equals(outputPlotFileName))
- try {
- aida.saveAs(outputPlotFileName);
- } catch (IOException ex) {
- Logger.getLogger(TrigRateDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);
- }
-
- }
-
-
- private void makePlots() {
-
- int nbins = 50;
- double bins_resz_min = -0.2;
- double bins_resz_max = 0.2;
- double bins_resy_min = -0.4;
- double bins_resy_max = 0.4;
-
- _plotter_resz_top = aida.analysisFactory().createPlotterFactory().create();
- _plotter_resz_top.setTitle("res z top");
- _plotter_resz_top.createRegions(5,1);
- _plotter_resy_top = aida.analysisFactory().createPlotterFactory().create();
- _plotter_resy_top.setTitle("res y top");
- _plotter_resy_top.createRegions(5,1);
- for(int i=1;i<6;++i) {
-// _h_resz_track_top[i-1] = aida.cloud1D("h_resz_track_top_layer"+i);
-// _h_resy_track_top[i-1] = aida.cloud1D("h_resy_track_top_layer"+i);
- _h_resz_track_top[i-1] = aida.histogram1D("h_resz_track_top_layer"+i,nbins,bins_resz_min,bins_resz_max);
- _h_resy_track_top[i-1] = aida.histogram1D("h_resy_track_top_layer"+i,nbins,bins_resy_min,bins_resy_max);
- _plotter_resz_top.region(i-1).plot(_h_resz_track_top[i-1]);
- _plotter_resy_top.region(i-1).plot(_h_resy_track_top[i-1]);
- }
-
- _plotter_resz_bottom = aida.analysisFactory().createPlotterFactory().create();
- _plotter_resz_bottom.setTitle("res z bottom");
- _plotter_resz_bottom.createRegions(5,1);
- _plotter_resy_bottom = aida.analysisFactory().createPlotterFactory().create();
- _plotter_resy_bottom.setTitle("res y bottom");
- _plotter_resy_bottom.createRegions(5,1);
- for(int i=1;i<6;++i) {
- _h_resz_track_bottom[i-1] = aida.histogram1D("h_resz_track_bottom_layer"+i,nbins,bins_resz_min,bins_resz_max);
- _h_resy_track_bottom[i-1] = aida.histogram1D("h_resy_track_bottom_layer"+i,nbins,bins_resy_min,bins_resy_max);
-// _h_resz_track_bottom[i-1] = aida.cloud1D("h_resz_track_bottom_layer"+i);
-// _h_resy_track_bottom[i-1] = aida.cloud1D("h_resy_track_bottom_layer"+i);
- _plotter_resz_bottom.region(i-1).plot(_h_resz_track_bottom[i-1]);
- _plotter_resy_bottom.region(i-1).plot(_h_resy_track_bottom[i-1]);
-
- }
-
- _plotter_mean_res = aida.analysisFactory().createPlotterFactory().create();
- _plotter_mean_res.setTitle("Mean res y");
- _plotter_mean_res.createRegions(2,2);
-
- IDataPointSetFactory dpsf = aida.analysisFactory().createDataPointSetFactory(null);
-
- dps_hth_y_b = dpsf.create("dps_hth_y_b", "Mean of y residual bottom",2);
- dps_hth_y_t = dpsf.create("dps_hth_y_t", "Mean of y residual top",2);
- _plotter_mean_res.region(1).plot(dps_hth_y_b);
- _plotter_mean_res.region(0).plot(dps_hth_y_t);
-
- dps_hth_z_b = dpsf.create("dps_hth_z_b", "Mean of z residual bottom",2);
- dps_hth_z_t = dpsf.create("dps_hth_z_t", "Mean of z residual top",2);
- _plotter_mean_res.region(3).plot(dps_hth_z_b);
- _plotter_mean_res.region(2).plot(dps_hth_z_t);
-
- ((PlotterRegion)_plotter_mean_res.region(0)).getPlot().setAllowUserInteraction(true);
- ((PlotterRegion) _plotter_mean_res.region(0)).getPlot().setAllowPopupMenus(true);
- ((PlotterRegion)_plotter_mean_res.region(1)).getPlot().setAllowUserInteraction(true);
- ((PlotterRegion) _plotter_mean_res.region(1)).getPlot().setAllowPopupMenus(true);
-
-
- plotterFrame = new AIDAFrame();
- plotterFrame.setTitle("HTH Residuals");
- plotterFrame.addPlotter(_plotter_resz_top);
- plotterFrame.addPlotter(_plotter_resz_bottom);
- plotterFrame.addPlotter(_plotter_resy_top);
- plotterFrame.addPlotter(_plotter_resy_bottom);
- plotterFrame.addPlotter(_plotter_mean_res);
- plotterFrame.pack();
- plotterFrame.setVisible(!hideFrame);
-
- }
-
-
- void updatePlots() {
- dps_hth_y_t.clear();
- dps_hth_z_t.clear();
- dps_hth_z_b.clear();
- dps_hth_z_t.clear();
-
- for(int i=1;i<6;++i) {
-
- double mean = this._h_resy_track_bottom[i-1].mean();
- double stddev = this._h_resy_track_bottom[i-1].rms();
- double N = this._h_resy_track_bottom[i-1].entries();
- double error = N >0 ? stddev/Math.sqrt(N) : 0;
- dps_hth_y_b.addPoint();
- dps_hth_y_b.point(i-1).coordinate(1).setValue(mean);
- dps_hth_y_b.point(i-1).coordinate(1).setErrorPlus(error);
- dps_hth_y_b.point(i-1).coordinate(0).setValue(i);
- dps_hth_y_b.point(i-1).coordinate(0).setErrorPlus(0);
-
- mean = this._h_resy_track_top[i-1].mean();
- stddev = this._h_resy_track_top[i-1].rms();
- N = this._h_resy_track_top[i-1].entries();
- error = N >0 ? stddev/Math.sqrt(N) : 0;
- dps_hth_y_t.addPoint();
- dps_hth_y_t.point(i-1).coordinate(1).setValue(mean);
- dps_hth_y_t.point(i-1).coordinate(1).setErrorPlus(error);
- dps_hth_y_t.point(i-1).coordinate(0).setValue(i);
- dps_hth_y_t.point(i-1).coordinate(0).setErrorPlus(0);
-
- mean = this._h_resz_track_top[i-1].mean();
- stddev = this._h_resz_track_top[i-1].rms();
- N = this._h_resz_track_top[i-1].entries();
- error = N >0 ? stddev/Math.sqrt(N) : 0;
- dps_hth_z_t.addPoint();
- dps_hth_z_t.point(i-1).coordinate(1).setValue(mean);
- dps_hth_z_t.point(i-1).coordinate(1).setErrorPlus(error);
- dps_hth_z_t.point(i-1).coordinate(0).setValue(i);
- dps_hth_z_t.point(i-1).coordinate(0).setErrorPlus(0);
-
- mean = this._h_resz_track_bottom[i-1].mean();
- stddev = this._h_resz_track_bottom[i-1].rms();
- N = this._h_resz_track_bottom[i-1].entries();
- error = N >0 ? stddev/Math.sqrt(N) : 0;
- dps_hth_z_b.addPoint();
- dps_hth_z_b.point(i-1).coordinate(1).setValue(mean);
- dps_hth_z_b.point(i-1).coordinate(1).setErrorPlus(error);
- dps_hth_z_b.point(i-1).coordinate(0).setValue(i);
- dps_hth_z_b.point(i-1).coordinate(0).setErrorPlus(0);
-
-
-
-
-
- }
-
-
-
-
- }
-
-
-}
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]
diff -N WTrackUtils.java --- WTrackUtils.java 15 Nov 2012 02:46:15 -0000 1.4 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,259 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-package org.lcsim.hps.users.phansson;
-
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.VecOp;
-import java.util.ArrayList;
-import java.util.List;
-
-/**
- *
- * @author phansson
- */
-public class WTrackUtils {
- private int _debug = 0;
- List<Double> _delta_point = new ArrayList<Double>();
-
- public WTrackUtils() {
- }
-
- public WTrackUtils(boolean debug) {
- _debug = debug?1:0;
- }
-
- public WTrackUtils(int debug) {
- _debug = debug;
- }
-
-
- public void setDebug(boolean debug) {
- _debug = debug?1:0;
- }
-
- public boolean getDebug() {
- return _debug>0?true:false;
- }
-
- public List<Double> getDeltas() {
- return _delta_point;
- }
-
- private Hep3Vector getMomentumOnHelix(WTrack track, double s) {
- double a = track.getBfieldConstant();
- Hep3Vector p0 = track.getP0();
- double rho = a / p0.magnitude();
- double px = p0.x()*Math.cos(rho*s) - p0.y()*Math.sin(rho*s);
- double py = p0.y()*Math.cos(rho*s) + p0.x()*Math.sin(rho*s);
- double pz = p0.z();
- return (new BasicHep3Vector(px,py,pz));
- }
-
- private Hep3Vector getPointOnHelix(WTrack track, double s) {
- double a = track.getBfieldConstant();
- Hep3Vector p0 = track.getP0();
- Hep3Vector x0 = track.getX0();
- double rho = a / p0.magnitude();
- double x = x0.x() + p0.x()/a*Math.sin(rho*s) - p0.y()/a*(1-Math.cos(rho*s));
- double y = x0.y() + p0.y()/a*Math.sin(rho*s) + p0.x()/a*(1-Math.cos(rho*s));
- double z = x0.z() + p0.z()/p0.magnitude()*s;
- return (new BasicHep3Vector(x,y,z));
- }
-
- public double getPathLengthToPlaneApprox(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
- /*
- * Find the approximate path length to the point xp
- * in arbitrary oriented, constant magnetic field with unit vector h
- */
- double a = track.getBfieldConstant();
- Hep3Vector p0 = track.getP0();
- Hep3Vector x0 = track.getX0();
- double p = p0.magnitude();
- double rho = a / p;
- double A = VecOp.dot(eta,VecOp.cross(p0, h))/p*0.5*rho;
- double B = VecOp.dot(p0,eta)/p;
- double C = VecOp.dot(VecOp.sub(x0,xp),eta);
- double t = B*B-4*A*C;
- if(t<0) {
- System.out.println(this.getClass().getSimpleName() + ": getPathLengthToPlaneApprox ERROR t is negative (" + t + ")!" );
- System.out.println(this.getClass().getSimpleName() + ": p " + p + " rho " + rho + " A " + A + " B " + B + " C " + C);
- System.out.println(this.getClass().getSimpleName() + ": track params\n" + track.toString());
- System.out.println(this.getClass().getSimpleName() + ": xp " + xp.toString());
- System.out.println(this.getClass().getSimpleName() + ": eta " + eta.toString());
- System.out.println(this.getClass().getSimpleName() + ": h " + h.toString());
- System.exit(1);
- }
- double root1 = (-B + Math.sqrt(t)) /(2*A);
- double root2 = (-B - Math.sqrt(t)) /(2*A);
-
- // choose the smallest positive solution
- // if both negative choose the smallest negative ???
- //if(root1==0 || root2==0) root=0;
- double root = Math.abs(root1) <= Math.abs(root2) ? root1 : root2;
-// else if(Math.signum(root1)>0 && Math.signum(root2)<0) root = root1;
-// else if(Math.signum(root2)>0 && Math.signum(root1)<0) root = root2;
-// else if(Math.signum(root1)>0 && Math.signum(root2)>0) root = root1 > root2 ? root2 : root1;
-// else if(Math.signum(root1)<0 && Math.signum(root2)<0) root = root1 < root2 ? root2 : root1;
-// else {
-// System.out.println(this.getClass().getSimpleName() + ": I should never get here! (root1=" + root1 + " root2=" + root2+")");
-// System.exit(1);
-// }
- if((_debug>1)) {
- System.out.println(this.getClass().getSimpleName() + ": getPathLengthToPlaneApprox ");
- System.out.println(this.getClass().getSimpleName() + ": " + track.toString());
- System.out.println(this.getClass().getSimpleName() + ": xp " + xp.toString());
- System.out.println(this.getClass().getSimpleName() + ": eta " + eta.toString());
- System.out.println(this.getClass().getSimpleName() + ": h " + h.toString());
- System.out.println(this.getClass().getSimpleName() + ": p " + p + " rho " + rho + " t " + t + " A " + A + " B " + B + " C " + C);
- System.out.println(this.getClass().getSimpleName() + ": root1 " + root1 + " root2 " + root2 + " -> root " + root);
- }
- return root;
-
- }
-
-
- private Hep3Vector getPointOnHelix(WTrack track, double s, Hep3Vector h) {
- /*
- * Get point on helix at path lenght s
- * in arbitrary oriented, constant magnetic field with unit vector h
- */
-
- double a = track.getBfieldConstant();
- Hep3Vector p0 = track.getP0();
- double p = p0.magnitude();
- Hep3Vector x0 = track.getX0();
- double rho = a / p0.magnitude();
- double srho = s*rho;
- Hep3Vector a1 = VecOp.mult(1/a*Math.sin(srho), p0);
- Hep3Vector a2 = VecOp.mult(1/a*(1-Math.cos(srho)),VecOp.cross(p0,h));
- Hep3Vector a3 = VecOp.mult(VecOp.dot(p0, h)/p,h);
- Hep3Vector a4 = VecOp.mult(s-Math.sin(srho)/rho, a3);
- Hep3Vector x = VecOp.add(x0,a1);
- x = VecOp.sub(x,a2);
- x = VecOp.add(x,a4);
- return x;
- }
-
- private Hep3Vector getMomentumOnHelix(WTrack track, double s, Hep3Vector h) {
- /*
- * Get point on helix at path lenght s
- * in arbitrary oriented, constant magnetic field with unit vector h
- */
-
- double a = track.getBfieldConstant();
- Hep3Vector p0 = track.getP0();
- double rho = a / p0.magnitude();
- double srho = s*rho;
- Hep3Vector a1 = VecOp.mult(Math.cos(srho), p0);
- Hep3Vector a2 = VecOp.cross(p0, VecOp.mult(Math.sin(srho),h));
- Hep3Vector a3 = VecOp.mult(VecOp.dot(p0,h),VecOp.mult(1-Math.cos(srho),h));
- Hep3Vector p = VecOp.sub(a1,a2);
- p = VecOp.add(p,a3);
- return p;
- }
-
-
- public double[] getHelixParametersAtPathLength(WTrack track, double s, Hep3Vector h) {
- /*
- * Calculate the exact position of the new helix parameters at path length s
- * in an arbitrarily oriented, constant magnetic field
- *
- * point xp is the point
- * h is a unit vector in the direction of the magnetic field
- */
-
- // Find track parameters at that path length
- Hep3Vector p = this.getMomentumOnHelix(track, s, h);
- Hep3Vector x = this.getPointOnHelix(track, s, h);
-
- Hep3Vector p_tmp = this.getMomentumOnHelix(track, s);
- Hep3Vector x_tmp = this.getPointOnHelix(track, s);
-
- if((_debug>1)) {
- System.out.println(this.getClass().getSimpleName() + ": point on helix at s");
- System.out.println(this.getClass().getSimpleName() + ": p " + p.toString() + " p_tmp " + p_tmp.toString());
- System.out.println(this.getClass().getSimpleName() + ": x " + x.toString() + " x_tmp " + x_tmp.toString());
- }
-
-
- //Create the new parameter array
- double [] pars = new double[7];
- pars[0] = p.x();
- pars[1] = p.y();
- pars[2] = p.z();
- pars[3] = track.getParameters()[3]; //E is unchanged
- pars[4] = x.x();
- pars[5] = x.y();
- pars[6] = x.z();
- return pars;
- }
-
-
-
- public Hep3Vector getHelixAndPlaneIntercept(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
-
- /*
- * Find the interception point between the helix and plane
- * track: the helix in W representation
- * xp: point on the plane
- * eta: unit vector of the plane
- * h: unit vector of magnetic field
- */
-
- _delta_point.clear();
-
- int max_it = 10;
- double epsilon = 1e-4;
- int iteration = 1;
- double delta_s = 9999;
- double s_prev = 999999;
- List<WTrack> tracks = new ArrayList<WTrack>();
- WTrack trk = new WTrack(track);
- while(iteration<=max_it && Math.abs(s_prev)>epsilon) {
- if((_debug>1)) {
- System.out.println(this.getClass().getSimpleName() + ": iteration " + iteration + " --- ");
- System.out.println(this.getClass().getSimpleName() + ": s_prev " + s_prev + " delta_s " + delta_s);
- System.out.println(this.getClass().getSimpleName() + ": " + trk.toString());
- }
- // Start by approximating the path length to the point
- double s = this.getPathLengthToPlaneApprox(trk, xp, eta, h);
-
- if((_debug>1)) {
- System.out.println(this.getClass().getSimpleName() + ": Approx. path s " + s);
- }
-
- // Find the track parameters at this point
- double[] params = this.getHelixParametersAtPathLength(trk, s, h);
- // update the track parameters
- trk.setTrackParameters(params);
-
- if((_debug>1)) {
- System.out.println(this.getClass().getSimpleName() + ": updated params " + trk.toString());
- }
-
- tracks.add(trk);
- iteration++;
- delta_s = Math.abs(s-s_prev);
- s_prev = s;
-
- //Save distance between point and estimate
- Hep3Vector dpoint = VecOp.sub(xp, trk.getX0());
- _delta_point.add(Math.abs(VecOp.dot(eta, dpoint)));
-
-
- }
-
- if((_debug>1)) {
- System.out.println(this.getClass().getSimpleName() + ": Final s " + s_prev + " delta_s " + delta_s);
- System.out.println(this.getClass().getSimpleName() + ": Final track params " + trk.toString());
- }
-
- Hep3Vector point = trk.getX0();
- return point;
-
- }
-
-}
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