Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
RunMPAlignment.java | +76 | -78 | 1.4 -> 1.5 |
MPAlignmentParameters.java | +1202 | -507 | 1.8 -> 1.9 |
+1278 | -585 |
Large commit of alignment classes: added debug histograms, debug output, new local and global der calculations, debug for simulation.
diff -u -r1.4 -r1.5 --- RunMPAlignment.java 29 Aug 2012 02:06:22 -0000 1.4 +++ RunMPAlignment.java 28 Sep 2012 22:12:04 -0000 1.5 @@ -1,20 +1,17 @@
package org.lcsim.hps.users.phansson;
-import hep.aida.*;
import java.io.BufferedReader; import java.io.FileNotFoundException; import java.io.FileReader;
-import org.lcsim.hps.users.phansson.MPAlignmentParameters;
import java.io.IOException; import java.util.ArrayList; import java.util.List; import java.util.logging.Level; import java.util.logging.Logger;
- -import org.lcsim.event.EventHeader; -import org.lcsim.event.Track; -import org.lcsim.event.TrackerHit;
+import org.lcsim.event.*; +import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
import org.lcsim.geometry.Detector;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
import org.lcsim.util.Driver; import org.lcsim.util.aida.AIDA;
@@ -35,6 +32,8 @@
int totalTracks=0; int totalTracksProcessed=0; private String _resLimitFileName="";
+ private String outputPlotFileName=""; + private boolean hideFrame = false;
// flipSign is a kludge... // HelicalTrackFitter doesn't deal with B-fields in -ive Z correctly // so we set the B-field in +iveZ and flip signs of fitted tracks
@@ -43,22 +42,22 @@
int flipSign = 1; private boolean _debug = false; double _resLayer1MinY;
-// _resLayer2MinY,_resLayer3MinY,_resLayer4MinY,_resLayer5MinY,_resLayer6MinY,_resLayer7MinY,_resLayer8MinY,_resLayer9MinY,_resLayer10MinY; -// double _resLayer1MinX,_resLayer2MinX,_resLayer3MinX,_resLayer4MinX,_resLayer5MinX,_resLayer6MinX,_resLayer7MinX,_resLayer8MinX,_resLayer9MinX,_resLayer10MinX; -// double _resLayer1MaxY,_resLayer2MaxY,_resLayer3MaxY,_resLayer4MaxY,_resLayer5MaxY,_resLayer6MaxY,_resLayer7MaxY,_resLayer8MaxY,_resLayer9MaxY,_resLayer10MaxY; -// double _resLayer1MaxX,_resLayer2MaxX,_resLayer3MaxX,_resLayer4MaxX,_resLayer5MaxX,_resLayer6MaxX,_resLayer7MaxX,_resLayer8MaxX,_resLayer9MaxX,_resLayer10MaxX; -// - -
+ String simTrackerHitCollectionName = "TrackerHits";
public void setDebug(boolean v) { this._debug = v; }
+ public void setOutputPlotFileName(String filename) { + outputPlotFileName = filename; + }
+ public void setHideFrame(boolean hide) { + hideFrame = hide; + }
public RunMPAlignment() { nlayers[0] = 10;
@@ -81,9 +80,12 @@
}
+ @Override
public void detectorChanged(Detector detector) {
- ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug);
+ //ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug,hideFrame); + String outputMilleTxtFile = "alignMP.txt"; + ap = new MPAlignmentParameters(outputMilleTxtFile,_debug,hideFrame);
//Initialize the res limits for(int i=1;i<=10;++i) {
@@ -96,30 +98,6 @@
} } loadResidualLimits();
- /* - //ap.setResLimits(top/bottom,layer,direction: u/v/w, min, max); - ap.setResLimits(0,1,0, -0.3, 0.3); - ap.setResLimits(0,2,0, -0.3, 0.3); - ap.setResLimits(0,3,0, -0.4, 0.1); - ap.setResLimits(0,4,0, -0.2, 0.3); - ap.setResLimits(0,5,0, -0.1, 0.4); - ap.setResLimits(0,6,0, -0.3,0.1); - ap.setResLimits(0,7,0, -0.4, 0.45); - ap.setResLimits(0,8,0, -0.5, 0.3); - ap.setResLimits(0,9,0, -1, 1); - ap.setResLimits(0,10,0, -1, 1); - - ap.setResLimits(1,1,0, -0.2,0.2); - ap.setResLimits(1,2,0, -0.2, 0.2); - ap.setResLimits(1,3,0, -0.2, 0.2); - ap.setResLimits(1,4,0, -0.3, 0.3); - ap.setResLimits(1,5,0, -0.2, 0.4); - ap.setResLimits(1,6,0, -0.5, 0.); - ap.setResLimits(1,7,0, -0.7, 0.15); - ap.setResLimits(1,8,0, -0.1, 0.7); - ap.setResLimits(1,9,0, -1, 0.7); - ap.setResLimits(1,10,0, -0.5, 1); - */
@@ -127,15 +105,51 @@
+ @Override
public void process( EventHeader event) {
- // Create a map between tracks and the associated MCParticle - List<Track> tracklist = event.get(Track.class, "MatchedTracks"); -// System.out.println("Number of Tracks = " + tracklist.size()); - double duRange=0.1; - for (Track trk : tracklist) {
+ List<Track> tracklist = null; + if(event.hasCollection(Track.class,"MatchedTracks")) { + tracklist = event.get(Track.class, "MatchedTracks"); + if(_debug) { + System.out.println("Number of Tracks = " + tracklist.size()); + } + } + + List<SimTrackerHit> simTrackerHits = new ArrayList<SimTrackerHit>(); + if(event.hasCollection(SimTrackerHit.class, simTrackerHitCollectionName)){ + simTrackerHits.addAll(event.get(SimTrackerHit.class, simTrackerHitCollectionName)); + if(_debug){ + System.out.println(this.getClass().getSimpleName() + ": Event: " + event.getEventNumber() + " Number of SimTrackerHits: " + simTrackerHits.size()); + System.out.print(this.getClass().getSimpleName() + ": MC Particles: "); + for(MCParticle mcParticle : event.getMCParticles()){ + System.out.print(mcParticle.getPDGID() + " "); + } + System.out.print("\n"); + } + } + List<HelicalTrackStrip> strips = null; + if(event.hasCollection(HelicalTrackStrip.class, "HelicalTrackStrips")) { + strips = event.get(HelicalTrackStrip.class, "HelicalTrackStrips"); + if(_debug) System.out.println("Event has " + strips.size() + "HelicalTrackStrips"); + } + + List<SiTrackerHit> trackerHits = null; + + if(event.hasCollection(SiTrackerHit.class, "StripClusterer_SiTrackerHitStrip1D")) { + trackerHits = event.get(SiTrackerHit.class, "StripClusterer_SiTrackerHitStrip1D"); + if(_debug) System.out.println("Event has " + trackerHits.size() + "SiTrackerHit"); + } + + + + for (Track trk : tracklist) { + + //if(trk.getCharge()<0) continue; + if(trk.getTrackStates().get(0).getMomentum()[0]<0.8) continue; +
totalTracks++; if(hitsOnBothSides(trk)) continue;
@@ -143,40 +157,22 @@
totalTracksProcessed++;
- ap.PrintResidualsAndDerivatives(trk,totalTracks); - - if(1==1){ - aida.histogram1D("Track d0",50,-0.5,0.5).fill(trk.getTrackParameter(0)); - aida.histogram1D("Track sin(phi0)",50,-0.5,0.5).fill(Math.sin(trk.getTrackParameter(1))); - aida.histogram1D("Track z0",50,-0.1,0.1).fill(Math.sin(trk.getTrackParameter(3))); - aida.histogram1D("Track chi^2",50,0,25).fill(trk.getChi2()); - for (int i = 1; i < 11; i++) { - double[] res = ap.getResidual(trk, i); - int mylayer=(int)res[6]; - if(mylayer<11){ - aida.histogram1D("Track chi^2 Positive Side",50,0,25).fill(trk.getChi2()); - }else{ - aida.histogram1D("Track chi^2 Negative Side",50,0,25).fill(trk.getChi2()); - }
+ ap.PrintResidualsAndDerivatives(trk,totalTracks,"LOCAL",simTrackerHits); + + +
- aida.histogram1D("deltaU -- Layer " + mylayer,50,-duRange,duRange).fill(res[0]); - aida.histogram1D("deltaU Pull-- Layer " + mylayer,50,-3,3).fill(res[0]/res[3]); - if(i==3&&Math.sin(trk.getTrackParameter(1))>0){ - aida.histogram1D("Positive phi0 deltaU -- Layer " + mylayer,50,-duRange,duRange).fill(res[0]); - aida.histogram1D("Positive phi0 deltaU Pull-- Layer " + mylayer,50,-3,3).fill(res[0]/res[3]); - } - if(i==3&&Math.sin(trk.getTrackParameter(1))<0){ - aida.histogram1D("Negative phi0 deltaU -- Layer " + mylayer,50,-duRange,duRange).fill(res[0]); - aida.histogram1D("Negative phi0 deltaU Pull-- Layer " + mylayer,50,-3,3).fill(res[0]/res[3]); - } - - } - } - }
+ }
+ ap.CalculateResidualSimHits(strips,trackerHits,simTrackerHits); + + + +
} public void endOfData() {
+ ap.updatePlots();
try { System.out.println("Total Number of Tracks Found = "+totalTracks); System.out.println("Total Number of Tracks Processed = "+totalTracksProcessed);
@@ -184,6 +180,13 @@
} catch (IOException ex) { Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE, null, ex); }
+ if (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); + } +
}
@@ -247,16 +250,11 @@
}
- ap.setResLimits(0,1,0, -0.3, 0.3);
+ ap.setResLimits(0,1,0, -0.3, 0.3);
}
- - public void setResLayer1MinY(double val) { - this._resLayer1MinY = val; - } -
public void setResidualLimitFileName(String fileName) { this._resLimitFileName = fileName; }
diff -u -r1.8 -r1.9 --- MPAlignmentParameters.java 4 Sep 2012 03:01:13 -0000 1.8 +++ MPAlignmentParameters.java 28 Sep 2012 22:12:04 -0000 1.9 @@ -21,12 +21,14 @@
import org.lcsim.detector.tracker.silicon.ChargeCarrier; import org.lcsim.detector.tracker.silicon.SiSensor; import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
-import org.lcsim.event.RawTrackerHit; -import org.lcsim.event.Track; -import org.lcsim.event.TrackerHit;
+import org.lcsim.event.*;
import org.lcsim.fit.helicaltrack.*; import org.lcsim.hps.event.HPSTransformations; import org.lcsim.hps.monitoring.AIDAFrame;
+import org.lcsim.hps.recon.tracking.TrackUtils; +import org.lcsim.hps.recon.tracking.TrackerHitUtils; +import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit; +import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate; import org.lcsim.recon.tracking.seedtracker.SeedTrack; import org.lcsim.util.aida.AIDA;
@@ -62,38 +64,47 @@
private double[] _error = new double[3]; FileWriter fWriter; PrintWriter pWriter;
- Set<SiSensor> _process_sensors = new HashSet<SiSensor>();
private HPSTransformations _detToTrk; boolean _DEBUG=false; double smax = 1e3; ResLimit _resLimits = new ResLimit();
- AlignmentUtils _alignUtils = new AlignmentUtils(_DEBUG);
+ AlignmentUtils _alignUtils; + TrackUtils trackUtil; + TrackerHitUtils trackerHitUtil;
private AIDAFrame plotterFrame;
+ private AIDAFrame plotterFrameHTH; + private AIDAFrame plotterFrameSimHit; + private AIDAFrame plotterFrameSummary; + private boolean hideFrame = false;
private AIDA aida = AIDA.defaultInstance(); private IAnalysisFactory af = aida.analysisFactory();
- //private List< List<IHistogram1D> > hres = new ArrayList<List<IHistogram1D> >(); -// private List< List<IHistogram1D> > resy_org_layallhit = new ArrayList<List<IHistogram1D> >(); -// private List< List<IHistogram1D> > resy_org_layallsinglehit = new ArrayList<List<IHistogram1D> >(); -// private List< List<IHistogram1D> > resy_org_lay1hit = new ArrayList<List<IHistogram1D> >(); -// private List< IHistogram1D > nhits_tracker = new ArrayList<IHistogram1D>(); -// private List< List<IHistogram1D> > nhits_tracker_layer = new ArrayList<List<IHistogram1D> >(); -// private List< IHistogram1D > ncl_ecal = new ArrayList<IHistogram1D>(); -// private List< IHistogram1D > selcl_ecal_e = new ArrayList<IHistogram1D>(); -// private List< IHistogram1D > cl_ecal_e = new ArrayList<IHistogram1D>(); -// private List< IHistogram2D > ncl_ecal_map = new ArrayList<IHistogram2D>(); -// private List< IHistogram2D > nselcl_ecal_map = new ArrayList<IHistogram2D>(); -// private List< IHistogram2D> resy_2d_org_layallhit = new ArrayList<IHistogram2D>(); -// private List< IHistogram2D> resy_2d_org_layallsinglehit = new ArrayList<IHistogram2D>(); -//
+ IDataPointSet dps_t; + IDataPointSet dps_b; + IDataPointSet dps_hth_y_b; + IDataPointSet dps_hth_y_t; + IDataPointSet dps_hth_z_b; + IDataPointSet dps_hth_z_t; + IDataPointSet dps_pull_t; + IDataPointSet dps_pull_b; + IDataPointSet dps_pull_hth_y_b; + IDataPointSet dps_pull_hth_y_t; + IDataPointSet dps_pull_hth_z_b; + IDataPointSet dps_pull_hth_z_t;
+
- public MPAlignmentParameters(String outfile, boolean debug) {
+ + public MPAlignmentParameters(String outfile, boolean debug, boolean hidePlots) {
_DEBUG = debug; _detToTrk = new HPSTransformations();
+ _alignUtils = new AlignmentUtils(_DEBUG); + trackUtil = new TrackUtils(); + trackerHitUtil = new TrackerHitUtils(_DEBUG); + hideFrame = hidePlots;
try { //open things up fWriter = new FileWriter(outfile);
@@ -106,174 +117,7 @@
//if(_DEBUG) _globalParameterList.print();
- - - - - - aida.tree().cd("/"); - plotterFrame = new AIDAFrame(); - plotterFrame.setTitle("HPS Alignment Plots"); - - - - //IHistogramFactory hf = aida.histogramFactory(); - String[] side = {"top","bottom"}; - String[] direction = {"u","v","w"}; - - for (int d=0;d<3;++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[d] + " " + side[s]); - IPlotter plotter_reserr = af.createPlotterFactory().create(); - plotter_reserr.createRegions(5,2,0); - plotter_reserr.setTitle("reserr_" + direction[d] + " " + side[s]); - IPlotter plotter_respull = af.createPlotterFactory().create(); - plotter_respull.createRegions(5,2,0); - plotter_respull.setTitle("respull_" + direction[d] + " " + side[s]); - - for (int iLayer=1;iLayer<11;++iLayer) { - double xmin,xmax; - if(direction[d]=="v") { - xmin = -100; - xmax = 100; - } else { - xmin = -2; - xmax = 2; - - } - // IHistogram h = aida.histogram1D("res_"+direction[d]+"_layer" + (iLayer) + "_" + side[iSide] , 50, -10, 10); - IHistogram h = aida.histogram1D("res_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, xmin, xmax); - IHistogram h1 = aida.histogram1D("reserr_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); - IHistogram h2 = aida.histogram1D("respull_"+direction[d]+"_layer" + (iLayer) + "_" + side[s] , 50, -20, 20); - - int region = (iLayer-1);//*2+iSide; - - plotter_res.region(region).plot(h); - plotter_reserr.region(region).plot(h1); - plotter_respull.region(region).plot(h2); - - ((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); - - - } - - //plotter_res.show(); - plotterFrame.addPlotter(plotter_res); - plotterFrame.addPlotter(plotter_reserr); - plotterFrame.addPlotter(plotter_respull); - } - } - - - for (int s=0;s<2;++s) { - - IPlotter plotter_hthresy = af.createPlotterFactory().create(); - plotter_hthresy.createRegions(5,1,0); - plotter_hthresy.setTitle("hthres_y " + side[s]); - IPlotter plotter_hthresypull = af.createPlotterFactory().create(); - plotter_hthresypull.createRegions(5,1,0); - plotter_hthresypull.setTitle("hthrespull_y " + side[s]); - IPlotter plotter_hthresyerr = af.createPlotterFactory().create(); - plotter_hthresyerr.createRegions(5,1,0); - plotter_hthresyerr.setTitle("hthreserr_y " + side[s]); - IPlotter plotter_hthresz = af.createPlotterFactory().create(); - plotter_hthresz.createRegions(5,1,0); - plotter_hthresz.setTitle("hthres_z " + side[s]); - IPlotter plotter_hthreszerr = af.createPlotterFactory().create(); - plotter_hthreszerr.createRegions(5,1,0); - plotter_hthreszerr.setTitle("hthreserr_z " + side[s]); - IPlotter plotter_hthreszpull = af.createPlotterFactory().create(); - plotter_hthreszpull.createRegions(5,1,0); - plotter_hthreszpull.setTitle("hthrespull_z " + side[s]); - IPlotter plotter_hthreszerrms = af.createPlotterFactory().create(); - plotter_hthreszerrms.createRegions(5,1,0); - plotter_hthreszerrms.setTitle("hthreserrms_z " + side[s]); - IPlotter plotter_hthreszerrres = af.createPlotterFactory().create(); - plotter_hthreszerrres.createRegions(5,1,0); - plotter_hthreszerrres.setTitle("hthreserrres_z " + side[s]); - - - - - - for (int iLayer=1;iLayer<6;++iLayer) { - - - IHistogram h = aida.histogram1D("hthres_y_layer" + (iLayer) + "_" + side[s] , 50, -2, 2); - IHistogram h1 = aida.histogram1D("hthreserr_y_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); - IHistogram h2 = aida.histogram1D("hthrespull_y_layer" + (iLayer) + "_" + side[s] , 50, -5, 5); - - IHistogram hz = aida.histogram1D("hthres_z_layer" + (iLayer) + "_" + side[s] , 50, -2, 2); - IHistogram hz1 = aida.histogram1D("hthreserr_z_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); - IHistogram hz2 = aida.histogram1D("hthrespull_z_layer" + (iLayer) + "_" + side[s] , 50, -5, 5); - IHistogram hz11 = aida.histogram1D("hthreserrms_z_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); - IHistogram hz111 = aida.histogram1D("hthreserrres_z_layer" + (iLayer) + "_" + side[s] , 50, 0, 1); - - - - - - int region = (iLayer-1);//*2+iSide; - - plotter_hthresy.region(region).plot(h); - plotter_hthresyerr.region(region).plot(h1); - plotter_hthresypull.region(region).plot(h2); - - plotter_hthresz.region(region).plot(hz); - plotter_hthreszerr.region(region).plot(hz1); - plotter_hthreszpull.region(region).plot(hz2); - - plotter_hthreszerrms.region(region).plot(hz11); - plotter_hthreszerrres.region(region).plot(hz111); - - - ((PlotterRegion) plotter_hthresy.region(region)).getPlot().setAllowUserInteraction(true); - ((PlotterRegion) plotter_hthresy.region(region)).getPlot().setAllowPopupMenus(true); - ((PlotterRegion) plotter_hthresyerr.region(region)).getPlot().setAllowUserInteraction(true); - ((PlotterRegion) plotter_hthresyerr.region(region)).getPlot().setAllowPopupMenus(true); - ((PlotterRegion) plotter_hthresypull.region(region)).getPlot().setAllowUserInteraction(true); - ((PlotterRegion) plotter_hthresypull.region(region)).getPlot().setAllowPopupMenus(true); - - ((PlotterRegion) plotter_hthresz.region(region)).getPlot().setAllowUserInteraction(true); - ((PlotterRegion) plotter_hthresz.region(region)).getPlot().setAllowPopupMenus(true); - ((PlotterRegion) plotter_hthreszerr.region(region)).getPlot().setAllowUserInteraction(true); - ((PlotterRegion) plotter_hthreszerr.region(region)).getPlot().setAllowPopupMenus(true); - ((PlotterRegion) plotter_hthreszpull.region(region)).getPlot().setAllowUserInteraction(true); - ((PlotterRegion) plotter_hthreszpull.region(region)).getPlot().setAllowPopupMenus(true); - - ((PlotterRegion) plotter_hthreszerrms.region(region)).getPlot().setAllowUserInteraction(true); - ((PlotterRegion) plotter_hthreszerrms.region(region)).getPlot().setAllowPopupMenus(true); - ((PlotterRegion) plotter_hthreszerrres.region(region)).getPlot().setAllowUserInteraction(true); - ((PlotterRegion) plotter_hthreszerrres.region(region)).getPlot().setAllowPopupMenus(true); - - - - } - - plotterFrame.addPlotter(plotter_hthresy); - plotterFrame.addPlotter(plotter_hthresyerr); - plotterFrame.addPlotter(plotter_hthresypull); - plotterFrame.addPlotter(plotter_hthresz); - plotterFrame.addPlotter(plotter_hthreszerr); - plotterFrame.addPlotter(plotter_hthreszpull); - plotterFrame.addPlotter(plotter_hthreszerrms); - plotterFrame.addPlotter(plotter_hthreszerrres); - } - - - plotterFrame.pack(); - plotterFrame.setVisible(true); - -
+ makeAlignmentPlots();
@@ -289,7 +133,7 @@
_resLimits.add(s,l,d, low,high); }
- public void PrintResidualsAndDerivatives(Track track, int itrack) {
+ public void PrintResidualsAndDerivatives(Track track, int itrack, String type, List<SimTrackerHit> simhits) {
SeedTrack st = (SeedTrack) track; SeedCandidate seed = st.getSeedCandidate();
@@ -298,6 +142,7 @@
List<TrackerHit> hitsOnTrack = track.getTrackerHits(); String half = hitsOnTrack.get(0).getPosition()[2]>0 ? "top" : "bottom"; pWriter.printf("TRACK %s (%d)\n",half,itrack);
+ aida.histogram1D("Track Chi2 "+ half).fill(track.getChi2());
for (TrackerHit hit : hitsOnTrack) { HelicalTrackHit htc = (HelicalTrackHit) hit; double msdrphi = msmap.get(htc).drphi();
@@ -309,20 +154,32 @@
TrackDirection trkdir = HelixUtils.CalculateTrackDirection(_trk, sHit); cross.setTrackDirection(trkdir, _trk.covariance()); for (HelicalTrackStrip cl : clusterlist) {
- CalculateLocalDerivatives(cl); - CalculateGlobalDerivatives(cl); - CalculateResidual(cl, msdrphi, msdz); -// CalculateResidual(cl, 0,0); - //PrintStripResiduals(cl); - PrintStripResidualsNew(cl); -
+ if(type=="GLOBAL") { + + //CalculateLocalDerivativesGLOBAL(cl); + CalculateGlobalDerivatives(cl); + CalculateResidualGLOBAL(cl, msdrphi, msdz); + + } else { + CalculateResidual(cl, msdrphi, msdz); + //CalculateResidualSim(cl, msdrphi, msdz,simhits); + + CalculateLocalDerivatives(cl); + CalculateGlobalDerivatives(cl); + + // CalculateResidual(cl, 0,0); + //PrintStripResiduals(cl); + PrintStripResidualsNew(cl); + }
} }
+ + if(itrack%50==0) this.updatePlots();
//AddTarget(0.1, 0.02); }
- private void CalculateLocalDerivatives(HelicalTrackStrip strip) { - //get track parameters.
+ private BasicMatrix CalculateLocalDerivativesOld(HelicalTrackStrip strip) { + //get track parameters.
double d0 = _trk.dca(); double z0 = _trk.z0(); double slope = _trk.slope();
@@ -348,15 +205,33 @@
} BasicMatrix dfdqGlobal = FillMatrix(dfdq, 3, 5);
- BasicMatrix dfdqGlobalNew = _alignUtils.calculateLocalHelixDerivatives(_trk, strip, smax, _nlc); - Hep3Matrix trkToStrip = getTrackToStripRotation(strip);
+ return dfdqGlobal; + } + + private void CalculateLocalDerivatives(HelicalTrackStrip strip) { + + BasicMatrix dfdqGlobalOld = CalculateLocalDerivativesOld(strip); + + BasicMatrix dfdqGlobal = _alignUtils.calculateLocalHelixDerivatives(_trk, strip, smax, _nlc); + Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip);
//_dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal);
- _dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobalNew); -
+ _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(); + //strip origin is defined in the tracking coordinate system (x=beamline) + double xint = strip.origin().x(); + double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0); + double phi = -s/R + phi0;
double[] trackpars = {d0, z0, slope, phi0, R, s, xint};
- System.out.println(" d0 z0 slope phi0 R"); - System.out.printf("Values %5.5f %5.5f %5.5f %5.5f %5.5f\n", d0, z0, slope, phi0, R);
+ System.out.printf("%10s%10s%10s%10s%10s\n","d0","z0","slope","phi0","R"); + System.out.printf("%10.5f%10.5f%10.5f%10.5f%10.5f\n", d0, z0, slope, phi0, R);
System.out.println("Strip Origin: "); System.out.println(strip.origin()); System.out.println("s xint");
@@ -364,12 +239,17 @@
System.out.println("Compare local derivatives "); System.out.println("dfdqGlobal:"); System.out.println(dfdqGlobal.toString());
- System.out.println("dfdqGlobalNew:"); - System.out.println(dfdqGlobalNew.toString());
+ System.out.println("dfdqGlobalOld:"); + System.out.println(dfdqGlobalOld.toString());
} }
+ + + + +
private void CalculateGlobalDerivatives(HelicalTrackStrip strip) {
@@ -383,6 +263,10 @@
double z0 = _trk.z0(); double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0); double phi = -s/R + phi0;
+ double umeas = strip.umeas(); + Hep3Vector corigin = strip.origin(); + double vmeas = corigin.y(); //THis is wrong +
if(_DEBUG) { System.out.println("Calculate global derivatives for this strip hit in layer " + strip.layer());
@@ -390,8 +274,8 @@
System.out.printf(" %10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n", d0, z0, slope, phi0, R,xint,phi); }
- -
+ if(1==1) + return;
//1st index = alignment parameter (only u so far) //2nd index = residual coordinate (on du so far)
@@ -421,109 +305,287 @@
int side = strip.origin().z()>0 ? 10000 : 20000; int layer = strip.layer();
- //Flag to tell if this hit is affected by the given global parameter
+ //Rotation matrix from the tracking fram to the sensor/strip frame + Hep3Matrix T = this.trackerHitUtil.getTrackToStripRotation(strip); + + //Derivatives of the rotation matrix deltaR' w.r.t. rotations a,b,c around axis x,y,z + Hep3Matrix ddeltaRprime_da = new BasicHep3Matrix(0,0,0,0,0,1,0,-1,0); + Hep3Matrix ddeltaRprime_db = new BasicHep3Matrix(0,0,-1,0,0,0,1,0,0); + Hep3Matrix ddeltaRprime_dc = new BasicHep3Matrix(0,1,0,-1,0,0,0,0,0); + + + //**************************************************************************** + //Calculate jacobian da/db + + // Upper left 3x3 + Hep3Vector deltaX_gl = new BasicHep3Vector(1,0,0); + Hep3Vector deltaY_gl = new BasicHep3Vector(0,1,0); + Hep3Vector deltaZ_gl = new BasicHep3Vector(0,0,1); + Hep3Vector dq_dx = VecOp.mult(T, deltaX_gl); + Hep3Vector dq_dy = VecOp.mult(T, deltaY_gl); + Hep3Vector dq_dz = VecOp.mult(T, deltaZ_gl); + + // Upper right 3x3 + Hep3Vector x_vec = new BasicHep3Vector(xint,corigin.y(),corigin.z()); //THIS IS WRONG!!!! + Hep3Vector x_vec_tmp = VecOp.mult(ddeltaRprime_da, x_vec); + Hep3Vector dq_da = VecOp.mult(T,x_vec_tmp); + x_vec_tmp = VecOp.mult(ddeltaRprime_db, x_vec); + Hep3Vector dq_db = VecOp.mult(T,x_vec_tmp); + x_vec_tmp = VecOp.mult(ddeltaRprime_dc, x_vec); + Hep3Vector dq_dc = VecOp.mult(T,x_vec_tmp); + + // Lower left 3x3 + //deltaR = (alpha,beta,gamma) the rotation alignment parameters in the local frame + Hep3Vector ddeltaR_dx = new BasicHep3Vector(0,0,0); + Hep3Vector ddeltaR_dy = new BasicHep3Vector(0,0,0); + Hep3Vector ddeltaR_dz = new BasicHep3Vector(0,0,0); + + // Lower right 3x3 + //deltaR = (alpha,beta,gamma) the rotation alignment parameters in the local frame + //Expressing T in Euler angles i,j,k + double i = 0; + double j = 0; + double k = 0; + double dalpha_da = Math.cos(k)*Math.sin(i)*Math.sin(j) + Math.cos(i)*Math.sin(k); + double dbeta_da = Math.cos(i)*Math.cos(k) - Math.sin(i)*Math.sin(j)*Math.sin(k); + double dgamma_da = -Math.cos(j)*Math.sin(i); + double dalpha_db = -Math.cos(i)*Math.cos(k)*Math.sin(j) + Math.sin(i)*Math.sin(k); + double dbeta_db = Math.cos(k)*Math.sin(i) + Math.cos(i)*Math.sin(j)*Math.sin(k); + double dgamma_db = Math.cos(i)*Math.cos(j); + double dalpha_dc = -Math.cos(i)*Math.cos(k)*Math.sin(j) + Math.sin(i)*Math.sin(k); + double dbeta_dc = Math.cos(k)*Math.sin(i) + Math.cos(i)*Math.sin(j)*Math.sin(k); + double dgamma_dc = Math.sin(i)*Math.sin(j); + + Hep3Vector ddeltaR_da = new BasicHep3Vector(dalpha_da,dbeta_da,dgamma_da); + Hep3Vector ddeltaR_db = new BasicHep3Vector(dalpha_db,dbeta_db,dgamma_db); + Hep3Vector ddeltaR_dc = new BasicHep3Vector(dalpha_dc,dbeta_dc,dgamma_dc); + + //Now fill the Jacobian 6x6 matrix + BasicMatrix da_db = new BasicMatrix(6,6); + //upper left 3x3 + da_db.setElement(0,0,dq_dx.x()); + da_db.setElement(1,0,dq_dx.y()); + da_db.setElement(2,0,dq_dx.z()); + da_db.setElement(0,1,dq_dy.x()); + da_db.setElement(1,1,dq_dy.y()); + da_db.setElement(2,1,dq_dy.z()); + da_db.setElement(0,2,dq_dz.x()); + da_db.setElement(1,2,dq_dz.y()); + da_db.setElement(2,2,dq_dz.z()); + //upper right 3x3 + da_db.setElement(0,3,dq_da.x()); + da_db.setElement(1,3,dq_da.y()); + da_db.setElement(2,3,dq_da.z()); + da_db.setElement(0,4,dq_db.x()); + da_db.setElement(1,4,dq_db.y()); + da_db.setElement(2,4,dq_db.z()); + da_db.setElement(0,5,dq_dc.x()); + da_db.setElement(1,5,dq_dc.y()); + da_db.setElement(2,5,dq_dc.z()); + //lower right 3x3 + da_db.setElement(3,3,ddeltaR_da.x()); + da_db.setElement(4,3,ddeltaR_da.y()); + da_db.setElement(5,3,ddeltaR_da.z()); + da_db.setElement(3,4,ddeltaR_db.x()); + da_db.setElement(4,4,ddeltaR_db.y()); + da_db.setElement(5,4,ddeltaR_db.z()); + da_db.setElement(3,5,ddeltaR_dc.x()); + da_db.setElement(4,5,ddeltaR_dc.y()); + da_db.setElement(5,5,ddeltaR_dc.z()); + //lower left 3x3 + da_db.setElement(3,0,ddeltaR_dx.x()); + da_db.setElement(4,0,ddeltaR_dx.y()); + da_db.setElement(5,0,ddeltaR_dx.z()); + da_db.setElement(3,1,ddeltaR_dy.x()); + da_db.setElement(4,1,ddeltaR_dy.y()); + da_db.setElement(5,1,ddeltaR_dy.z()); + da_db.setElement(3,2,ddeltaR_dz.x()); + da_db.setElement(4,2,ddeltaR_dz.y()); + da_db.setElement(5,2,ddeltaR_dz.z()); + + if(_DEBUG) { + System.out.println("da_db: \n" + da_db.toString()); + } + + + BasicMatrix db_da = da_db; + db_da.invert(); + + + if(_DEBUG) { + System.out.println("db_da: \n" + db_da.toString()); + } + + //**************************************************************************** + // Calcualte the global derivatives in the global frame dz/db + // where z is the residual in the global frame and b are the alignment + // parameters in the global frame + + + + //Flag to tell if this hit is affected by the given global parameter
boolean useGL = false; //Clear the old parameter list _glp.clear();
-
+
+
+
+
//****************************************************************************
- //Calculate derivatives for a residual in x,y,z-direction for translation in x
+ //Derivatives of the predicted hit position qp for a translation of in x
- BasicMatrix dfdpTRACK = new BasicMatrix(3,1); - dfdpTRACK.setElement(0, 0, _alignUtils.dx_dx()); - dfdpTRACK.setElement(1, 0, _alignUtils.dy_dx(xint,xr,yr,d0,phi0,R,phi)); - dfdpTRACK.setElement(2, 0, _alignUtils.dz_dx(xint, xr, yr, d0, phi0, slope, R, phi));
+ BasicMatrix dqp_da_TRACK = new BasicMatrix(3,1); + dqp_da_TRACK.setElement(0, 0, _alignUtils.dx_dx()); + dqp_da_TRACK.setElement(1, 0, _alignUtils.dy_dx(xint, d0, phi0, R, phi)); + dqp_da_TRACK.setElement(2, 0, _alignUtils.dz_dx(xint, xr, yr, d0, phi0, slope, R, phi));
//Put it into a matrix to be able to transform easily
- //BasicMatrix _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1);
+ //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 = getTrackToStripRotation(strip);
+ Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip);
if(_DEBUG) { System.out.println("Final transformation from tracking frame to strip frame:\n" + trkToStrip.toString()); }
- - //Transform derivatives to sensor frame! - BasicMatrix dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdpTRACK);
+ //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(dfdp);
+ gp_tx.setDfDp(dqp_da);
_glp.add(gp_tx); if (_DEBUG) { gp_tx.print();
- System.out.println("Track frame dfdp: " + dfdpTRACK); - //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());
+ 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());
} //****************************************************************************
- //Calculate derivatives for a residual in x,y,z-direction for translation in y - - dfdpTRACK.setElement(0, 0, _alignUtils.dx_dy(xint,xr,yr,d0,phi0,R,phi)); - dfdpTRACK.setElement(1, 0, _alignUtils.dy_dy()); - dfdpTRACK.setElement(2, 0, _alignUtils.dz_dy(xint, xr, d0, phi0, slope, R, phi));
+ //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(xint-(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 _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1);
+ //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1);
//Transform derivatives to sensor frame!
- dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdpTRACK);
+ 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(dfdp);
+ gp_ty.setDfDp(dqp_da);
_glp.add(gp_ty); if (_DEBUG) { gp_ty.print();
- System.out.println("Track frame dfdp: " + dfdpTRACK);
+ 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()); } //****************************************************************************
- //Calculate derivatives for a residual in x,y,z-direction for translation in z
+ //Derivatives of the predicted hit position qp for a translation of in z
- dfdpTRACK.setElement(0, 0, _alignUtils.dx_dz(xint, xr, yr, d0, phi0, slope, R, phi)); - dfdpTRACK.setElement(1, 0, _alignUtils.dy_dz(xint, xr, yr, d0, phi0, slope, R, phi)); - dfdpTRACK.setElement(2, 0, _alignUtils.dz_dz());
+ 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 _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1);
+ //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1);
//Transform derivatives to sensor frame!
- dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdpTRACK);
+ 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(dfdp);
+ gp_tz.setDfDp(dqp_da);
_glp.add(gp_tz); if (_DEBUG) { gp_tz.print();
- System.out.println("Track frame dfdp: " + dfdpTRACK);
+ 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()); }
+ //I need to convert the dqp/da in the global frame to dqp_da in the sensor fram + + //place holders + double dup_dw = 0; + double dvp_dw = 0; + double um = umeas; + double vm = vmeas; //should always be zero, I think + + BasicMatrix dres_ddu = new BasicMatrix(3,1); + dres_ddu.setElement(0, 0, 1); + dres_ddu.setElement(1, 0, 0); + dres_ddu.setElement(2, 0, 0); + BasicMatrix dres_ddv = new BasicMatrix(3,1); + dres_ddv.setElement(0, 0, 0); + dres_ddv.setElement(1, 0, 1); + dres_ddv.setElement(2, 0, 0); + BasicMatrix dres_ddw = new BasicMatrix(3,1); + dres_ddw.setElement(0, 0, -dup_dw); + dres_ddw.setElement(1, 0, -dvp_dw); + dres_ddw.setElement(2, 0, 0); + BasicMatrix dres_dalpha = new BasicMatrix(3,1); + dres_dalpha.setElement(0, 0, vm*dup_dw); + dres_dalpha.setElement(1, 0, vm*dvp_dw); + dres_dalpha.setElement(2, 0, 0); + BasicMatrix dres_dbeta = new BasicMatrix(3,1); + dres_dbeta.setElement(0, 0, -um*dup_dw); + dres_dbeta.setElement(1, 0, -um*dvp_dw); + dres_dbeta.setElement(2, 0, 0); + BasicMatrix dres_dgamma = new BasicMatrix(3,1); + dres_dgamma.setElement(0, 0, vm); + dres_dgamma.setElement(1, 0, -um); + dres_dgamma.setElement(2, 0, 0); + + //compact version + BasicMatrix dres_da = new BasicMatrix(3,6); + for(int icol=0;icol<6;++icol) { + for(int irow=0;irow<3;++irow) { + if(icol==0) dres_da.setElement(irow, icol, dres_ddu.e(irow, 0)); + if(icol==1) dres_da.setElement(irow, icol, dres_ddv.e(irow, 0)); + if(icol==2) dres_da.setElement(irow, icol, dres_ddw.e(irow, 0)); + if(icol==3) dres_da.setElement(irow, icol, dres_dalpha.e(irow, 0)); + if(icol==4) dres_da.setElement(irow, icol, dres_dbeta.e(irow, 0)); + if(icol==5) dres_da.setElement(irow, icol, dres_dgamma.e(irow, 0)); + } + } + + + + + + + + + + + /* + +
//****************************************************************************
- //Calculate derivatives for a residual in x,y,z-direction for rotation alpha (rotation around beamline direction x )
+ //Derivatives of the predicted hit position qp for a rotation alpha around the x-axis
- dfdpTRACK.setElement(0, 0, _alignUtils.dx_dalpha()); - dfdpTRACK.setElement(1, 0, _alignUtils.dy_dalpha()); - dfdpTRACK.setElement(2, 0, _alignUtils.dz_dalpha());
+ dqp_da_TRACK.setElement(0, 0, _alignUtils.dx_dalpha()); + dqp_da_TRACK.setElement(1, 0, _alignUtils.dy_dalpha()); + dqp_da_TRACK.setElement(2, 0, _alignUtils.dz_dalpha());
//Put it into a matrix to be able to transform easily
- //BasicMatrix _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1);
+ //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1);
//Transform derivatives to sensor frame!
- dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdpTRACK);
+ dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dqp_da_TRACK);
//Add it to the global parameter object GlobalParameter gp_ra = new GlobalParameter("Rotation alpha",side,layer,2000,100,true); gp_ra.setDfDp(dfdp); _glp.add(gp_ra); if (_DEBUG) { gp_ra.print();
- System.out.println("Track frame dfdp: " + dfdpTRACK);
+ System.out.println("Track frame dfdp: " + 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()); }
@@ -531,22 +593,22 @@
//**************************************************************************** //Calculate derivatives for a residual in x,y,z-direction for rotation beta (rotation around beamline direction y )
- dfdpTRACK.setElement(0, 0, _alignUtils.dx_dbeta()); - dfdpTRACK.setElement(1, 0, _alignUtils.dy_dbeta()); - dfdpTRACK.setElement(2, 0, _alignUtils.dz_dbeta());
+ dqp_da_TRACK.setElement(0, 0, _alignUtils.dx_dbeta()); + dqp_da_TRACK.setElement(1, 0, _alignUtils.dy_dbeta()); + dqp_da_TRACK.setElement(2, 0, _alignUtils.dz_dbeta());
//Put it into a matrix to be able to transform easily
- //BasicMatrix _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1);
+ //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1);
//Transform derivatives to sensor frame!
- dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdpTRACK);
+ dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dqp_da_TRACK);
//Add it to the global parameter object GlobalParameter gp_rb = new GlobalParameter("Rotation beta",side,layer,2000,200,true); gp_rb.setDfDp(dfdp); _glp.add(gp_rb); if (_DEBUG) { gp_rb.print();
- System.out.println("Track frame dfdp: " + dfdpTRACK);
+ System.out.println("Track frame dfdp: " + 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()); }
@@ -554,137 +616,38 @@
//**************************************************************************** //Calculate derivatives for a residual in x,y,z-direction for rotation gamma (rotation around beamline direction z )
- dfdpTRACK.setElement(0, 0, _alignUtils.dx_dgamma()); - dfdpTRACK.setElement(1, 0, _alignUtils.dy_dgamma()); - dfdpTRACK.setElement(2, 0, _alignUtils.dz_dgamma());
+ dqp_da_TRACK.setElement(0, 0, _alignUtils.dx_dgamma()); + dqp_da_TRACK.setElement(1, 0, _alignUtils.dy_dgamma()); + dqp_da_TRACK.setElement(2, 0, _alignUtils.dz_dgamma());
//Put it into a matrix to be able to transform easily
- //BasicMatrix _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1);
+ //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1);
//Transform derivatives to sensor frame!
- dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdpTRACK);
+ dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, dqp_da_TRACK);
//Add it to the global parameter object GlobalParameter gp_rg = new GlobalParameter("Rotation gamma",side,layer,2000,300,true); gp_rg.setDfDp(dfdp); _glp.add(gp_rg); if (_DEBUG) { gp_rg.print();
- System.out.println("Track frame dfdp: " + dfdpTRACK);
+ System.out.println("Track frame dfdp: " + 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()); }
+ */
- - -// //Go through each global parameter and see if it has non-zero contribution -// for (GlobalParameter gp : _globalParameterList.getList()) { -// -// -// useGL = false; -// -//// dfdpTRACK.setElement(0, 0, 0); -//// dfdpTRACK.setElement(1, 0, 0); -//// dfdpTRACK.setElement(2, 0, 0); -//// -// -// //Is this global parameter on the same side -// if(gp.getSide()==side) { -// -// -// // correct side -// if(gp.getLayer()==layer) { -// //correct layer -// -// //Calcuate dfdp derivatives depending on type of global parameter -// if(gp.getType()==1000) { -// -// //This residual is affected by a global parameter -// useGL = true; -// -// //translation -// if(gp.getDirection()==100) { -// //x - tracking frame --> beamline direction -// dfdpTRACK[0][0] = 1; //df/dx -// dfdpTRACK[1][0] = 0; //df/dy -// dfdpTRACK[2][0] = 0; //df/dz -// } -// else if(gp.getDirection()==200) { -// //y - tracking frame --> almost unmeasured coordinate direction -// dfdpTRACK[0][0] = 0; //df/dx -// dfdpTRACK[1][0] = 1; //df/dy -// dfdpTRACK[2][0] = 0; //df/dz -// -// } -// else { -// //z - tracking frame --> almost measured coordinate direction -// dfdpTRACK[0][0] = 0; //df/dx -// dfdpTRACK[1][0] = 0; //df/dy -// dfdpTRACK[2][0] = 1; //df/dz -// } -// -// }//type -// else if(gp.getType()==2000) { -// //rotation -// -// //This residual is affected by a global parameter -// useGL = true; -// -// -// if(gp.getDirection()==100) { -// //x - tracking frame --> beamline direction -// dfdpTRACK[0][0] = 0; //df_x/d(alpha) -// dfdpTRACK[1][0] = -1; //df_y/d(alpha) -// dfdpTRACK[2][0] = 1; //df_z/d(alpha) -// } -// else if(gp.getDirection()==200) { -// //y - tracking frame --> almost unmeasured coordinate direction -// dfdpTRACK[0][0] = 1; //df_x/d(beta) -// dfdpTRACK[1][0] = 0; //df_y/d(beta) -// dfdpTRACK[2][0] = -1; //df_z/d(beta) -// } -// else { -// // Rotation around z -// dfdpTRACK[0][0] = -1; //df_x/d(gamma) -// dfdpTRACK[1][0] = 1; //df_y/d(gamma) -// dfdpTRACK[2][0] = 0; //df_z/d(gamma) -// -// } -// -// }//type -// }//layer -// }//side -// -// if(useGL) { -// -// //Put it into a matrix to be able to transform easily -// BasicMatrix _dfdpTRACK = FillMatrix(dfdpTRACK, 3, 1); -// //Get transformation matrix from tracking frame to sensor frame where the residuals are calculated -// Hep3Matrix trkToStrip = getTrackToStripRotation(strip); -// //Transform derivatives to sensor frame! -// BasicMatrix dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, _dfdpTRACK); -// //Add it to the global parameter object -// gp.setDfDp(dfdp); -// _glp.add(gp); -// if (_DEBUG) { -// 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()); -// } -// } else { -// if (_DEBUG) { -// System.out.printf("Global parameters %s name %s didn't affect this strip on side %d and layer %d\n",gp.getLabel(),gp.getName(),side,layer); -// } -// } -// } //gps -// -// -// - -
}
+ + + + +
private void CalculateResidual(HelicalTrackStrip strip, double msdrdphi, double msdz) { Hep3Vector u = strip.u();
@@ -693,7 +656,8 @@
Hep3Vector corigin = strip.origin(); double phi0 = _trk.phi0(); double R = _trk.R();
- double xint = strip.origin().x();
+ //double xint = strip.origin().x();
+ double xint = trackUtil.calculateHelixInterceptXPlane(_trk, strip);
double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0); double phi = -s/R + phi0; Hep3Vector trkpos = HelixUtils.PointOnHelix(_trk, s);
@@ -703,7 +667,7 @@
Hep3Vector mserr = new BasicHep3Vector(msdrdphi * Math.sin(phi), msdrdphi * Math.sin(phi), msdz); Hep3Vector vdiffTrk = VecOp.sub(trkpos, corigin);
- Hep3Matrix trkToStrip = getTrackToStripRotation(strip);
+ Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip);
Hep3Vector vdiff = VecOp.mult(trkToStrip, vdiffTrk); //Hep3Vector mserrrot = VecOp.mult(trkToStrip, mserr); double umc = vdiff.x();
@@ -724,6 +688,7 @@
_resid[2] = wmeas - wmc; _error[2] = wError; if (_DEBUG) {
+ System.out.println("---- " + this.getClass().getSimpleName() + " CalculateResidual ----");
System.out.println("Strip Origin: "); System.out.println(corigin.toString()); System.out.println("Position on the track (tracking coordinates) at the strip:");
@@ -746,6 +711,271 @@
} }
+ + + private void CalculateResidualSim(HelicalTrackStrip strip, double msdrdphi, double msdz, List<SimTrackerHit> simTrackerHits) { + + + Hep3Vector u = strip.u(); + Hep3Vector v = strip.v(); + Hep3Vector w = strip.w(); + Hep3Vector corigin = strip.origin(); + 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, smax, _nlc).get(0); + double phi = -s/R + phi0; + Hep3Vector trkpos = HelixUtils.PointOnHelix(_trk, s); + + //System.out.println("trkpos = "+trkpos.toString()); + //System.out.println("origin = "+corigin.toString()); + + 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); + //Hep3Vector mserrrot = VecOp.mult(trkToStrip, mserr); + 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; +// //System.out.println("strip error="+uError+"; ms error ="+msuError); +// _resid[0] = umeas - umc; +// _error[0] = Math.sqrt(uError * uError + msuError * msuError); +// _resid[1] = vmeas - vmc; +// _error[1] = vError; +// _resid[2] = wmeas - wmc; +// _error[2] = wError; + if (_DEBUG) { + System.out.println("---- " + this.getClass().getSimpleName() + " CalculateResidualSim ----"); + System.out.println("Strip Origin: "); + System.out.println(corigin.toString()); + System.out.println("Position on the track (tracking coordinates) at the strip:"); + System.out.println(trkpos.toString()); + System.out.println("Difference between the track position and strip origin (tracking coordinates) at the strip:"); + System.out.println("vdiffTrk :"); + System.out.println(vdiffTrk.toString()); + System.out.println("Difference between the track position and strip origin (\"strip\" coordinates) at the strip:"); + System.out.println("vdiff :"); + System.out.println(vdiff.toString()); + System.out.println("u :"); + System.out.println(u.toString()); + System.out.println("umeas = " + umeas + "; umc = " + umc + " ( prediction)"); + System.out.println("udiff = " + _resid[0] + " +/- " + _error[0] + " ( uError=" + uError + ", msuError=" + msuError + ")"); + System.out.println("MS: drdphi=" + msdrdphi + ", dz=" + msdz); + System.out.println("MS: phi=" + phi + " => msvec=" + mserr.toString());[truncated at 1000 lines; 998 more skipped]
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