Commit in hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl on MAIN | |||
GBLFileIO.java | +32 | -43 | 1.6 -> 1.7 |
GBLOutput.java | +67 | -33 | 1.9 -> 1.10 |
GBLOutputDriver.java | -1 | 1.4 -> 1.5 | |
+99 | -77 |
Added iter prop, fixed issues in prop in output and fixed prec issues
diff -u -r1.6 -r1.7 --- GBLFileIO.java 29 Aug 2013 21:10:43 -0000 1.6 +++ GBLFileIO.java 12 Sep 2013 17:57:38 -0000 1.7 @@ -21,6 +21,7 @@
/** * * @author phansson
+ *
*/ public class GBLFileIO {
@@ -63,11 +64,11 @@
} void printOldPerTrackParam(HelicalTrackFit htf) {
- addLine(String.format("Track perPar (R phi0 slope d0 z0) %f %f %f %f %f",htf.R(),htf.phi0(),htf.slope(),htf.dca(),htf.z0()));
+ addLine(String.format("Track perPar (R phi0 slope d0 z0) %.12f %.12f %.12f %.12f %.12f",htf.R(),htf.phi0(),htf.slope(),htf.dca(),htf.z0()));
} String getPerTrackParamStr(PerigeeParams perPar) {
- return String.format("Track perPar (kappa theta phi d0 z0) %f %f %f %f %f",perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
+ return String.format("Track perPar (R theta phi d0 z0) %.12f %.12f %.12f %.12f %.12f",1.0/perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
} void printPerTrackParam(PerigeeParams perPar) {
@@ -75,7 +76,7 @@
} String getPerTrackParamTruthStr(PerigeeParams perPar) {
- return String.format("Truth perPar (kappa theta phi d0 z0) %f %f %f %f %f",perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
+ return String.format("Truth perPar (kappa theta phi d0 z0) %.12f %.12f %.12f %.12f %.12f",perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
} void printPerTrackParamTruth(PerigeeParams perPar) {
@@ -83,7 +84,7 @@
} String getClTrackParamTruthStr(ClParams perPar) {
- return String.format("Truth clPar (q/p lambda phi xT yT) %f %f %f %f %f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
+ return String.format("Truth clPar (q/p lambda phi xT yT) %.12f %.12f %.12f %.12f %.12f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
} void printClTrackParamTruth(ClParams perPar) {
@@ -91,7 +92,7 @@
} String getClTrackParamStr(ClParams perPar) {
- return String.format("Track clPar (q/p lambda phi xT yT) %f %f %f %f %f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
+ return String.format("Track clPar (q/p lambda phi xT yT) %.12f %.12f %.12f %.12f %.12f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
} void printClTrackParam(ClParams perPar) { addLine(String.format("%s",this.getClTrackParamStr(perPar)));
@@ -100,24 +101,12 @@
void printNrHits(int n) { addLine(String.format("Track nr hits <%d>",n)); }
- -// void printStrip(int id, int layer, double s, Hep3Vector trkpos, Hep3Vector trkposlocal, Hep3Vector m, Hep3Vector res_local,Hep3Vector res_local_error, Hep3Vector u, Hep3Vector origin) { -// addLine(String.format("Strip id <%d> layer <%d> org <%f %f %f> s <%f> tp <%f %f %f> tplocal <%f %f %f> m <%f %f %f> res <%f %f %f> reserr <%f %f %f> u <%f %f %f> ",id,layer, -// origin.x(),origin.y(),origin.z(), -// s, -// trkpos.x(),trkpos.y(),trkpos.z(), -// trkposlocal.x(),trkposlocal.y(),trkposlocal.z(), -// m.x(),m.y(),m.z(), -// res_local.x(),res_local.y(),res_local.z(), -// res_local_error.x(),res_local_error.y(),res_local_error.z(), -// u.x(),u.y(),u.z())); -// }
void printStripJacPointToPoint(int id, int layer, double s, BasicMatrix jac) {
- String str = String.format("Strip id <%d> layer <%d> s <%f> jac <",id,layer,s);
+ String str = String.format("Strip id <%d> layer <%d> s <%.10f> jac <",id,layer,s);
for(int r=0;r<jac.getNRows();++r) { for(int c=0;c<jac.getNColumns();++c) {
- str += String.format("%f ", jac.e(r, c));
+ str += String.format("%.10f ", jac.e(r, c));
} } str += ">";
@@ -126,10 +115,10 @@
void printStripScatJacPointToPoint(int id, int layer, double s, BasicMatrix jac) {
- String str = String.format("Strip scat id <%d> layer <%d> s <%f> jac <",id,layer,s);
+ String str = String.format("Strip scat id <%d> layer <%d> s <%.10f> jac <",id,layer,s);
for(int r=0;r<jac.getNRows();++r) { for(int c=0;c<jac.getNColumns();++c) {
- str += String.format("%f ", jac.e(r, c));
+ str += String.format("%.10f ", jac.e(r, c));
} } str += ">";
@@ -137,10 +126,10 @@
} void printStripL2m(int id, int layer, double s, BasicMatrix jac) {
- String str = String.format("Strip LocalToMeas id <%d> layer <%d> s <%f> L2m <",id,layer,s);
+ String str = String.format("Strip LocalToMeas id <%d> layer <%d> s <%.10f> L2m <",id,layer,s);
for(int r=0;r<jac.getNRows();++r) { for(int c=0;c<jac.getNColumns();++c) {
- str += String.format("%f ", jac.e(r, c));
+ str += String.format("%.10f ", jac.e(r, c));
} } str += ">";
@@ -162,7 +151,7 @@
String str = "Track clCov "; for(int irow=0;irow<cov.getNRows();++irow) { for(int icol=0;icol<cov.getNColumns();++icol) {
- str += String.format("%f ", cov.e(irow, icol));
+ str += String.format("%.10f ", cov.e(irow, icol));
} } addLine(str);
@@ -170,16 +159,16 @@
void printStripTrackDir(double sinPhi, double sinLambda) {
- String str = String.format("Strip sinPhi sinLambda %f %f",sinPhi,sinLambda);
+ String str = String.format("Strip sinPhi sinLambda %.10f %.10f",sinPhi,sinLambda);
addLine(str); } void printStripTrackDirFull(Hep3Vector dir) {
- addLine(String.format("Strip track dir %f %f %f",dir.x(),dir.y(),dir.z()));
+ addLine(String.format("Strip track dir %.10f %.10f %.10f",dir.x(),dir.y(),dir.z()));
} void printStripTrackPos(Hep3Vector pos) {
- addLine(String.format("Strip track pos %f %f %f",pos.x(),pos.y(),pos.z()));
+ addLine(String.format("Strip track pos %.10f %.10f %.10f",pos.x(),pos.y(),pos.z()));
} void printStrip(int id, int layer) {
@@ -187,74 +176,74 @@
} void printStripPathLen(double s) {
- addLine(String.format("Strip pathLen %f", s));
+ addLine(String.format("Strip pathLen %.10f", s));
} void printStripPathLen3D(double s) {
- addLine(String.format("Strip pathLen3D %f", s));
+ addLine(String.format("Strip pathLen3D %.10f", s));
} void printStereoAngle(double stereoAngle) {
- addLine(String.format("Strip stereo angle %f", stereoAngle));
+ addLine(String.format("Strip stereo angle %.10f", stereoAngle));
} void printStripMeas(double u) {
- addLine(String.format("Strip u %f", u));
+ addLine(String.format("Strip u %.10f", u));
} void printStripMeasRes(double ures, double uresErr) {
- addLine(String.format("Strip ures %f %f", ures, uresErr));
+ addLine(String.format("Strip ures %.10f %.10f", ures, uresErr));
} void printStripMeasResTruth(double ures, double uresErr) {
- addLine(String.format("Strip truth ures %f %f", ures, uresErr));
+ addLine(String.format("Strip truth ures %.10f %.10f", ures, uresErr));
} void printStripMeasResSimHit(double ures, double uresErr) {
- addLine(String.format("Strip simhit ures %f %f", ures, uresErr));
+ addLine(String.format("Strip simhit ures %.10f %.10f", ures, uresErr));
} void printStripScat(double scatAngle) {
- addLine(String.format("Strip scatangle %f",scatAngle));
+ addLine(String.format("Strip scatangle %.10f",scatAngle));
} void printMeasDir(Hep3Vector u) {
- addLine(String.format("Strip meas dir %f %f %f", u.x(), u.y(), u.z()));
+ addLine(String.format("Strip meas dir %.10f %.10f %.10f", u.x(), u.y(), u.z()));
} void printNonMeasDir(Hep3Vector u) {
- addLine(String.format("Strip non-meas dir %f %f %f", u.x(), u.y(), u.z()));
+ addLine(String.format("Strip non-meas dir %.10f %.10f %.10f", u.x(), u.y(), u.z()));
} void printNormalDir(Hep3Vector u) {
- addLine(String.format("Strip normal dir %f %f %f", u.x(), u.y(), u.z()));
+ addLine(String.format("Strip normal dir %.10f %.10f %.10f", u.x(), u.y(), u.z()));
} void printMomentum(double p, double p_truth) {
- addLine(String.format("Track mom %f %f", p, p_truth));
+ addLine(String.format("Track mom %.10f %.10f", p, p_truth));
} void printPerToClPrj(Hep3Matrix perToClPrj) { String str = "Track perToClPrj "; for(int irow=0;irow<perToClPrj.getNRows();++irow) { for(int icol=0;icol<perToClPrj.getNColumns();++icol) {
- str += String.format("%f ", perToClPrj.e(irow, icol));
+ str += String.format("%.10f ", perToClPrj.e(irow, icol));
} } addLine(str); } void printChi2(double[] chisq, int[] ndf) {
- addLine(String.format("Track chi2/ndf (circle,zfit) %f %d %f %d",chisq[0],ndf[0],chisq[1],ndf[1]));
+ addLine(String.format("Track chi2/ndf (circle,zfit) %.10f %d %.10f %d",chisq[0],ndf[0],chisq[1],ndf[1]));
} void printOrigin(Hep3Vector pos) {
- addLine(String.format("Strip origin pos %f %f %f",pos.x(),pos.y(),pos.z()));
+ addLine(String.format("Strip origin pos %.10f %.10f %.10f",pos.x(),pos.y(),pos.z()));
} void printHitPos3D(Hep3Vector pos) {
- addLine(String.format("Strip 3D hit pos %f %f %f",pos.x(),pos.y(),pos.z()));
+ addLine(String.format("Strip 3D hit pos %.10f %.10f %.10f",pos.x(),pos.y(),pos.z()));
}
diff -u -r1.9 -r1.10 --- GBLOutput.java 5 Sep 2013 00:05:14 -0000 1.9 +++ GBLOutput.java 12 Sep 2013 17:57:38 -0000 1.10 @@ -36,16 +36,12 @@
*/ public class GBLOutput {
- private FileWriter _fWriter; - private PrintWriter _pWriter;
private int _debug;
- private boolean _hideFrame = false;
private GBLFileIO file; private Hep3Vector _B; private TrackerHitUtils _trackerHitUtils = new TrackerHitUtils(); private MaterialSupervisor _materialmanager; private MultipleScattering _scattering;
- private Detector detector;
private HPSTransformations _hpstrans = new HPSTransformations(); private double _beamEnergy = 2.2; //GeV
@@ -67,9 +63,6 @@
public void setDebug(int debug) { _debug = debug; }
- public void setHideFrame(boolean hide) { - _hideFrame = hide; - }
public void buildModel(Detector detector) { _materialmanager.buildModel(detector); }
@@ -249,22 +242,71 @@
//Find intercept point with sensor in tracking frame Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B.z())); Hep3Vector trkposXPlane = TrackUtils.getHelixXPlaneIntercept(htf, strip.w(), origin);
- Hep3Vector trkposXPlaneDCAFlip = TrackUtils.getHelixXPlaneIntercept(htf_dca_flip, strip.w(), origin);
+ Hep3Vector trkposXPlaneIter = TrackUtils.getHelixPlanePositionIter(htf, strip.origin(), strip.w(), 1.0e-8); + Hep3Vector trkposDiff = VecOp.sub(trkposXPlaneIter, trkpos); + if(trkposDiff.magnitude() > 1.0e-7) { + System.out.printf("WARNING trkposDiff mag = %.10f [%.10f %.10f %.10f]\n",trkposDiff.magnitude(),trkposDiff.x(),trkposDiff.y(),trkposDiff.z()); + System.exit(1); + }
if(_debug>0) {
- System.out.printf("trkpos at intercept %s\n",trkpos.toString()); - System.out.printf("trkposXPlane at intercept %s\n",trkposXPlane.toString()); - System.out.printf("trkposXPlaneDCAFlip at intercept %s\n",trkposXPlaneDCAFlip.toString());
+ System.out.printf("trkpos at intercept [%.10f %.10f %.10f]\n",trkpos.x(),trkpos.y(),trkpos.z()); + //System.out.printf("trkposXPlane at intercept %s\n",trkposXPlane.toString()); + System.out.printf("trkposXPlaneIter at intercept [%.10f %.10f %.10f]\n",trkposXPlaneIter.x(),trkposXPlaneIter.y(),trkposXPlaneIter.z());
// is it on the plane?
- double trkpos_d = VecOp.dot(strip.w(),VecOp.unit(VecOp.sub(trkpos, strip.origin()))); - double trkposXPlane_d = VecOp.dot(strip.w(),VecOp.unit(VecOp.sub(trkposXPlane, strip.origin()))); - System.out.printf("trkpos_d %f with diff %s trkpos %s org %s \n",trkpos_d,VecOp.unit(VecOp.sub(trkpos, strip.origin())), trkpos.toString(),strip.origin().toString()); - System.out.printf("trkposXPlane_d %f with diff %s trkpos %s org %s \n",trkposXPlane_d,VecOp.unit(VecOp.sub(trkposXPlane, strip.origin())), trkposXPlane.toString(),strip.origin().toString());
+ //double trkpos_d = VecOp.dot(strip.w(),VecOp.unit(VecOp.sub(trkpos, strip.origin()))); + //double trkposXPlaneIter_d = VecOp.dot(strip.w(),VecOp.unit(VecOp.sub(trkposXPlaneIter, strip.origin()))); + //System.out.printf("trkpos_d %f with diff %s trkpos %s org %s \n",trkpos_d,VecOp.unit(VecOp.sub(trkpos, strip.origin())), trkpos.toString(),strip.origin().toString()); + //System.out.printf("trkposXPlane_d %f with diff %s trkpos %s org %s \n",trkposXPlaneIter_d,VecOp.unit(VecOp.sub(trkposXPlaneIter, strip.origin())), trkposXPlaneIter.toString(),strip.origin().toString());
} //trkpos = trkposXPlane; //trkpos = trkposXPlaneDCAFlip;
+// if(strip.layer()==1) trkpos = new BasicHep3Vector(88.8078567539 ,0.213035827282, 1.5292560945); +// if(strip.layer()==2) trkpos = new BasicHep3Vector(96.2602646374, 0.378914765887 , 1.686096828); +// if(strip.layer()==3) trkpos = new BasicHep3Vector(188.751065257, 3.59937204783, 3.63336621837); +// if(strip.layer()==4) trkpos = new BasicHep3Vector(196.197763346 ,3.95227914386, 3.7902238256); +// if(strip.layer()==5) trkpos = new BasicHep3Vector(288.617639967, 9.49588626783 ,5.7383167412); +// if(strip.layer()==6) trkpos = new BasicHep3Vector(296.05862198 ,10.0360590169 ,5.89529022351 ); +// if(strip.layer()==7) trkpos = new BasicHep3Vector(488.64433778 ,28.9100728781, 9.96718984117); +// if(strip.layer()==8) trkpos = new BasicHep3Vector(496.071890128, 29.8276792157, 10.1246568417); +// if(strip.layer()==9) trkpos = new BasicHep3Vector(687.835519286, 58.4555308451, 14.2045621589); +// if(strip.layer()==10) trkpos = new BasicHep3Vector(695.251412131, 59.7550500287, 14.3629733178); +// if(strip.layer()==11) trkpos = new BasicHep3Vector(886.710664575, 98.3529112595, 18.472815436); +// if(strip.layer()==12) trkpos = new BasicHep3Vector(894.114639235, 100.042816527 ,18.6326045072 ); +// + + + if(_debug>0) { + //calculate track position for a given number of x-plane positions + double[] x_tmp = {88.8078567539,96.2602646374,188.751065257,196.197763346,288.617639967,296.05862198,488.64433778,496.071890128,687.835519286,695.251412131,886.710664575,894.114639235}; + for(int i=0;i<x_tmp.length;++i) { + double s_tmp = HelixUtils.PathToXPlane(htf, x_tmp[i], 0., 0).get(0); + Hep3Vector trkpos_x_tmp = HelixUtils.PointOnHelix(htf, s_tmp); + + //Calculate path length by hand for x plane + double dca = htf.dca(); + double xc = (htf.R() - dca) * Math.sin(htf.phi0()); + double sinPhi = (xc - x_tmp[i])/htf.R(); + double phi_at_x = Math.asin(sinPhi); + //double phi_at_x = Math.asin( 1/htf.R()* ( -1*dca*Math.sin(htf.phi0()) + htf.R()*Math.sin(htf.phi0()) - trkpos.x() ) ); + double dphi_at_x = phi_at_x - htf.phi0(); + if(dphi_at_x > Math.PI) dphi_at_x -= 2.0*Math.PI; + if(dphi_at_x < -Math.PI) dphi_at_x += 2.0*Math.PI; + double s_at_x = -1.0*dphi_at_x*htf.R(); + double y_at_x = dca*Math.cos(htf.phi0()) - htf.R()*Math.cos(htf.phi0()) + htf.R()*Math.cos(phi_at_x); + double z_at_x = htf.z0() + s_at_x*htf.slope(); + Hep3Vector trkposmanual_x_tmp = new BasicHep3Vector(x_tmp[i],y_at_x,z_at_x); + //System.out.printf("s = %f trkpos %s (initial manual xplane)\n", s_at_x , trkposmanual_at_x.toString()); + //System.out.printf("params %f %f %f %f %f\n", dca,htf.z0(),htf.phi0(),htf.slope(),htf.curvature()); + System.out.printf("x %f gives trkpos at %s for s %f s3D %f\n", x_tmp[i],trkpos_x_tmp.toString(),s_tmp,s_tmp/Math.cos(Math.atan(htf.slope()))); + System.out.printf("x %f gives trkpos at %s for s %f s3D %f (manual)\n", x_tmp[i],trkposmanual_x_tmp.toString(),s_at_x,s_at_x/Math.cos(Math.atan(htf.slope()))); + + } + } + +
Hep3Vector trkposTruth = TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z()));
@@ -313,21 +355,9 @@
double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0); double s3D = s / Math.cos(Math.atan(htf.slope()));
- if(_debug>0) { - - //Calculate path length by hand for x plane - { - double phi_at_x = Math.asin( 1/htf.R()* ( -1*htf.dca()*Math.sin(htf.phi0()) + htf.R()*Math.sin(htf.phi0()) - trkpos.x() ) ); - double dphi_at_x = phi_at_x - htf.phi0(); - //if(dphi_at_x > Math.PI) dphi_at_x -= 2.0*Math.PI; - double s_at_x = -1.0*dphi_at_x*htf.R(); - double y_at_x = htf.dca()*Math.cos(htf.phi0()) - htf.R()*Math.cos(htf.phi0()) + htf.R()*Math.cos(phi_at_x); - double z_at_x = htf.z0() + s_at_x*htf.slope(); - Hep3Vector trkposmanual_at_x = new BasicHep3Vector(trkpos.x(),y_at_x,z_at_x); - System.out.printf("s = %f trkpos %s (initial manual xplane)\n", s_at_x , trkposmanual_at_x.toString()); - - }
+
+ if(_debug>0) {
double s_XPlane = HelixUtils.PathToXPlane(htf,trkposXPlane.x(),0,0).get(0);
@@ -342,7 +372,7 @@
System.out.printf("s = %f trkpos %s (initial)\n", s_test, trkpos_test.toString()); s_test = 88.8103791522 * Math.cos(Math.atan(htfTruth.slope())); trkpos_test = HelixUtils.PointOnHelix(htfTruth, s_test);
- System.out.printf("s = %f trkpos %s (truth)\n", s_test, trkpos_test.toString());
+ System.out.printf("s = %f trkpos %s (truth)\n", s_test, trkpos_test.toString());
} //print the path length
@@ -370,6 +400,7 @@
// start by find the distance vector between the center and the track position Hep3Vector vdiffTrk = VecOp.sub(trkpos, origin); Hep3Vector vdiffTrkTruth = VecOp.sub(trkposTruth, origin);
+ Hep3Vector vdiffTrkIter = VecOp.sub(trkposXPlaneIter, origin);
// then find the rotation from tracking to measurement frame Hep3Matrix trkToStripRot = _trackerHitUtils.getTrackToStripRotation(strip);
@@ -396,7 +427,12 @@
// residual in measurement frame Hep3Vector res_meas = VecOp.sub(m_meas, trkpos_meas); Hep3Vector resTruth_meas = VecOp.sub(m_meas, trkposTruth_meas);
+ double pred_meas_iter = VecOp.dot(strip.u(),vdiffTrkIter); + double res_meas_iter = umeas - pred_meas_iter;
+ if(_debug>0) { + System.out.printf("layer %d uRes %.10f uRes_alt %.10f\n",strip.layer(),res_meas.x(),res_meas_iter); + }
// residual uncertainty in measurement frame Hep3Vector res_err_meas = new BasicHep3Vector(uError,vError,wError);
@@ -434,15 +470,13 @@
System.out.printf("%s: WARNING cannot find scatter for detector %s with strip cluster at %s\n",this.getClass(),((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement().getName(),strip.origin().toString()); //can be edge case where helix is outside, but close to sensor, so use hack with the sensor origin closest to hit -> FIX THIS!
- MaterialManager matMan = ((MaterialManager)_scattering.getMaterialManager());
DetectorPlane closest = null; double dx = 999999.9; if(MaterialSupervisor.class.isInstance(_scattering.getMaterialManager())) { MaterialSupervisor matSup = (MaterialSupervisor)_scattering.getMaterialManager(); for(ScatteringDetectorVolume vol : matSup.getMaterialVolumes()) { DetectorPlane plane = (DetectorPlane)vol;
- double x_origin = plane.origin().x(); - double dx_loop = strip.origin().x();
+ double dx_loop = Math.abs(strip.origin().x() - plane.origin().x());
if(dx_loop<dx) { dx = dx_loop; closest = plane;
diff -u -r1.4 -r1.5 --- GBLOutputDriver.java 30 Aug 2013 01:23:22 -0000 1.4 +++ GBLOutputDriver.java 12 Sep 2013 17:57:38 -0000 1.5 @@ -63,7 +63,6 @@
System.out.printf("%s: B-field in z %s\n",this.getClass().getSimpleName(),bfield.toString()); gbl = new GBLOutput(gblFile,bfield); gbl.setDebug(_debug);
- gbl.setHideFrame(hideFrame);
gbl.buildModel(detector); truthRes = new TruthResiduals(bfield); truthRes.setDebug(_debug);
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