6 added + 4 modified, total 10 files
java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment/AlignmentParameters.java (rev 0)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment/AlignmentParameters.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -0,0 +1,559 @@
+package org.lcsim.hps.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 java.io.FileWriter;
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+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.RawTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.fit.helicaltrack.MultipleScatter;
+import org.lcsim.fit.helicaltrack.TrackDirection;
+import org.lcsim.hps.event.HPSTransformations;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+
+/**
+ * Class to calculate and print the residuals and derivatives
+ * of the alignment parameters...used as input for MillePede
+ * Notation follows the MillePede manual:
+ * http://www.desy.de/~blobel/Mptwo.pdf
+ *
+ * the track is measured in the HelicalTrackFit frame
+ * and residuals are in the sensor frame (u,v,w)
+ *
+ * ordering of track parameters is
+ * double d0 = _trk.dca();
+ * double z0 = _trk.z0();
+ * double slope = _trk.slope();
+ * double phi0 = _trk.phi0();
+ * double R = _trk.R();
+ *
+ * @author mgraham
+ */
+public class AlignmentParameters {
+
+ private int _nlc = 5; //the five track parameters
+ private int _ngl = 1; //delta(u) and delta(gamma) for each plane
+ private BasicMatrix _dfdq;
+ private BasicMatrix _dfdp;
+ private HelicalTrackFit _trk;
+ private double[] _resid = new double[3];
+ private double[] _error = new double[3];
+ private int[] _globalLabel = new int[1];
+ FileWriter fWriter;
+ PrintWriter pWriter;
+ Set<SiSensor> _process_sensors = new HashSet<SiSensor>();
+ boolean _DEBUG = false;
+ double smax = 1e3;
+
+ public AlignmentParameters(String outfile) {
+ 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);
+ }
+
+ }
+
+ public void PrintResidualsAndDerivatives(Track track) {
+ SeedTrack st = (SeedTrack) track;
+ SeedCandidate seed = st.getSeedCandidate();
+ Map<HelicalTrackHit, MultipleScatter> msmap = seed.getMSMap();
+ _trk = seed.getHelix();
+ List<TrackerHit> hitsOnTrack = track.getTrackerHits();
+ for (TrackerHit hit : hitsOnTrack) {
+ HelicalTrackHit htc = (HelicalTrackHit) hit;
+ double msdrphi = msmap.get(htc).drphi();
+ double msdz = msmap.get(htc).dz();
+ double sHit = _trk.PathMap().get(htc);
+ HelicalTrackCross cross = (HelicalTrackCross) htc;
+ List<HelicalTrackStrip> clusterlist = cross.getStrips();
+ 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);
+ }
+ }
+ AddTarget(0.1, 0.02);
+ }
+
+ private void CalculateLocalDerivatives(HelicalTrackStrip strip) {
+ //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[][] dfdq = new double[3][5];
+ //dx/dq
+ //these are wrong for X, but for now it doesn't matter
+ dfdq[0][0] = Math.sin(phi0);
+ dfdq[0][1] = 0;
+ dfdq[0][2] = 0;
+ dfdq[0][3] = d0 * Math.cos(phi0) + R * Math.sin(phi0) - s * Math.cos(phi0);
+ dfdq[0][4] = (phi - phi0) * Math.cos(phi0);
+ double[] mydydq = dydq(R, d0, phi0, xint, s);
+ double[] mydzdq = dzdq(R, d0, phi0, xint, slope, s);
+ for (int i = 0; i < 5; i++) {
+ dfdq[1][i] = mydydq[i];
+ dfdq[2][i] = mydzdq[i];
+ }
+
+ BasicMatrix dfdqGlobal = FillMatrix(dfdq, 3, 5);
+ Hep3Matrix trkToStrip = getTrackToStripRotation(strip);
+ _dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal);
+
+ if (_DEBUG) {
+ double[] trackpars = {d0, z0, slope, phi0, R, s, xint};
+ System.out.println("Strip Origin: ");
+ System.out.println(strip.origin());
+ System.out.println("trkToStrip Rotation:");
+ System.out.println(trkToStrip.toString());
+ printDerivatives(trackpars, dfdq);
+ }
+ }
+
+ private void CalculateGlobalDerivatives(HelicalTrackStrip strip) {
+ //1st index = alignment parameter (only u so far)
+ //2nd index = residual coordinate (on du so far)
+
+ double[][] dfdpLab = new double[3][1];
+ dfdpLab[0][0] = 0; //df/dx
+ dfdpLab[1][0] = 0; //df/dy
+ dfdpLab[2][0] = 1; //df/dz
+ BasicMatrix _dfdpLab = FillMatrix(dfdpLab, 3, 1);
+ Hep3Matrix trkToStrip = getTrackToStripRotation(strip);
+ _dfdp = (BasicMatrix) MatrixOp.mult(trkToStrip, _dfdpLab);
+ if (_DEBUG) {
+ System.out.printf("dfdz = %5.5f %5.5f %5.5f\n", _dfdp.e(0, 0), _dfdp.e(1, 0), _dfdp.e(2, 0));
+ }
+ _globalLabel[0] = GetIdentifier(strip);
+// _globalLabel[0] = GetIdentifierModule(strip);
+
+ }
+
+ private void CalculateResidual(HelicalTrackStrip strip, double msdrdphi, double msdz) {
+
+ 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 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 = getTrackToStripRotation(strip);
+ Hep3Vector vdiff = VecOp.mult(trkToStrip, vdiffTrk);
+ 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 = 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("Strip Origin: ");
+ System.out.println(corigin.toString());
+ System.out.println("Position on Track:");
+ System.out.println(trkpos.toString());
+ 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);
+ System.out.println("udiff = " + _resid[0] + " +/- " + _error[0]);
+
+ }
+
+ }
+
+ public double[] getResidual(Track track, int layer) {
+ double[] res = new double[7];
+ SeedTrack st = (SeedTrack) track;
+ SeedCandidate seed = st.getSeedCandidate();
+ Map<HelicalTrackHit, MultipleScatter> msmap = seed.getMSMap();
+ _trk = seed.getHelix();
+ List<TrackerHit> hitsOnTrack = track.getTrackerHits();
+ for (TrackerHit hit : hitsOnTrack) {
+ HelicalTrackHit htc = (HelicalTrackHit) hit;
+ double sHit = _trk.PathMap().get(htc);
+ HelicalTrackCross cross = (HelicalTrackCross) htc;
+ List<HelicalTrackStrip> clusterlist = cross.getStrips();
+ TrackDirection trkdir = HelixUtils.CalculateTrackDirection(_trk, sHit);
+ double msdrphi = msmap.get(htc).drphi();
+ double msdz = msmap.get(htc).dz();
+ cross.setTrackDirection(trkdir, _trk.covariance());
+ for (HelicalTrackStrip cl : clusterlist) {
+ if (cl.layer() == layer) {
+ CalculateResidual(cl, msdrphi, msdz);
+ res[0] = _resid[0];
+ res[1] = _resid[1];
+ res[2] = _resid[2];
+ res[3] = _error[0];
+ res[4] = _error[1];
+ res[5] = _error[2];
+ res[6] = layer;
+ if(hit.getPosition()[2]<0)res[6]=layer+10;
+ }
+ }
+ }
+ return res;
+
+ }
+
+ public void AddTarget(double beamdy, double beamdz) {
+ double[][] dfdp = new double[3][1];
+ double d0 = _trk.dca();
+ double z0 = _trk.z0();
+ double slope = _trk.slope();
+ double phi0 = _trk.phi0();
+ double R = _trk.R();
+ double xint = 0; //target
+ double s = HelixUtils.PathToXPlane(_trk, xint, smax, _nlc).get(0);
+ Hep3Vector ptAtTarget = HelixUtils.PointOnHelix(_trk, s);
+ double[] mydydq = dydq(R, d0, phi0, xint, s);
+ double[] mydzdq = dzdq(R, d0, phi0, xint, slope, s);
+ _resid[0] = ptAtTarget.z();
+ _resid[1] = ptAtTarget.y();
+ _resid[2] = ptAtTarget.x();
+ _error[0] = beamdz;
+ _error[1] = beamdy;
+ _error[2] = 666;
+ dfdp[0][0] = 1;
+ dfdp[1][0] = 0;
+ dfdp[2][0] = 0;
+ _dfdp = FillMatrix(dfdp, 3, 1);
+ _globalLabel[0] = 666;
+ pWriter.printf("%4d\n", 666);
+ pWriter.printf("%5.5e %5.5e %5.5e\n", _resid[0], _resid[1], _resid[2]);
+ pWriter.printf("%5.5e %5.5e %5.5e\n", _error[0], _error[1], _error[2]);
+ for (int i = 0; i < _nlc; i++) {
+ pWriter.printf("%5.5e %5.5e -1.0\n", mydzdq[i], mydydq[i]);
+ }
+ for (int j = 0; j < _ngl; j++) {
+ pWriter.printf("%5.5e %5.5e %5.5e %5d\n", _dfdp.e(0, j), _dfdp.e(1, j), _dfdp.e(2, j), _globalLabel[j]);
+ }
+
+ }
+
+ private void PrintStripResiduals(HelicalTrackStrip strip) {
+ if (_DEBUG) {
+ System.out.printf("Strip Layer = %4d\n", strip.layer());
+ 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");
+ for (int i = 0; i < _nlc; i++) {
+ System.out.printf("%s %5.5e %5.5e %5.5e\n", q[i], _dfdq.e(0, i), _dfdq.e(1, i), _dfdq.e(2, i));
+ }
+ String[] p = {"u-displacement"};
+ System.out.println("global parameter derivatives");
+ for (int j = 0; j < _ngl; j++) {
+ 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]);
+ }
+
+ }
+ pWriter.printf("%4d\n", strip.layer());
+ pWriter.printf("%5.5e %5.5e %5.5e\n", _resid[0], _resid[1], _resid[2]);
+ pWriter.printf("%5.5e %5.5e %5.5e\n", _error[0], _error[1], _error[2]);
+ for (int i = 0; i < _nlc; i++) {
+ pWriter.printf("%5.5e %5.5e %5.5e\n", _dfdq.e(0, i), _dfdq.e(1, i), _dfdq.e(2, i));
+ }
+ for (int j = 0; j < _ngl; j++) {
+ pWriter.printf("%5.5e %5.5e %5.5e %5d\n", _dfdp.e(0, j), _dfdp.e(1, j), _dfdp.e(2, j), _globalLabel[j]);
+ }
+ }
+
+ private Hep3Matrix getTrackToStripRotation(HelicalTrackStrip strip) {
+ ITransform3D detToStrip = GetGlobalToLocal(strip);
+ Hep3Matrix detToStripMatrix = (BasicHep3Matrix) detToStrip.getRotation().getRotationMatrix();
+ Hep3Matrix detToTrackMatrix = (BasicHep3Matrix) HPSTransformations.getMatrix();
+
+ if (_DEBUG) {
+ System.out.println("gblToLoc translation:");
+ System.out.println(detToStrip.getTranslation().toString());
+ System.out.println("gblToLoc Rotation:");
+ System.out.println(detToStrip.getRotation().toString());
+ System.out.println("detToTrack Rotation:");
+ System.out.println(detToTrackMatrix.toString());
+ }
+
+ return (Hep3Matrix) VecOp.mult(detToStripMatrix, VecOp.inverse(detToTrackMatrix));
+ }
+
+ private ITransform3D GetGlobalToLocal(HelicalTrackStrip strip) {
+ RawTrackerHit rth = (RawTrackerHit) strip.rawhits().get(0);
+ IDetectorElement ide = rth.getDetectorElement();
+ SiSensor sensor = ide.findDescendants(SiSensor.class).get(0);
+ SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
+ return electrodes.getGlobalToLocal();
+ }
+
+ private int GetIdentifier(HelicalTrackStrip strip) {
+ RawTrackerHit rth = (RawTrackerHit) strip.rawhits().get(0);
+ IDetectorElement ide = rth.getDetectorElement();
+ SiSensor sensor = ide.findDescendants(SiSensor.class).get(0);
+ // return rth.getIdentifierFieldValue(sensor.getName());
+ return sensor.getSensorID(); //individual sensor positions
+// int sid=sensor.getSensorID();
+// int global=1;
+// if(sid>10)global=2;
+// return global; //return top/bottom plates
+ }
+
+ private int GetIdentifierModule(HelicalTrackStrip strip) {
+ RawTrackerHit rth = (RawTrackerHit) strip.rawhits().get(0);
+ IDetectorElement ide = rth.getDetectorElement();
+ SiSensor sensor = ide.findDescendants(SiSensor.class).get(0);
+ // return rth.getIdentifierFieldValue(sensor.getName());
+// return sensor.getSensorID(); //individual sensor positions
+ int sid = sensor.getSensorID();
+ int gid = -1;
+ switch (sid) {
+ case 1:
+ gid = 1; break;
+ case 2:
+ gid = 1;break;
+ case 3:
+ gid = 2;break;
+ case 4:
+ gid = 2;break;
+ case 5:
+ gid = 3;break;
+ case 6:
+ gid = 3;break;
+ case 7:
+ gid = 4;break;
+ case 8:
+ gid = 4;break;
+ case 9:
+ gid = 5;break;
+ case 10:
+ gid = 5;break;
+ case 11:
+ gid = 11;break;
+ case 12:
+ gid = 11;break;
+ case 13:
+ gid = 12;break;
+ case 14:
+ gid = 12;break;
+ case 15:
+ gid = 13;break;
+ case 16:
+ gid = 13;break;
+ case 17:
+ gid = 14;break;
+ case 18:
+ gid = 14;break;
+ case 19:
+ gid = 15;break;
+ case 20:
+ gid = 15;break;
+ }
+
+ return gid; //return top/bottom plates
+ }
+
+ private BasicMatrix FillMatrix(double[][] array, int nrow, int ncol) {
+ BasicMatrix retMat = new BasicMatrix(nrow, ncol);
+ for (int i = 0; i < nrow; i++) {
+ for (int j = 0; j < ncol; j++) {
+ retMat.setElement(i, j, array[i][j]);
+ }
+ }
+ return retMat;
+ }
+
+ public void closeFile() throws IOException {
+ pWriter.close();
+ fWriter.close();
+ }
+
+ private double dsdR(double R, double d0, double phi0, double xint) {
+ double sqrtTerm = Sqrt(R * R - Math.pow(((d0 - R) * Sin(phi0) + xint), 2));
+
+ double rsign = Math.signum(R);
+ double dsdr = (1 / sqrtTerm) * ((-rsign * xint) + (-rsign) * d0 * Sin(phi0)
+ + ArcTan(R * Cos(phi0), (-R) * Sin(phi0))
+ * sqrtTerm
+ - ArcTan(rsign * sqrtTerm, xint + (d0 - R) * Sin(phi0))
+ * sqrtTerm);
+
+
+ if (_DEBUG)
+ System.out.println("xint = " + xint + "; dsdr = " + dsdr);
+ return dsdr;
+
+ }
+
+ private double dsdphi(double R, double d0, double phi0, double xint) {
+ double sqrtTerm = Sqrt(R * R - Math.pow(((d0 - R) * Sin(phi0) + xint), 2));
+ double rsign = Math.signum(R);
+ double dsdphi = R * (sqrtTerm + rsign * d0 * Cos(phi0) - rsign * R * Cos(phi0)) / sqrtTerm;
+ if (_DEBUG)
+ System.out.println("xint = " + xint + "; dsdphi = " + dsdphi);
+ return dsdphi;
+ }
+
+ private double dsdd0(double R, double d0, double phi0, double xint) {
+ double sqrtTerm = Sqrt(R * R - Math.pow(((d0 - R) * Sin(phi0) + xint), 2));
+ double rsign = Math.signum(R);
+ double dsdd0 = rsign * (R * Sin(phi0)) / sqrtTerm;
+ if (_DEBUG)
+ System.out.println("xint = " + xint + "; dsdd0 = " + dsdd0);
+ return dsdd0;
+ }
+
+ private double[] dydq(double R, double d0, double phi0, double xint, double s) {
+ double[] dy = new double[5];
+// dy[0] = Cos(phi0) + Cot(phi0 - s / R) * Csc(phi0 - s / R) * dsdd0(R, d0, phi0, xint);
+ dy[0] = Cos(phi0) - Sec(phi0 - s / R) * Tan(phi0 - s / R) * dsdd0(R, d0, phi0, xint);
+ dy[1] = 0;
+ dy[2] = 0;
+// dy[3] = (-(d0 - R)) * Sin(phi0) - R * Cot(phi0 - s / R) * Csc(phi0 - s / R) * (1 - dsdphi(R, d0, phi0, xint) / R);
+ dy[3] = (-(d0 - R)) * Sin(phi0) + Sec(phi0 - s / R) * Tan(phi0 - s / R) * (R - dsdphi(R, d0, phi0, xint));
+ // dy[4] = -Cos(phi0) + Csc(phi0 - s / R) - R * Cot(phi0 - s / R) * Csc(phi0 - s / R) * (s / (R * R) - dsdR(R, d0, phi0, xint) / R);
+ dy[4] = -Cos(phi0) + Sec(phi0 - s / R) + (1 / R) * Sec(phi0 - s / R) * Tan(phi0 - s / R) * (s - R * dsdR(R, d0, phi0, xint));
+ return dy;
+ }
+
+ private double[] dzdq(double R, double d0, double phi0, double xint, double slope, double s) {
+ double[] dz = new double[5];
+ dz[0] = slope * dsdd0(R, d0, phi0, xint);
+ dz[1] = 1;
+ dz[2] = s;
+ dz[3] = slope * dsdphi(R, d0, phi0, xint);
+ dz[4] = slope * dsdR(R, d0, phi0, xint);
+ return dz;
+ }
+
+ private double Csc(double val) {
+ return 1 / Math.sin(val);
+ }
+
+ private double Cot(double val) {
+ return 1 / Math.tan(val);
+ }
+
+ private double Sec(double val) {
+ return 1 / Math.cos(val);
+ }
+
+ private double Sin(double val) {
+ return Math.sin(val);
+ }
+
+ private double Cos(double val) {
+ return Math.cos(val);
+ }
+
+ private double Tan(double val) {
+ return Math.tan(val);
+ }
+
+ private double ArcTan(double val1, double val2) {
+ return Math.atan2(val1, val2);
+ }
+
+ private double Sign(double val) {
+ return Math.signum(val);
+ }
+
+ private double Sqrt(double val) {
+ return Math.sqrt(val);
+ }
+
+ private void printDerivatives(double[] trackpars, double[][] dfdq) {
+ System.out.println("======================================================");
+ System.out.println("s xint");
+ System.out.printf("%5.5f %5.5f\n", trackpars[5], trackpars[6]);
+ System.out.println(" d0 z0 slope phi0 R");
+ System.out.printf("Values %5.5f %5.5f %5.5f %5.5f %5.5f\n", trackpars[0], trackpars[1], trackpars[2], trackpars[3], trackpars[4]);
+ System.out.printf("dzdq ");
+ for (int i = 0; i < 5; i++) {
+ System.out.printf("%5.3e ", dfdq[2][i]);
+ }
+ System.out.println();
+ System.out.printf("dudq ");
+ for (int i = 0; i < _nlc; i++) {
+ System.out.printf("%5.3e ", _dfdq.e(0, i));
+ }
+ System.out.println();
+ System.out.println();
+ System.out.printf("dydq ");
+ for (int i = 0; i < 5; i++) {
+ System.out.printf("%5.3e ", dfdq[1][i]);
+ }
+ System.out.println();
+ System.out.printf("dvdq ");
+ for (int i = 0; i < _nlc; i++) {
+ System.out.printf("%5.3e ", _dfdq.e(1, i));
+ }
+ System.out.println();
+ System.out.println();
+ System.out.printf("dxdq ");
+ for (int i = 0; i < 5; i++) {
+ System.out.printf("%5.3e ", dfdq[0][i]);
+ }
+ System.out.println();
+ System.out.printf("dwdq ");
+ for (int i = 0; i < _nlc; i++) {
+ System.out.printf("%5.3e ", _dfdq.e(2, i));
+ }
+ System.out.println();
+ // System.out.println( _trk.xc()+ "; "+_trk.yc());
+// System.out.println( _trk.x0()+ "; "+_trk.y0());
+ }
+}
java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment/HPSStrips.java (rev 0)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment/HPSStrips.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -0,0 +1,482 @@
+package org.lcsim.hps.alignment;
+/*
+ * SiStrips.java
+ *
+ * Created on July 22, 2005, 4:07 PM
+ *
+ * To change this template, choose Tools | Options and locate the template under
+ * the Source Creation and Management node. Right-click the template and choose
+ * Open. You can then make changes to the template in the Source Editor.
+ */
+
+//import static org.lcsim.units.clhep.SystemOfUnits.*;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.ITransform3D;
+import org.lcsim.detector.Transform3D;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.SortedMap;
+import java.util.TreeMap;
+import java.util.Set;
+import java.util.HashSet;
+import java.util.List;
+import org.lcsim.detector.solids.GeomOp2D;
+import org.lcsim.detector.solids.GeomOp3D;
+import org.lcsim.detector.solids.Line3D;
+import org.lcsim.detector.solids.LineSegment3D;
+import org.lcsim.detector.solids.Point3D;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.ChargeDistribution;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.detector.tracker.silicon.SiStrips;
+
+/**
+ *
+ * @author tknelson
+ */
+public class HPSStrips extends SiStrips
+{
+
+ // Fields
+
+ // Object definition
+ private ChargeCarrier _carrier; // charge carrier collected
+ private int _nstrips; // number of strips
+ private double _pitch; // sense pitch
+ private IDetectorElement _detector; // associated detector element
+ private ITransform3D _parent_to_local; // parent to local transform
+ private ITransform3D _local_to_global; // transformation to global coordinates
+ private ITransform3D _global_to_local; // transformation from global coordinates
+ private Polygon3D _geometry; // region in which strips are defined
+ private double _capacitance_intercept = 10.; // fixed capacitance independent of strip length
+ private double _capacitance_slope = 0.1; // capacitance per unit length of strip
+
+ // Cached for convenience
+ private double _strip_offset;
+
+
+
+ public HPSStrips(ChargeCarrier carrier, double pitch, IDetectorElement detector, ITransform3D parent_to_local,ITransform3D misalignment)
+ {
+
+// System.out.println("Plane of polygon in sensor coordinates has... ");
+// System.out.println(" normal: "+((SiSensor)detector).getBiasSurface(carrier).getNormal());
+// System.out.println(" distance: "+((SiSensor)detector).getBiasSurface(carrier).getDistance());
+
+ setCarrier(carrier);
+ setPitch(pitch);
+ setGeometry(((SiSensor)detector).getBiasSurface(carrier).transformed(parent_to_local));
+ setStripNumbering();
+ setDetectorElement(detector);
+ setParentToLocal(parent_to_local);
+ setGlobalToLocal(Transform3D.multiply(Transform3D.multiply(parent_to_local,detector.getGeometry().getGlobalToLocal()),misalignment));
+ setLocalToGlobal(getGlobalToLocal().inverse());
+ }
+
+ public HPSStrips(ChargeCarrier carrier, double pitch, int nstrips, IDetectorElement detector, ITransform3D parent_to_local,ITransform3D misalignment)
+ {
+ setCarrier(carrier);
+ setPitch(pitch);
+ setGeometry(((SiSensor)detector).getBiasSurface(carrier).transformed(parent_to_local));
+ setNStrips(nstrips);
+ setDetectorElement(detector);
+ setParentToLocal(parent_to_local);
+ setGlobalToLocal(Transform3D.multiply(Transform3D.multiply(parent_to_local,detector.getGeometry().getGlobalToLocal()),misalignment));
+ setLocalToGlobal(getGlobalToLocal().inverse());
+ }
+
+
+ // SiSensorElectrodes interface
+ //=============================
+
+ // Mechanical properties
+ public int getNAxes()
+ {
+ return 1;
+ }
+
+ public IDetectorElement getDetectorElement()
+ {
+ return _detector;
+ }
+
+ public ITransform3D getParentToLocal()
+ {
+ return _parent_to_local;
+ }
+
+ public ITransform3D getLocalToGlobal()
+ {
+ return _local_to_global;
+ }
+
+ public ITransform3D getGlobalToLocal()
+ {
+ return _global_to_local;
+ }
+
+ public Polygon3D getGeometry()
+ {
+ return _geometry;
+ }
+
+ public Hep3Vector getMeasuredCoordinate(int axis)
+ {
+ if (axis == 0) return new BasicHep3Vector(1.0,0.0,0.0);
+ else return null;
+ }
+
+ public Hep3Vector getUnmeasuredCoordinate(int axis)
+ {
+ if (axis == 0) return new BasicHep3Vector(0.0,1.0,0.0);
+ else return null;
+ }
+
+ public int getNeighborCell(int cell, int ncells_0, int ncells_1)
+ {
+ int neighbor_cell = cell + ncells_0;
+ if (isValidCell(neighbor_cell)) return neighbor_cell;
+ else return -1;
+ }
+
+ public Set<Integer> getNearestNeighborCells(int cell)
+ {
+ Set<Integer> neighbors = new HashSet<Integer>();
+ for (int ineigh = -1 ; ineigh <= 1; ineigh=ineigh+2)
+ {
+ int neighbor_cell = getNeighborCell(cell,ineigh,0);
+ if (isValidCell(neighbor_cell)) neighbors.add(neighbor_cell);
+ }
+ return neighbors;
+ }
+
+ public boolean isValidCell(int cell)
+ {
+ return (cell >= 0 && cell < getNCells());
+ }
+
+ public int getNCells()
+ {
+ return _nstrips;
+ }
+
+ public int getNCells(int axis)
+ {
+ if (axis == 0)
+ {
+ return _nstrips;
+ }
+ else return 1;
+ }
+
+ public double getPitch(int axis)
+ {
+ if (axis == 0)
+ {
+ return _pitch;
+ }
+ else return 0;
+ }
+ public int getCellID(Hep3Vector position)
+ {
+ return (int)Math.round((position.x()+_strip_offset)/_pitch);
+ }
+
+
+ public int getRowNumber(Hep3Vector position)
+ {
+ return 0;
+ }
+
+ public int getColumnNumber(Hep3Vector position)
+ {
+ return getCellID(position);
+ }
+
+ public int getCellID(int row_number, int column_number)
+ {
+ return column_number;
+ }
+
+ public int getRowNumber(int cell_id)
+ {
+ return 0;
+ }
+
+ public int getColumnNumber(int cell_id)
+ {
+ return cell_id;
+ }
+
+ public Hep3Vector getPositionInCell(Hep3Vector position)
+ {
+ return VecOp.sub(position,getCellPosition(getCellID(position)));
+ }
+
+ public Hep3Vector getCellPosition(int strip_number)
+ {
+ return new BasicHep3Vector(strip_number*_pitch-_strip_offset,0.0,0.0);
+ }
+
+ // Electrical properties
+
+ /**
+ * Capacitance intercept parameter. Units are pF.
+ *
+ * Capacitance is calculated as:
+ * C = capacitance_intercept + strip_length * capacitance slope
+ *
+ * @param capacitance_intercept
+ */
+ public void setCapacitanceIntercept(double capacitance_intercept) {
+ _capacitance_intercept = capacitance_intercept;
+ }
+
+ /**
+ * Capacitance per unit strip length. Units are pF / mm.
+ *
+ * @param capacitance_slope
+ */
+ public void setCapacitanceSlope(double capacitance_slope) {
+ _capacitance_slope = capacitance_slope;
+ }
+
+ public ChargeCarrier getChargeCarrier()
+ {
+ return _carrier;
+ }
+
+ /**
+ * Capacitance for a particular cell. Units are pF.
+ *
+ * @param cell_id
+ * @return
+ */
+ public double getCapacitance(int cell_id) // capacitance in pF
+ {
+ return _capacitance_intercept + _capacitance_slope*getStripLength(cell_id);
+ }
+
+ /**
+ * Nominal capacitance used for throwing random noise in the sensor.
+ * Calculated using middle strip. Units are pF.
+ *
+ * @return
+ */
+ public double getCapacitance() {
+ return getCapacitance(getNCells(0) / 2);
+ }
+
+ public SortedMap<Integer,Integer> computeElectrodeData(ChargeDistribution distribution)
+ {
+ SortedMap<Integer,Integer> electrode_data = new TreeMap<Integer,Integer>();
+
+ int base_strip = getCellID(distribution.getMean());
+
+ // put charge on strips in window 3-sigma strips on each side of base strip
+ int axis = 0;
+ int window_size = (int)Math.ceil(3.0*distribution.sigma1D(getMeasuredCoordinate(axis))/getPitch(axis));
+
+ double integral_lower = distribution.getNormalization();
+ double integral_upper = distribution.getNormalization();
+
+ for (int istrip = base_strip-window_size; istrip <= base_strip+window_size; istrip++)
+ {
+ double cell_edge_upper = getCellPosition(istrip).x() + getPitch(axis)/2.0;
+
+// System.out.println("cell_edge_upper: "+cell_edge_upper);
+
+ double integration_limit = cell_edge_upper; //cell_edge_upper-distribution.mean().x();
+
+// System.out.println("integration_limit: "+integration_limit);
+
+ integral_upper = distribution.upperIntegral1D(getMeasuredCoordinate(axis),integration_limit);
+
+// System.out.println("integral_upper: "+integral_upper);
+
+ if (integral_lower<integral_upper)
+ {
+ throw new RuntimeException("Error in integrating Gaussian charge distribution!");
+ }
+
+ int strip_charge = (int)Math.round(integral_lower-integral_upper);
+
+// System.out.println("strip_charge: "+strip_charge);
+
+ if (strip_charge != 0)
+ {
+ electrode_data.put(istrip,strip_charge);
+ }
+
+ integral_lower = integral_upper;
+ }
+
+ return electrode_data;
+
+ }
+
+ // Strip specific methods
+
+ // length of strip
+ public double getStripLength(int cell_id)
+ {
+// System.out.println("strip_length: "+getStrip(cell_id).getLength());
+ return getStrip(cell_id).getLength();
+ }
+
+ // center of strip
+ public Hep3Vector getStripCenter(int cell_id)
+ {
+ LineSegment3D strip = getStrip(cell_id);
+ return strip.getEndPoint(strip.getLength()/2);
+ }
+
+ // line segment for strip
+ public LineSegment3D getStrip(int cell_id)
+ {
+ Line3D strip_line = new Line3D(new Point3D(getCellPosition(cell_id)),getUnmeasuredCoordinate(0));
+
+// System.out.println("Number of strips: "+this._nstrips);
+// System.out.println("Strip offset: "+this._strip_offset);
+// System.out.println("Pitch: "+this._pitch);
+// System.out.println("cell_id: "+cell_id);
+// System.out.println("strip_line start point: "+strip_line.getStartPoint());
+// System.out.println("strip_line direction: "+strip_line.getDirection());
+
+ List<Point3D> intersections = new ArrayList<Point3D>();
+
+ // Get intersections between strip line and edges of electrode polygon
+ for (LineSegment3D edge : _geometry.getEdges())
+ {
+// System.out.println("edge start point: "+edge.getStartPoint());
+// System.out.println("edge end point: "+edge.getEndPoint());
+
+ if (GeomOp2D.intersects(strip_line,edge))
+ {
+ intersections.add(GeomOp3D.lineBetween(strip_line,edge).getStartPoint());
+ }
+ }
+
+ // Check for rare occurrence of duplicates (can happen at corners of polygon)
+ List<Point3D> strip_ends = new ArrayList<Point3D>(intersections);
+ if (intersections.size() > 2)
+ {
+ for (int ipoint1 = 0; ipoint1 < intersections.size(); ipoint1++)
+ {
+ Point3D point1 = intersections.get(ipoint1);
+ for (int ipoint2 = ipoint1+1; ipoint2 < intersections.size(); ipoint2++)
+ {
+ Point3D point2 = intersections.get(ipoint2);
+ if (GeomOp3D.intersects(point1,point2))
+ {
+ strip_ends.remove(point2);
+ if (strip_ends.size() == 2) break;
+ }
+ }
+ }
+ }
+
+ return new LineSegment3D(strip_ends.get(0),strip_ends.get(1));
+ }
+
+ // Private setters
+ //==================
+ public void setCarrier(ChargeCarrier carrier)
+ {
+ _carrier = carrier;
+ }
+
+ public void setGeometry(Polygon3D geometry)
+ {
+// System.out.println("Plane of polygon has... ");
+// System.out.println(" normal: "+geometry.getNormal());
+// System.out.println(" distance: "+geometry.getDistance());
+//
+// System.out.println("Working plane has... ");
+// System.out.println(" normal: "+GeomOp2D.PLANE.getNormal());
+// System.out.println(" distance: "+GeomOp2D.PLANE.getDistance());
+
+ if (GeomOp3D.equals(geometry.getPlane(),GeomOp2D.PLANE))
+ {
+ _geometry = geometry;
+ }
+ else
+ {
+ throw new RuntimeException("Electrode geometry must be defined in x-y plane!!");
+ }
+ }
+
+ private void setStripNumbering()
+ {
+ double xmin = Double.MAX_VALUE;
+ double xmax = Double.MIN_VALUE;
+ for (Point3D vertex : _geometry.getVertices())
+ {
+ xmin = Math.min(xmin,vertex.x());
+ xmax = Math.max(xmax,vertex.x());
+ }
+
+// System.out.println("xmin: " + xmin);
+// System.out.println("xmax: " + xmax);
+//
+//
+// System.out.println("# strips: " + (int)Math.ceil((xmax-xmin)/getPitch(0)) ) ;
+
+ setNStrips( (int)Math.ceil((xmax-xmin)/getPitch(0)) ) ;
+ }
+
+ private void setNStrips(int nstrips)
+ {
+ _nstrips = nstrips;
+ setStripOffset();
+// _strip_offset = (_nstrips-1)*_pitch/2.;
+ }
+
+ private void setStripOffset()
+ {
+ double xmin = Double.MAX_VALUE;
+ double xmax = Double.MIN_VALUE;
+ for (Point3D vertex : _geometry.getVertices())
+ {
+ xmin = Math.min(xmin,vertex.x());
+ xmax = Math.max(xmax,vertex.x());
+ }
+
+ double strips_center = (xmin+xmax)/2;
+
+ _strip_offset = ((_nstrips-1)*_pitch)/2 - strips_center;
+
+ }
+
+ private void setPitch(double pitch)
+ {
+ _pitch = pitch;
+ }
+
+ private void setDetectorElement(IDetectorElement detector)
+ {
+ _detector = detector;
+ }
+
+ private void setParentToLocal(ITransform3D parent_to_local)
+ {
+ _parent_to_local = parent_to_local;
+ }
+
+ private void setLocalToGlobal(ITransform3D local_to_global)
+ {
+ _local_to_global = local_to_global;
+ }
+
+ private void setGlobalToLocal(ITransform3D global_to_local)
+ {
+ _global_to_local = global_to_local;
+ }
+
+}
+
+
java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment/RunAlignment.java (rev 0)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment/RunAlignment.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -0,0 +1,97 @@
+package org.lcsim.hps.alignment;
+
+import hep.aida.IAnalysisFactory;
+import java.io.IOException;
+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.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author mgraham
+ */
+public class RunAlignment extends Driver {
+
+ 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 = "1pt8";
+ AlignmentParameters ap;
+ int totalTracks=0;
+// 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
+// note: this should be -1 for Test configurations and +1 for Full (v3.X and lower) configurations
+// this is set by the _config variable (detType in HeavyPhotonDriver)
+ int flipSign = 1;
+
+ public RunAlignment(int trackerLayers, int mintrkLayers, String config) {
+ nlayers[0] = trackerLayers;
+ _minLayers = mintrkLayers;
+ _config = config;
+ if (_config.contains("Test"))
+ flipSign = -1;
+ ap = new AlignmentParameters("/Users/mgraham/HPS/align.txt");
+
+ }
+
+ 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) {
+ totalTracks++;
+ ap.PrintResidualsAndDerivatives(trk);
+
+ 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());
+ }
+
+ 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]);
+ }
+
+ }
+ }
+ }
+
+ }
+
+ public void endOfData() {
+ try {
+ System.out.println("Total Number of Tracks Found = "+totalTracks);
+ ap.closeFile();
+ } catch (IOException ex) {
+ Logger.getLogger(RunAlignment.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+}
java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment/SiTrackerSpectrometerSensorSetup.java (rev 0)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/alignment/SiTrackerSpectrometerSensorSetup.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -0,0 +1,265 @@
+package org.lcsim.hps.alignment;
+
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.List;
+
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IRotation3D;
+import org.lcsim.detector.ITranslation3D;
+import org.lcsim.detector.RotationPassiveXYZ;
+import org.lcsim.detector.Transform3D;
+import org.lcsim.detector.Translation3D;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.solids.Trd;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.subdetector.SiTrackerSpectrometer;
+import org.lcsim.util.Driver;
+
+public class SiTrackerSpectrometerSensorSetup extends Driver {
+
+ String subdetectorName;
+
+ public SiTrackerSpectrometerSensorSetup() {
+ }
+
+ public SiTrackerSpectrometerSensorSetup(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void setSubdetectorName(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void detectorChanged(Detector detector) {
+ if (subdetectorName == null) {
+ throw new RuntimeException("The subdetectorName was not set.");
+ }
+
+ Subdetector subdetector = detector.getSubdetector(subdetectorName);
+ if (subdetector instanceof SiTrackerSpectrometer) {
+ setupSensorDetectorElements(subdetector);
+ } else {
+ throw new RuntimeException("The subdetector " + subdetectorName + " is not an instance of SiTrackerSpectrometer.");
+ }
+ }
+
+ private void setupSensorDetectorElements(Subdetector subdet) {
+ System.out.println(this.getClass().getCanonicalName() + " - Setting up sensors for " + subdet.getName() + " ...");
+ int sensorId = 0;
+
+ for (IDetectorElement endcap : subdet.getDetectorElement().getChildren()) {
+ for (IDetectorElement layer : endcap.getChildren()) {
+ //int nwedges = layer.getChildren().size();
+ for (IDetectorElement wedge : layer.getChildren()) {
+ for (IDetectorElement module : wedge.getChildren()) {
+ List<SiSensor> sensors = module.findDescendants(SiSensor.class);
+
+ if (sensors.size() == 0) {
+ throw new RuntimeException("No sensors found in module.");
+ }
+
+ for (SiSensor sensor : sensors) {
+ Trd sensor_solid = (Trd) sensor.getGeometry().getLogicalVolume().getSolid();
+
+ Polygon3D n_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, -1, 0)).get(0);
+ Polygon3D p_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, 1, 0)).get(0);
+
+ // Bias the sensor
+// sensor.setBiasSurface(ChargeCarrier.ELECTRON, p_side);
+// sensor.setBiasSurface(ChargeCarrier.HOLE, n_side);
+
+ sensor.setBiasSurface(ChargeCarrier.HOLE, p_side);
+ sensor.setBiasSurface(ChargeCarrier.ELECTRON, n_side);
+
+// double strip_angle = Math.atan2(sensor_solid.getXHalfLength2() - sensor_solid.getXHalfLength1(), sensor_solid.getZHalfLength() * 2);
+ double strip_angle = 0.00;
+ ITranslation3D electrodes_position = new Translation3D(VecOp.mult(-p_side.getDistance(), new BasicHep3Vector(0, 0, 1))); // translate to outside of polygon
+ //ITranslation3D electrodes_position = new Translation3D(VecOp.mult(n_side.getDistance(), new BasicHep3Vector(0, 0, 1))); // translate to outside of polygon
+ // System.out.println("SensorID = " + sensorId + " " + electrodes_position.toString());
+ IRotation3D electrodes_rotation = new RotationPassiveXYZ(-Math.PI / 2, 0, strip_angle);
+ Transform3D electrodes_transform = new Transform3D(electrodes_position, electrodes_rotation);
+
+ // Free calculation of readout electrodes, sense electrodes determined thereon
+// SiStrips readout_electrodes = new SiStrips(ChargeCarrier.HOLE, 0.060, sensor, electrodes_transform);
+// SiStrips sense_electrodes = new SiStrips(ChargeCarrier.HOLE,0.030,(readout_electrodes.getNCells()*2-1),sensor,electrodes_transform);
+ ITranslation3D misalign_position;
+// System.out.println(layer.getName());
+// if (layer.getName().contains("3")) {
+// if (layer.getName().contains("3")&&layer.getName().contains("positive")) {
+/*
+ if (layer.getName().contains("positive")) {
+ System.out.println("Putting in a misalignment for layer "+layer.getName());
+ misalign_position = new Translation3D(0, 0.05, 0.0); // translate to outside of polygon
+ } else {
+ // misalign_position = new Translation3D(0, 0.0, 0.0);
+ misalign_position = new Translation3D(0, -0.05, 0.0);
+ }
+ */
+
+// if ((layer.getName().contains("3")||layer.getName().contains("4"))&&layer.getName().contains("positive")) {
+ // if (layer.getName().contains("positive")) {
+// System.out.println("Putting in a misalignment for layer "+layer.getName());
+// misalign_position = new Translation3D(0, 0.010, 0.0);
+// } else {
+// misalign_position = new Translation3D(0, 0.0, 0.0);
+ // misalign_position = new Translation3D(0, -0.05, 0.0);
+// }
+/*
+ if (layer.getName().contains("positive")) {
+ int gid = GetIdentifierModule(layer.getName());
+ misalign_position = new Translation3D(0, 0.0, 0.0);
+ if (gid == 1)
+ misalign_position = new Translation3D(0, -0.0144, 0.0);
+ if (gid == 2)
+ misalign_position = new Translation3D(0, 0.05-0.0297, 0.0);
+ if (gid == 3)
+ misalign_position = new Translation3D(0, -0.0253, 0.0);
+ if (gid == 4)
+ misalign_position = new Translation3D(0, -0.0346, 0.0);
+ if (gid == 5)
+ misalign_position = new Translation3D(0, -0.0433, 0.0);
+
+ } else {
+ misalign_position = new Translation3D(0, 0.0, 0.0);
+ }
+ */
+
+ int gid = GetIdentifierLayer(layer.getName());
+ misalign_position = new Translation3D(0, 0.0, 0.0);
+ /*
+ if (layer.getName().contains("positive")) {
+ if (gid == 1)
+ misalign_position = new Translation3D(0, -0.00144, 0.0);
+ if (gid == 2)
+ misalign_position = new Translation3D(0, 0.005, 0.0);
+ if (gid == 3)
+ misalign_position = new Translation3D(0, -0.00253, 0.0);
+ if (gid == 4)
+ misalign_position = new Translation3D(0, -0.00346, 0.0);
+ if (gid == 5)
+ misalign_position = new Translation3D(0, 0.00433, 0.0);
+ if (gid == 6)
+ misalign_position = new Translation3D(0, 0.0002, 0.0);
+ if (gid == 7)
+ misalign_position = new Translation3D(0, 0.002, 0.0);
+ if (gid == 8)
+ misalign_position = new Translation3D(0, -0.004, 0.0);
+ if (gid == 9)
+ misalign_position = new Translation3D(0, 0.006, 0.0);
+ if (gid == 10)
+ misalign_position = new Translation3D(0, -0.001, 0.0);
+ } else {
+ if (gid == 1)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 2)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 3)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 4)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 5)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 6)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 7)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 8)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 9)
+ misalign_position = new Translation3D(0, 0.00, 0.0);
+ if (gid == 10)
+ misalign_position = new Translation3D(0, 0.01, 0.0);
+ }
+
+*/
+ IRotation3D misalign_rotation = new RotationPassiveXYZ(0, 0, 0);
+ Transform3D misalign_transform = new Transform3D(misalign_position, misalign_rotation);
+
+ HPSStrips readout_electrodes = new HPSStrips(ChargeCarrier.HOLE, 0.060, sensor, electrodes_transform, misalign_transform);
+ HPSStrips sense_electrodes = new HPSStrips(ChargeCarrier.HOLE, 0.030, (readout_electrodes.getNCells() * 2 - 1), sensor, electrodes_transform, misalign_transform);
+
+
+// SiStrips readout_electrodes = new SiStrips(ChargeCarrier.ELECTRON, 0.060, sensor, electrodes_transform);
+// SiStrips sense_electrodes = new SiStrips(ChargeCarrier.ELECTRON, 0.030, (readout_electrodes.getNCells() * 2 - 1), sensor, electrodes_transform);
+
+ // SiSensorElectrodes sense_electrodes = new SiStrips(ChargeCarrier.HOLE, 0.060, sensor, electrodes_transform);
+
+//pristine conditions
+/*
+ readout_electrodes.setCapacitanceIntercept(0);
+ readout_electrodes.setCapacitanceSlope(0.12);
+ sense_electrodes.setCapacitanceIntercept(0);
+ sense_electrodes.setCapacitanceSlope(0.12);
+ */
+
+ readout_electrodes.setCapacitanceIntercept(0);
+ readout_electrodes.setCapacitanceSlope(0.16);
+ sense_electrodes.setCapacitanceIntercept(0);
+ sense_electrodes.setCapacitanceSlope(0.16);
+
+ sensor.setSenseElectrodes(sense_electrodes);
+ sensor.setReadoutElectrodes(readout_electrodes);
+//
+
+// double[][] transfer_efficiencies = {{1.0}};
+ double[][] transfer_efficiencies = {{0.986, 0.419}};
+ sensor.setTransferEfficiencies(ChargeCarrier.HOLE, new BasicMatrix(transfer_efficiencies));
+// sensor.setTransferEfficiencies(ChargeCarrier.ELECTRON, new BasicMatrix(transfer_efficiencies));
+ // here
+
+ sensor.setSensorID(++sensorId);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ private int GetIdentifierModule(String mylayer) {
+ int gid = -1;
+ if (mylayer.contains("1") || mylayer.contains("2"))
+ gid = 1;
+ if (mylayer.contains("3") || mylayer.contains("4"))
+ gid = 2;
+ if (mylayer.contains("5") || mylayer.contains("6"))
+ gid = 3;
+ if (mylayer.contains("7") || mylayer.contains("8"))
+ gid = 4;
+ if (mylayer.contains("9") || mylayer.contains("10"))
+ gid = 5;
+
+ return gid; //return top/bottom plates
+ }
+
+ private int GetIdentifierLayer(String mylayer) {
+ int gid = -1;
+ if (mylayer.contains("1"))
+ gid = 1;
+ if (mylayer.contains("2"))
+ gid = 2;
+ if (mylayer.contains("3"))
+ gid = 3;
+ if (mylayer.contains("4"))
+ gid = 4;
+ if (mylayer.contains("5"))
+ gid = 5;
+ if (mylayer.contains("6"))
+ gid = 6;
+ if (mylayer.contains("7"))
+ gid = 7;
+ if (mylayer.contains("8"))
+ gid = 8;
+ if (mylayer.contains("9"))
+ gid = 9;
+ if (mylayer.contains("10"))
+ gid = 10;
+
+ return gid; //return top/bottom plates
+ }
+}
java/trunk/hps-java/src/main/java/org/lcsim/hps/examples
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/examples/StarterAnalysisDriver.java 2013-12-03 16:59:02 UTC (rev 52)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/examples/StarterAnalysisDriver.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -14,9 +14,9 @@
import org.lcsim.event.base.BaseTrackState;
import org.lcsim.geometry.Detector;
import org.lcsim.geometry.compact.Field;
-import org.lcsim.hps.users.meeg.LCIOTrackAnalysis;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
+import org.lcsim.hps.recon.tracking.LCIOTrackAnalysis;
/*
* Example analysis driver.
java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/DataTrackerFakeHitDriver.java 2013-12-03 16:59:02 UTC (rev 52)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/DataTrackerFakeHitDriver.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -28,7 +28,6 @@
import org.lcsim.fit.helicaltrack.HitIdentifier;
import org.lcsim.geometry.Detector;
import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
-import org.lcsim.hps.users.phansson.WTrack;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType;
java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/LCIOTrackAnalysis.java (rev 0)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/LCIOTrackAnalysis.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -0,0 +1,167 @@
+package org.lcsim.hps.recon.tracking;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.identifier.Identifier;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.hps.recon.tracking.SvtUtils;
+
+/**
+ *
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: LCIOTrackAnalysis.java,v 1.3 2013/10/24 18:11:43 meeg Exp $
+ */
+public class LCIOTrackAnalysis {
+
+ protected Track track;
+ protected MCParticle _mcp = null;
+ protected double _purity;
+ protected int _nhits;
+ protected int _nbadhits;
+ private int _nAxialhits;
+ private int _nZhits;
+ protected boolean _hasLayerOne;
+ private int[] _nStripHitsPerLayer = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ protected Map<Integer, Hep3Vector> _hitLocationPerLayer = new HashMap<Integer, Hep3Vector>();
+ protected int _nhitsNew;
+
+ public Track getTrack() {
+ return track;
+ }
+
+ public LCIOTrackAnalysis(Track trk, RelationalTable hittomc, RelationalTable hittostrip, RelationalTable hittorotated) {
+ track = trk;
+
+ // Get the number of hits on the track
+ _nhits = trk.getTrackerHits().size();
+
+ // Create a map containing the number of hits for each MCParticle associated with the track
+ Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
+ _hasLayerOne = false;
+ // Loop over the hits on the track (HelicalTrackHits)
+ for (TrackerHit rotatedHit : trk.getTrackerHits()) {
+ TrackerHit hit = (TrackerHit) hittorotated.from(rotatedHit);
+ // get the set of MCParticles associated with this hit and update the hit count for each MCParticle
+ Set<MCParticle> mclist = hittomc.allFrom(hit);
+// System.out.println("MCParticle count: " + mclist.size());
+ for (MCParticle mcp : mclist) {
+ if (mcp != null) {
+// System.out.println(mcp.getOrigin());
+ Integer mchits = 0;
+ if (mcmap.containsKey(mcp)) {
+ mchits = mcmap.get(mcp);
+ }
+ mchits++;
+ mcmap.put(mcp, mchits);
+ }
+ }
+
+ Set<TrackerHit> hitlist = hittostrip.allFrom(hit);
+ for (TrackerHit cl : hitlist) {
+ int layer = -1;
+ int module = -1;
+ List<RawTrackerHit> rawHits = cl.getRawHits();
+// System.out.println("RawHits: " + rawHits.size());
+ for (RawTrackerHit rawHit : rawHits) {
+// System.out.println(rawHit.getCellID());
+ IIdentifier id = new Identifier(rawHit.getCellID());
+ int newLayer = SvtUtils.getInstance().getHelper().getValue(id, "layer");
+ if (layer != -1 && layer != newLayer) {
+ System.out.format("TrackerHit has hits from multiple layers: %d and %d\n", layer, newLayer);
+ }
+ layer = newLayer;
+ int newModule = SvtUtils.getInstance().getHelper().getValue(id, "module");
+ if (module != -1 && module != newModule) {
+ System.out.format("TrackerHit has hits from multiple modules: %d and %d\n", module, newModule);
+ }
+ module = newModule;
+// System.out.println(SvtUtils.getInstance().getHelper().getValue(id, "strip"));
+ }
+// System.out.format("layer %d, module %d\n", layer, module);
+ if (layer == 1) {
+ _hasLayerOne = true;
+ }
+
+
+ _nStripHitsPerLayer[layer - 1] = rawHits.size();
+ _hitLocationPerLayer.put(layer, new BasicHep3Vector(cl.getPosition()));
+ _nhitsNew++;
+
+ boolean isAxial = SvtUtils.getInstance().isAxial(SvtUtils.getInstance().getSensor(module, layer - 1));
+ if (isAxial) {
+ _nAxialhits++;
+ } else {
+ _nZhits++;
+
+ }
+ }
+ }
+
+ // Find the MCParticle that has the most hits on the track
+
+ int nbest = 0;
+ MCParticle mcbest = null;
+ for (MCParticle mcp : mcmap.keySet()) {
+ int count = mcmap.get(mcp);
+ if (count > nbest) {
+ nbest = count;
+ mcbest = mcp;
+ }
+ }
+
+ if (nbest > 0) {
+ _mcp = mcbest;
+ }
+ _purity = (double) nbest / (double) _nhits;
+ _nbadhits = _nhits - nbest;
+ }
+
+ public MCParticle getMCParticle() {
+ return _mcp;
+ }
+
+ public int getNHits() {
+ return _nhits;
+ }
+
+ public int getNBadHits() {
+ return _nbadhits;
+ }
+
+ public double getPurity() {
+ return _purity;
+ }
+
+ public int getNHitsNew() {
+ return _nhitsNew;
+ }
+
+ public int getNAxialHits() {
+ return _nAxialhits;
+ }
+
+ public int getNZHits() {
+ return _nZhits;
+ }
+
+ public boolean hasLayerOne() {
+ return _hasLayerOne;
+ }
+
+ public Hep3Vector getClusterPosition(Integer layer) {
+ return _hitLocationPerLayer.get(layer);
+ }
+
+ public int getNumberOfStripHits(int layer) {
+ return _nStripHitsPerLayer[layer - 1];
+ }
+}
java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/TrackUtils.java 2013-12-03 16:59:02 UTC (rev 52)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/TrackUtils.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -34,7 +34,6 @@
import org.lcsim.hps.event.HPSTransformations;
import org.lcsim.hps.recon.vertexing.HelixConverter;
import org.lcsim.hps.recon.vertexing.StraightLineTrack;
-import org.lcsim.hps.users.phansson.WTrack;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
import org.lcsim.util.swim.Helix;
java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/WTrack.java (rev 0)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/WTrack.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -0,0 +1,339 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.hps.recon.tracking;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.constants.Constants;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+
+/**
+ *
+ * @author phansson
+ */
+public class WTrack {
+
+ private boolean _debug = false;
+ public enum PARAM{TEST;}
+ private double[] _parameters = new double[7];
+ public HelicalTrackFit _htf = null;
+ private double _bfield;
+ private double _a;
+
+ static int max_iterations_intercept = 10;
+ static double epsilon_intercept = 1e-4;
+
+ public WTrack(WTrack trk) {
+ _bfield = trk._bfield;
+ _a = trk._a;
+ _parameters = trk._parameters;
+ _htf = trk._htf;
+ _debug = trk._debug;
+ }
+
+ public WTrack(HelicalTrackFit track, double bfield) {
+ initWithTrack(track, bfield, false);
+ }
+
+ public WTrack(HelicalTrackFit track, double bfield, boolean flip) {
+ initWithTrack(track, bfield, flip);
+ }
+ public void initWithTrack(HelicalTrackFit track, double bfield, boolean flip) {
+ _htf = track;
+ _bfield = flip ? -1.0*bfield : bfield; // flip if needed
+ _a = -1*Constants.fieldConversion*_bfield*Math.signum(track.R());
+ double p = track.p(Math.abs(_bfield));
+ double theta = Math.PI/2.0 - Math.atan(track.slope());
+ double phi = track.phi0();
+ _parameters[0] = p*Math.cos(phi)*Math.sin(theta);
+ _parameters[1] = p*Math.sin(phi)*Math.sin(theta);
+ _parameters[2] = p*Math.cos(theta);
+ _parameters[3] = Math.sqrt(_parameters[0]*_parameters[0]+_parameters[1]*_parameters[1]+_parameters[2]*_parameters[2]);
+ _parameters[4] = -1*track.dca()*Math.sin(phi); //x0
+ _parameters[5] = track.dca()*Math.cos(phi); //y0
+ _parameters[6] = track.z0(); //z0
+ if(_debug) {
+ System.out.printf("%s: WTrack initialized (p=%f,bfield=%f,theta=%f,phi=%f) from HelicalTrackFit:\n%s:%s\n",this.getClass().getSimpleName(),
+ p,_bfield,theta,phi,
+ this.getClass().getSimpleName(),this.toString());
+ }
+ }
+ public void setTrackParameters(double [] params) {
+ _parameters = params;
+ }
+
+ public double[] getParameters() {
+ return _parameters;
+ }
+
+ private boolean goingForward() {
+ // assuming the track should go in the x-direction --> not very general -> FIX THIS!?
+ return getP0().x()>0 ? true : false;
+ }
+
+
+ public double a() {
+ return _a;
+
+ }
+
+ private int getCharge() {
+ return (int) Math.signum(_htf.R());
+ }
+
+ public Hep3Vector getP0() {
+ return ( new BasicHep3Vector(_parameters[0],_parameters[1],_parameters[2]));
+ }
+
+ public Hep3Vector getX0() {
+ return ( new BasicHep3Vector(_parameters[4],_parameters[5],_parameters[6]));
+ }
+
+ public String paramsToString() {
+ String str = "";
+ for(int i=0;i<7;++i) str += _parameters[i] + ", ";
+ return str;
+ }
+ public String toString() {
+
+ String str = "WTrack params [" + paramsToString() + "]";
+ if(this._htf!=null) {
+ str += "\n with corresponding HelicalTrackFit:\n";
+ str += this._htf.toString();
+ }
+ return str;
+ }
+
+
+
+
+
+
+ private Hep3Vector getMomentumOnHelix(double s) {
+ WTrack track = this;
+ double a = track.a();
+ Hep3Vector p0 = track.getP0();
+ double rho = a / p0.magnitude();
+ double px = p0.x()*Math.cos(rho*s) - p0.y()*Math.sin(rho*s);
+ double py = p0.y()*Math.cos(rho*s) + p0.x()*Math.sin(rho*s);
+ double pz = p0.z();
+ return (new BasicHep3Vector(px,py,pz));
+ }
+
+ private Hep3Vector getPointOnHelix(double s) {
+ WTrack track = this;
+ double a = track.a();
+ Hep3Vector p0 = track.getP0();
+ Hep3Vector x0 = track.getX0();
+ double rho = a / p0.magnitude();
+ double x = x0.x() + p0.x()/a*Math.sin(rho*s) - p0.y()/a*(1-Math.cos(rho*s));
+ double y = x0.y() + p0.y()/a*Math.sin(rho*s) + p0.x()/a*(1-Math.cos(rho*s));
+ double z = x0.z() + p0.z()/p0.magnitude()*s;
+ return (new BasicHep3Vector(x,y,z));
+ }
+
+ private double getPathLengthToPlaneApprox(Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
+ /*
+ * Find the approximate path length to the point xp
+ * in arbitrary oriented, constant magnetic field with unit vector h
+ */
+ WTrack track = this;
+ double a = track.a();
+ Hep3Vector p0 = track.getP0();
+ Hep3Vector x0 = track.getX0();
+ double p = p0.magnitude();
+ double rho = a / p;
+ double A = VecOp.dot(eta,VecOp.cross(p0, h))/p*0.5*rho;
+ double B = VecOp.dot(p0,eta)/p;
+ double C = VecOp.dot(VecOp.sub(x0,xp),eta);
+ double t = B*B-4*A*C;
+ if(t<0) {
+ System.out.println(" getPathLengthToPlaneApprox ERROR t is negative (" + t + ")!" );
+ System.out.println(" p " + p + " rho " + rho + " a " + a + " A " + A + " B " + B + " C " + C);
+ System.out.println(" track params: " + track.paramsToString());
+ System.out.println(" xp " + xp.toString());
+ System.out.println(" eta " + eta.toString());
+ System.out.println(" h " + h.toString());
+ System.exit(1);
+ }
+ double root1 = (-B + Math.sqrt(t)) /(2*A);
+ double root2 = (-B - Math.sqrt(t)) /(2*A);
+
+ // choose the smallest positive solution
+ // if both negative choose the smallest negative ???
+ //if(root1==0 || root2==0) root=0;
+ double root = Math.abs(root1) <= Math.abs(root2) ? root1 : root2;
+// else if(Math.signum(root1)>0 && Math.signum(root2)<0) root = root1;
+// else if(Math.signum(root2)>0 && Math.signum(root1)<0) root = root2;
+// else if(Math.signum(root1)>0 && Math.signum(root2)>0) root = root1 > root2 ? root2 : root1;
+// else if(Math.signum(root1)<0 && Math.signum(root2)<0) root = root1 < root2 ? root2 : root1;
+// else {
+// System.out.println(" I should never get here! (root1=" + root1 + " root2=" + root2+")");
+// System.exit(1);
+// }
+ if(_debug) {
+ System.out.println(" getPathLengthToPlaneApprox ");
+ System.out.println(" " + track.paramsToString());
+ System.out.println(" xp " + xp.toString());
+ System.out.println(" eta " + eta.toString());
+ System.out.println(" h " + h.toString());
+ System.out.println(" p " + p + " rho " + rho + " t " + t + " A " + A + " B " + B + " C " + C);
+ System.out.println(" root1 " + root1 + " root2 " + root2 + " -> root " + root);
+ }
+ return root;
+
+ }
+
+
+ private Hep3Vector getPointOnHelix(double s, Hep3Vector h) {
+ /*
+ * Get point on helix at path lenght s
+ * in arbitrary oriented, constant magnetic field with unit vector h
+ */
+ WTrack track = this;
+ double a = track.a();
+ Hep3Vector p0 = track.getP0();
+ double p = p0.magnitude();
+ Hep3Vector x0 = track.getX0();
+ double rho = a / p0.magnitude();
+ double srho = s*rho;
+ Hep3Vector a1 = VecOp.mult(1/a*Math.sin(srho), p0);
+ Hep3Vector a2 = VecOp.mult(1/a*(1-Math.cos(srho)),VecOp.cross(p0,h));
+ Hep3Vector a3 = VecOp.mult(VecOp.dot(p0, h)/p,h);
+ Hep3Vector a4 = VecOp.mult(s-Math.sin(srho)/rho, a3);
+ Hep3Vector x = VecOp.add(x0,a1);
+ x = VecOp.sub(x,a2);
+ x = VecOp.add(x,a4);
+ return x;
+ }
+
+ private Hep3Vector getMomentumOnHelix(double s, Hep3Vector h) {
+ /*
+ * Get point on helix at path lenght s
+ * in arbitrary oriented, constant magnetic field with unit vector h
+ */
+ WTrack track = this;
+ double a = track.a();
+ Hep3Vector p0 = track.getP0();
+ double rho = a / p0.magnitude();
+ double srho = s*rho;
+ Hep3Vector a1 = VecOp.mult(Math.cos(srho), p0);
+ Hep3Vector a2 = VecOp.cross(p0, VecOp.mult(Math.sin(srho),h));
+ Hep3Vector a3 = VecOp.mult(VecOp.dot(p0,h),VecOp.mult(1-Math.cos(srho),h));
+ Hep3Vector p = VecOp.sub(a1,a2);
+ p = VecOp.add(p,a3);
+ return p;
+ }
+
+
+ private double[] getHelixParametersAtPathLength(double s, Hep3Vector h) {
+ /*
+ * Calculate the exact position of the new helix parameters at path length s
+ * in an arbitrarily oriented, constant magnetic field
+ *
+ * point xp is the point
+ * h is a unit vector in the direction of the magnetic field
+ */
+
+ // Find track parameters at that path length
+ Hep3Vector p = getMomentumOnHelix(s, h);
+ Hep3Vector x = getPointOnHelix(s, h);
+
+ Hep3Vector p_tmp = getMomentumOnHelix(s);
+ Hep3Vector x_tmp = getPointOnHelix(s);
+
+ if(_debug) {
+ System.out.println(" point on helix at s");
+ System.out.println(" p " + p.toString() + " p_tmp " + p_tmp.toString());
+ System.out.println(" x " + x.toString() + " x_tmp " + x_tmp.toString());
+ }
+
+
+ //Create the new parameter array
+ double [] pars = new double[7];
+ pars[0] = p.x();
+ pars[1] = p.y();
+ pars[2] = p.z();
+ pars[3] = getParameters()[3]; //E is unchanged
+ pars[4] = x.x();
+ pars[5] = x.y();
+ pars[6] = x.z();
+ return pars;
+ }
+
+
+ public Hep3Vector getHelixAndPlaneIntercept(Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
+
+ /*
+ * Find the interception point between the helix and plane
+ * xp: point on the plane
+ * eta: unit vector of the plane
+ * h: unit vector of magnetic field
+ */
+
+
+ int iteration = 1;
+ double s_total = 0.;
+ double step = 9999999.9;
+ //List<WTrack> tracks = new ArrayList<WTrack>();
+ WTrack trk = this;
+ while(iteration<=max_iterations_intercept && Math.abs(step)>epsilon_intercept) {
+
+ if(_debug) {
+ System.out.printf("%s: Iteration %d\n", this.getClass().getSimpleName(),iteration);
+ System.out.printf("%s: s_total %f prev_step %.3f current trk params: %s \n",
+ this.getClass().getSimpleName(),s_total,step,trk.paramsToString());
+ }
+
+ // check that the track is not looping
+
+ if(trk.goingForward()) {
+
+
+ // Start by approximating the path length to the point
+ step = getPathLengthToPlaneApprox(xp, eta, h);
+
+ if(_debug) System.out.printf("%s: path length step s=%.3f\n",this.getClass().getSimpleName(),step);
+
+ // Find the track parameters at this point
+ double[] params = getHelixParametersAtPathLength(step, h);
+
+ // update the track parameters
+ trk.setTrackParameters(params);
+
+ if(_debug) System.out.printf("%s: updated track params: [%s]\n",this.getClass().getSimpleName(),trk.paramsToString());
+
+ //tracks.add(trk);
+ iteration++;
+ s_total += step;
+
+ //Save distance between point and estimate
+ //Hep3Vector dpoint = VecOp.sub(xp, trk.getX0());
+
+ } else {
+ //if(_debug)
+ System.out.printf("%s: this track started to go backwards?! params [%s]\n",this.getClass().getSimpleName(),trk.toString());
+ return null;
+ }
+
+
+ }
+
+ if(_debug) System.out.printf("%s: final total_s=%f with final step %f after %d iterations gave track params: %s\n",
+ this.getClass().getSimpleName(),s_total,step,iteration,trk.paramsToString());
+
+ return trk.getX0();
+
+ }
+
+
+
+
+
+
+}
java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
--- java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl/GBLFileIO.java 2013-12-03 16:59:02 UTC (rev 52)
+++ java/trunk/hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl/GBLFileIO.java 2013-12-03 17:08:28 UTC (rev 53)
@@ -16,7 +16,7 @@
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.hps.recon.tracking.gbl.GBLOutput.ClParams;
import org.lcsim.hps.recon.tracking.gbl.GBLOutput.PerigeeParams;
-import org.lcsim.hps.users.mgraham.alignment.RunAlignment;
+import org.lcsim.hps.alignment.RunAlignment;
/**
* Handles text file printing for the GBL
SVNspam 0.1