Author: [log in to unmask] Date: Thu Oct 1 19:49:39 2015 New Revision: 3750 Log: start work on simple GBL refitter Added: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLRefitterDriver.java - copied, changed from r3749, java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/MergeTrackCollections.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/MergeTrackCollections.java ============================================================================= --- java/trunk/tracking/src/main/java/org/hps/recon/tracking/MergeTrackCollections.java (original) +++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/MergeTrackCollections.java Thu Oct 1 19:49:39 2015 @@ -7,7 +7,6 @@ import java.util.Map; import java.util.Set; import org.lcsim.event.EventHeader; -import org.lcsim.event.RelationalTable; import org.lcsim.event.Track; import org.lcsim.event.TrackerHit; import org.lcsim.lcio.LCIOConstants; @@ -68,15 +67,13 @@ @Override public void process(EventHeader event) { - RelationalTable hitToStrips = TrackUtils.getHitToStripsTable(event); - RelationalTable hitToRotated = TrackUtils.getHitToRotatedTable(event); List<List<Track>> trackCollections = event.get(Track.class); // Loop over each of the track collections retrieved from the event Map<Set<TrackerHit>, List<Track>> hitsToTracksMap = new HashMap<Set<TrackerHit>, List<Track>>(); for (List<Track> tracklist : trackCollections) { for (Track trk : tracklist) { - Set<TrackerHit> trackHits = new HashSet<TrackerHit>(TrackUtils.getStripHits(trk, hitToStrips, hitToRotated)); + Set<TrackerHit> trackHits = new HashSet<TrackerHit>(trk.getTrackerHits()); List<Track> matchingTracks = hitsToTracksMap.get(trackHits); if (matchingTracks == null) { matchingTracks = new ArrayList<Track>(); Copied: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLRefitterDriver.java (from r3749, 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/GBLRefitterDriver.java Thu Oct 1 19:49:39 2015 @@ -1,754 +1,108 @@ package org.hps.recon.tracking.gbl; -import hep.physics.vec.BasicHep3Vector; -import hep.physics.vec.Hep3Vector; -import hep.physics.vec.VecOp; -import static java.lang.Math.abs; -import static java.lang.Math.sin; -import static java.lang.Math.sqrt; import java.util.ArrayList; import java.util.HashMap; +import java.util.HashSet; import java.util.List; import java.util.Map; -import java.util.logging.Formatter; -import java.util.logging.Level; -import java.util.logging.Logger; +import java.util.Set; +import org.hps.recon.tracking.MaterialSupervisor; +import org.hps.recon.tracking.MultipleScattering; 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; -import org.hps.util.BasicLogFormatter; -import org.lcsim.constants.Constants; import org.lcsim.event.EventHeader; -import org.lcsim.event.GenericObject; -import org.lcsim.event.LCRelation; +import org.lcsim.event.RawTrackerHit; +import org.lcsim.event.RelationalTable; import org.lcsim.event.Track; +import org.lcsim.event.TrackerHit; +import org.lcsim.fit.helicaltrack.HelicalTrackFit; import org.lcsim.geometry.Detector; -import org.lcsim.geometry.compact.converter.MilleParameter; +import org.lcsim.lcio.LCIOConstants; import org.lcsim.util.Driver; -import org.lcsim.util.log.LogUtil; /** - * A Driver which refits tracks using GBL. Modeled on the hps-dst code written - * by Per Hansson and Omar Moreno. Requires the GBL Collections and Relations to + * A Driver which refits tracks using GBL. Does not require GBL collections to * be present in the event. - * - * @author Norman A Graf - * @author Per Hansson Adrian <[log in to unmask]> - * - * @version $Id: HpsGblRefitter.java 3460 2015-08-29 01:45:39Z - * [log in to unmask] $ */ -public class HpsGblRefitter extends Driver { +public class GBLRefitterDriver extends Driver { - static Formatter f = new BasicLogFormatter(); - private static Logger logger = LogUtil.create(HpsGblRefitter.class.getSimpleName(), f, Level.WARNING); - //private static final Logger logger = Logger.getLogger(HpsGblRefitter.class.getName()); - private boolean _debug = false; - private final String trackCollectionName = "MatchedTracks"; - private final String track2GblTrackRelationName = "TrackToGBLTrack"; - private final String gblTrack2StripRelationName = "GBLTrackToStripData"; + private String inputCollectionName = "MatchedTracks"; + private String outputCollectionName = "GBLTracks"; - private MilleBinary mille; - private String milleBinaryFileName = MilleBinary.DEFAULT_OUTPUT_FILE_NAME; - private boolean writeMilleBinary = false; + private double bfield; + private final MultipleScattering _scattering = new MultipleScattering(new MaterialSupervisor()); - private final MakeGblTracks _makeTracks; - - public void setDebug(boolean debug) { - _debug = debug; - _makeTracks.setDebug(debug); + public void setInputCollectionName(String inputCollectionName) { + this.inputCollectionName = inputCollectionName; } - public void setMilleBinaryFileName(String filename) { - milleBinaryFileName = filename; - } - - public void setWriteMilleBinary(boolean writeMillepedeFile) { - writeMilleBinary = writeMillepedeFile; - } - - public HpsGblRefitter() { - _makeTracks = new MakeGblTracks(); - _makeTracks.setDebug(_debug); - logger.setLevel(Level.WARNING); - System.out.println("level " + logger.getLevel().toString()); - } - - //@Override - //public void setLogLevel(String logLevel) { - // logger.setLevel(Level.parse(logLevel)); - //} - @Override - protected void startOfData() { - if (writeMilleBinary) { - mille = new MilleBinary(milleBinaryFileName); - } + public void setOutputCollectionName(String outputCollectionName) { + this.outputCollectionName = outputCollectionName; } @Override - protected void endOfData() { - if (writeMilleBinary) { - mille.close(); - } + protected void detectorChanged(Detector detector) { + bfield = Math.abs(TrackUtils.getBField(detector).magnitude()); + _scattering.getMaterialManager().buildModel(detector); + _scattering.setBField(bfield); // only absolute of B is needed as it's used for momentum calculation only } @Override protected void process(EventHeader event) { - double bfield = TrackUtils.getBField(event.getDetector()).y(); - //double bfac = 0.0002998 * bfield; - double bfac = Constants.fieldConversion * bfield; - - // get the tracks - if (!event.hasCollection(Track.class, trackCollectionName)) { - if (_debug) { - System.out.printf("%s: No tracks in Event %d \n", this.getClass().getSimpleName(), event.getEventNumber()); - } + if (!event.hasCollection(Track.class, inputCollectionName)) { return; } + List<Track> tracks = event.get(Track.class, inputCollectionName); + RelationalTable hitToStrips = TrackUtils.getHitToStripsTable(event); + RelationalTable hitToRotated = TrackUtils.getHitToRotatedTable(event); - //get the relations to the GBLtracks - if (!event.hasItem(track2GblTrackRelationName)) { - System.out.println("Need Relations " + track2GblTrackRelationName); - return; - } - // and strips - if (!event.hasItem(gblTrack2StripRelationName)) { - System.out.println("Need Relations " + gblTrack2StripRelationName); - return; + List<Track> refittedTracks = new ArrayList<Track>(); + + Map<Track, Track> inputToRefitted = new HashMap<Track, Track>(); + for (Track track : tracks) { + HelicalTrackFit helix = TrackUtils.getHTF(track); + FittedGblTrajectory fit = GblUtils.doGBLFit(helix, TrackUtils.getStripHits(track, hitToStrips, hitToRotated), _scattering, bfield, 0); + + Track newTrack = MakeGblTracks.makeCorrectedTrack(fit, helix, track.getTrackerHits(), track.getType(), bfield); + refittedTracks.add(newTrack); + inputToRefitted.put(track, newTrack); } - List<LCRelation> track2GblTrackRelations = event.get(LCRelation.class, track2GblTrackRelationName); - //need a map of GBLTrackData keyed on the Generic object from which it created - Map<GenericObject, GBLTrackData> gblObjMap = new HashMap<GenericObject, GBLTrackData>(); - //need a map of SeedTrack to GBLTrackData keyed on the track object from which it created - Map<GBLTrackData, Track> gblToSeedMap = new HashMap<GBLTrackData, Track>(); + Map<Set<TrackerHit>, List<Track>> hitSetToTrackList = new HashMap<Set<TrackerHit>, List<Track>>(); - // loop over the relations - for (LCRelation relation : track2GblTrackRelations) { - Track t = (Track) relation.getFrom(); - GenericObject gblTrackObject = (GenericObject) relation.getTo(); - GBLTrackData gblT = new GBLTrackData(gblTrackObject); - gblObjMap.put(gblTrackObject, gblT); - gblToSeedMap.put(gblT, t); - } //end of loop over tracks + for (Track track : tracks) { + Set<TrackerHit> trackHth = new HashSet<TrackerHit>(track.getTrackerHits()); + for (Track otherTrack : tracks) { + Set<TrackerHit> allHth = new HashSet<TrackerHit>(otherTrack.getTrackerHits()); + allHth.addAll(trackHth); + List<TrackerHit> hthList = new ArrayList<TrackerHit>(allHth); + if (hthList.size() == trackHth.size()) { + continue; + } - //get the strip hit relations - List<LCRelation> gblTrack2StripRelations = event.get(LCRelation.class, gblTrack2StripRelationName); + boolean[] hasHit = new boolean[6]; + boolean isGood = true; - // need a map of lists of strip data keyed by the gblTrack to which they correspond - Map<GBLTrackData, List<GBLStripClusterData>> stripsGblMap = new HashMap<GBLTrackData, List<GBLStripClusterData>>(); - for (LCRelation relation : gblTrack2StripRelations) { - //from GBLTrackData to GBLStripClusterData - GenericObject gblTrackObject = (GenericObject) relation.getFrom(); - //Let's get the GBLTrackData that corresponds to this object... - GBLTrackData gblT = gblObjMap.get(gblTrackObject); - GBLStripClusterData sd = new GBLStripClusterData((GenericObject) relation.getTo()); - if (stripsGblMap.containsKey(gblT)) { - stripsGblMap.get(gblT).add(sd); - } else { - stripsGblMap.put(gblT, new ArrayList<GBLStripClusterData>()); - stripsGblMap.get(gblT).add(sd); - } - } + for (TrackerHit hit : hthList) { + int layer = (TrackUtils.getLayer(hit) - 1) / 2; + if (hasHit[layer]) { + isGood = false; + break; + } + hasHit[layer] = true; + } + if (isGood) { + HelicalTrackFit helix = TrackUtils.getHTF(track); + Set<TrackerHit> allStripHits = new HashSet<TrackerHit>(TrackUtils.getStripHits(track, hitToStrips, hitToRotated)); + allStripHits.addAll(TrackUtils.getStripHits(otherTrack, hitToStrips, hitToRotated)); - // loop over the tracks and do the GBL fit - List<FittedGblTrajectory> trackFits = new ArrayList<FittedGblTrajectory>(); - logger.info("Trying to fit " + stripsGblMap.size() + " tracks"); - for (GBLTrackData t : stripsGblMap.keySet()) { - FittedGblTrajectory traj = fit(stripsGblMap.get(t), bfac, _debug); - if (traj != null) { - logger.info("GBL fit successful"); - if (_debug) { - System.out.printf("%s: GBL fit successful.\n", getClass().getSimpleName()); - } - // write to MP binary file - if (writeMilleBinary) { - traj.get_traj().milleOut(mille); - } - traj.set_seed(gblToSeedMap.get(t)); - trackFits.add(traj); - } else { - logger.info("GBL fit failed"); - if (_debug) { - System.out.printf("%s: GBL fit failed.\n", getClass().getSimpleName()); + FittedGblTrajectory fit = GblUtils.doGBLFit(helix, new ArrayList<TrackerHit>(allStripHits), _scattering, bfield, 0); + Track newTrack = MakeGblTracks.makeCorrectedTrack(fit, helix, new ArrayList<TrackerHit>(allHth), 0, bfield); + System.out.format("%f %f %f\n", fit.get_chi2(), inputToRefitted.get(track).getChi2(), inputToRefitted.get(otherTrack).getChi2()); } } } - - logger.info(event.get(Track.class, trackCollectionName).size() + " tracks in collection \"" + trackCollectionName + "\""); - logger.info(gblObjMap.size() + " tracks in gblObjMap"); - logger.info(gblToSeedMap.size() + " tracks in gblToSeedMap"); - logger.info(stripsGblMap.size() + " tracks in stripsGblMap"); - logger.info(trackFits.size() + " fitted GBL tracks before adding to event"); - - _makeTracks.Process(event, trackFits, bfield); - - if (_debug) { - System.out.printf("%s: Done.\n", getClass().getSimpleName()); - } - - } - - public static FittedGblTrajectory fit(List<GBLStripClusterData> hits, double bfac, boolean debug) { - // path length along trajectory - double s = 0.; - - // jacobian to transport errors between points along the path - Matrix jacPointToPoint = new Matrix(5, 5); - jacPointToPoint.UnitMatrix(); - // Vector of the strip clusters used for the GBL fit - List<GblPoint> listOfPoints = new ArrayList<GblPoint>(); - - // Store the projection from local to measurement frame for each strip cluster - Map< Integer, Matrix> proL2m_list = new HashMap<Integer, Matrix>(); - // Save the association between strip cluster and label - - //start trajectory at refence point (s=0) - this point has no measurement - GblPoint ref_point = new GblPoint(jacPointToPoint); - listOfPoints.add(ref_point); - - // Loop over strips - int n_strips = hits.size(); - for (int istrip = 0; istrip != n_strips; ++istrip) { - GBLStripClusterData strip = hits.get(istrip); - //MG--9/18/2015--beamspot has Id=666/667...don't include it in the GBL fit - if (strip.getId() > 99) { - continue; - } - if (debug) { - System.out.println("HpsGblFitter: Processing strip " + istrip + " with id/layer " + strip.getId()); - } - // Path length step for this cluster - double step = strip.getPath3D() - s; - if (debug) { - System.out.println("HpsGblFitter: " + "Path length step " + step + " from " + s + " to " + strip.getPath3D()); - } - - // get measurement frame unit vectors - Hep3Vector u = strip.getU(); - Hep3Vector v = strip.getV(); - Hep3Vector w = strip.getW(); - - // Measurement direction (perpendicular and parallel to strip direction) - Matrix mDir = new Matrix(2, 3); - mDir.set(0, 0, u.x()); - mDir.set(0, 1, u.y()); - mDir.set(0, 2, u.z()); - mDir.set(1, 0, v.x()); - mDir.set(1, 1, v.y()); - mDir.set(1, 2, v.z()); - if (debug) { - System.out.println("HpsGblFitter: " + "mDir"); - mDir.print(4, 6); - } - Matrix mDirT = mDir.copy().transpose(); - - if (debug) { - System.out.println("HpsGblFitter: " + "mDirT"); - mDirT.print(4, 6); - } - - // Track direction - double sinLambda = sin(strip.getTrackLambda());//->GetLambda()); - double cosLambda = sqrt(1.0 - sinLambda * sinLambda); - double sinPhi = sin(strip.getTrackPhi());//->GetPhi()); - double cosPhi = sqrt(1.0 - sinPhi * sinPhi); - - if (debug) { - System.out.println("HpsGblFitter: " + "Track direction sinLambda=" + sinLambda + " sinPhi=" + sinPhi); - } - - // Track direction in curvilinear frame (U,V,T) - // U = Z x T / |Z x T|, V = T x U - Matrix uvDir = new Matrix(2, 3); - uvDir.set(0, 0, -sinPhi); - uvDir.set(0, 1, cosPhi); - uvDir.set(0, 2, 0.); - uvDir.set(1, 0, -sinLambda * cosPhi); - uvDir.set(1, 1, -sinLambda * sinPhi); - uvDir.set(1, 2, cosLambda); - - if (debug) { - System.out.println("HpsGblFitter: " + "uvDir"); - uvDir.print(6, 4); - } - - // projection from measurement to local (curvilinear uv) directions (duv/dm) - Matrix proM2l = uvDir.times(mDirT); - - //projection from local (uv) to measurement directions (dm/duv) - Matrix proL2m = proM2l.copy(); - proL2m = proL2m.inverse(); - proL2m_list.put(strip.getId(), proL2m.copy()); // is a copy needed or is that just a C++/root thing? - - if (debug) { - System.out.println("HpsGblFitter: " + "proM2l:"); - proM2l.print(4, 6); - System.out.println("HpsGblFitter: " + "proL2m:"); - proL2m.print(4, 6); - System.out.println("HpsGblFitter: " + "proM2l*proL2m (should be unit matrix):"); - (proM2l.times(proL2m)).print(4, 6); - } - - // measurement/residual in the measurement system - // only 1D measurement in u-direction, set strip measurement direction to zero - Vector meas = new Vector(2); -// double uRes = strip->GetUmeas() - strip->GetTrackPos().x(); // how can this be correct? - double uRes = strip.getMeas() - strip.getTrackPos().x(); - meas.set(0, uRes); - meas.set(1, 0.); -// //meas[0][0] += deltaU[iLayer] # misalignment - Vector measErr = new Vector(2); - measErr.set(0, strip.getMeasErr()); - measErr.set(1, 0.); - Vector measPrec = new Vector(2); - measPrec.set(0, 1.0 / (measErr.get(0) * measErr.get(0))); - measPrec.set(1, 0.); - - if (debug) { - System.out.println("HpsGblFitter: " + "meas: "); - meas.print(4, 6); - System.out.println("HpsGblFitter: " + "measErr:"); - measErr.print(4, 6); - System.out.println("HpsGblFitter: " + "measPrec:"); - measPrec.print(4, 6); - } - - //Find the Jacobian to be able to propagate the covariance matrix to this strip position - jacPointToPoint = gblSimpleJacobianLambdaPhi(step, cosLambda, abs(bfac)); - - if (debug) { - System.out.println("HpsGblFitter: " + "jacPointToPoint to extrapolate to this point:"); - jacPointToPoint.print(4, 6); - } - - // Get the transpose of the Jacobian - Matrix jacPointToPointTransposed = jacPointToPoint.copy().transpose(); - - // Option to use uncorrelated MS errors - // This is similar to what is done in lcsim seedtracker - // The msCov below holds the MS errors - // This is for testing purposes only. - boolean useUncorrMS = false; - Matrix msCov = new Matrix(5, 5); - - // Propagate the MS covariance matrix (in the curvilinear frame) to this strip position - msCov = msCov.times(jacPointToPointTransposed); - msCov = jacPointToPoint.times(msCov); - - // Get the MS covariance for the measurements in the measurement frame - Matrix proL2mTransposed = proL2m.copy().transpose(); - - Matrix measMsCov = proL2m.times((msCov.getMatrix(3, 4, 3, 4)).times(proL2mTransposed)); - - if (debug) { - System.out.println("HpsGblFitter: " + " msCov at this point:"); - msCov.print(4, 6); - System.out.println("HpsGblFitter: " + "measMsCov at this point:"); - measMsCov.print(4, 6); - } - - GblPoint point = new GblPoint(jacPointToPoint); - //Add measurement to the point - point.addMeasurement(proL2m, meas, measPrec, 0.); - //Add scatterer in curvilinear frame to the point - // no direction in this frame - Vector scat = new Vector(2); - - // Scattering angle in the curvilinear frame - //Note the cosLambda to correct for the projection in the phi direction - Vector scatErr = new Vector(2); - scatErr.set(0, strip.getScatterAngle()); - scatErr.set(1, strip.getScatterAngle() / cosLambda); - Vector scatPrec = new Vector(2); - scatPrec.set(0, 1.0 / (scatErr.get(0) * scatErr.get(0))); - scatPrec.set(1, 1.0 / (scatErr.get(1) * scatErr.get(1))); - - // add scatterer if not using the uncorrelated MS covariances for testing - if (!useUncorrMS) { - point.addScatterer(scat, scatPrec); - if (debug) { - System.out.println("HpsGblFitter: " + "adding scatError to this point:"); - scatErr.print(4, 6); - } - } - - // Add this GBL point to list that will be used in fit - listOfPoints.add(point); - int iLabel = listOfPoints.size(); - - // Update MS covariance matrix - msCov.set(1, 1, msCov.get(1, 1) + scatErr.get(0) * scatErr.get(0)); - msCov.set(2, 2, msCov.get(2, 2) + scatErr.get(1) * scatErr.get(1)); - - //// Calculate global derivatives for this point - // track direction in tracking/global frame - Hep3Vector tDirGlobal = new BasicHep3Vector(cosPhi * cosLambda, sinPhi * cosLambda, sinLambda); - - // Cross-check that the input is consistent - if (VecOp.sub(tDirGlobal, strip.getTrackDirection()).magnitude() > 0.00001) { - throw new RuntimeException("track directions are inconsistent: " + tDirGlobal.toString() + " and " + strip.getTrackDirection().toString()); - } - // rotate track direction to measurement frame - Hep3Vector tDirMeas = new BasicHep3Vector(VecOp.dot(tDirGlobal, u), VecOp.dot(tDirGlobal, v), VecOp.dot(tDirGlobal, w)); - // TODO this is a trivial one. Fix it. - Hep3Vector normalMeas = new BasicHep3Vector(VecOp.dot(w, u), VecOp.dot(w, v), VecOp.dot(w, w)); - - // vector coplanar with measurement plane from origin to prediction - Hep3Vector tPosMeas = strip.getTrackPos(); - - // measurements: non-measured directions - double vmeas = 0.; - double wmeas = 0.; - - // calculate and add derivatives to point - GlobalDers glDers = new GlobalDers(strip.getId(), meas.get(0), vmeas, wmeas, tDirMeas, tPosMeas, normalMeas); - - //TODO find a more robust way to get half. - boolean isTop = Math.sin(strip.getTrackLambda()) > 0; - - // Get the list of millepede parameters - List<MilleParameter> milleParameters = glDers.getDers(isTop); - - // need to make vector and matrices for interface - List<Integer> labGlobal = new ArrayList<Integer>(); - Matrix addDer = new Matrix(1, milleParameters.size()); - for (int i = 0; i < milleParameters.size(); ++i) { - labGlobal.add(milleParameters.get(i).getId()); - addDer.set(0, i, milleParameters.get(i).getValue()); - } - point.addGlobals(labGlobal, addDer); - String logders = ""; - for (int i = 0; i < milleParameters.size(); ++i) { - logders += labGlobal.get(i) + "\t" + addDer.get(0, i) + "\n"; - } - logger.info("\n" + logders); - - logger.info("uRes " + strip.getId() + " uRes " + uRes + " pred (" + strip.getTrackPos().x() + "," + strip.getTrackPos().y() + "," + strip.getTrackPos().z() + ") s(3D) " + strip.getPath3D()); - - //go to next point - s += step; - - } //strips - - //create the trajectory - GblTrajectory traj = new GblTrajectory(listOfPoints); //,seedLabel, clSeed); - - if (!traj.isValid()) { - System.out.println("HpsGblFitter: " + " Invalid GblTrajectory -> skip"); - return null; //1;//INVALIDTRAJ; - } - - // print the trajectory - if (debug) { - System.out.println("%%%% Gbl Trajectory "); - traj.printTrajectory(1); - traj.printData(); - traj.printPoints(4); - } - // fit trajectory - double[] dVals = new double[2]; - int[] iVals = new int[1]; - traj.fit(dVals, iVals, ""); - logger.info("fit result: Chi2=" + dVals[0] + " Ndf=" + iVals[0] + " Lost=" + dVals[1]); - Vector aCorrection = new Vector(5); - SymMatrix aCovariance = new SymMatrix(5); - traj.getResults(1, aCorrection, aCovariance); - if (debug) { - System.out.println(" cor "); - aCorrection.print(6, 4); - System.out.println(" cov "); - aCovariance.print(6, 4); - } - - logger.fine("locPar " + aCorrection.toString()); - -// - return new FittedGblTrajectory(traj, dVals[0], iVals[0], dVals[1]); - } - - @Override - protected void detectorChanged(Detector detector) { - } - - private static Matrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac) { - /** - * Simple jacobian: quadratic in arc length difference. using lambda phi - * as directions - * - * @param ds: arc length difference - * @type ds: float - * @param cosl: cos(lambda) - * @type cosl: float - * @param bfac: Bz*c - * @type bfac: float - * @return: jacobian to move by 'ds' on trajectory - * @rtype: matrix(float) ajac(1,1)= 1.0D0 ajac(2,2)= 1.0D0 - * ajac(3,1)=-DBLE(bfac*ds) ajac(3,3)= 1.0D0 - * ajac(4,1)=-DBLE(0.5*bfac*ds*ds*cosl) ajac(4,3)= DBLE(ds*cosl) - * ajac(4,4)= 1.0D0 ajac(5,2)= DBLE(ds) ajac(5,5)= 1.0D0 ''' jac = - * np.eye(5) jac[2, 0] = -bfac * ds jac[3, 0] = -0.5 * bfac * ds * ds * - * cosl jac[3, 2] = ds * cosl jac[4, 1] = ds return jac - */ - Matrix jac = new Matrix(5, 5); - jac.UnitMatrix(); - jac.set(2, 0, -bfac * ds); - jac.set(3, 0, -0.5 * bfac * ds * ds * cosl); - jac.set(3, 2, ds * cosl); - jac.set(4, 1, ds); - return jac; - } - - private static class GlobalDers { - - private final int _layer; - private final double _umeas; // measurement direction - private final double _vmeas; // unmeasured direction - private final double _wmeas; // normal to plane - private final Hep3Vector _t; // track direction - private final Hep3Vector _p; // track prediction - private final Hep3Vector _n; // normal to plane - private final Matrix _dm_dg; //Global derivaties of the local measurements - private final Matrix _dr_dm; //Derivatives of residuals w.r.t. measurement - private final Matrix _dr_dg; //Derivatives of residuals w.r.t. global parameters - - public GlobalDers(int layer, double umeas, double vmeas, double wmeas, Hep3Vector tDir, Hep3Vector tPred, Hep3Vector normal) { - _layer = layer; - _umeas = umeas; - _vmeas = vmeas; - _wmeas = wmeas; - _t = tDir; - _p = tPred; - _n = normal; - // Derivatives of residuals w.r.t. perturbed measurement - _dr_dm = getResDers(); - // Derivatives of perturbed measurement w.r.t. global parameters - _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' - //print dr_dm - //print 'dr_dg' - //print self.dr_dg - } - - /** - * Derivative of mt, the perturbed measured coordinate vector m w.r.t. - * to global parameters: u,v,w,alpha,beta,gamma - */ - private Matrix getMeasDers() { - - //Derivative of the local measurement for a translation in u - double dmu_du = 1.; - double dmv_du = 0.; - double dmw_du = 0.; - // Derivative of the local measurement for a translation in v - double dmu_dv = 0.; - double dmv_dv = 1.; - double dmw_dv = 0.; - // Derivative of the local measurement for a translation in w - double dmu_dw = 0.; - double dmv_dw = 0.; - double dmw_dw = 1.; - // Derivative of the local measurement for a rotation around u-axis (alpha) - double dmu_dalpha = 0.; - double dmv_dalpha = _p.z(); // self.wmeas - double dmw_dalpha = -1.0 * _p.y(); // -1.0 * self.vmeas - // Derivative of the local measurement for a rotation around v-axis (beta) - double dmu_dbeta = -1.0 * _p.z(); //-1.0 * self.wmeas - double dmv_dbeta = 0.; - double dmw_dbeta = _p.x(); //self.umeas - // Derivative of the local measurement for a rotation around w-axis (gamma) - double dmu_dgamma = _p.y(); // self.vmeas - double dmv_dgamma = -1.0 * _p.x(); // -1.0 * self.umeas - double dmw_dgamma = 0.; - // put into matrix - Matrix dm_dg = new Matrix(3, 6); - dm_dg.set(0, 0, dmu_du); - dm_dg.set(0, 1, dmu_dv); - dm_dg.set(0, 2, dmu_dw); - dm_dg.set(0, 3, dmu_dalpha); - dm_dg.set(0, 4, dmu_dbeta); - dm_dg.set(0, 5, dmu_dgamma); - dm_dg.set(1, 0, dmv_du); - dm_dg.set(1, 1, dmv_dv); - dm_dg.set(1, 2, dmv_dw); - dm_dg.set(1, 3, dmv_dalpha); - dm_dg.set(1, 4, dmv_dbeta); - dm_dg.set(1, 5, dmv_dgamma); - dm_dg.set(2, 0, dmw_du); - dm_dg.set(2, 1, dmw_dv); - dm_dg.set(2, 2, dmw_dw); - dm_dg.set(2, 3, dmw_dalpha); - dm_dg.set(2, 4, dmw_dbeta); - dm_dg.set(2, 5, dmw_dgamma); - //dmdg = np.array([[dmu_du, dmu_dv, dmu_dw, dmu_dalpha, dmu_dbeta, dmu_dgamma],[dmv_du, dmv_dv, dmv_dw, dmv_dalpha, dmv_dbeta, dmv_dgamma],[dmw_du, dmw_dv, dmw_dw, dmw_dalpha, dmw_dbeta, dmw_dgamma]]) - return dm_dg; - } - - /** - * Derivatives of the local perturbed residual w.r.t. the measurements m - * (u,v,w)' - */ - private Matrix getResDers() { - 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) { - delta = i == j ? 1. : 0.; - val = delta - _t.v()[i] * _n.v()[j] / tdotn; - dr_dm.set(i, j, val); - } - } - return dr_dm; - } - - /** - * Turn derivative matrix into @Milleparameter - * - * @param isTop - top or bottom track - * @return list of @Milleparameters - */ - private List<MilleParameter> getDers(boolean isTop) { - int transRot; - int direction; - int label; - double value; - List<MilleParameter> milleParameters = new ArrayList<MilleParameter>(); - int topBot = isTop == true ? 1 : 2; - for (int ip = 1; ip < 7; ++ip) { - if (ip > 3) { - transRot = 2; - direction = ((ip - 1) % 3) + 1; - } else { - transRot = 1; - direction = ip; - } - label = topBot * MilleParameter.half_offset + transRot * MilleParameter.type_offset + direction * MilleParameter.dimension_offset + _layer; - value = _dr_dg.get(0, ip - 1); - MilleParameter milleParameter = new MilleParameter(label, value, 0.0); - milleParameters.add(milleParameter); - } - return milleParameters; - } - - - /* - - class globalDers: - def __init__(self,layer,umeas,vmeas,wmeas, tDir, tPred, normal): - self.layer = layer # measurement direction - self.umeas = umeas # measurement direction - self.vmeas = vmeas # unmeasured direction - self.wmeas = wmeas # normal to plane - self.t = tDir # track direction - self.p = tPred # track prediction - self.n = normal # normal to plane - # Global derivaties of the local measurements - self.dm_dg = self.getMeasDers() - # Derivatives of residuals w.r.t. measurement - self.dr_dm = self.getResDers() - # Derivatives of residuals w.r.t. global parameters - self.dr_dg = np.dot(self.dr_dm, self.dm_dg) - #print 'dm_dg' - #print dm_dg - #print 'dr_dm' - #print dr_dm - #print 'dr_dg' - #print self.dr_dg - - def dump(self): - print 'globalDers:' - print 'layer ', self.layer - print 'umeas ', self.umeas, ' vmeas ', self.vmeas, ' wmeas ', self.wmeas - print 't ', self.t, ' p ', self.p, ' n ', self.n - print 'dm_dg\n',self.dm_dg, '\ndr_dm\n',self.dr_dm,'\ndr_dg\n',self.dr_dg - - def getDers(self,isTop): - half_offset = 10000 - translation_offset = 1000 - direction_offset = 100 - topBot = 1 - transRot = 1 - direction = 1 - if not isTop: - topBot = 2 - res = {} - labels = [] - ders = [] - for ip, name in global_params.iteritems(): - if ip > 3: - transRot = 2 - direction = ((ip-1) % 3) + 1 - else: - direction = ip - label = (int)(topBot * half_offset + transRot * translation_offset + direction * direction_offset + self.layer) - labels.append(label) - ders.append(self.dr_dg[0,ip-1]) - - return {'labels':np.array([labels]) , 'ders':np.array([ders])} - - - - - def getResDers(self): - # Derivatives of the local perturbed residual w.r.t. the measurements m (u,v,w)' - tdotn = np.dot(self.t.T,self.n) - drdg = np.eye(3) - #print 't ', self.t, ' n ', self.n, ' dot(t,n) ', tdotn - for i in range(3): - for j in range(3): - delta = 0. - if i==j: - delta = 1. - drdg[i][j] = delta - self.t[i]*self.n[j]/tdotn[0] - return drdg - - - - - def getMeasDers(self): - # Derivative of mt, the perturbed measured coordinate vector m - # w.r.t. to global parameters: u,v,w,alpha,beta,gamma - - # Derivative of the local measurement for a translation in u - dmu_du = 1. - dmv_du = 0. - dmw_du = 0. - # Derivative of the local measurement for a translation in v - dmu_dv = 0. - dmv_dv = 1. - dmw_dv = 0. - # Derivative of the local measurement for a translation in w - dmu_dw = 0. - dmv_dw = 0. - dmw_dw = 1. - # Derivative of the local measurement for a rotation around u-axis (alpha) - dmu_dalpha = 0. - dmv_dalpha = self.p[2] # self.wmeas - dmw_dalpha = -1.0 * self.p[1] # -1.0 * self.vmeas - # Derivative of the local measurement for a rotation around v-axis (beta) - dmu_dbeta = -1.0 * self.p[2] #-1.0 * self.wmeas - dmv_dbeta = 0. - dmw_dbeta = self.p[0] #self.umeas - # Derivative of the local measurement for a rotation around w-axis (gamma) - dmu_dgamma = self.p[1] # self.vmeas - dmv_dgamma = -1.0 * self.p[0] # -1.0 * self.umeas - dmw_dgamma = 0. - # put into matrix - dmdg = np.array([[dmu_du, dmu_dv, dmu_dw, dmu_dalpha, dmu_dbeta, dmu_dgamma],[dmv_du, dmv_dv, dmv_dw, dmv_dalpha, dmv_dbeta, dmv_dgamma],[dmw_du, dmw_dv, dmw_dw, dmw_dalpha, dmw_dbeta, dmw_dgamma]]) - #print dmw_dbeta - #dmdg = np.array([[dmu_du, dmu_dv],[dmu_dw, dmu_dalpha], [dmw_dbeta, dmw_dgamma]]) - return dmdg - */ + // Put the tracks back into the event and exit + int flag = 1 << LCIOConstants.TRBIT_HITS; + event.put(outputCollectionName, refittedTracks, Track.class, flag); } } 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 19:49:39 2015 @@ -26,7 +26,6 @@ import org.lcsim.event.base.BaseTrack; import org.lcsim.event.base.BaseTrackState; import org.lcsim.fit.helicaltrack.HelicalTrackFit; -import org.lcsim.fit.helicaltrack.HelicalTrackHit; import org.lcsim.lcio.LCIOConstants; import org.lcsim.recon.tracking.seedtracker.SeedCandidate; import org.lcsim.recon.tracking.seedtracker.SeedTrack; @@ -42,7 +41,7 @@ public class MakeGblTracks { private String _TrkCollectionName = "GBLTracks"; - private static Logger logger = LogUtil.create(MakeGblTracks.class, new BasicLogFormatter(), Level.INFO); + private static Logger logger = LogUtil.create(MakeGblTracks.class, new BasicLogFormatter(), Level.OFF); /** * Creates a new instance of MakeTracks.