Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
RunMPAlignment.java | +80 | -118 | 1.8 -> 1.9 |
MPAlignmentParameters.java | +277 | -567 | 1.13 -> 1.14 |
+357 | -685 |
Restructures initialization of parameter class. Fixed global derivatives jacobian.
diff -u -r1.8 -r1.9 --- RunMPAlignment.java 12 Oct 2012 05:59:36 -0000 1.8 +++ RunMPAlignment.java 23 Oct 2012 01:31:17 -0000 1.9 @@ -25,8 +25,6 @@
private AIDA aida = AIDA.defaultInstance(); String[] detNames = {"Tracker"};
- Integer _minLayers = 8; - Integer[] nlayers = {8};
int nevt = 0; double[] beamsize = {0.001, 0.02, 0.02}; String _config = "";
@@ -43,16 +41,17 @@
// this is set by the _config variable (detType in HeavyPhotonDriver) int flipSign = 1; private boolean _debug = false;
- double _resLayer1MinY; -
+ private String _type = "LOCAL"; //GLOBAL OR LOCAL RESIDUALS
- String simTrackerHitCollectionName = "TrackerHits"; -
+ private String simTrackerHitCollectionName = "TrackerHits"; + private String outputMilleFile = "alignMP.txt";
public void setDebug(boolean v) { this._debug = v; }
-
+ public void setMilleFile(String filename) { + outputMilleFile = filename; + }
public void setOutputPlotFileName(String filename) { outputPlotFileName = filename; }
@@ -61,60 +60,35 @@
hideFrame = hide; }
- public RunMPAlignment() { - nlayers[0] = 10; - _minLayers = 8; -//// if (_config.contains("Test")) -//// flipSign = -1; - - -
+ public void setResidualLimitFileName(String fileName) { + this._resLimitFileName = fileName;
}
- - public RunMPAlignment(int trackerLayers, int mintrkLayers, String config) { - nlayers[0] = trackerLayers; - _minLayers = mintrkLayers; - _config = config; - if (_config.contains("Test")) - flipSign = -1; -
+ + public RunMPAlignment() {
} @Override public void detectorChanged(Detector detector) {
- //ap = new MPAlignmentParameters("/Users/phansson/work/HPS/software/reco/run/alignMP.txt",_debug,hideFrame); - String outputMilleTxtFile = "alignMP.txt"; - Hep3Vector IP = new BasicHep3Vector(0., 0., 1.); - Hep3Vector _bfield = new BasicHep3Vector(0,0,detector.getFieldMap().getField(IP).y()); - - ap = new MPAlignmentParameters(outputMilleTxtFile,_bfield,_debug,hideFrame); - - //Initialize the res limits - for(int i=1;i<=10;++i) { - for(int j=0;j<3;++j) { - double xmin = -50; - double xmax = 50; - for(int side=0;side<2;++side) { - ap.setResLimits(side,i,j, xmin, xmax); - } - } - }
+ ap = new MPAlignmentParameters(this.outputMilleFile); + ap._DEBUG = _debug; + ap.setType(this._type); + ap.setHideFrame(hideFrame); + double bfield = detector.getFieldMap().getField(new BasicHep3Vector(0., 0., 1.)).y(); + System.out.printf("%s: B-field in z %.3f\n",this.getClass().getSimpleName(),bfield); + ap.setBfield(bfield);
loadResidualLimits();
- -
} @Override
- public void process( - EventHeader event) { -
+ public void process(EventHeader event) {
+
List<Track> tracklist = null; if(event.hasCollection(Track.class,"MatchedTracks")) { tracklist = event.get(Track.class, "MatchedTracks");
@@ -123,33 +97,6 @@
} }
- 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(this.getClass().getSimpleName() + ": 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(this.getClass().getSimpleName() + ": Event has " + trackerHits.size() + " SiTrackerHit"); - } - - -
for (Track trk : tracklist) { //if(trk.getCharge()>0) continue;
@@ -159,23 +106,28 @@
if(hitsOnBothSides(trk)) continue;
- totalTracksProcessed++;
+ double Px = trk.getTrackStates().get(0).getMomentum()[0]; + if(Px < 0.1) { + System.out.printf("%s: Trk p = [%.3f,%.3f,%.3f] is low skip!?\n",this.getClass().getSimpleName(), + trk.getTrackStates().get(0).getMomentum()[0],trk.getTrackStates().get(0).getMomentum()[1],trk.getTrackStates().get(0).getMomentum()[2]); + continue; + }
+ totalTracksProcessed++;
- ap.PrintResidualsAndDerivatives(trk,totalTracks,"LOCAL",simTrackerHits);
+ ap.PrintResidualsAndDerivatives(trk, totalTracks);
}
- - ap.CalculateResidualSimHits(strips,trackerHits,simTrackerHits);
}
+ @Override
public void endOfData() { ap.updatePlots(); try {
@@ -185,11 +137,12 @@
} 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);
+ if (!"".equals(outputPlotFileName)) { + try { + aida.saveAs(outputPlotFileName); + } catch (IOException ex) { + Logger.getLogger(TrigRateDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex); + }
} }
@@ -220,49 +173,58 @@
private void loadResidualLimits() {
- FileReader fReader = null; - BufferedReader bReader = null; - try { - fReader = new FileReader(this._resLimitFileName); - bReader = new BufferedReader(fReader); - - String line; - while( (line = bReader.readLine()) != null) { - if (line.contains("#")) continue; - String[] vec = line.split("\\s+"); - if(vec.length!=5) { - System.out.println(this.getClass().getSimpleName() + ": Error: residual limits line has wrong format -> " + line); - System.exit(1); - } - try { - int side = Integer.parseInt(vec[0]); - int layer = Integer.parseInt(vec[1]); - int direction = Integer.parseInt(vec[2]); - double min = Double.parseDouble(vec[3]); - double max = Double.parseDouble(vec[4]); - ap.setResLimits(side, layer, direction, min, max); - } catch(NumberFormatException e) { - Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE,null,e); - } - } - bReader.close(); - fReader.close(); - } - catch (FileNotFoundException ex) { - Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE, null, ex); - } catch (IOException e) { - Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE,null,e); - }
+ //Initialize the res limits + for(int i=1;i<=10;++i) { + for(int j=0;j<3;++j) { + double xmin = -50; + double xmax = 50; + for(int side=0;side<2;++side) { + ap.setResLimits(side,i,j, xmin, xmax); + } + } + }
+ if(!"".equals(this._resLimitFileName)) { + FileReader fReader; + BufferedReader bReader; + try { + fReader = new FileReader(this._resLimitFileName); + bReader = new BufferedReader(fReader); + + String line; + while( (line = bReader.readLine()) != null) { + if (line.contains("#")) continue; + String[] vec = line.split("\\s+"); + if(vec.length!=5) { + System.out.println(this.getClass().getSimpleName() + ": Error: residual limits line has wrong format -> " + line); + System.exit(1); + } + try { + int side = Integer.parseInt(vec[0]); + int layer = Integer.parseInt(vec[1]); + int direction = Integer.parseInt(vec[2]); + double min = Double.parseDouble(vec[3]); + double max = Double.parseDouble(vec[4]); + ap.setResLimits(side, layer, direction, min, max); + } catch(NumberFormatException e) { + Logger.getLogger(RunMPAlignment.class.getName()).log(Level.SEVERE,null,e); + } + } + bReader.close(); + fReader.close(); + } + catch (FileNotFoundException ex) { + Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE, null, ex); + } catch (IOException e) { + Logger.getLogger(TrigRateAna.class.getName()).log(Level.SEVERE,null,e); + } + }
- ap.setResLimits(0,1,0, -0.3, 0.3);
+
}
- public void setResidualLimitFileName(String fileName) { - this._resLimitFileName = fileName; - }
diff -u -r1.13 -r1.14 --- MPAlignmentParameters.java 17 Oct 2012 23:00:38 -0000 1.13 +++ MPAlignmentParameters.java 23 Oct 2012 01:31:17 -0000 1.14 @@ -2,34 +2,29 @@
import hep.aida.*; import hep.aida.ref.plotter.PlotterRegion;
-import org.lcsim.hps.users.mgraham.alignment.*;
import hep.physics.matrix.BasicMatrix; import hep.physics.matrix.MatrixOp;
-import hep.physics.vec.BasicHep3Matrix; -import hep.physics.vec.BasicHep3Vector; -import hep.physics.vec.Hep3Matrix; -import hep.physics.vec.Hep3Vector; -import hep.physics.vec.VecOp;
+import hep.physics.vec.*;
import java.io.FileWriter; import java.io.IOException; import java.io.PrintWriter;
-import java.util.*;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Map;
import java.util.logging.Level; import java.util.logging.Logger;
-import org.lcsim.detector.IDetectorElement; -import org.lcsim.detector.ITransform3D; -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.*;
+import org.lcsim.event.RawTrackerHit; +import org.lcsim.event.SimTrackerHit; +import org.lcsim.event.Track; +import org.lcsim.event.TrackerHit;
import org.lcsim.fit.helicaltrack.*; import org.lcsim.hps.event.HPSTransformations; import org.lcsim.hps.monitoring.AIDAFrame; import org.lcsim.hps.recon.tracking.SvtUtils; import org.lcsim.hps.recon.tracking.TrackUtils; import org.lcsim.hps.recon.tracking.TrackerHitUtils;
-import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit; -import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.hps.users.mgraham.alignment.RunAlignment;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate; import org.lcsim.recon.tracking.seedtracker.SeedTrack; import org.lcsim.util.aida.AIDA;
@@ -56,30 +51,27 @@
private int _nlc = 5; //the five track parameters private int _ngl = 1; //delta(u) and delta(gamma) for each plane
- //private GlobalParameters _globalParameterList;
private List<GlobalParameter> _glp = new ArrayList<GlobalParameter>(); private BasicMatrix _dfdq;
- //private BasicMatrix _dfdp;
private HelicalTrackFit _trk; private double[] _resid = new double[3]; private double[] _error = new double[3]; FileWriter fWriter; PrintWriter pWriter;
- private HPSTransformations _detToTrk;
+ private String _type = "LOCAL";
boolean _DEBUG=false; double smax = 1e3;
+ Hep3Vector _bfield = null;
ResLimit _resLimits = new ResLimit();
+ private HPSTransformations _detToTrk;
AlignmentUtils _alignUtils; TrackUtils trackUtil; TrackerHitUtils trackerHitUtil;
- Hep3Vector _bfield = null;
WTrackUtils wutils = new WTrackUtils();
- boolean _doSimHitsResiduals = false; - boolean _includeMS = false;
+ private boolean _includeMS = true;
private AIDAFrame plotterFrame;
- private AIDAFrame plotterFrameSimHit;
private AIDAFrame plotterFrameSummary; private boolean hideFrame = false; private AIDA aida = AIDA.defaultInstance();
@@ -99,27 +91,21 @@
- public MPAlignmentParameters(String outfile, Hep3Vector bfield_vec, boolean debug, boolean hidePlots) { - _DEBUG = debug;
+ public MPAlignmentParameters(String outfile) { +
_detToTrk = new HPSTransformations(); _alignUtils = new AlignmentUtils(_DEBUG); trackUtil = new TrackUtils(); trackerHitUtil = new TrackerHitUtils(_DEBUG);
- hideFrame = hidePlots; - _bfield = bfield_vec; - this.wutils.setDebug(debug);
+ _bfield = new BasicHep3Vector(0., 0., 1.); + this.wutils.setDebug(this._DEBUG);
try {
-//open things up
fWriter = new FileWriter(outfile); pWriter = new PrintWriter(fWriter); } catch (IOException ex) { Logger.getLogger(RunAlignment.class.getName()).log(Level.SEVERE, null, ex); }
- //_globalParameterList = new GlobalParameters(); - //if(_DEBUG) _globalParameterList.print(); - -
makeAlignmentPlots();
@@ -136,7 +122,29 @@
_resLimits.add(s,l,d, low,high); }
- public void PrintResidualsAndDerivatives(Track track, int itrack, String type, List<SimTrackerHit> simhits) {
+ public void setHideFrame(boolean hide) { + this.hideFrame = hide; + } + + public void setDebug(boolean debug) { + this._DEBUG = debug; + } + + public void setBfield(double bfield) { + _bfield = new BasicHep3Vector(0., 0., bfield); + } + + public void setType(String t) { + this._type = t; + } + + public void setIncludeMS(boolean include) { + this._includeMS = include; + } + + + + public void PrintResidualsAndDerivatives(Track track, int itrack) {
SeedTrack st = (SeedTrack) track; SeedCandidate seed = st.getSeedCandidate();
@@ -145,8 +153,8 @@
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()); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": Loop over " + hitsOnTrack.size() + " hitsOnTrack");
+ aida.cloud1D("Track Chi2 "+ half).fill(track.getChi2()); + if(_DEBUG) System.out.printf("%s: track %d has %d hits\n",this.getClass().getSimpleName(),itrack,hitsOnTrack.size());
for (TrackerHit hit : hitsOnTrack) { HelicalTrackHit htc = (HelicalTrackHit) hit;
@@ -164,36 +172,29 @@
List<HelicalTrackStrip> clusterlist = cross.getStrips();
- - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": Loop over " + clusterlist.size() + " clusterlist for this hitontrack"); - - for (HelicalTrackStrip cl : clusterlist) {
- if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": cluster size " + cl.rawhits().size()); - - if(type=="GLOBAL") { - - //CalculateLocalDerivativesGLOBAL(cl); - CalculateGlobalDerivatives(cl); - CalculateResidualGLOBAL(cl, msdrphi, msdz);
+ if(_DEBUG) { + System.out.printf("%s: This hit has %d clusters msdrphi=%.4f msdz=%.4f\n",this.getClass().getSimpleName(),clusterlist.size(),msdrphi,msdz); + }
- } else {
+ for (HelicalTrackStrip cl : clusterlist) { + if(!"GLOBAL".equals(_type)) {
CalculateResidual(cl, msdrphi, msdz);
- - if(_doSimHitsResiduals) CalculateResidualSim(cl, msdrphi, msdz,simhits); -
CalculateLocalDerivatives(cl); CalculateGlobalDerivatives(cl);
- PrintStripResidualsNew(cl);
+ PrintStripResiduals(cl); + } else { + //CalculateLocalDerivativesGLOBAL(cl); + CalculateGlobalDerivatives(cl); + CalculateResidualGLOBAL(cl, msdrphi, msdz);
} } } if(itrack%50==0) this.updatePlots();
- //AddTarget(0.1, 0.02);
}
- private BasicMatrix CalculateLocalDerivativesOld(HelicalTrackStrip strip) {
+ private BasicMatrix CalculateLocalDerivativesOld(double xint) {
//get track parameters. double d0 = _trk.dca(); double z0 = _trk.z0();
@@ -201,7 +202,6 @@
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[][] dfdq = new double[3][5];
@@ -224,12 +224,10 @@
} 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);
+ double xint = strip.origin().x(); + BasicMatrix dfdqGlobalOld = CalculateLocalDerivativesOld(xint); + BasicMatrix dfdqGlobal = _alignUtils.calculateLocalHelixDerivatives(_trk, xint); + Hep3Matrix trkToStrip = trackerHitUtil.getTrackToStripRotation(strip);
_dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal); if (_DEBUG) {
@@ -240,22 +238,17 @@
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.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"); - System.out.printf("%5.5f %5.5f\n", s, xint); - System.out.println("Compare local derivatives "); - System.out.println("dfdqGlobal:"); - System.out.println(dfdqGlobal.toString()); - System.out.println("dfdqGlobalOld:"); - System.out.println(dfdqGlobalOld.toString());
+ System.out.printf("%s: --- CalculateLocalDerivatives Result --- \n",this.getClass().getSimpleName()); + System.out.printf("%s: Strip Origin %s \n",this.getClass().getSimpleName(),strip.origin().toString()); + System.out.printf("%s: %10s%10s%10s%10s%10s%10s%10s\n",this.getClass().getSimpleName(),"d0","z0","slope","phi0","R", "xint", "s"); + System.out.printf("%s: %10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n",this.getClass().getSimpleName(), d0, z0, slope, phi0, R,xint,s); + System.out.printf("%s: Local derivatives:\n",this.getClass().getSimpleName()); + System.out.printf("%s: \n%s\n",this.getClass().getSimpleName(),dfdqGlobal.toString()); + System.out.printf("%s: OLD Local derivatives:\n",this.getClass().getSimpleName()); + System.out.printf("%s: \n%s\n",this.getClass().getSimpleName(),dfdqGlobalOld.toString());
} }
@@ -267,30 +260,38 @@
private void CalculateGlobalDerivatives(HelicalTrackStrip strip) {
+ //Find interception with plane that the strips belongs to + Hep3Vector trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z())); + + if(Double.isNaN(trkpos.x())) { + System.out.printf("%s: error this trkpos is wrong %s\n",this.getClass().getSimpleName(),trkpos.toString()); + System.out.printf("%s: origin %s trk \n%s\n",this.getClass().getSimpleName(),strip.origin(),_trk.toString()); + trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()),true); + System.exit(1); + }
double slope = _trk.slope();
- double xint = strip.origin().x();
+ double xint = trkpos.x();
double xr = 0.0; double yr = 0.0; double d0 = _trk.dca(); double phi0 = _trk.phi0(); double R = _trk.R(); double z0 = _trk.z0();
- double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
+ double s = HelixUtils.PathToXPlane(_trk, xint, 0, 0).get(0);
double phi = -s/R + phi0; double umeas = strip.umeas(); Hep3Vector corigin = strip.origin(); double vmeas = corigin.y(); //THis is wrong
-
+ int layer = strip.layer(); + int side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? 10000 : 20000;
if(_DEBUG) {
- System.out.println("Calculate global derivatives for this strip hit in layer " + strip.layer()); - System.out.printf(" %10s%10s%10s%10s%10s%10s%10s\n","d0","z0","slope","phi0","R","xint","phi"); - 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); -
+ System.out.printf("%s: --- Calculate global derivatives ---\n",this.getClass().getSimpleName()); + System.out.printf("%s: Side %d, layer %d, strip origin %s\n",this.getClass().getSimpleName(),side,layer,corigin.toString()); + System.out.printf("%s: %10s%10s%10s%10s%10s%10s%10s%10s%10s\n",this.getClass().getSimpleName(),"d0","z0","slope","phi0","R","xint","phi", "xint","s"); + System.out.printf("%s: %10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f%10.4f\n",this.getClass().getSimpleName(), d0, z0, slope, phi0, R,xint,phi,xint,s);
}
- if(1==1) - return;
//1st index = alignment parameter (only u so far) //2nd index = residual coordinate (on du so far)
@@ -316,18 +317,25 @@
//into the sensor frame because that is where the residuals are calculated
- //Which half of the detecotor? - int side = strip.origin().z()>0 ? 10000 : 20000; - int layer = strip.layer(); -
+
//Rotation matrix from the tracking fram to the sensor/strip frame
- Hep3Matrix T = this.trackerHitUtil.getTrackToStripRotation(strip);
+ Hep3Matrix T = 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);
+ if(_DEBUG) { + System.out.printf("%s: Rotation matrx from tracking/global to local frame T:\n %s\n",this.getClass().getSimpleName(),T.toString()); + System.out.printf("%s: Derivatives of the rotation matrix deltaR' w.r.t. rotation a,b,c around x,y,z axis:\n",this.getClass().getSimpleName()); + System.out.printf("%s: ddeltaRprime_da:\n %s\n",this.getClass().getSimpleName(),ddeltaRprime_da.toString()); + System.out.printf("%s: ddeltaRprime_db:\n %s\n",this.getClass().getSimpleName(),ddeltaRprime_db.toString()); + System.out.printf("%s: ddeltaRprime_dc:\n %s\n",this.getClass().getSimpleName(),ddeltaRprime_dc.toString()); + } + + +
//**************************************************************************** //Calculate jacobian da/db
@@ -340,14 +348,60 @@
Hep3Vector dq_dy = VecOp.mult(T, deltaY_gl); Hep3Vector dq_dz = VecOp.mult(T, deltaZ_gl);
+ if(_DEBUG) { + System.out.printf("%s: - Upper left 3x3 of Jacobian da/db: dq_dx,dq_dy,dq_dz\n",this.getClass().getSimpleName()); + System.out.printf("%s: dq_dx: %s\n",this.getClass().getSimpleName(),dq_dx.toString()); + System.out.printf("%s: dq_dy: %s\n",this.getClass().getSimpleName(),dq_dy.toString()); + System.out.printf("%s: dq_dz: %s\n",this.getClass().getSimpleName(),dq_dz.toString()); + } + + + + if(_DEBUG) { + System.out.printf("%s: - Upper right 3x3 of Jacobian da/db: dq_dx,dq_dy,dq_dz\n",this.getClass().getSimpleName()); + } + +
// 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);
+ Hep3Vector x_vec = trkpos; //position around where the small rotation happens + Hep3Vector x_vec_tmp = VecOp.mult(ddeltaRprime_da, x_vec); // derivative of the position + Hep3Vector x_vec_tmp2 = VecOp.sub(x_vec_tmp, x_vec); // subtract position + Hep3Vector dq_da = VecOp.mult(T,x_vec_tmp2); //rotated into local frame + + if(_DEBUG) { + System.out.printf("%s: position %s rotation derivative w.r.t. a ddeltaR'/da %s\n",this.getClass().getSimpleName(),x_vec.toString(),x_vec_tmp.toString()); + System.out.printf("%s: subtracted %s and rotated to local %s \n",this.getClass().getSimpleName(),x_vec_tmp2.toString(),dq_da.toString()); + } + +
x_vec_tmp = VecOp.mult(ddeltaRprime_db, x_vec);
- Hep3Vector dq_db = VecOp.mult(T,x_vec_tmp);
+ x_vec_tmp2 = VecOp.sub(x_vec_tmp, x_vec); + Hep3Vector dq_db = VecOp.mult(T,x_vec_tmp2); + + if(_DEBUG) { + System.out.printf("%s: position %s rotation derivative w.r.t. a ddeltaR'/db %s\n",this.getClass().getSimpleName(),x_vec.toString(),x_vec_tmp.toString()); + System.out.printf("%s: subtracted %s and rotated to local %s \n",this.getClass().getSimpleName(),x_vec_tmp2.toString(),dq_db.toString()); + } +
x_vec_tmp = VecOp.mult(ddeltaRprime_dc, x_vec);
- Hep3Vector dq_dc = VecOp.mult(T,x_vec_tmp);
+ x_vec_tmp2 = VecOp.sub(x_vec_tmp, x_vec); + Hep3Vector dq_dc = VecOp.mult(T,x_vec_tmp2); + + if(_DEBUG) { + System.out.printf("%s: position %s rotation derivative w.r.t. a ddeltaR'/dc %s\n",this.getClass().getSimpleName(),x_vec.toString(),x_vec_tmp.toString()); + System.out.printf("%s: subtracted %s and rotated to local %s \n",this.getClass().getSimpleName(),x_vec_tmp2.toString(),dq_dc.toString()); + } + + + if(_DEBUG) { + System.out.printf("%s: Summary:\n",this.getClass().getSimpleName()); + System.out.printf("%s: dq_da: %s\n",this.getClass().getSimpleName(),dq_da.toString()); + System.out.printf("%s: dq_db: %s\n",this.getClass().getSimpleName(),dq_db.toString()); + System.out.printf("%s: dq_dc: %s\n",this.getClass().getSimpleName(),dq_dc.toString()); + } + + +
// Lower left 3x3 //deltaR = (alpha,beta,gamma) the rotation alignment parameters in the local frame
@@ -355,7 +409,19 @@
Hep3Vector ddeltaR_dy = new BasicHep3Vector(0,0,0); Hep3Vector ddeltaR_dz = new BasicHep3Vector(0,0,0);
+ + if(_DEBUG) { + System.out.printf("%s: - Lower left 3x3 of Jacobian ddeltaR/d(x,y,z): dDeltaR_dx,dDeltaR_dy,dDeltaR_dz\n",this.getClass().getSimpleName()); + System.out.printf("%s: ddeltaR_dx: %s\n",this.getClass().getSimpleName(),ddeltaR_dx.toString()); + System.out.printf("%s: ddeltaR_dy: %s\n",this.getClass().getSimpleName(),ddeltaR_dy.toString()); + System.out.printf("%s: ddeltaR_dz: %s\n",this.getClass().getSimpleName(),ddeltaR_dz.toString()); + } + + + +
// 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;
@@ -374,7 +440,15 @@
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);
+*/ + + if(_DEBUG) { + System.out.printf("%s: - Lower right 3x3 of Jacobian ddeltaR/d(a,b,c): \n",this.getClass().getSimpleName()); + System.out.printf("%s: T: %s\n",this.getClass().getSimpleName(),T.toString()); + }
+ +
//Now fill the Jacobian 6x6 matrix BasicMatrix da_db = new BasicMatrix(6,6); //upper left 3x3
@@ -398,15 +472,23 @@
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());
+ da_db.setElement(3,3,T.e(0, 0)); + da_db.setElement(4,3,T.e(1, 0)); + da_db.setElement(5,3,T.e(2, 0)); + da_db.setElement(3,4,T.e(0, 1)); + da_db.setElement(4,4,T.e(1, 1)); + da_db.setElement(5,4,T.e(2, 1)); + da_db.setElement(3,5,T.e(0, 2)); + da_db.setElement(4,5,T.e(1, 2)); + da_db.setElement(5,5,T.e(2, 2)); +// 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());
@@ -419,18 +501,30 @@
da_db.setElement(5,2,ddeltaR_dz.z()); if(_DEBUG) {
- System.out.println("da_db: \n" + da_db.toString());
+ System.out.printf("%s: da_db:\n%s \n",this.getClass().getSimpleName(),da_db.toString()); + System.out.printf("%s: det(da_db) = %.3f \n",this.getClass().getSimpleName(),da_db.det());
}
- BasicMatrix db_da = da_db; - db_da.invert();
+ BasicMatrix db_da = (BasicMatrix) MatrixOp.inverse(da_db); + //db_da.invert();
if(_DEBUG) {
- System.out.println("db_da: \n" + db_da.toString());
+ System.out.printf("%s: db_da: \n%s \n",this.getClass().getSimpleName(),db_da.toString()); + BasicMatrix prod = (BasicMatrix) MatrixOp.mult(db_da,da_db); + System.out.printf("%s: db_da*da_db: \n%s \n",this.getClass().getSimpleName(),prod.toString());
}
+ + + if(1==1) + return; + + + + +
//**************************************************************************** // Calcualte the global derivatives in the global frame dz/db // where z is the residual in the global frame and b are the alignment
@@ -664,51 +758,54 @@
private void CalculateResidual(HelicalTrackStrip strip, double msdrdphi, double msdz) {
- if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": CalculateResidual");
+ if(_DEBUG) System.out.printf("%s: --- CalculateResidual ---\n",this.getClass().getSimpleName());
Hep3Vector u = strip.u(); Hep3Vector v = strip.v(); Hep3Vector w = strip.w(); Hep3Vector eta = VecOp.cross(u,v); Hep3Vector corigin = strip.origin();
+ String side = SvtUtils.getInstance().isTopLayer((SiSensor)((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom"; + + double bfield = Math.abs(this._bfield.z()); + if(_DEBUG) { + System.out.printf("%s: Finding interception point for residual calculation (B=%s)\n",this.getClass().getSimpleName(),this._bfield.toString()); + double xint_wrong = trackUtil.calculateHelixInterceptXPlane(_trk, strip); + double s_wrong = HelixUtils.PathToXPlane(_trk, xint_wrong, smax, _nlc).get(0); + Hep3Vector trkpos_wrong= HelixUtils.PointOnHelix(_trk, s_wrong); + System.out.printf("%s: using wrong method xint_wrong=%.3f s_wrong=%.3f -> trkpos_wrong=%s\n",this.getClass().getSimpleName(),xint_wrong,s_wrong,trkpos_wrong.toString()); + } + //Find interception with plane that the strips belongs to + Hep3Vector trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, bfield); + + if(_DEBUG) { + System.out.printf("%s: found interception point at %s \n",this.getClass().getSimpleName(),trkpos.toString()); + } + + + if(Double.isNaN(trkpos.x()) || Double.isNaN(trkpos.y()) || Double.isNaN(trkpos.z())) { + System.out.printf("%s: failed to get interception point (%s) \n",this.getClass().getSimpleName(),trkpos.toString()); + System.out.printf("%s: track params\n%s\n",this.getClass().getSimpleName(),_trk.toString()); + System.out.printf("%s: track pT=%.3f chi2=[%.3f][%.3f] \n",this.getClass().getSimpleName(),_trk.pT(bfield),_trk.chisq()[0],_trk.chisq()[1]); + trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, bfield,true); + System.exit(1); + }
- if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": Finding interception point for residual calculation");
- WTrack wtrack = new WTrack(_trk,Math.abs(_bfield.z()),true); //flip the charge since I'm looking at fitted tracks with B negative?
- Hep3Vector trkpos = wutils.getHelixAndPlaneIntercept(wtrack, strip.origin(), eta, new BasicHep3Vector(0,0,1));
double xint = trkpos.x(); double phi0 = _trk.phi0(); double R = _trk.R();
- //double xint = strip.origin().x();
double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0); double phi = -s/R + phi0;
- if(_DEBUG) { - double xint_wrong = trackUtil.calculateHelixInterceptXPlane(_trk, strip); - double s_wrong = HelixUtils.PathToXPlane(_trk, xint_wrong, smax, _nlc).get(0); - Hep3Vector trkpos_wrong= HelixUtils.PointOnHelix(_trk, s_wrong); - System.out.println(this.getClass().getSimpleName() + ": trkposdiff " + VecOp.sub(trkpos, trkpos_wrong).toString() + " from trkpos_cor " + trkpos.toString() + " trkpos_wrong " + trkpos_wrong.toString()); - }
- //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); - //int idiffbin = 0;
- //double idiffbin = (vdiffTrk.y())/10; -// if(idiffbin<-2) { -// idiffbin = -2; -// System.out.println("WARNING vdiffTrk.y() = " + vdiffTrk.y() + " merge to -> idiffbin= " + idiffbin); -// } -// if(idiffbin>2) { -// idiffbin = 2; -// System.out.println("WARNING vdiffTrk.y() = " + vdiffTrk.y() + " merge to -> idiffbin= " + idiffbin); -// }
double umc = vdiff.x();
@@ -721,288 +818,38 @@
double vError = (strip.vmax() - strip.vmin()) / Math.sqrt(12); double wmeas = 0; double wError = 10.0/Math.sqrt(12); //0.001;
- //if(idiffbin==0) { - _resid[0] = umeas - umc; - _error[0] = this._includeMS ? Math.sqrt(uError * uError + msuError * msuError) : uError; - _resid[1] = vmeas - vmc; - _error[1] = vError; - _resid[2] = wmeas - wmc; - _error[2] = wError; - //if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": idiffbin= " + idiffbin + " resid[0]=" + _resid[0]); - //} else { - // _resid[0] = 9999999.9; - //} - - - //Calcualte the distance from the center of the sensor - String side = corigin.z()>0. ? "top" : "bottom"; - //if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": idiffbin= " + idiffbin + " filling ydiff bist with " + _resid[0]); - //aida.histogram1D("res_u_vs_ydiff_"+idiffbin+"_layer_" + strip.layer() + "_" + side).fill(_resid[0]); - aida.profile1D("res_u_vs_ydiff_layer_" + strip.layer() + "_" + side).fill(vdiffTrk.y(),_resid[0]); - - 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:"); - 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()); - System.out.println("MS: msuError = " + msuError + " (msvec*u = " + mserr.toString() + " * " + u.toString() + ")" ); - //System.out.println("MS: msvec rotated to strip? = " + mserrrot.toString()); - - }
+ _resid[0] = umeas - umc; + _error[0] = this._includeMS ? Math.sqrt(uError * uError + msuError * msuError) : uError; + _resid[1] = vmeas - vmc; + _error[1] = vError; + _resid[2] = wmeas - wmc; + _error[2] = wError;
- } - - - 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()); - System.out.println("MS: msuError = " + msuError + " (msvec*u = " + mserr.toString() + " * " + u.toString() + ")" ); - - } - - if(simTrackerHits!=null && simTrackerHits.size()>0) { - List<SimTrackerHit> simHits = trackerHitUtil.stripClusterToSimHits(strip, simTrackerHits,true); - if(simHits.size()>0) { - //Should only be one! - if(simHits.size()>1) { - System.out.println(this.getClass().getSimpleName() + " ERROR multiple simhits for this strip"); - System.exit(1); - } - if(_DEBUG) System.out.println( simHits.size() + " sim hits for this raw hit at: "); - SimTrackerHit simHit =simHits.get(0); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in JLab frame" + simHit.getPositionVec().toString()); - //Transformation between the JLAB and tracking coordinate systems - Hep3Matrix detToTrackMatrix = trackerHitUtil.detToTrackRotationMatrix(); - Hep3Vector simHitPos = VecOp.mult(detToTrackMatrix, simHit.getPositionVec()); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in tracking frame" + simHitPos.toString()); - - String side = strip.origin().z()>0 ? "top" : "bottom"; - Hep3Vector stripPosition = trackerHitUtil.getClusterPosition(strip,true); - aida.histogram1D("simhitstripres_z_layer" + strip.layer() + "_" + side).fill(simHitPos.z()-stripPosition.z()); - aida.histogram1D("simhitstripres_y_layer" + strip.layer() + "_" + side).fill(simHitPos.y()-stripPosition.y()); - aida.histogram1D("simhitstripres_x_layer" + strip.layer() + "_" + side).fill(simHitPos.x()-stripPosition.x()); - - //Transform the sim hit to do a u residual - Hep3Vector vDiffSimTrk = VecOp.sub(simHitPos, corigin); - //Rotate this into the local sensor frame - Hep3Vector vdiffSim = VecOp.mult(trkToStrip, vDiffSimTrk); - double umc_sim = vdiffSim.x(); - double vmc_sim = vdiffSim.y(); - double wmc_sim = vdiffSim.z(); - aida.histogram1D("simhitstripres_u_layer" + strip.layer() + "_" + side).fill(umeas - umc_sim); - aida.histogram1D("simhitstripres_v_layer" + strip.layer() + "_" + side).fill(vmeas - vmc_sim); - aida.histogram1D("simhitstripres_w_layer" + strip.layer() + "_" + side).fill(wmeas - wmc_sim); - - } - - } else { - if(_DEBUG) System.out.println("No simTrackerHit collection"); - } - - -
+ //Fill residuals vs distrance from center of plane in the v directions + aida.profile1D("res_u_vs_ydiff_layer_" + strip.layer() + "_" + side).fill(vdiffTrk.y(),_resid[0]);
- } - - - - public void CalculateResidualSimHits(List<HelicalTrackStrip> strips, List<SiTrackerHit> strip1dhits, List<SimTrackerHit> simTrackerHits) { -// SeedTrack st = (SeedTrack) track; -// SeedCandidate seed = st.getSeedCandidate(); -// Map<HelicalTrackHit, MultipleScatter> msmap = seed.getMSMap(); -// HelicalTrackFit trk = seed.getHelix(); -// double msdrdphi = 0; - double msdz = 0; - - if (_DEBUG) { - System.out.println("---- " + this.getClass().getSimpleName() + " CalculateResidualSimHits ----"); - } - if(strips==null) { - if (_DEBUG) { - System.out.println(this.getClass().getSimpleName() + ": list of strips is null"); - } - return; - } -
if (_DEBUG) {
- System.out.println(this.getClass().getSimpleName() + ": " + strips.size() + " strips");
+ System.out.printf("%s: ---- CalculateResidual ----\n",this.getClass().getSimpleName()); + System.out.printf("%s: Strip Origin: %s\n",this.getClass(),corigin.toString()); + System.out.printf("%s: Position on the track (tracking coordinates) at the strip: %s\n",this.getClass(),trkpos.toString()); + System.out.printf("%s: vdiffTrk %s\n",this.getClass(),vdiffTrk.toString()); + System.out.printf("%s: vdiff %s\n",this.getClass(),vdiff.toString()); + System.out.printf("%s: u %s\n",this.getClass(),u.toString()); + System.out.printf("%s: umeas = %.4f umc = %.4f\n",this.getClass(),umeas,umc); + System.out.printf("%s: udiff = %.4f +/- %.4f ( uError=%.4f , msuError=%.4f\n",this.getClass(),_resid[0],_error[0],uError,msuError); + System.out.printf("%s: MS: drdphi=%.4f, dz=%.4f\n",this.getClass(),msdrdphi,msdz); + System.out.printf("%s: MS: phi=%.4f => msvec=%s\n",this.getClass(),phi,mserr.toString()); + System.out.printf("%s: MS: msuError = %.4f (msvec*u = %s * %s\n",this.getClass(),msuError,mserr.toString(),u.toString());
}
- for(HelicalTrackStrip strip : strips) { - //Rotate into tracking frame - Hep3Matrix detToTrackMatrix = this.trackerHitUtil.detToTrackRotationMatrix(); - - Hep3Vector u = VecOp.mult(detToTrackMatrix,strip.u()); - Hep3Vector v = VecOp.mult(detToTrackMatrix,strip.v()); - Hep3Vector w = VecOp.mult(detToTrackMatrix,strip.w()); - Hep3Vector corigin = VecOp.mult(detToTrackMatrix,strip.origin()); - double umeas = strip.umeas(); - double vmeas = 0; - double wmeas = 0; - if (_DEBUG) { - System.out.println("----"); - System.out.println("Strip Origin: "); - System.out.println(corigin.toString()); - System.out.println("SIMHITS"); - } - if(simTrackerHits!=null && simTrackerHits.size()>0) { - List<SimTrackerHit> simHits = this.trackerHitUtil.stripClusterToSimHits(strip, simTrackerHits,false); - if(simHits.size()>0) { - //Should only be one! - if(simHits.size()>1) { - System.out.println(this.getClass().getSimpleName() + " ERROR multiple simhits for this strip"); - System.exit(1); - } - - if(_DEBUG) System.out.println( simHits.size() + " sim hits for this raw hit at: "); - SimTrackerHit simHit =simHits.get(0); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in JLab frame" + simHit.getPositionVec().toString()); - //Transformation between the JLAB and tracking coordinate systems - Hep3Vector simHitPos = VecOp.mult(detToTrackMatrix, simHit.getPositionVec()); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in tracking frame" + simHitPos.toString()); - - String side = corigin.z()>0 ? "top" : "bottom"; - Hep3Vector stripPosition = this.trackerHitUtil.getClusterPosition(strip,false); - aida.histogram1D("simhitallstripres_z_layer" + strip.layer() + "_" + side).fill(simHitPos.z()-stripPosition.z()); - aida.histogram1D("simhitallstripres_y_layer" + strip.layer() + "_" + side).fill(simHitPos.y()-stripPosition.y()); - aida.histogram1D("simhitallstripres_x_layer" + strip.layer() + "_" + side).fill(simHitPos.x()-stripPosition.x()); - - // //Transform the sim hit to do a u residual - // Hep3Vector vDiffSimTrk = VecOp.sub(simHitPos, corigin); - // //Rotate this into the local sensor frame - // Hep3Vector vdiffSim = VecOp.mult(trkToStrip, vDiffSimTrk); - // double umc_sim = vdiffSim.x(); - // double vmc_sim = vdiffSim.y(); - // double wmc_sim = vdiffSim.z(); - // aida.histogram1D("simallhitstripres_u_layer" + strip.layer() + "_" + side).fill(umeas - umc_sim); - // aida.histogram1D("simallhitstripres_v_layer" + strip.layer() + "_" + side).fill(vmeas - vmc_sim); - // aida.histogram1D("simallhitstripres_w_layer" + strip.layer() + "_" + side).fill(wmeas - wmc_sim); - } - - - } else { - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " No simTrackerHit collection"); - } - - /* - if(strip1dhits!=null && strip1dhits.size()>0) { - List<SiTrackerHit> siHits = this.trackerHitUtil.stripClusterToSiHits(strip, strip1dhits, false); - if(siHits.size()>0) { - //Should only be one! - if(siHits.size()>1) { - System.out.println(this.getClass().getSimpleName() + " ERROR multiple simhits for this strip"); - System.exit(1); - } - SiTrackerHitStrip1D h = (SiTrackerHitStrip1D) siHits.get(0); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " SiTrackerHitStrip1D position in JLab frame" + h.getPositionAsVector().toString()); - //Transformation between the JLAB and tracking coordinate systems - Hep3Vector siHitPos = VecOp.mult(detToTrackMatrix, h.getPositionAsVector()); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + " sim hit position in tracking frame" + siHitPos.toString()); - -// String side = corigin.z()>0 ? "top" : "bottom"; -// Hep3Vector stripPosition = this.getClusterPosition(strip,false); -// aida.histogram1D("simhitallstripres_z_layer" + strip.layer() + "_" + side).fill(simHitPos.z()-stripPosition.z()); -// aida.histogram1D("simhitallstripres_y_layer" + strip.layer() + "_" + side).fill(simHitPos.y()-stripPosition.y()); -// aida.histogram1D("simhitallstripres_x_layer" + strip.layer() + "_" + side).fill(simHitPos.x()-stripPosition.x()); - - // //Transform the sim hit to do a u residual - // Hep3Vector vDiffSimTrk = VecOp.sub(simHitPos, corigin); - // //Rotate this into the local sensor frame - // Hep3Vector vdiffSim = VecOp.mult(trkToStrip, vDiffSimTrk); - // double umc_sim = vdiffSim.x(); - // double vmc_sim = vdiffSim.y(); - // double wmc_sim = vdiffSim.z(); - // aida.histogram1D("simallhitstripres_u_layer" + strip.layer() + "_" + side).fill(umeas - umc_sim); - // aida.histogram1D("simallhitstripres_v_layer" + strip.layer() + "_" + side).fill(vmeas - vmc_sim); - // aida.histogram1D("simallhitstripres_w_layer" + strip.layer() + "_" + side).fill(wmeas - wmc_sim); - } else { - System.out.println(this.getClass().getSimpleName() + " ERROR no SiTrackerHitStrip1D for this strip"); - //System.exit(1); - - } - - - } else { - if(_DEBUG) System.out.println("No SiTrackerHitStrip1D matched to this strip"); - } - - */ - - } - -
}
+ +
@@ -1042,20 +889,12 @@
- private void PrintStripResidualsNew(HelicalTrackStrip strip) {
+ private void PrintStripResiduals(HelicalTrackStrip strip) {
if (_DEBUG) {
- System.out.println(this.getClass().getSimpleName() + ": PrintStripResidualsNew"); - int s = 1; - if(strip.origin().z()>0) s = 0; - System.out.printf("Strip Layer = %4d\n", strip.layer()); - System.out.printf("Residuals (u,v,w) : %5.5e %5.5e %5.5e (limits: [%.1e,%.1e] [%.1e,%.1e] [%.1e,%.1e])\n" - , _resid[0], _resid[1], _resid[2], - _resLimits.getMin(s,strip.layer(),0), - _resLimits.getMax(s,strip.layer(),0), - _resLimits.getMin(s,strip.layer(),1), - _resLimits.getMax(s,strip.layer(),1), - _resLimits.getMin(s,strip.layer(),2), - _resLimits.getMax(s,strip.layer(),2));
+ System.out.println(this.getClass().getSimpleName() + ": --- PrintStripResiduals ---"); + String side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom"; + System.out.printf("%s: Strip layer %4d in %s at origin %s\n",this.getClass().getSimpleName(), strip.layer(), side, strip.origin().toString()); + System.out.printf("Residuals (u,v,w) : %5.5e %5.5e %5.5e\n", _resid[0], _resid[1], _resid[2]);
System.out.printf("Errors (u,v,w) : %5.5e %5.5e %5.5e\n", _error[0], _error[1], _error[2]); String[] q = {"d0", "z0", "slope", "phi0", "R"}; System.out.println("track parameter derivatives");
@@ -1068,7 +907,7 @@
System.out.printf("%s %5.5e %5.5e %5.5e %5d %10s\n", "", gl.dfdp(0), gl.dfdp(1), gl.dfdp(2), gl.getLabel(),gl.getName()); //System.out.printf("%s %5.5e %5.5e %5.5e %5d\n", p[j], _dfdp.e(0, j), _dfdp.e(1, j), _dfdp.e(2, j), _globalLabel[j]); }
-
+ System.out.println(this.getClass().getSimpleName() + ": --- END PrintStripResiduals ---");
} pWriter.printf("%d\n", strip.layer()); // Loop over the three directions u,v,w and print residuals and derivatives
@@ -1257,50 +1096,46 @@
- aida.tree().cd("/");
+ aida.tree().cd("/");
plotterFrame = new AIDAFrame(); plotterFrame.setTitle("Residuals");
-
plotterFrameSummary = new AIDAFrame(); plotterFrameSummary.setTitle("Summary");
- plotterFrameSimHit = new AIDAFrame(); - plotterFrameSimHit.setTitle("Residuals SimHits"); -
String[] side = {"top","bottom"}; IPlotter plotter_chi2 = af.createPlotterFactory().create(); plotter_chi2.createRegions(2,1,0); plotter_chi2.setTitle("Track Chi2"); plotterFrame.addPlotter(plotter_chi2);
- IHistogram hchi2_top = aida.histogram1D("Track Chi2 top",50,0,15); - IHistogram hchi2_bot = aida.histogram1D("Track Chi2 bottom",50,0,3);
+ ICloud1D hchi2_top = aida.cloud1D("Track Chi2 top"); + ICloud1D hchi2_bot = aida.cloud1D("Track Chi2 bottom");
plotter_chi2.region(0).plot(hchi2_top); plotter_chi2.region(1).plot(hchi2_bot); String[] direction = {"u"}; double xbins_u_res[][] = {
- {-0.01,0.01}, - {-0.01,0.01}, - {-0.01,0.01}, - {-0.01,0.01}, - {-0.01,0.01}, - {-0.01,0.01}, - {-0.01,0.01}, - {-0.01,0.01}, - {-0.01,0.01}, - {-0.01,0.01} -// {-0.3,0.3}, -// {-0.7,0.7}, -// {-0.7,0.7}, -// {-0.5,0.5}, -// {-0.5,0.5}, -// {-0.5,0.5}, -// {-2,2}, -// {-2,2}, -// {-2,2}, -// {-2,2}
+// {-0.01,0.01}, +// {-0.01,0.01}, +// {-0.01,0.01}, +// {-0.01,0.01}, +// {-0.01,0.01}, +// {-0.01,0.01}, +// {-0.01,0.01}, +// {-0.01,0.01}, +// {-0.01,0.01}, +// {-0.01,0.01} + {-0.15,0.15}, + {-0.4,0.4}, + {-0.7,0.7}, + {-0.5,0.5}, + {-0.5,0.5}, + {-0.5,0.5},[truncated at 1000 lines; 150 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