hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
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()));
}
hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
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;