7 added files
lcsim/src/org/lcsim/recon/tracking/seedtracker/diagnostic
diff -N AbstractSeedTrackerDiagnostics.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AbstractSeedTrackerDiagnostics.java 27 Aug 2008 17:57:46 -0000 1.1
@@ -0,0 +1,45 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.diagnostic;
+
+import org.lcsim.contrib.seedtracker.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.recon.tracking.seedtracker.HitManager;
+import org.lcsim.recon.tracking.seedtracker.MaterialManager;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+
+/**
+ *This class defines some common functionality that most classes implementing ISeedTrackerDiagnostics would want
+ * @author cozzy
+ */
+public abstract class AbstractSeedTrackerDiagnostics {
+
+ protected SeedStrategy currentStrategy;
+ protected double bField;
+ protected EventHeader event;
+ protected HitManager hitManager;
+ protected MaterialManager materialManager;
+
+ public void setBField(double bField) {
+ this.bField=bField;
+ }
+
+ public void setEvent(EventHeader event) {
+ this.event = event;
+ }
+
+ public void setHitManager(HitManager hm){
+ this.hitManager = hm;
+ }
+
+ public void setMaterialManager(MaterialManager mm){
+ this.materialManager=mm;
+ }
+
+ public void fireStrategyChanged(SeedStrategy strategy){
+ this.currentStrategy = strategy;
+ }
+}
lcsim/src/org/lcsim/recon/tracking/seedtracker/diagnostic
diff -N EmptySeedTrackerDiagnostics.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EmptySeedTrackerDiagnostics.java 27 Aug 2008 17:57:46 -0000 1.1
@@ -0,0 +1,41 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.diagnostic;
+
+import java.util.List;
+import org.lcsim.recon.tracking.seedtracker.ConfirmerExtender.Task;
+import org.lcsim.recon.tracking.seedtracker.HelixFitter;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFitter.FitStatus;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.line.SlopeInterceptLineFit;
+import org.lcsim.fit.zsegment.ZSegmentFit;
+
+/**
+ * This is the class you probably want to extend if you want to create your own diagnostics.
+ * As it extends from AbstractSeedTrackerDiagnostics, it already implements some common methods,
+ * so it is different from NullDiagnostics.
+ *
+ * @author cozzy
+ */
+public class EmptySeedTrackerDiagnostics extends AbstractSeedTrackerDiagnostics implements ISeedTrackerDiagnostics {
+
+ public void fireFinalDiagnostics(List<SeedCandidate> finalSeeds) {return;}
+ public void fireMergeIsBetterDiagnostics(SeedCandidate newSeed, SeedCandidate oldSeed, boolean better) {return;}
+ public void fireMergeIsDuplicateDiagnostics(SeedCandidate newSeed, SeedCandidate oldseed, boolean duplicate) {return;}
+ public void fireMergeStartDiagnostics(List<SeedCandidate> newSeedList) {return;}
+ public void fireFinderDPhiCut(int hitNumber, double dphi, boolean pass, HelicalTrackHit hit) {return;}
+ public void fireFinderDone(int maxseeds, int nseed, int nfit, int nconfirm, List<SeedCandidate> confirmedSeeds) {return;}
+ public void fireFitterFitMade(SeedCandidate seed, HelicalTrackFit hfit, CircleFit cfit, SlopeInterceptLineFit lfit, ZSegmentFit zfit, FitStatus fitStatus, boolean success) {return;}
+ public void fireConfirmerExtenderWorkingSeedInfo(Task task, SeedCandidate seed, List<HelicalTrackHit> hitlist){return;}
+ public void fireConfirmerExtenderFitNoSuccess(Task task, SeedCandidate Seed, HelicalTrackHit hit, HelixFitter fitter, boolean optimize) {return;}
+ public void fireConfirmerExtenderFitSuccess(Task task, SeedCandidate Seed, HelicalTrackHit hit, HelixFitter fitter, double chisqbest, boolean optimize) {return;}
+ public void fireConfirmerExtenderAllHitsPotentiallyBad(Task task, SeedCandidate workingSeed, List<HelicalTrackHit> hitlist, double chisqbest, double oldchisq) {return;}
+ public void fireConfirmerExtenderLayerDone(Task task, int numSeeds, List<SeedCandidate> seedlist) {return;}
+
+}
lcsim/src/org/lcsim/recon/tracking/seedtracker/diagnostic
diff -N FakeHelixFit.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FakeHelixFit.java 27 Aug 2008 17:57:46 -0000 1.1
@@ -0,0 +1,85 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.diagnostic;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.MCParticle;
+
+/**
+ *
+ * @author cozzy
+ */
+
+ //this is used to determine fit-like parameters for an MCParticle
+ public class FakeHelixFit{
+
+ private double z0;
+ private double dca;
+ private double phi0;
+ private double curv;
+ private double slope;
+ private Hep3Vector ref = new BasicHep3Vector(0.,0.,0.);
+
+ public FakeHelixFit(MCParticle part, double bfield){
+
+ Hep3Vector p = part.getMomentum();
+
+ double pt = Math.sqrt(p.x()*p.x() + p.y()*p.y());
+
+ //see (12) in http://www-flc.desy.de/lcnotes/notes/LC-DET-2006-004.pdf
+ curv = (bfield * 2.9979258e-4) / pt; //the magic number is c * unit conversion factor...
+
+ //change sign if necessary (i.e. sign of bfield and charge are different)
+ if (Math.signum(bfield) != Math.signum(part.getCharge()))
+ curv *=-1;
+
+ double r = 1/curv;
+
+ slope = p.z()/pt;
+
+ Hep3Vector o = part.getOrigin();
+
+ //find center of the circle...
+ double phi = Math.atan2(p.y(), p.x());
+ double xc = o.x() + Math.sin(phi)*r;
+ double yc = o.y() - Math.cos(phi)*r;
+
+ //calculate dca... should be the x-y distance from the reference point to the center of the circle minus the radius...
+ dca = Math.sqrt(xc*xc + yc*yc) - Math.abs(r);
+ dca *= Math.signum(curv); //i think that's the correct sign convention in this case...
+ phi0 = Math.atan2(xc/(r-dca),-yc/(r-dca));
+
+ //calculate z0... first we need the angle (alpha) between the lines connecting the reference point to the center of the circle and the
+ //line connecting the track origin to the center of the circle...
+
+ //x-y distance from origin of the particle to the reference point...
+ double s = (phi0 - phi)*r;
+
+ z0 = o.z() - s*slope;
+ }
+
+ public double z0() {
+ return z0;
+ }
+
+ public double dca() {
+ return dca;
+ }
+
+ public double phi0() {
+ return phi0;
+ }
+
+ public double curvature() {
+ return curv;
+ }
+
+ public double slope() {
+ return slope;
+ }
+
+ }
lcsim/src/org/lcsim/recon/tracking/seedtracker/diagnostic
diff -N FindableTracks.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FindableTracks.java 27 Aug 2008 17:57:46 -0000 1.1
@@ -0,0 +1,207 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.diagnostic;
+
+import org.lcsim.contrib.seedtracker.*;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.SpacePoint;
+import java.util.ArrayList;
+import java.util.Collection;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import org.lcsim.event.MCParticle;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.recon.tracking.seedtracker.HitManager;
+import org.lcsim.recon.tracking.seedtracker.SeedLayer;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+import org.lcsim.util.swim.Helix;
+import org.lcsim.util.swim.HelixSwimmer;
+
+/**
+ *
+ * @author cozzy
+ */
+public class FindableTracks {
+
+ private Set<MCParticle> goodMCs;
+ private double minPt;
+ private double maxDCA;
+ private double b_field = 5;
+ private double maxZ0;
+
+ private Hep3Vector ip = new BasicHep3Vector(0,0,0);
+
+ public FindableTracks(SeedStrategy strategy, HitManager hm){
+
+ goodMCs = determineGoodMCs(strategy, hm);
+ minPt = strategy.getMinPT();
+ maxDCA = strategy.getMaxDCA();
+ maxZ0 = strategy.getMaxZ0();
+ }
+ /**
+ * Checks whether the given MCParticle is findable according
+ * to the strategy given to this class. It checks the following:
+ * <ul>
+ * <li>That the particle has hits in all seeed layers
+ * <li>That the particle has sufficient confirm hits
+ * <li>That the particle has sufficient extension hits
+ * <li>That the particle's DCA is not too great
+ * <li>That the particle's Pt is not too small
+ * <li>That the particle's Z0 is not too great
+ * </ul>
+ * @param p MCParticle to check
+ * @return wheter or not the MCParticle is findable
+ */
+ public boolean isFindable(MCParticle p){
+
+ if (!goodMCs.contains(p)) return false;
+
+ double px = p.getMomentum().x();
+ double py = p.getMomentum().y();
+
+ double pt = Math.sqrt(px*px+py*py);
+
+ if (pt < minPt)
+ return false;
+
+
+ FakeHelixFit fhf = new FakeHelixFit(p, b_field);
+
+ if (Math.abs(fhf.z0()) > maxZ0 || Math.abs(fhf.dca()) > maxDCA)
+ return false;
+
+// HelixSwimmer h = new HelixSwimmer(b_field);
+//
+// h.setTrack(p.getMomentum(), new SpacePoint(p.getOrigin()),(int)p.getCharge());
+//
+// Helix hx = (Helix) h.getTrajectory();
+// double dca = hx.getSignedClosestDifferenceToPoint(ip);
+// double z0 = h.getPointAtLength(h.getTrackLengthToPoint(ip)).z() - ip.z();
+//
+
+// if (Math.abs(z0) > maxZ0 || Math.abs(dca) > maxDCA)
+// return false;
+
+ return true;
+
+ }
+
+ public void setBField(double b){
+ this.b_field = b;
+ }
+
+ //this method does the fun job of figuring out if the particle has sufficient hits in each type of layer
+ private Set<MCParticle> determineGoodMCs(SeedStrategy strat, HitManager hm){
+
+ //first, we determine all the MCParticles that go through all the seed layers
+ List<SeedLayer> seedLayers = strat.getLayers(SeedLayer.SeedType.Seed);
+ List<Set<MCParticle>> seedLayerMCs = new ArrayList<Set<MCParticle>>();
+
+ //make a set of MCParticles existing in each seed layer
+ for (SeedLayer sl: seedLayers) {
+ Set<MCParticle> thisLayerMCs = new HashSet<MCParticle>();
+ List<HelicalTrackHit> th = hm.getTrackerHits(sl);
+ for (HelicalTrackHit h : th)
+ thisLayerMCs.addAll(h.getMCParticles());
+ seedLayerMCs.add(thisLayerMCs);
+ }
+
+ if(seedLayerMCs.isEmpty()) return Collections.EMPTY_SET;
+
+ //take the first set of MCs and intersect it with all the sets
+ Set<MCParticle> allSeedMCs = seedLayerMCs.remove(0);
+ for (Set<MCParticle> otherSeedLayerMC : seedLayerMCs)
+ allSeedMCs.retainAll(otherSeedLayerMC);
+
+ if(allSeedMCs.isEmpty()) return Collections.EMPTY_SET; //shortcut out
+
+ //now we ensure we have the minimum number of confirmation hits
+ int minConfirm = strat.getMinConfirm();
+ List<SeedLayer> confirmLayers = strat.getLayers(SeedLayer.SeedType.Confirm);
+ OccurrenceCounter<MCParticle> confirmCounter = getMCCounter(confirmLayers, hm);
+
+ //intersect MCs in all the seeds with MCs that have enough confirm hits
+ allSeedMCs.retainAll(confirmCounter.getAllWithCountGreater(minConfirm));
+
+ if(allSeedMCs.isEmpty()) return Collections.EMPTY_SET; //shortcut out
+
+ //finally, check for enough extension hits...
+ int minExtend = strat.getMinHits() - minConfirm - seedLayers.size();
+ List<SeedLayer> extendLayers = strat.getLayers(SeedLayer.SeedType.Extend);
+ OccurrenceCounter<MCParticle> extendCounter = getMCCounter(extendLayers, hm);
+
+ //intersect MCs with all the seeds with MCs that have enough extend hits
+ allSeedMCs.retainAll(extendCounter.getAllWithCountGreater(minExtend));
+
+ return allSeedMCs;
+ }
+
+ private OccurrenceCounter<MCParticle> getMCCounter(List<SeedLayer> layers, HitManager hm) {
+
+ OccurrenceCounter<MCParticle> MCCounter = new OccurrenceCounter<MCParticle>();
+ for (SeedLayer l : layers) {
+
+ List<HelicalTrackHit> th = hm.getTrackerHits(l.getDetName(), l.getLayer(), l.getBarrelEndcapFlag());
+ for (HelicalTrackHit h : th) {
+ MCCounter.addAll(h.getMCParticles());
+ }
+ }
+
+ return MCCounter;
+ }
+
+ //utility class to count occurrences of stuff... A map of the object to
+ // number of occurrences
+ class OccurrenceCounter<T> {
+
+ Map<T,Integer> map;
+
+ public OccurrenceCounter(){
+ map = new HashMap<T,Integer>();
+ }
+
+ public void add(T obj){
+
+ if(map.containsKey(obj)){
+ int curr = map.get(obj);
+ map.put(obj, curr+1);
+ }
+ else map.put(obj, 1);
+ }
+
+ public void addAll(Collection<T> coll){
+ for (T t : coll)
+ add(t);
+ }
+
+ public int getCount(T obj) {
+
+ if (map.containsKey(obj))
+ return map.get(obj).intValue();
+ else return 0;
+ }
+
+ public Set<T> keySet() {
+ return map.keySet();
+ }
+
+ public Set<T> getAllWithCountGreater(int min) {
+
+ Set<T> ret = new HashSet<T>();
+ for (T t : map.keySet()) {
+ if (map.get(t).intValue() >= min)
+ ret.add(t);
+ }
+ return ret;
+ }
+ }
+
+}
lcsim/src/org/lcsim/recon/tracking/seedtracker/diagnostic
diff -N ISeedTrackerDiagnostics.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ISeedTrackerDiagnostics.java 27 Aug 2008 17:57:46 -0000 1.1
@@ -0,0 +1,205 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.diagnostic;
+
+import org.lcsim.contrib.seedtracker.*;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackFitter;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.line.SlopeInterceptLineFit;
+import org.lcsim.fit.zsegment.ZSegmentFit;
+import org.lcsim.recon.tracking.seedtracker.ConfirmerExtender;
+import org.lcsim.recon.tracking.seedtracker.HelixFitter;
+import org.lcsim.recon.tracking.seedtracker.HitManager;
+import org.lcsim.recon.tracking.seedtracker.MaterialManager;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+
+/**
+ *
+ * @author cozzy
+ */
+public interface ISeedTrackerDiagnostics {
+
+
+
+ /**
+ * SeedTracker will call this method with the event's bField
+ * @param bField the magnetic field, in Teslas
+ */
+ public void setBField(double bField);
+
+ /**
+ * SeedTracker will call this method with the current event
+ * @param event the current event
+ */
+ public void setEvent(EventHeader event);
+
+ /**
+ * SeedTracker will call this method with t
+ * @param hm the HitManager
+ */
+ public void setHitManager(HitManager hm);
+
+ /**SeedTracker will pass the MaterialManager with this method
+ *
+ * @param mm the MaterialManager
+ */
+ public void setMaterialManager(MaterialManager mm);
+
+ /**
+ * Fired whenever a new strategy is employed and passes the strategy.
+ * @param strategy
+ */
+ public void fireStrategyChanged(SeedStrategy strategy);
+
+ /**
+ * This is fired after merging with the final track lists. Any purity or efficiency diagnostics should go here.
+ * This is basically equivalent to extracting the list of tracks from the event.
+ *
+ * @param finalSeeds The final list of SeedCanditates determined by this strategy
+ */
+ public void fireFinalDiagnostics(List<SeedCandidate> finalSeeds);
+
+ /**
+ * This is fired whenever isBetter() is called in MergeSeedLists. The isBetter() function compares two
+ * SeedCandidates as follows:
+ *
+ * <ul>
+ * <li> if newSeed has 2 or more hits than oldSeed, return true;
+ * <li> if newSeed has 1 more hit than oldSeed, return true if the chi-squared
+ * difference is less than the badHitChisq, otherwise false
+ * <li>if newSeed has the same number of hits as oldSeed, return true
+ * if newSeed's fit has a better chi-squared, otherwise false
+ * <li>if newSeed has 1 hit less than oldSeed, return true if chi-squared
+ * difference is less than negative badHitChisq, otherwise false
+ * </ul>
+ *
+ * If returnValue is true, then newSeed should be better than oldSeed according to the strategy.
+ * This call gives the ability to check the performance of the isBetter algorithm.
+ * @param newSeed the new seed being compared
+ * @param oldSeed the old seed being compared against
+ * @param better the returnValue of the isBetter() function
+ */
+ public void fireMergeIsBetterDiagnostics(SeedCandidate newSeed, SeedCandidate oldSeed, boolean better);
+
+
+ /**
+ * This is fired when isDuplicated is called in MergeSeedLists. The
+ * isDuplicate() function returns true if it thinks that the two seedCandidates
+ * refer to the same track.
+ *
+ * @param newSeed the newSeed being compared
+ * @param oldSeed the oldSeed being compared
+ * @param duplicate the return value of isDuplicate
+ */
+ public void fireMergeIsDuplicateDiagnostics(SeedCandidate newSeed, SeedCandidate olDseed, boolean duplicate);
+
+ /**
+ * This is fired at the beginning of MergeSeedLists. It passes a list of the
+ * unmerged SeedCandidates (which means that many are likely duplicates).
+ *
+ * @param newSeedList the unmerged SeedCandidates
+ */
+ public void fireMergeStartDiagnostics(List<SeedCandidate> newSeedList);
+
+
+ /**
+ * This is fired by the finder whenever a dPhi cut is made.
+ *
+ * @param hitNumber Either 2 or 3, depending on whether it's the second
+ * or third hit being added
+ * @param dphi the difference in phi
+ * @param pass whether or not the hit passes
+ * @param hit the hit in question
+ */
+ public void fireFinderDPhiCut(int hitNumber, double dphi, boolean pass, HelicalTrackHit hit);
+
+ /**
+ * This is fired at the end of the finding routine
+ * @param maxseeds
+ * @param nseed
+ * @param nfit
+ * @param nconfirm
+ * @param confirmedSeeds
+ */
+ public void fireFinderDone(int maxseeds, int nseed, int nfit, int nconfirm, List<SeedCandidate> confirmedSeeds);
+
+
+ /**
+ * This is fired whenever a fit is made.
+ * @param seed The seed being fitted
+ * @param hfit The helicalTrack fit
+ * @param cfit The circle fit
+ * @param lfit The line fit
+ * @param zfit The ZSegment fit
+ * @param fitStatus The fit status
+ * @param success Whether or not the fit was a success
+ */
+
+ public void fireFitterFitMade(SeedCandidate seed, HelicalTrackFit hfit, CircleFit cfit, SlopeInterceptLineFit lfit, ZSegmentFit zfit, HelicalTrackFitter.FitStatus fitStatus, boolean success);
+
+
+ /**
+ * This is fired for each working seed that is being confirmed or extended at the beginning of confirmation / extension
+ * @param task Confirm or Extend
+ * @param seed The working seed
+ * @param hitlist The list of hits tob be used.
+ */
+ public void fireConfirmerExtenderWorkingSeedInfo(ConfirmerExtender.Task task, SeedCandidate seed, List<HelicalTrackHit> hitlist);
+
+ /**
+ * This is fired whenever the fit with a new hit fails
+ * @param task Confirm or extend
+ * @param seed The original working seed before the hit is added
+ * @param hit The hit added
+ * @param fitter The HelixFitter object
+ * @param optimize Whether the optimize bit is enabled or not
+
+
+ */
+ public void fireConfirmerExtenderFitNoSuccess(ConfirmerExtender.Task task, SeedCandidate seed, HelicalTrackHit hit, HelixFitter fitter, boolean optimize);
+
+
+
+
+ /**
+ * This is fired whenever the fit with a new hit succeeds
+ * @param task Confirm or extend
+ * @param seed The original working seed before the hit is added
+ * @param hit The hit added
+ * @param fitter The HelixFitter object
+ * @param chisqbest The current best chisq
+ * @param optimize Whether the optimize bit is enabled or not
+
+ */
+ public void fireConfirmerExtenderFitSuccess(ConfirmerExtender.Task task, SeedCandidate seed, HelicalTrackHit hit, HelixFitter fitter, double chisqbest, boolean optimize);
+
+
+
+
+ /**
+ * This is fired when all fit tries for a particular layer are potentially bad hits.
+ * @param task Confirm or extend
+ * @param workingSeed the current working seed
+ * @param hitlist the list of hits
+ * @param chisqbest the best chisq
+ * @param oldchisq the old chisq
+ */
+ public void fireConfirmerExtenderAllHitsPotentiallyBad(ConfirmerExtender.Task task, SeedCandidate workingSeed, List<HelicalTrackHit> hitlist, double chisqbest, double oldchisq);
+
+
+ /**
+ * This is fired at the end of confirmation.
+ * @param task Confirm or extend
+ * @param numSeeds the numer of seed confirmed or extended
+ * @param seedlist the final seed list
+ */
+ public void fireConfirmerExtenderLayerDone(ConfirmerExtender.Task task, int numSeeds, List<SeedCandidate> seedlist);
+}
lcsim/src/org/lcsim/recon/tracking/seedtracker/diagnostic
diff -N SeedTrackerDiagnostics.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SeedTrackerDiagnostics.java 27 Aug 2008 17:57:46 -0000 1.1
@@ -0,0 +1,247 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.diagnostic;
+
+import org.lcsim.contrib.seedtracker.*;
+import hep.physics.vec.Hep3Vector;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.Iterator;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import org.lcsim.recon.tracking.seedtracker.ConfirmerExtender.Task;
+import org.lcsim.event.MCParticle;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.math.chisq.ChisqProb;
+import org.lcsim.recon.tracking.seedtracker.HelixFitter;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author cozzy
+ */
+public class SeedTrackerDiagnostics extends EmptySeedTrackerDiagnostics implements ISeedTrackerDiagnostics{
+
+ AIDA aida = AIDA.defaultInstance();
+
+ private double purity_cutoff = 0.0;
+ private String prefix = "";
+ private boolean debugOut;
+
+ public SeedTrackerDiagnostics(){
+
+ }
+
+ /**
+ * Generates diagnostic plots given a list of seed candidates and a strategy.
+ * @param seedlist
+ * @param strategy
+ */
+ @Override
+ public void fireFinalDiagnostics(List<SeedCandidate> seedlist){
+
+ System.out.println(" After merging: "+seedlist.size());
+ Map<SeedCandidate,SeedValidator> vmap = generateValidators(seedlist);
+ makePurityPlots(vmap);
+ Set<MCParticle> likelyMCs = generateLikelyMCSet(vmap);
+ makeEfficiencyPlots(likelyMCs,currentStrategy);
+ }
+
+ /**
+ * Sets the purity cutoff for a track's MC to be considered found when
+ * calculating efficiency. Default is 0.0.
+ * @param new_purity_cutoff
+ */
+ public void setPurityCutoff(double new_purity_cutoff){
+ purity_cutoff = new_purity_cutoff;
+ }
+
+ /**
+ * Sets the prefix to be appended to the beginning of each plot (for example, a strategy name
+ * @param str the new prefix
+ */
+ public void setPrefix(String str){
+ this.prefix = str;
+ }
+
+
+ @Override
+ public void fireStrategyChanged(SeedStrategy strategy) {
+ setPrefix("chisq cut = "+strategy.getMaxChisq()+"/");
+ super.fireStrategyChanged(strategy);
+ }
+
+
+ @Override
+ public void fireMergeIsBetterDiagnostics(SeedCandidate newSeed, SeedCandidate oldSeed, boolean returnValue) {
+ if (returnValue){
+
+ SeedValidator oldV = new SeedValidator(oldSeed);
+ SeedValidator newV = new SeedValidator(newSeed);
+
+ if (oldV.getPurity() > newV.getPurity()){
+ makeBadDecisionPlots(newSeed,oldSeed, oldV,newV);
+ }
+ }
+ }
+
+ @Override
+ public void fireFinderDone(int maxseeds, int nseed, int nfit, int nconfirm, List<SeedCandidate> confirmedSeeds) {
+ System.out.println("Possible seeds: "+maxseeds);
+ System.out.print(" Trial seeds: "+nseed);
+ System.out.println(" Fitted seeds: "+nfit);
+ System.out.println(" Confirmed seeds: "+nconfirm);
+ System.out.println(" Confirmed seed candidates: "+confirmedSeeds.size());
+ }
+
+ @Override
+ public void fireConfirmerExtenderWorkingSeedInfo(Task task, SeedCandidate seed, List<HelicalTrackHit> hitlist) {
+ if (debugOut) {
+ System.out.println("Oldcirclechisq: "+seed.getHelix().chisq()[0]);
+ System.out.println("Oldchisq: "+seed.getHelix().chisqtot());
+ System.out.println("Oldhits: "+seed.getHits().size());
+ System.out.println("Old Helix: "+seed.getHelix().toString());
+ }
+ return;
+ }
+
+ @Override
+ public void fireConfirmerExtenderFitNoSuccess(Task task, SeedCandidate Seed, HelicalTrackHit hit, HelixFitter fitter, boolean optimize) {
+ if (debugOut) System.out.println("Circlechissq" + fitter.getCircleFit().chisq());
+ }
+
+ @Override
+ public void fireConfirmerExtenderFitSuccess(Task task, SeedCandidate Seed, HelicalTrackHit hit, HelixFitter fitter, double chisqbest, boolean optimize) {
+ if (debugOut) {
+ System.out.println("Good fit");
+ System.out.println("Chisq: "+fitter.getHelix().chisqtot());
+ System.out.println("Circle Chisq: "+fitter.getHelix().chisq()[0]);
+ System.out.println("New Helix: "+fitter.getHelix().toString());
+ }
+ }
+
+ @Override
+ public void fireConfirmerExtenderLayerDone(Task task, int numSeeds, List<SeedCandidate> seedlist) {
+ if (debugOut){
+ System.out.println(" "+task.toString()+" seeds: "+numSeeds);
+ System.out.println(" "+task.toString()+" seed candidates: "+seedlist.size());
+ }
+ }
+
+ private void makeBadDecisionPlots(SeedCandidate newseed, SeedCandidate seed, SeedValidator oldV, SeedValidator newV) {
+
+ aida.cloud1D(prefix+"Bad Decision newseed chisq").fill(newseed.getHelix().chisqtot());
+ aida.cloud1D(prefix+"Bad Decision oldseed chisq").fill(seed.getHelix().chisqtot());
+ aida.cloud2D(prefix+"Bad Decision newseed vs. oldseed chisq").fill(newseed.getHelix().chisqtot(), seed.getHelix().chisqtot());
+ aida.cloud1D(prefix+"Bad Decision old seedpurity").fill(oldV.getPurity());
+ aida.cloud1D(prefix+"Bad Decision new seedpurity").fill(newV.getPurity());
+
+
+ if(oldV.getVerdict().isGoodValue()){
+ Hep3Vector p = oldV.getLikelyMC().getMomentum();
+ double pt = Math.sqrt(p.x() * p.x() + p.y() * p.y());
+ aida.cloud1D(prefix+"Bad Decision real particle transverse momentum").fill(pt);
+
+ Hep3Vector r = oldV.getLikelyMC().getOrigin();
+ double rc = Math.sqrt(r.x()*r.x()+r.y()*r.y());
+ aida.cloud1D(prefix + "Bad Decision real particle start radius (cylindrical)").fill(rc);
+ }
+
+ }
+
+
+ //makes purity plots
+ private void makePurityPlots(Map<SeedCandidate, SeedValidator> vmap) {
+
+ for (SeedCandidate s : vmap.keySet()) {
+
+ SeedValidator v = vmap.get(s);
+
+ if (v.getVerdict().isGoodValue()) {
+ aida.cloud1D(prefix+"Good seeds chisq").fill(s.getHelix().chisqtot());
+ aida.cloud1D(prefix+"Good seeds numhits").fill(s.getHits().size());
+ aida.cloud1D(prefix+"Good seeds pt").fill(s.getHelix().pT(5.0));
+ } else {
+ aida.cloud1D(prefix+"Bad seeds chisq").fill(s.getHelix().chisqtot());
+ aida.cloud1D(prefix+"Bad seeds numhits").fill(s.getHits().size());
+ aida.cloud1D(prefix+"Bad seeds pt").fill(s.getHelix().pT(5.0));
+ }
+
+ aida.cloud1D(prefix+"purity").fill(v.getPurity());
+ aida.cloud2D(prefix+"purity vs. numHits").fill(v.getPurity(), s.getHits().size());
+ aida.cloud2D(prefix+"purity vs. pt").fill(v.getPurity(), s.getHelix().pT(5.0));
+ aida.cloud2D(prefix+"purity vs. cth").fill(v.getPurity(), s.getHelix().cth());
+ aida.cloud2D(prefix+"purity vs. chisq").fill(v.getPurity(), s.getHelix().chisqtot());
+
+ double chisqdof = ChisqProb.gammq(s.getHelix().chisqtot(), s.getHits().size()-1);
+
+ aida.cloud2D(prefix+"purity vs. (1 - chisq_cdf(chisq,ndof)").fill(v.getPurity(), chisqdof);
+ aida.cloud2D(prefix+"purity cs. dca").fill(v.getPurity(), s.getHelix().dca());
+ }
+ }
+
+ private void makeEfficiencyPlots(Set<MCParticle> likelyMCs, SeedStrategy strat){
+
+ FindableTracks findable = new FindableTracks(strat, hitManager);
+ findable.setBField(bField);
+
+ List<MCParticle> MCs = new ArrayList<MCParticle>();
+ MCs.addAll(event.getMCParticles());
+
+ Iterator iter = MCs.iterator();
+
+ //remove all none findable MCs
+ while(iter.hasNext()){
+ if(!findable.isFindable((MCParticle)iter.next()))
+ iter.remove();
+ }
+
+ int numFindable = MCs.size();
+ aida.cloud1D(prefix+"Number findable").fill(numFindable);
+
+
+ // remove all the MC's we've found... this gives us findable - found
+ MCs.removeAll(likelyMCs);
+
+ int numFound = numFindable - MCs.size(); // findable - (findable - found) = found
+
+ aida.cloud1D(prefix+"Num found (purity cutoff="+purity_cutoff+") ").fill(numFound);
+
+ if(numFindable > 0){
+ double eff = (double) numFound / (double) numFindable;
+ aida.cloud1D(prefix+"Efficiency").fill(eff);
+ }
+
+ }
+
+
+ private Map<SeedCandidate,SeedValidator> generateValidators(List<SeedCandidate> seedlist) {
+
+ Map<SeedCandidate,SeedValidator> ret = new HashMap<SeedCandidate,SeedValidator>();
+ for (SeedCandidate s : seedlist) {
+ ret.put(s, new SeedValidator(s));
+ }
+
+ return ret;
+ }
+
+ private Set<MCParticle> generateLikelyMCSet(Map<SeedCandidate, SeedValidator> vmap) {
+
+ Set<MCParticle> set = new HashSet<MCParticle>();
+
+ for (SeedValidator v : vmap.values()) {
+ if(v.getPurity()<purity_cutoff) continue;
+ set.add(v.getLikelyMC());
+ }
+ return set;
+ }
+}
+
\ No newline at end of file
lcsim/src/org/lcsim/recon/tracking/seedtracker/diagnostic
diff -N SeedValidator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SeedValidator.java 27 Aug 2008 17:57:46 -0000 1.1
@@ -0,0 +1,201 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.recon.tracking.seedtracker.diagnostic;
+
+import org.lcsim.contrib.seedtracker.*;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.event.MCParticle;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+
+/**
+ *
+ * @author cozzy
+ */
+public class SeedValidator {
+
+ private SeedCandidate candidate;
+ private MCParticle likelyMC;
+ private Map<MCParticle,Integer> occurs;
+ private SeedVerdict verdict;
+ private double purity=0;
+
+ public enum SeedVerdict{
+ SEED_VALID_BUT_MULTIPLE_MCS_IN_HITS,
+ SEED_VALID,
+ SEED_INVALID,
+ MULTIPLE_POSSIBLE_MCS,
+ GET_MCPARTICLES_FAILED;
+
+ public boolean isGoodValue(){
+ return (this==SEED_VALID || this == SEED_VALID_BUT_MULTIPLE_MCS_IN_HITS);
+ }
+ }
+
+
+ public SeedValidator(SeedCandidate s){
+ candidate = s;
+ occurs = getOccurrences(s);
+ evaluateOccurrences();
+ }
+
+ private void evaluateOccurrences() {
+
+ //this will happen if getMCParticles fails for whatever reason;
+ if(occurs == null) {
+ verdict = SeedVerdict.GET_MCPARTICLES_FAILED;
+ return;
+ }
+
+ //if there's only one MCParticle in the map, then that's the one we want.
+ if(occurs.size()==1) {
+ likelyMC = (MCParticle) occurs.keySet().toArray()[0];
+ verdict = SeedVerdict.SEED_VALID;
+ purity = 1;
+ return;
+ }
+
+
+ //get a list of the MCParticles that appear in all the hits
+ List<MCParticle> mcs_in_all_hits = new ArrayList<MCParticle>();
+ int size = candidate.getHits().size();
+
+ int max_hits = 0;
+ for (MCParticle p : occurs.keySet()){
+ int num_hits = occurs.get(p).intValue();
+ if(num_hits==size) mcs_in_all_hits.add(p);
+ if(num_hits>=max_hits){
+ max_hits=num_hits;
+ likelyMC = p;
+ }
+ }
+
+ purity = ((double)max_hits)/((double)size);
+
+ //if no MCParticles meet this criteria, we have a baad seed.
+ if (mcs_in_all_hits.size()==0){
+ verdict = SeedVerdict.SEED_INVALID;
+ return;
+ }
+
+ //if there are multiple MCs meeting this criteria, that's also problematic
+ if (mcs_in_all_hits.size()>1){
+ verdict = SeedVerdict.MULTIPLE_POSSIBLE_MCS;
+ return;
+ }
+
+ verdict = SeedVerdict.SEED_VALID;
+
+ }
+
+
+ /**
+ * Returns an enum of type SeedVerdict encoding the status of the MC...
+ * @return
+ */
+ public SeedVerdict getVerdict(){
+ return verdict;
+ }
+
+ /**
+ * Returns the purity of the candidate. The purity is defined as the number
+ * of hits containing the likelyMC (see the {@link #getLikelyMC() getLikelyMC}
+ * method divided by the total number of hits.
+ * @return
+ */
+ public double getPurity(){
+ return purity;
+ }
+ /**
+ * Returns the most likely MC belonging to the seed. If there is an MC that
+ * all the hits contain, that MC is returned. Otherwise, the MC with the
+ * most hits containing it is returned. If more than one MC has the most hits,
+ * only one will be returned, though it's not easy to predict whigh.
+ * @return the most likely MC
+ */
+ public MCParticle getLikelyMC(){
+ return likelyMC;
+ }
+
+
+ public enum NewAdditionVerdict{
+ ADDITION_IS_GOOD,
+ ADDITION_IS_GOOD_BUT_MULTIPLE_MCS,
+ ADDITION_IS_BAD,
+ SEED_IS_BAD_OR_CANNOT_DETERMINE;
+
+ public boolean isGoodValue(){
+ return (this==ADDITION_IS_GOOD ||
+ this==ADDITION_IS_GOOD_BUT_MULTIPLE_MCS);
+ }
+ }
+
+
+ /**
+ * Determines whether or not this hit is a good hit to add.
+ *
+ * @param A helical track hit.
+ * @return a NewAdditionVerdict enum
+ */
+
+ public NewAdditionVerdict isGoodAddition(HelicalTrackHit h){
+ if (likelyMC==null)
+ return NewAdditionVerdict.SEED_IS_BAD_OR_CANNOT_DETERMINE;
+
+ List<MCParticle> parts = h.getMCParticles();
+
+ if (parts==null)
+ return NewAdditionVerdict.SEED_IS_BAD_OR_CANNOT_DETERMINE;
+
+ if(!parts.contains(likelyMC))
+ return NewAdditionVerdict.ADDITION_IS_BAD;
+
+ if(parts.size()==1)
+ return NewAdditionVerdict.ADDITION_IS_GOOD;
+
+ return NewAdditionVerdict.ADDITION_IS_GOOD_BUT_MULTIPLE_MCS;
+ }
+
+
+
+ /**
+ * Given a seed candidate, this returns a map of MCParticles to the number
+ * of hits in the seed candidate containing them.
+ *
+ * Returns null if getMCParticles fails on a hit.
+ *
+ * This is static because it might be useful for some other purpose...
+ *
+ * @param s a SeedCandidate
+ * @return a map of MCParticles to number of hits containing that MCParticle,
+ * or null if getMCParticles fails.
+ */
+ public static Map<MCParticle,Integer> getOccurrences(SeedCandidate s){
+
+ Map<MCParticle,Integer> occurrences = new HashMap<MCParticle,Integer>();
+
+ for (HelicalTrackHit h : s.getHits()){
+
+ List<MCParticle> mcs = h.getMCParticles();
+ if (mcs == null) return null;
+
+ for (MCParticle p : mcs) {
+ if(!occurrences.containsKey(p)){
+ occurrences.put(p, 1);
+ }
+ else {
+ int currVal = occurrences.get(p);
+ occurrences.put(p,currVal+1);
+ }
+ }
+ }
+
+ return occurrences;
+ }
+}
CVSspam 0.2.8