Author: [log in to unmask]
Date: Thu Oct 1 18:39:46 2015
New Revision: 3749
Log:
more cleanup to GBL
Modified:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.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/GBLOutputDriver.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.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/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.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 Oct 1 18:39:46 2015
@@ -994,16 +994,12 @@
}
public static double getTrackTime(Track track, RelationalTable hitToStrips, RelationalTable hitToRotated) {
- int nStrips = 0;
double meanTime = 0;
- for (TrackerHit hit : track.getTrackerHits()) {
- Set<TrackerHit> htsList = hitToStrips.allFrom(hitToRotated.from(hit));
- for (TrackerHit hts : htsList) {
- nStrips++;
- meanTime += hts.getTime();
- }
- }
- meanTime /= nStrips;
+ List<TrackerHit> stripHits = getStripHits(track, hitToStrips, hitToRotated);
+ for (TrackerHit hit : stripHits) {
+ meanTime += hit.getTime();
+ }
+ meanTime /= stripHits.size();
return meanTime;
}
@@ -1280,4 +1276,8 @@
}
return null;
}
+
+ public static Hep3Vector getBField(Detector detector) {
+ return detector.getFieldMap().getField(new BasicHep3Vector(0., 0., 500.0));
+ }
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/FittedGblTrajectory.java Thu Oct 1 18:39:46 2015
@@ -56,12 +56,6 @@
_ndf = ndf;
_lost = lost;
}
- public void set_track_data(GBLTrackData t) {
- _t = t;
- }
- public GBLTrackData get_track_data() {
- return _t;
- }
public void set_seed(Track seed) {
_seed = seed;
}
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 Oct 1 18:39:46 2015
@@ -13,7 +13,6 @@
import java.util.HashMap;
import java.util.List;
import java.util.Map;
-import org.apache.commons.math3.util.Pair;
import org.hps.recon.tracking.CoordinateTransformations;
import org.hps.recon.tracking.MaterialSupervisor;
import org.hps.recon.tracking.MultipleScattering;
@@ -21,15 +20,8 @@
import org.hps.recon.tracking.MultipleScattering.ScatterPoints;
import org.hps.recon.tracking.TrackUtils;
import org.hps.recon.tracking.TrackerHitUtils;
-import org.lcsim.detector.DetectorIdentifierHelper;
import org.lcsim.detector.IDetectorElement;
-import org.lcsim.detector.ITransform3D;
-import org.lcsim.detector.identifier.IIdentifier;
-import org.lcsim.detector.identifier.IIdentifierHelper;
-import org.lcsim.detector.tracker.silicon.ChargeCarrier;
import org.lcsim.detector.tracker.silicon.HpsSiSensor;
-import org.lcsim.detector.tracker.silicon.SiSensor;
-import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
import org.lcsim.event.MCParticle;
import org.lcsim.event.RawTrackerHit;
import org.lcsim.event.SimTrackerHit;
@@ -468,7 +460,7 @@
if (_debug > 0) {
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());
}
- scatAngle = GblUtils.estimateScatter(strip, htf, _scattering, _B.magnitude());
+ scatAngle = GblUtils.estimateScatter(sensor, htf, _scattering, _B.magnitude());
}
//GBLDATA
@@ -486,223 +478,6 @@
addBeamspotToHitList(htf, stripClusterDataList, istrip, htfTruth);
}
- }
-
- static List<GBLStripClusterData> printGBL(HelicalTrackFit htf, List<SiTrackerHitStrip1D> hitsOnTrack, MultipleScattering _scattering, double _B, int _debug) {
- List<GBLStripClusterData> stripClusterDataList = new ArrayList<GBLStripClusterData>();
-
- // Find scatter points along the path
- ScatterPoints scatters = _scattering.FindHPSScatterPoints(htf);
-
- if (_debug > 0) {
- System.out.printf("perPar covariance matrix\n%s\n", htf.covariance().toString());
- }
-
- for (SiTrackerHitStrip1D hitOnTrack : hitsOnTrack) {
- HelicalTrackStripGbl strip = new HelicalTrackStripGbl(makeDigiStrip(hitOnTrack), true);
-
- // find Millepede layer definition from DetectorElement
- IDetectorElement de = ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement();
- HpsSiSensor sensor;
- if (de instanceof HpsSiSensor) {
- sensor = (HpsSiSensor) de;
- } else {
- throw new ClassCastException("Detector element " + de.getName() + " couldn't be casted to " + HpsSiSensor.class.getName());
- }
-
- int millepedeId = sensor.getMillepedeId();
-
- if (_debug > 0) {
- System.out.printf("layer %d millepede %d (DE=\"%s\", origin %s) \n", strip.layer(), millepedeId, sensor.getName(), strip.origin().toString());
- }
-
- //Center of the sensor
- Hep3Vector origin = strip.origin();
-
- //Find intercept point with sensor in tracking frame
- Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B));
- if (trkpos == null) {
- if (_debug > 0) {
- System.out.println("Can't find track intercept; use sensor origin");
- }
- trkpos = strip.origin();
- }
- if (_debug > 0) {
- System.out.printf("trkpos at intercept [%.10f %.10f %.10f]\n", trkpos.x(), trkpos.y(), trkpos.z());
- }
-
- //GBLDATA
- GBLStripClusterData stripData = new GBLStripClusterData(millepedeId);
- //Add to output list
- stripClusterDataList.add(stripData);
-
- //path length to intercept
- double s = HelixUtils.PathToXPlane(htf, trkpos.x(), 0, 0).get(0);
- double s3D = s / Math.cos(Math.atan(htf.slope()));
-
- //GBLDATA
- stripData.setPath(s);
- stripData.setPath3D(s3D);
-
- //GBLDATA
- stripData.setU(strip.u());
- stripData.setV(strip.v());
- stripData.setW(strip.w());
-
- //Print track direction at intercept
- Hep3Vector tDir = HelixUtils.Direction(htf, s);
- double phi = htf.phi0() - s / htf.R();
- double lambda = Math.atan(htf.slope());
-
- //GBLDATA
- stripData.setTrackDir(tDir);
- stripData.setTrackPhi(phi);
- stripData.setTrackLambda(lambda);
-
- //Print residual in measurement system
- // start by find the distance vector between the center and the track position
- Hep3Vector vdiffTrk = VecOp.sub(trkpos, origin);
-
- // then find the rotation from tracking to measurement frame
- Hep3Matrix trkToStripRot = getTrackToStripRotation(strip.getStrip());
-
- // then rotate that vector into the measurement frame to get the predicted measurement position
- Hep3Vector trkpos_meas = VecOp.mult(trkToStripRot, vdiffTrk);
-
- //GBLDATA
- stripData.setMeas(strip.umeas());
- stripData.setTrackPos(trkpos_meas);
- stripData.setMeasErr(strip.du());
-
- if (_debug > 1) {
- System.out.printf("rotation matrix to meas frame\n%s\n", VecOp.toString(trkToStripRot));
- System.out.printf("tPosGlobal %s origin %s\n", trkpos.toString(), origin.toString());
- System.out.printf("tDiff %s\n", vdiffTrk.toString());
- System.out.printf("tPosMeas %s\n", trkpos_meas.toString());
- }
-
- if (_debug > 0) {
- System.out.printf("layer %d millePedeId %d uRes %.10f\n", strip.layer(), millepedeId, stripData.getMeas() - stripData.getTrackPos().x());
- }
-
- // find scattering angle
- ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement());
- double scatAngle;
-
- if (scatter != null) {
- scatAngle = scatter.getScatterAngle().Angle();
- } else {
- if (_debug > 0) {
- System.out.printf("WARNING cannot find scatter for detector %s with strip cluster at %s\n", ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement().getName(), strip.origin().toString());
- }
- scatAngle = GblUtils.estimateScatter(strip, htf, _scattering, _B);
- }
-
- //GBLDATA
- stripData.setScatterAngle(scatAngle);
- }
- return stripClusterDataList;
- }
-
- private static Hep3Matrix getTrackToStripRotation(HelicalTrackStrip strip) {
- // This function transforms the hit to the sensor coordinates
-
- // Transform from JLab frame to sensor frame (done through the RawTrackerHit)
- 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);
- ITransform3D detToStrip = electrodes.getGlobalToLocal();
- // Get rotation matrix
- Hep3Matrix detToStripMatrix = detToStrip.getRotation().getRotationMatrix();
- // Transformation between the JLAB and tracking coordinate systems
- Hep3Matrix detToTrackMatrix = CoordinateTransformations.getMatrix();
-
- return VecOp.mult(detToStripMatrix, VecOp.inverse(detToTrackMatrix));
- }
-
- private static String getName(IDetectorElement de) {
- // Find the first level down from the top of the de tree
- while (de.getParent().getParent() != null) {
- de = de.getParent();
- }
- // Find the name of this detector
- String detname = de.getName();
- return detname;
- }
-
- private static int getLayer(IDetectorElement de) {
- int layer = -1;
- IIdentifierHelper hlp = de.getIdentifierHelper();
- if (hlp instanceof DetectorIdentifierHelper) {
- DetectorIdentifierHelper dehlp = (DetectorIdentifierHelper) hlp;
- // Get the identifier
- IIdentifier id = de.getIdentifier();
- // Get the layer number
- layer = dehlp.getLayerValue(id);
- }
- return layer;
- }
-
- private static BarrelEndcapFlag getBarrelEndcapFlag(IDetectorElement de) {
- BarrelEndcapFlag beflag = BarrelEndcapFlag.UNKNOWN;
- // Find the second level down from the top of the de tree
- while (de.getParent().getParent().getParent() != null) {
- de = de.getParent();
- }
- // Get the DetectorIdentifierHelper
- IIdentifierHelper hlp = de.getIdentifierHelper();
- if (hlp instanceof DetectorIdentifierHelper) {
- DetectorIdentifierHelper dehlp = (DetectorIdentifierHelper) hlp;
- // Get the identifier
- IIdentifier id = de.getIdentifier();
- // Get the BarrelEndcapFlag
- if (dehlp.isBarrel(id)) {
- beflag = BarrelEndcapFlag.BARREL;
- } else if (dehlp.isEndcapPositive(id)) {
- beflag = BarrelEndcapFlag.ENDCAP_NORTH;
- } else if (dehlp.isEndcapNegative(id)) {
- beflag = BarrelEndcapFlag.ENDCAP_SOUTH;
- }
- }
- return beflag;
- }
-
- private static HelicalTrackStrip makeDigiStrip(SiTrackerHitStrip1D h) {
- SiTrackerHitStrip1D local = h.getTransformedHit(TrackerHitType.CoordinateSystem.SENSOR);
- SiTrackerHitStrip1D global = h.getTransformedHit(TrackerHitType.CoordinateSystem.GLOBAL);
-
- ITransform3D trans = local.getLocalToGlobal();
- Hep3Vector org = trans.transformed(new BasicHep3Vector(0., 0., 0.));
- Hep3Vector u = global.getMeasuredCoordinate();
- Hep3Vector v = global.getUnmeasuredCoordinate();
-
- double umeas = local.getPosition()[0];
- double vmin = VecOp.dot(local.getUnmeasuredCoordinate(), local.getHitSegment().getStartPoint());
- double vmax = VecOp.dot(local.getUnmeasuredCoordinate(), local.getHitSegment().getEndPoint());
- double du = Math.sqrt(local.getCovarianceAsMatrix().diagonal(0));
-
- IDetectorElement de = h.getSensor();
- String det = getName(de);
- int lyr = getLayer(de);
- BarrelEndcapFlag be = getBarrelEndcapFlag(de);
-
- double dEdx = h.getdEdx();
- double time = h.getTime();
- List<RawTrackerHit> rawhits = h.getRawHits();
- HelicalTrackStrip strip = new HelicalTrackStrip(org, u, v, umeas, du, vmin, vmax, dEdx, time, rawhits, det, lyr, be);
-
- try {
- if (h.getMCParticles() != null) {
- for (MCParticle p : h.getMCParticles()) {
- strip.addMCParticle(p);
- }
- }
- } catch (RuntimeException e) {
- // Okay when MC info not present.
- }
-
- return strip;
}
/**
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java Thu Oct 1 18:39:46 2015
@@ -65,7 +65,7 @@
@Override
public void detectorChanged(Detector detector) {
- Hep3Vector bfield = detector.getFieldMap().getField(new BasicHep3Vector(0., 0., 500.0));
+ Hep3Vector bfield = TrackUtils.getBField(detector);
// Create the class that handles all the GBL output
gbl = new GBLOutput(gblFileName, bfield); // if filename is empty no text file is written
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java Thu Oct 1 18:39:46 2015
@@ -1,12 +1,30 @@
package org.hps.recon.tracking.gbl;
import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Matrix;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+import org.hps.recon.tracking.CoordinateTransformations;
import org.hps.recon.tracking.MaterialSupervisor;
import org.hps.recon.tracking.MultipleScattering;
+import org.hps.recon.tracking.TrackUtils;
+import org.lcsim.constants.Constants;
import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.ITransform3D;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.HpsSiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.TrackerHit;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType;
import org.lcsim.recon.tracking.seedtracker.ScatterAngle;
/**
@@ -82,12 +100,11 @@
return mat;
}
- public static double estimateScatter(HelicalTrackStripGbl strip, HelicalTrackFit htf, MultipleScattering scattering, double _B) {
+ public static double estimateScatter(IDetectorElement hitElement, HelicalTrackFit htf, MultipleScattering scattering, double _B) {
//can be edge case where helix is outside, but close to sensor, so make a new scatter point assuming the helix does pass through the sensor
MaterialSupervisor.DetectorPlane hitPlane = null;
if (MaterialSupervisor.class.isInstance(scattering.getMaterialManager())) {
MaterialSupervisor matSup = (MaterialSupervisor) scattering.getMaterialManager();
- IDetectorElement hitElement = ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement();
for (MaterialSupervisor.ScatteringDetectorVolume vol : matSup.getMaterialVolumes()) {
if (vol.getDetectorElement() == hitElement) {
hitPlane = (MaterialSupervisor.DetectorPlane) vol;
@@ -107,4 +124,174 @@
throw new UnsupportedOperationException("Should not happen. This problem is only solved with the MaterialSupervisor.");
}
}
+
+ public static FittedGblTrajectory doGBLFit(HelicalTrackFit htf, List<TrackerHit> stripHits, MultipleScattering _scattering, double bfield, int debug) {
+ List<GBLStripClusterData> stripData = makeStripData(htf, stripHits, _scattering, bfield, debug);
+ double bfac = Constants.fieldConversion * bfield;
+
+ FittedGblTrajectory fit = HpsGblRefitter.fit(stripData, bfac, debug > 0);
+ return fit;
+ }
+
+ public static List<GBLStripClusterData> makeStripData(HelicalTrackFit htf, List<TrackerHit> stripHits, MultipleScattering _scattering, double _B, int _debug) {
+ List<GBLStripClusterData> stripClusterDataList = new ArrayList<GBLStripClusterData>();
+
+ // Find scatter points along the path
+ MultipleScattering.ScatterPoints scatters = _scattering.FindHPSScatterPoints(htf);
+
+ if (_debug > 0) {
+ System.out.printf("perPar covariance matrix\n%s\n", htf.covariance().toString());
+ }
+
+ for (TrackerHit stripHit : stripHits) {
+ HelicalTrackStripGbl strip;
+ if (stripHit instanceof SiTrackerHitStrip1D) {
+ strip = new HelicalTrackStripGbl(makeDigiStrip((SiTrackerHitStrip1D) stripHit), true);
+ } else {
+ SiTrackerHitStrip1D newHit = new SiTrackerHitStrip1D(stripHit);
+ strip = new HelicalTrackStripGbl(makeDigiStrip(newHit), true);
+ }
+
+ // find Millepede layer definition from DetectorElement
+ HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit) stripHit.getRawHits().get(0)).getDetectorElement();
+
+ int millepedeId = sensor.getMillepedeId();
+
+ if (_debug > 0) {
+ System.out.printf("layer %d millepede %d (DE=\"%s\", origin %s) \n", strip.layer(), millepedeId, sensor.getName(), strip.origin().toString());
+ }
+
+ //Center of the sensor
+ Hep3Vector origin = strip.origin();
+
+ //Find intercept point with sensor in tracking frame
+ Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B));
+ if (trkpos == null) {
+ if (_debug > 0) {
+ System.out.println("Can't find track intercept; use sensor origin");
+ }
+ trkpos = strip.origin();
+ }
+ if (_debug > 0) {
+ System.out.printf("trkpos at intercept [%.10f %.10f %.10f]\n", trkpos.x(), trkpos.y(), trkpos.z());
+ }
+
+ //GBLDATA
+ GBLStripClusterData stripData = new GBLStripClusterData(millepedeId);
+ //Add to output list
+ stripClusterDataList.add(stripData);
+
+ //path length to intercept
+ double s = HelixUtils.PathToXPlane(htf, trkpos.x(), 0, 0).get(0);
+ double s3D = s / Math.cos(Math.atan(htf.slope()));
+
+ //GBLDATA
+ stripData.setPath(s);
+ stripData.setPath3D(s3D);
+
+ //GBLDATA
+ stripData.setU(strip.u());
+ stripData.setV(strip.v());
+ stripData.setW(strip.w());
+
+ //Print track direction at intercept
+ Hep3Vector tDir = HelixUtils.Direction(htf, s);
+ double phi = htf.phi0() - s / htf.R();
+ double lambda = Math.atan(htf.slope());
+
+ //GBLDATA
+ stripData.setTrackDir(tDir);
+ stripData.setTrackPhi(phi);
+ stripData.setTrackLambda(lambda);
+
+ //Print residual in measurement system
+ // start by find the distance vector between the center and the track position
+ Hep3Vector vdiffTrk = VecOp.sub(trkpos, origin);
+
+ // then find the rotation from tracking to measurement frame
+ Hep3Matrix trkToStripRot = getTrackToStripRotation(sensor);
+
+ // then rotate that vector into the measurement frame to get the predicted measurement position
+ Hep3Vector trkpos_meas = VecOp.mult(trkToStripRot, vdiffTrk);
+
+ //GBLDATA
+ stripData.setMeas(strip.umeas());
+ stripData.setTrackPos(trkpos_meas);
+ stripData.setMeasErr(strip.du());
+
+ if (_debug > 1) {
+ System.out.printf("rotation matrix to meas frame\n%s\n", VecOp.toString(trkToStripRot));
+ System.out.printf("tPosGlobal %s origin %s\n", trkpos.toString(), origin.toString());
+ System.out.printf("tDiff %s\n", vdiffTrk.toString());
+ System.out.printf("tPosMeas %s\n", trkpos_meas.toString());
+ }
+
+ if (_debug > 0) {
+ System.out.printf("layer %d millePedeId %d uRes %.10f\n", strip.layer(), millepedeId, stripData.getMeas() - stripData.getTrackPos().x());
+ }
+
+ // find scattering angle
+ MultipleScattering.ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement());
+ double scatAngle;
+
+ if (scatter != null) {
+ scatAngle = scatter.getScatterAngle().Angle();
+ } else {
+ if (_debug > 0) {
+ System.out.printf("WARNING cannot find scatter for detector %s with strip cluster at %s\n", ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement().getName(), strip.origin().toString());
+ }
+ scatAngle = GblUtils.estimateScatter(sensor, htf, _scattering, _B);
+ }
+
+ //GBLDATA
+ stripData.setScatterAngle(scatAngle);
+ }
+ return stripClusterDataList;
+ }
+
+ private static Hep3Matrix getTrackToStripRotation(SiSensor sensor) {
+ // This function transforms the hit to the sensor coordinates
+
+ // Transform from JLab frame to sensor frame (done through the RawTrackerHit)
+ SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
+ ITransform3D detToStrip = electrodes.getGlobalToLocal();
+ // Get rotation matrix
+ Hep3Matrix detToStripMatrix = detToStrip.getRotation().getRotationMatrix();
+ // Transformation between the JLAB and tracking coordinate systems
+ Hep3Matrix detToTrackMatrix = CoordinateTransformations.getMatrix();
+
+ return VecOp.mult(detToStripMatrix, VecOp.inverse(detToTrackMatrix));
+ }
+
+ private static HelicalTrackStrip makeDigiStrip(SiTrackerHitStrip1D h) {
+ SiTrackerHitStrip1D local = h.getTransformedHit(TrackerHitType.CoordinateSystem.SENSOR);
+ SiTrackerHitStrip1D global = h.getTransformedHit(TrackerHitType.CoordinateSystem.GLOBAL);
+
+ ITransform3D trans = local.getLocalToGlobal();
+ Hep3Vector org = trans.transformed(new BasicHep3Vector(0., 0., 0.));
+ Hep3Vector u = global.getMeasuredCoordinate();
+ Hep3Vector v = global.getUnmeasuredCoordinate();
+
+ //rotate to tracking frame
+ Hep3Vector neworigin = CoordinateTransformations.transformVectorToTracking(org);
+ Hep3Vector newu = CoordinateTransformations.transformVectorToTracking(u);
+ Hep3Vector newv = CoordinateTransformations.transformVectorToTracking(v);
+
+ double umeas = local.getPosition()[0];
+ double vmin = VecOp.dot(local.getUnmeasuredCoordinate(), local.getHitSegment().getStartPoint());
+ double vmax = VecOp.dot(local.getUnmeasuredCoordinate(), local.getHitSegment().getEndPoint());
+ double du = Math.sqrt(local.getCovarianceAsMatrix().diagonal(0));
+
+ //don't fill fields we don't use
+// IDetectorElement de = h.getSensor();
+// String det = getName(de);
+// int lyr = getLayer(de);
+// BarrelEndcapFlag be = getBarrelEndcapFlag(de);
+ double dEdx = h.getdEdx();
+ double time = h.getTime();
+ List<RawTrackerHit> rawhits = h.getRawHits();
+ HelicalTrackStrip strip = new HelicalTrackStrip(neworigin, newu, newv, umeas, du, vmin, vmax, dEdx, time, rawhits, null, -1, null);
+
+ return strip;
+ }
}
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 Oct 1 18:39:46 2015
@@ -303,7 +303,7 @@
if (_debug) {
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());
}
- scatAngle = GblUtils.estimateScatter(strip, htf, _scattering, _B);
+ scatAngle = GblUtils.estimateScatter(((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement(), htf, _scattering, _B);
}
// Scattering angle in the curvilinear frame
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 Oct 1 18:39:46 2015
@@ -13,6 +13,7 @@
import java.util.logging.Formatter;
import java.util.logging.Level;
import java.util.logging.Logger;
+import org.hps.recon.tracking.TrackUtils;
import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.SymMatrix;
import org.hps.recon.tracking.gbl.matrix.Vector;
@@ -94,9 +95,7 @@
@Override
protected void process(EventHeader event) {
-// Hep3Vector bfieldvec = event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0., 1.));
- Hep3Vector bfieldvec = event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0., 500.));//mg...get the b-field in the middle of the magnet
- double bfield = bfieldvec.y();
+ double bfield = TrackUtils.getBField(event.getDetector()).y();
//double bfac = 0.0002998 * bfield;
double bfac = Constants.fieldConversion * bfield;
@@ -155,11 +154,9 @@
// loop over the tracks and do the GBL fit
List<FittedGblTrajectory> trackFits = new ArrayList<FittedGblTrajectory>();
- int trackNum = 0;
logger.info("Trying to fit " + stripsGblMap.size() + " tracks");
for (GBLTrackData t : stripsGblMap.keySet()) {
FittedGblTrajectory traj = fit(stripsGblMap.get(t), bfac, _debug);
- ++trackNum;
if (traj != null) {
logger.info("GBL fit successful");
if (_debug) {
@@ -170,7 +167,6 @@
traj.get_traj().milleOut(mille);
}
traj.set_seed(gblToSeedMap.get(t));
- traj.set_track_data(t);
trackFits.add(traj);
} else {
logger.info("GBL fit failed");
@@ -194,7 +190,7 @@
}
- private static FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac, boolean debug) {
+ public static FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac, boolean debug) {
// path length along trajectory
double s = 0.;
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java Thu Oct 1 18:39:46 2015
@@ -13,6 +13,7 @@
import org.apache.commons.math3.util.Pair;
import org.hps.recon.tracking.TrackType;
+import static org.hps.recon.tracking.gbl.GBLOutput.getPerToClPrj;
import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.SymMatrix;
import org.hps.recon.tracking.gbl.matrix.Vector;
@@ -50,10 +51,11 @@
}
public void setDebug(boolean debug) {
- if (debug)
+ if (debug) {
logger.setLevel(Level.INFO);
- else
+ } else {
logger.setLevel(Level.OFF);
+ }
}
/**
@@ -72,50 +74,14 @@
for (FittedGblTrajectory fittedTraj : gblTrajectories) {
- // Initialize the reference point to the origin
- double[] ref = new double[]{0., 0., 0.};
SeedTrack seedTrack = (SeedTrack) fittedTraj.get_seed();
SeedCandidate trackseed = seedTrack.getSeedCandidate();
- // Create a new SeedTrack
- BaseTrack trk = new BaseTrack();
-
- // Add the hits to the track
- for (HelicalTrackHit hit : trackseed.getHits()) {
- trk.addHit((TrackerHit) hit);
- }
-
- // Retrieve the helix
- HelicalTrackFit helix = trackseed.getHelix();
- // Set state at vertex
- Pair<double[], SymmetricMatrix> correctedHelixParams = getGblCorrectedHelixParameters(helix, fittedTraj, bfield, FittedGblTrajectory.GBLPOINT.IP);
- trk.setTrackParameters(correctedHelixParams.getFirst(), bfield);// hack to set the track charge
- trk.getTrackStates().clear();
-
- TrackState stateVertex = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtIP, bfield);
- trk.getTrackStates().add(stateVertex);
-
- // Set state at last point
- Pair<double[], SymmetricMatrix> correctedHelixParamsLast = getGblCorrectedHelixParameters(helix, fittedTraj, bfield, FittedGblTrajectory.GBLPOINT.LAST);
- TrackState stateLast = new BaseTrackState(correctedHelixParamsLast.getFirst(), ref, correctedHelixParamsLast.getSecond().asPackedArray(true), TrackState.AtLastHit, bfield);
- trk.getTrackStates().add(stateLast);
-
- // Set other info needed
- trk.setChisq(fittedTraj.get_chi2());
- trk.setNDF(fittedTraj.get_ndf());
- trk.setFitSuccess(true);
- trk.setRefPointIsDCA(true);
- trk.setTrackType(TrackType.setGBL(seedTrack.getType(), true));
-
+ // Create a new Track
+ Track trk = makeCorrectedTrack(fittedTraj, trackseed.getHelix(), seedTrack.getTrackerHits(), seedTrack.getType(), bfield);
+
// Add the track to the list of tracks
tracks.add(trk);
- logger.info(String.format("helix chi2 %f ndf %d gbl chi2 %f ndf %d\n", helix.chisqtot(), helix.ndf()[0] + helix.ndf()[1], trk.getChi2(), trk.getNDF()));
- if (logger.getLevel().intValue() <= Level.INFO.intValue()) {
- for (int i = 0; i < 5; ++i) {
- logger.info(String.format("param %d: %.10f -> %.10f helix-gbl= %f", i, helix.parameters()[i], trk.getTrackParameter(i), helix.parameters()[i] - trk.getTrackParameter(i)));
- }
- }
-
}
logger.info("adding " + Integer.toString(tracks.size()) + " Gbl tracks to event with " + event.get(Track.class, "MatchedTracks").size() + " matched tracks");
@@ -125,16 +91,59 @@
event.put(_TrkCollectionName, tracks, Track.class, flag);
}
-
+ public static Track makeCorrectedTrack(FittedGblTrajectory fittedTraj, HelicalTrackFit helix, List<TrackerHit> trackHits, int trackType, double bfield) {
+ // Initialize the reference point to the origin
+ double[] ref = new double[]{0., 0., 0.};
+
+ // Create a new SeedTrack
+ BaseTrack trk = new BaseTrack();
+
+ // Add the hits to the track
+ for (TrackerHit hit : trackHits) {
+ trk.addHit(hit);
+ }
+
+ // Set state at vertex
+ Pair<double[], SymmetricMatrix> correctedHelixParams = getGblCorrectedHelixParameters(helix, fittedTraj.get_traj(), bfield, FittedGblTrajectory.GBLPOINT.IP);
+ trk.setTrackParameters(correctedHelixParams.getFirst(), bfield);// hack to set the track charge
+ trk.getTrackStates().clear();
+
+ TrackState stateVertex = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtIP, bfield);
+ trk.getTrackStates().add(stateVertex);
+
+ // Set state at last point
+ Pair<double[], SymmetricMatrix> correctedHelixParamsLast = getGblCorrectedHelixParameters(helix, fittedTraj.get_traj(), bfield, FittedGblTrajectory.GBLPOINT.LAST);
+ TrackState stateLast = new BaseTrackState(correctedHelixParamsLast.getFirst(), ref, correctedHelixParamsLast.getSecond().asPackedArray(true), TrackState.AtLastHit, bfield);
+ trk.getTrackStates().add(stateLast);
+
+ // Set other info needed
+ trk.setChisq(fittedTraj.get_chi2());
+ trk.setNDF(fittedTraj.get_ndf());
+ trk.setFitSuccess(true);
+ trk.setRefPointIsDCA(true);
+ trk.setTrackType(TrackType.setGBL(trackType, true));
+
+ // Add the track to the list of tracks
+// tracks.add(trk);
+ logger.info(String.format("helix chi2 %f ndf %d gbl chi2 %f ndf %d\n", helix.chisqtot(), helix.ndf()[0] + helix.ndf()[1], trk.getChi2(), trk.getNDF()));
+ if (logger.getLevel().intValue() <= Level.INFO.intValue()) {
+ for (int i = 0; i < 5; ++i) {
+ logger.info(String.format("param %d: %.10f -> %.10f helix-gbl= %f", i, helix.parameters()[i], trk.getTrackParameter(i), helix.parameters()[i] - trk.getTrackParameter(i)));
+ }
+ }
+ return trk;
+ }
+
/**
- * Compute the updated helix parameters and covariance matrix at a given point along the trajectory.
+ * Compute the updated helix parameters and covariance matrix at a given
+ * point along the trajectory.
*
* @param helix - original seed track
* @param traj - fitted GBL trajectory
* @param point - the point along the track where the result is computed.
* @return corrected parameters.
*/
- private Pair<double[], SymmetricMatrix> getGblCorrectedHelixParameters(HelicalTrackFit helix, FittedGblTrajectory traj, double bfield, FittedGblTrajectory.GBLPOINT point) {
+ public static Pair<double[], SymmetricMatrix> getGblCorrectedHelixParameters(HelicalTrackFit helix, GblTrajectory traj, double bfield, FittedGblTrajectory.GBLPOINT point) {
// get seed helix parameters
double qOverP = helix.curvature() / (Constants.fieldConversion * Math.abs(bfield) * Math.sqrt(1 + Math.pow(helix.slope(), 2)));
@@ -151,14 +160,15 @@
Vector locPar = new Vector(5);
SymMatrix locCov = new SymMatrix(5);
int pointIndex;
- if( point.compareTo( FittedGblTrajectory.GBLPOINT.IP) == 0 )
+ if (point.compareTo(FittedGblTrajectory.GBLPOINT.IP) == 0) {
pointIndex = 1;
- else if (point.compareTo(FittedGblTrajectory.GBLPOINT.LAST) == 0)
- pointIndex = traj.get_traj().getNumPoints();
- else
+ } else if (point.compareTo(FittedGblTrajectory.GBLPOINT.LAST) == 0) {
+ pointIndex = traj.getNumPoints();
+ } else {
throw new RuntimeException("This GBLPOINT " + point.toString() + "( " + point.name() + ") is not valid");
-
- traj.get_traj().getResults(pointIndex, locPar, locCov); // vertex point
+ }
+
+ traj.getResults(pointIndex, locPar, locCov); // vertex point
// locCov.print(10, 8);
double qOverPCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.QOVERP.getValue());
double xTPrimeCorr = locPar.get(FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue());
@@ -169,7 +179,9 @@
logger.info((helix.slope() > 0 ? "top: " : "bot ") + "qOverPCorr " + qOverPCorr + " xtPrimeCorr " + xTPrimeCorr + " yTPrimeCorr " + yTPrimeCorr + " xTCorr " + xTCorr + " yTCorr " + yTCorr);
// calculate new d0 and z0
- Hep3Matrix perToClPrj = traj.get_track_data().getPrjPerToCl();
+// Hep3Matrix perToClPrj = traj.get_track_data().getPrjPerToCl();
+ Hep3Matrix perToClPrj = getPerToClPrj(helix);
+
Hep3Matrix clToPerPrj = VecOp.inverse(perToClPrj);
Hep3Vector corrPer = VecOp.mult(clToPerPrj, new BasicHep3Vector(xTCorr, yTCorr, 0.0));
