Author: [log in to unmask] Date: Thu Feb 19 10:04:04 2015 New Revision: 2166 Log: Adding new HelicalTrackStrip class that encapsulates the lcsim version. This fixes a problem where the local oordinate system for the sensor is flipped to always have the normal of the plane to point along the track. Added: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HelicalTrackStripGbl.java Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java java/trunk/users/src/main/java/org/hps/users/phansson/StripMPAlignmentInput.java Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java (original) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java Thu Feb 19 10:04:04 2015 @@ -12,6 +12,7 @@ import java.util.List; import java.util.Map; +import org.hps.recon.tracking.gbl.HelicalTrackStripGbl; import org.lcsim.detector.ITransform3D; import org.lcsim.detector.solids.Box; import org.lcsim.detector.solids.Point3D; @@ -149,7 +150,7 @@ * @param bfield - magnetic field value * @return point at intercept */ - public static Hep3Vector getHelixPlaneIntercept(HelicalTrackFit helfit, HelicalTrackStrip strip, double bfield) { + public static Hep3Vector getHelixPlaneIntercept(HelicalTrackFit helfit, HelicalTrackStripGbl strip, double bfield) { Hep3Vector point_on_plane = strip.origin(); Hep3Vector unit_vec_normal_to_plane = VecOp.cross(strip.u(), strip.v());// strip.w(); Hep3Vector intercept_point = getHelixPlaneIntercept(helfit, unit_vec_normal_to_plane, point_on_plane, bfield); @@ -463,7 +464,7 @@ return residuals; } - public static Map<String, Double> calculateLocalTrackHitResiduals(Track track, HelicalTrackHit hth, HelicalTrackStrip strip, double bFieldInZ) { + public static Map<String, Double> calculateLocalTrackHitResiduals(Track track, HelicalTrackHit hth, HelicalTrackStripGbl strip, double bFieldInZ) { SeedTrack st = (SeedTrack) track; SeedCandidate seed = st.getSeedCandidate(); @@ -474,7 +475,7 @@ return calculateLocalTrackHitResiduals(_trk, strip, msdrdphi, msdz, bFieldInZ); } - public static Map<String, Double> calculateLocalTrackHitResiduals(HelicalTrackFit _trk, HelicalTrackStrip strip, double msdrdphi, double msdz, double bFieldInZ) { + public static Map<String, Double> calculateLocalTrackHitResiduals(HelicalTrackFit _trk, HelicalTrackStripGbl strip, double msdrdphi, double msdz, double bFieldInZ) { boolean debug = false; boolean includeMS = true; @@ -508,7 +509,7 @@ Hep3Vector vdiffTrk = VecOp.sub(trkpos, corigin); TrackerHitUtils thu = new TrackerHitUtils(debug); - Hep3Matrix trkToStrip = thu.getTrackToStripRotation(strip); + Hep3Matrix trkToStrip = thu.getTrackToStripRotation(strip.getStrip()); Hep3Vector vdiff = VecOp.mult(trkToStrip, vdiffTrk); double umc = vdiff.x(); Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java (original) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java Thu Feb 19 10:04:04 2015 @@ -256,8 +256,8 @@ HelicalTrackCross htc = (HelicalTrackCross) hit; List<HelicalTrackStrip> strips = htc.getStrips(); - for(HelicalTrackStrip strip : strips) { - + for(HelicalTrackStrip stripOld : strips) { + HelicalTrackStripGbl strip = new HelicalTrackStripGbl(stripOld, true); if(_debug>0) System.out.printf("%s: layer %d\n",this.getClass().getSimpleName(),strip.layer()); if(textFile != null) { @@ -379,7 +379,7 @@ Hep3Vector vdiffTrkTruth = VecOp.sub(trkposTruth, origin); // then find the rotation from tracking to measurement frame - Hep3Matrix trkToStripRot = _trackerHitUtils.getTrackToStripRotation(strip); + Hep3Matrix trkToStripRot = _trackerHitUtils.getTrackToStripRotation(strip.getStrip()); // then rotate that vector into the measurement frame to get the predicted measurement position Hep3Vector trkpos_meas = VecOp.mult(trkToStripRot, vdiffTrk); @@ -397,6 +397,15 @@ //GBLDATA stripData.setMeas(strip.umeas()); stripData.setTrackPos(trkpos_meas); + + if(_debug>1) { + System.out.printf("%s: rotation matrix to meas frame\n%s\n", getClass().getSimpleName(), VecOp.toString(trkToStripRot)); + System.out.printf("%s: tPosGlobal %s origin %s\n", getClass().getSimpleName(), trkpos.toString(), origin.toString()); + System.out.printf("%s: tDiff %s\n", getClass().getSimpleName(), vdiffTrk.toString()); + System.out.printf("%s: tPosMeas %s\n", getClass().getSimpleName(), trkpos_meas.toString()); + } + + // residual in measurement frame Hep3Vector res_meas = VecOp.sub(m_meas, trkpos_meas); @@ -430,14 +439,14 @@ } // find scattering angle - ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()); + ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit)strip.getStrip().rawhits().get(0)).getDetectorElement()); double scatAngle; if(scatter != null) { scatAngle = scatter.getScatterAngle().Angle(); } else { - 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()); + System.out.printf("%s: WARNING cannot find scatter for detector %s with strip cluster at %s\n",this.getClass(),((RawTrackerHit)strip.getStrip().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! DetectorPlane closest = null; double dx = 999999.9; @@ -937,7 +946,6 @@ } } - - + } Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java (original) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java Thu Feb 19 10:04:04 2015 @@ -22,7 +22,7 @@ public static final int ID = 0; public static final int BANK_INT_SIZE = 1; } - private static class GBLDOUBLE { + public static class GBLDOUBLE { public static final int PERKAPPA =0; public static final int PERTHETA = 1; public static final int PERPHI = 2; @@ -84,7 +84,8 @@ this.bank_double[GBLDOUBLE.PERD0] = perPar.getD0(); this.bank_double[GBLDOUBLE.PERZ0] = perPar.getZ0(); } - + + public void setPrjPerToCl(int row, int col, double val) { int idx = col + row*3; Added: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HelicalTrackStripGbl.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HelicalTrackStripGbl.java (added) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HelicalTrackStripGbl.java Thu Feb 19 10:04:04 2015 @@ -0,0 +1,100 @@ +package org.hps.recon.tracking.gbl; + +import hep.physics.vec.Hep3Matrix; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp; + +import org.hps.recon.tracking.CoordinateTransformations; +import org.lcsim.detector.IDetectorElement; +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.fit.helicaltrack.HelicalTrackStrip; + +/** + * Encapsulate @HelicalTrackStrip to make sure that unit vectors are based on geometry. + * TODO should extend as a permanent solution. + * + * @author Per Hansson Adrian <[log in to unmask]> + * + */ +public class HelicalTrackStripGbl { + private HelicalTrackStrip _strip; + private SiSensorElectrodes _electrodes = null; + private Hep3Matrix _electrodesToTracking = null; + private Hep3Vector _u = null; + private Hep3Vector _v = null; + private Hep3Vector _w = null; + private boolean _useGeomDef = false; + public HelicalTrackStripGbl(HelicalTrackStrip strip, boolean useGeomDef) { + _strip = strip; + _useGeomDef = useGeomDef; + } + public double du() { + return _strip.du(); + } + public double vmin() { + return _strip.vmin(); + } + public double vmax() { + return _strip.vmax(); + } + public double umeas() { + return _strip.umeas(); + } + public HelicalTrackStrip getStrip() { + return _strip; + } + public Hep3Vector origin() { + return _strip.origin(); + } + public int layer() { + return _strip.layer(); + } + private SiSensorElectrodes getElectrodes() { + if(_electrodes==null) { + RawTrackerHit rth = (RawTrackerHit) _strip.rawhits().get(0); + IDetectorElement ide = rth.getDetectorElement(); + SiSensor sensor = ide.findDescendants(SiSensor.class).get(0); + _electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE); + } + return _electrodes; + } + private Hep3Matrix getElectrodeToTrackingMatrix() { + if(_electrodesToTracking==null) { + SiSensorElectrodes electrodes = getElectrodes(); + _electrodesToTracking = VecOp.mult(CoordinateTransformations.getMatrix(), electrodes.getLocalToGlobal().getRotation().getRotationMatrix()); + } + return _electrodesToTracking; + } + public Hep3Vector u() { + if(_u == null) { + if(_useGeomDef) { + _u = VecOp.mult(getElectrodeToTrackingMatrix(), getElectrodes().getMeasuredCoordinate(0)); + } + else { + _u = _strip.u(); + } + } + return _u; + } + public Hep3Vector v() { + if(_v == null) { + if(_useGeomDef) { + _v = VecOp.mult(getElectrodeToTrackingMatrix(), getElectrodes().getUnmeasuredCoordinate(0)); + } else { + _v = _strip.v(); + } + } + return _v; + } + public Hep3Vector w() { + if(_w==null) { + _w = VecOp.cross(u(), v()); + } + return _w; + } + + +} Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java (original) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java Thu Feb 19 10:04:04 2015 @@ -160,7 +160,7 @@ } for (int istrip = 0; istrip < stripClusters.size(); ++istrip) { - HelicalTrackStrip strip = stripClusters.get(istrip); + HelicalTrackStripGbl strip = new HelicalTrackStripGbl(stripClusters.get(istrip), true); if (_debug) { System.out.printf(" layer %d origin %s\n", strip.layer(), strip.origin().toString()); @@ -235,7 +235,7 @@ Hep3Vector vdiffTrk = VecOp.sub(trkpos, strip.origin()); // then find the rotation from tracking to measurement frame - Hep3Matrix trkToStripRot = _trackHitUtils.getTrackToStripRotation(strip); + Hep3Matrix trkToStripRot = _trackHitUtils.getTrackToStripRotation(strip.getStrip()); // then rotate that vector into the measurement frame to get the predicted measurement position Hep3Vector trkpos_meas = VecOp.mult(trkToStripRot, vdiffTrk); @@ -300,13 +300,13 @@ BasicMatrix scat = GblUtils.getInstance().zeroMatrix(0, 2); // find scattering angle - ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit) strip.rawhits().get(0)).getDetectorElement()); + ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement()); double scatAngle; if (scatter != null) { scatAngle = scatter.getScatterAngle().Angle(); } else { - 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()); + System.out.printf("%s: WARNING cannot find scatter for detector %s with strip cluster at %s\n", this.getClass(), ((RawTrackerHit) strip.getStrip().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 //TODO check if this really makes sense to do DetectorPlane closest = null; Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java (original) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java Thu Feb 19 10:04:04 2015 @@ -11,8 +11,10 @@ import java.util.HashMap; import java.util.List; import java.util.Map; +import java.util.logging.Level; import java.util.logging.Logger; +import org.hps.recon.tracking.gbl.GBLTrackData.GBLDOUBLE; import org.hps.recon.tracking.gbl.matrix.Matrix; import org.hps.recon.tracking.gbl.matrix.SymMatrix; import org.hps.recon.tracking.gbl.matrix.Vector; @@ -37,7 +39,8 @@ */ public class HpsGblRefitter extends Driver { - private static Logger logger = LogUtil.create(HpsGblRefitter.class); + //private static Logger logger = LogUtil.create(HpsGblRefitter.class); + private Logger logger = getLogger(); private boolean _debug = false; private final String trackCollectionName = "MatchedTracks"; private final String track2GblTrackRelationName = "TrackToGBLTrack"; @@ -157,6 +160,18 @@ List<FittedGblTrajectory> trackFits = new ArrayList<FittedGblTrajectory>(); int trackNum = 0; for (GBLTrackData t : stripsGblMap.keySet()) { + + if(Math.cos(t.getDoubleVal(GBLDOUBLE.PERTHETA))>0) { + //if(_debug) + logger.info(" top track"); + continue; + + } else { + logger.info(" bot track"); + + } + + FittedGblTrajectory traj = fit(stripsGblMap.get(t), bfac); ++trackNum; if(traj!=null) { @@ -409,11 +424,11 @@ addDer.set(0, i, milleParameters.get(i).getValue()); } point.addGlobals(labGlobal, addDer); - + String logders = ""; for(int i=0; i < milleParameters.size(); ++i) { - logger.info(labGlobal.get(i) + "\t" + addDer.get(0, i)); + logders += labGlobal.get(i) + "\t" + addDer.get(0, i) + "\n"; } - + logger.info("\n"+ logders); /* @@ -556,6 +571,9 @@ _dm_dg = getMeasDers(); // Calculate, by chain rule, derivatives of residuals w.r.t. global parameters _dr_dg = _dr_dm.times(_dm_dg); + + //logger.log(Level.FINER," dr_dm\n"+ _dr_dm.toString() + "\ndm_dg\n" + _dm_dg.toString() + "\ndr_dg\n" +_dr_dg.toString()); + //logger.info("loglevel " + logger.getLevel().toString()); //print 'dm_dg' //print dm_dg //print 'dr_dm' @@ -625,6 +643,7 @@ double tdotn = VecOp.dot(_t, _n); Matrix dr_dm = Matrix.identity(3,3); //print 't ', self.t, ' n ', self.n, ' dot(t,n) ', tdotn + //logger.info("t " + _t.toString() +" n " + _n.toString() + " dot(t,n) " + tdotn); double delta, val; for(int i=0; i<3; ++i) { for(int j=0; j<3; ++j) { Modified: java/trunk/users/src/main/java/org/hps/users/phansson/StripMPAlignmentInput.java ============================================================================= --- java/trunk/users/src/main/java/org/hps/users/phansson/StripMPAlignmentInput.java (original) +++ java/trunk/users/src/main/java/org/hps/users/phansson/StripMPAlignmentInput.java Thu Feb 19 10:04:04 2015 @@ -20,8 +20,10 @@ import java.util.Map; + //===> import org.hps.conditions.deprecated.SvtUtils; import org.hps.recon.tracking.TrackUtils; +import org.hps.recon.tracking.gbl.HelicalTrackStripGbl; import org.lcsim.detector.tracker.silicon.HpsSiSensor; import org.lcsim.event.RawTrackerHit; import org.lcsim.event.Track; @@ -126,7 +128,8 @@ if(_DEBUG) { System.out.printf("%s: This hit has %d clusters msdrphi=%.4f msdz=%.4f\n",this.getClass().getSimpleName(),clusterlist.size(),msdrphi,msdz); } - for (HelicalTrackStrip cl : clusterlist) { + for (HelicalTrackStrip strip : clusterlist) { + HelicalTrackStripGbl cl = new HelicalTrackStripGbl(strip,false); CalculateResidual(cl, msdrphi, msdz); CalculateLocalDerivatives(cl); CalculateGlobalDerivatives(cl); @@ -140,12 +143,12 @@ } - private void CalculateLocalDerivatives(HelicalTrackStrip strip) { + private void CalculateLocalDerivatives(HelicalTrackStripGbl strip) { double xint = strip.origin().x(); BasicMatrix dfdqGlobalOld = _oldAlignUtils.calculateLocalHelixDerivatives(_trk, xint); BasicMatrix dfdqGlobal = _alignUtils.calculateLocalHelixDerivatives(_trk, xint); BasicMatrix dfdqGlobalNum = this._numDerivatives.calculateLocalHelixDerivatives(_trk, xint); - Hep3Matrix trkToStrip = trackerHitUtil.getTrackToStripRotation(strip); + Hep3Matrix trkToStrip = trackerHitUtil.getTrackToStripRotation(strip.getStrip()); _dfdq = (BasicMatrix) MatrixOp.mult(trkToStrip, dfdqGlobal); if (_DEBUG) { @@ -175,7 +178,7 @@ - private void CalculateGlobalDerivatives(HelicalTrackStrip strip) { + private void CalculateGlobalDerivatives(HelicalTrackStripGbl strip) { /* * Residual in local sensor frame is defined as r = m_a-p @@ -198,7 +201,7 @@ TrackDirection trkdir = HelixUtils.CalculateTrackDirection(_trk, pathLengthToInterception); Hep3Vector t_TRACK = trkdir.Direction(); Hep3Vector n_TRACK = strip.w(); - Hep3Matrix T = trackerHitUtil.getTrackToStripRotation(strip); + Hep3Matrix T = trackerHitUtil.getTrackToStripRotation(strip.getStrip()); Hep3Vector t = VecOp.mult(T, t_TRACK); Hep3Vector n = VecOp.mult(T, n_TRACK); @@ -287,7 +290,7 @@ int layer = strip.layer(); - HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) strip.rawhits().get(0)).getDetectorElement(); + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement(); //===> int side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? 10000 : 20000; int side = sensor.isTopLayer() ? 10000 : 20000; @@ -376,7 +379,7 @@ - private void CalculateResidual(HelicalTrackStrip strip, double msdrdphi, double msdz) { + private void CalculateResidual(HelicalTrackStripGbl strip, double msdrdphi, double msdz) { if(_DEBUG) System.out.printf("%s: --- CalculateResidual ---\n",this.getClass().getSimpleName()); double bfield = Math.abs(this._bfield.z()); @@ -391,7 +394,7 @@ _error[2] = res_local.get("wreserr"); double vdiffy = res_local.get("vdiffTrky"); - HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) strip.rawhits().get(0)).getDetectorElement(); + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement(); String side = sensor.isTopLayer() ? "top" : "bottom"; //===> String side = SvtUtils.getInstance().isTopLayer((SiSensor)((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom"; //Fill residuals vs distrance from center of plane in the v directions @@ -496,9 +499,9 @@ - private void PrintStripResiduals(HelicalTrackStrip strip) { - - HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) strip.rawhits().get(0)).getDetectorElement(); + private void PrintStripResiduals(HelicalTrackStripGbl strip) { + + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement(); String side = sensor.isTopLayer() ? "top" : "bottom"; //===> String side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom"; @@ -834,7 +837,7 @@ - private void CalculateGlobalDerivativesWRONG(HelicalTrackStrip strip) { + private void CalculateGlobalDerivativesWRONG(HelicalTrackStripGbl strip) { //Find interception with plane that the strips belongs to Hep3Vector x_vec = TrackUtils.getHelixPlaneIntercept(_trk, strip, Math.abs(this._bfield.z())); @@ -858,7 +861,7 @@ Hep3Vector corigin = strip.origin(); double vmeas = corigin.y(); //THis is wrong int layer = strip.layer(); - HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) strip.rawhits().get(0)).getDetectorElement(); + HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement(); int side = sensor.isTopLayer() ? 10000 : 20000; //===> int side = SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? 10000 : 20000; @@ -887,7 +890,7 @@ //Rotation matrix from the tracking fram to the sensor/strip frame - Hep3Matrix T = trackerHitUtil.getTrackToStripRotation(strip); + Hep3Matrix T = trackerHitUtil.getTrackToStripRotation(strip.getStrip()); @@ -1069,7 +1072,7 @@ //Put it into a matrix to be able to transform easily //BasicMatrix _dqp_da_TRACK = FillMatrix(dqp_da_TRACK, 3, 1); //Get transformation matrix from tracking frame to sensor frame where the residuals are calculated - Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip); + Hep3Matrix trkToStrip = this.trackerHitUtil.getTrackToStripRotation(strip.getStrip()); if(_DEBUG) { System.out.println("Final transformation from tracking frame to strip frame:\n" + trkToStrip.toString()); }