hps-java/src/main/java/org/lcsim/hps/users/phansson
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;
}
hps-java/src/main/java/org/lcsim/hps/users/phansson
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]