lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -u -r1.3 -r1.4
--- TrackerHitTrackFinder.java 25 Jul 2007 00:24:17 -0000 1.3
+++ TrackerHitTrackFinder.java 27 Jul 2007 06:19:37 -0000 1.4
@@ -1,4 +1,4 @@
-/*
+/*
* TrackerHitTrackFinder.java
*
* Created on July 18, 2007, 1:53 PM
@@ -22,6 +22,8 @@
import org.lcsim.event.Track;
import org.lcsim.event.base.BaseTrackerHitMC;
import org.lcsim.event.TrackerHit;
+import org.lcsim.fit.line.SlopeInterceptLineFit;
+import org.lcsim.fit.line.SlopeInterceptLineFitter;
import org.lcsim.spacegeom.SpacePoint;
import org.lcsim.util.Driver;
import java.util.ArrayList;
@@ -35,15 +37,22 @@
*
* Method:
* * Sorts all hits by distance from center
- * * Goes through all combinations of 3 hits (starting from outside)
- * * For each combination of 3 hits, a helix is fit, and a list of candidate
- * hits is generated that lie close to the helix but are not too close to
- * other hits
- * * If the chisq for the candidates + the 3 seed hits is good enough, they
- * form a track, otherwise the farthest hit from the helix is removed until
- * either the chisq is good enough or there are too few hits.
- * * Tracks must satisfy a minimum number of hits (and a maximum number of hits too, which can be set arbitrarily high)
- * * Optionally, hits may be required to be on different layers
+ * * Goes through all pairs of points, and then finds the best third point to
+ * form a line (in the s-z plane). These are the seed points.
+ * * A helical fit is made from the seed points
+ * * If the chisq of the seed helix is good enough, a list of candidates is made.
+ * The candidates must be well-separated, close to the seed helix, and each of
+ * their chisqs with the seed must be good.
+ * * All candidates are then added to the track, and if the chisq sucks,
+ * the farthest hit is removed until either there are too few points
+ * or the recalculated chisq is good enough.
+ * * All hits must be close enough to the new helix. If a candidate is too far,
+ * it is removed and the helix refit.
+ * * If a seed is too far from the helix after the previous step, we give up on the
+ * seed and try a new one.
+ * * Tracks must satisfy a minimum number of hits (and a maximum number of hits too,
+ * which can be set arbitrarily high)
+ *
*
* Usage:
* Generally, you'll want to use the public methods to give the driver an input collection of hits (from the event).
@@ -66,12 +75,13 @@
private ArrayList<TrackerHit> _used = new ArrayList<TrackerHit>();
private HelicalTrackFitter _fitter = new HelicalTrackFitter();
private HelicalTrackFit _fit = null;
+ private SlopeInterceptLineFitter _lfitter = new SlopeInterceptLineFitter();
private String _outputSimTrackerHitCollection = "";
private String _outputTrackerHitCollection = "";
private double _chisq_l = 10; //chisq cutoff for the linear fit portion of the helical tracker
- private double _chisq_c = 10; //chisq cutoff for the circular fit portion of the helical tracker
+ private double _chisq_c = 1; //chisq cutoff for the circular fit portion of the helical tracker
private int _min_hits = 4; //minimum number of hits for a track
private int _max_hits = 200; // max number of hits for a track
@@ -137,204 +147,214 @@
{
if (hit1==hit2) continue;
if (_require_different_layers && sameLayer(hit1,hit2)) continue;
-
- for (TrackerHit hit3 : _input_hit_collection)
+ if (_used.contains(hit1) || _used.contains(hit2)) continue;
+
+ //for hit 3, find the best linear chisq
+ TrackerHit hit3 = null;
+ double minchisq = 100000000;
+
+ for (TrackerHit hit : _input_hit_collection)
{
- if (hit1==hit3 || hit2==hit3) continue;
- if (_require_different_layers && (sameLayer(hit1,hit3) || sameLayer(hit2,hit3))) continue;
- if (_used.contains(hit1)||_used.contains(hit2)||_used.contains(hit3)) continue; //don't reuse hits
-
- _hitset.clear();
- _hitset.add(hit1);
- _hitset.add(hit2);
- _hitset.add(hit3);
-
- _fit = doFit(_hitset);
+ if (hit1==hit || hit2==hit) continue;
+ if (_require_different_layers && (sameLayer(hit1,hit) || sameLayer(hit2,hit))) continue;
+ if (_used.contains(hit)) continue; //don't reuse hits
- generateHistogram("chisql_1",_fit.chisq()[1]);
-
- //check line chisquare to see if it's reasonable'
- if (_fit.chisq()[1]/3.0>_chisq_l) continue;
- generateHistogram("chisql_2",_fit.chisq()[1]);
-
- ArrayList<TrackerHit> candidates = new ArrayList<TrackerHit>();
-
- //look for well separated new hits that are candidates to attach
- for (TrackerHit new_hit : _input_hit_collection)
+ _fit = doFit(new ArrayList<TrackerHit>(Arrays.asList(new TrackerHit[]{hit1,hit2,hit})));
+ if (_fit.chisq()[1] < minchisq)
{
- if (_used.contains(new_hit)) continue;
- if (_hitset.contains(new_hit)) continue; //make sure not already used
-
-
- //check separation/distance
-
- double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(new_hit.getPosition()))));
-
- if ( distance> _required_distance || _hit_separation.get(new_hit)<_required_separation) continue;
-
-
- //try adding it to hitset and assess chisq
-
- _hitset.add(new_hit);
- _fit = doFit(_hitset);
-
- //check new chisq's and see if they're reasonable'
- if (_fit.chisq()[1]/((double)_hitset.size())>_chisq_l || _fit.chisq()[0]/((double)_hitset.size())>_chisq_c)
- {
- _hitset.remove(new_hit);
- continue;
- }
-
- //if requiring different layers, don't add anything in a layer already used... this is probably not the best way to do this
- if (_require_different_layers)
- {
- boolean fail = false;
- if (sameLayer(new_hit,hit1) || sameLayer(new_hit,hit2) || sameLayer(new_hit,hit3)) fail=true;
+ minchisq = _fit.chisq()[1];
+ hit3 = hit;
+ }
+ }
+
+ _hitset.clear();
+ _hitset.add(hit1);
+ _hitset.add(hit2);
+ _hitset.add(hit3);
- for (TrackerHit i : candidates)
- {
- if (sameLayer(new_hit,i)) fail=true;
- }
- if (fail) continue;
- }
-
- //if tests so far passed, add new_hit to candidates list
- candidates.add(new_hit);
- _hitset.remove(new_hit);
- }
-
- generateHistogram("candidates_size",candidates.size());
-
- ////first, try adding all candidates
- _hitset.addAll(candidates);
-
- if (_hitset.size() < _min_hits || _hitset.size() > _max_hits) continue; //if too many or too few hits, skip
-
+ ArrayList<TrackerHit> candidates = new ArrayList<TrackerHit>();
+
+ //look for well separated new hits that are candidates to attach
+ for (TrackerHit new_hit : _input_hit_collection)
+ {
+ if (_used.contains(new_hit)) continue;
+ if (_hitset.contains(new_hit)) continue; //make sure not already used
+
+
+ //check separation/distance
+
+ double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(new_hit.getPosition()))));
+
+ if ( distance> _required_distance || _hit_separation.get(new_hit)<_required_separation) continue;
+
+
+ //try adding it to hitset and assess chisq
+
+ _hitset.add(new_hit);
_fit = doFit(_hitset);
-
- generateHistogram("chisql_3",_fit.chisq()[1]);
- generateHistogram("chisqc_1",_fit.chisq()[0]);
-
- //if chisquared is too high, remove hits farthest from helix until chisquared is acceptable
- while (_fit.chisq()[1]/((double)_hitset.size())>_chisq_l || _fit.chisq()[0]/((double)_hitset.size())>_chisq_c)
+
+ //check new chisq's and see if they're reasonable'
+ if (_fit.chisq()[1]/((double)_hitset.size())>_chisq_l || _fit.chisq()[0]/((double)_hitset.size())>_chisq_c)
{
- if (_hitset.size()<_min_hits) break;
-
- int size = _hitset.size();
- double farthest_distance = -1;
-
- TrackerHit removeMe=null;
-
+ _hitset.remove(new_hit);
+ continue;
+ }
+
+ //if requiring different layers, don't add anything in a layer already used... this is probably not the best way to do this
+ if (_require_different_layers)
+ {
+ boolean fail = false;
+ if (sameLayer(new_hit,hit1) || sameLayer(new_hit,hit2) || sameLayer(new_hit,hit3)) fail=true;
+
for (TrackerHit i : candidates)
{
- double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(i.getPosition()))));
+ if (sameLayer(new_hit,i)) fail=true;
+ }
+ if (fail) continue;
+ }
+
+ //if tests so far passed, add new_hit to candidates list
+ candidates.add(new_hit);
+ _hitset.remove(new_hit);
+ }
+
+ generateHistogram("candidates_size",candidates.size());
+
+ ////first, try adding all candidates
+ _hitset.addAll(candidates);
+
+ if (_hitset.size() < _min_hits || _hitset.size() > _max_hits) continue; //if too many or too few hits, skip
+
+ _fit = doFit(_hitset);
+
+ generateHistogram("chisql_3",_fit.chisq()[1]);
+ generateHistogram("chisqc_1",_fit.chisq()[0]);
- if (distance > farthest_distance)
- {
- farthest_distance = distance;
- removeMe = i;
- }
+ //if chisquared is too high, remove hits farthest from helix until chisquared is acceptable
+ while (_fit.chisq()[1]/((double)_hitset.size())>_chisq_l || _fit.chisq()[0]/((double)_hitset.size())>_chisq_c)
+ {
+ if (_hitset.size()<_min_hits) break;
+
+ int size = _hitset.size();
+ double farthest_distance = -1;
+
+ TrackerHit removeMe=null;
+
+ for (TrackerHit i : candidates)
+ {
+ double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(i.getPosition()))));
+
+ if (distance > farthest_distance)
+ {
+ farthest_distance = distance;
+ removeMe = i;
}
-
- if(!_hitset.remove(removeMe)) System.out.println("Infinite Loop! Diagnostic information: "+removeMe+" "+candidates.size());
- candidates.remove(removeMe);
-
- //fit again
- _fit = doFit(_hitset);
}
-
- //now check if any candidate hits are too far and remove them if they are, refitting after each removal...
-
- boolean keepgoing = true;
-
- while(keepgoing) //while loop so that we start over whenever we remove a candidate
+
+ if(!_hitset.remove(removeMe)) System.out.println("Infinite Loop! Diagnostic information: "+removeMe+" "+candidates.size());
+ candidates.remove(removeMe);
+
+ //fit again
+ _fit = doFit(_hitset);
+ }
+
+ //now check if any candidate hits are too far and remove them if they are, refitting after each removal...should be a better way to write this by sorting the points or something
+
+ boolean keepgoing = true;
+
+ while(keepgoing) //while loop so that we start over whenever we remove a candidate
+ {
+ keepgoing = false; //assume all candidates fall within critical distance
+ ArrayList<TrackerHit> removeFromCandidates = new ArrayList<TrackerHit>(); //since java won't let you remove an element from a List you're looping over...'
+
+ for (TrackerHit i : candidates)
{
- keepgoing = false; //assume all candidates fall within critical distance
- ArrayList<TrackerHit> removeFromCandidates = new ArrayList<TrackerHit>(); //since java won't let you remove an element from a List you're looping over...'
-
- for (TrackerHit i : candidates)
+ double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(i.getPosition()))));
+
+ if(distance>_required_distance_second_pass) //check distance to make sure it's reasonable'
{
- double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(i.getPosition()))));
+ _hitset.remove(i);
+ _fit = doFit(_hitset); //refit
+ removeFromCandidates.add(i);
+ keepgoing = true;
+ break;
+ }
+ }
+ candidates.removeAll(removeFromCandidates);
+ }
+
+ double max_distance = -1; //make sure that the maximum distance from the helix isn't too great...this is mostly for the seed points, I guess'
+
+ for (TrackerHit i : _hitset)
+ {
+ double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(i.getPosition()))));
+ if (distance>max_distance) max_distance=distance;
+ }
+
+ //now, I think we have a track
+ if (_hitset.size()>=_min_hits && max_distance<_required_distance_second_pass)
+ {
+ MCParticle majority_particle = findMajorityParticle(_hitset);
+ double purity = findPurity(_hitset,majority_particle);
+ double b_field = event.getDetector().getFieldMap().getField(_origin)[2];
+
+ StandaloneOuterTrack track = new StandaloneOuterTrack(b_field,_fit,_hitset,majority_particle,purity);
+
+ _tracks.add(track);
+ _used.addAll(_hitset);
+ purities.add(purity);
+
+ //We love plots!
+ generateHistogram("Number of hits per track",_hitset.size());
+ generateHistogram("purity",purity);
+ generateScatterPlot("purity vs. hits per track",purity,_hitset.size());
+ generateScatterPlot("purity vs. chisql per nhits",purity,_fit.chisq()[1]/_hitset.size());
+ generateScatterPlot("purity vs. chisqc per nhits",purity,_fit.chisq()[0]/_hitset.size());
+ generateScatterPlot("purity vs. majority_particle momentum",purity,majority_particle.getMomentum().magnitude());
+ generateScatterPlot("purity vs. z momentum",purity,majority_particle.getMomentum().z());
+ generateScatterPlot("purity vs. theta",purity,Math.atan(_fit.parameters()[4]));
+ generateScatterPlot("purity vs. phi",purity,_fit.parameters()[1]);
+ generateScatterPlot("purity vs. z",purity,_fit.parameters()[3]);
+ generateScatterPlot("purity vs. curvature",purity,_fit.parameters()[2]);
+ generateScatterPlot("purity vs. dca",purity,_fit.parameters()[0]);
+
+ HashMap<String,Integer> det_hits = new HashMap<String,Integer>();
+
+ double total_distance = 0;
+ double min_sep = 100000000;
+ double total_sep = 0;
+
+ for (TrackerHit i : _hitset) //calculate some info for plotting purposes
+ {
+ double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(i.getPosition()))));
+ total_distance = total_distance + distance;
+ if (distance>max_distance) max_distance=distance;
+
+ double separation = _hit_separation.get(i).doubleValue();
+ if (separation<min_sep) min_sep = separation;
+ total_sep = total_sep + separation;
- if(distance>_required_distance_second_pass) //check distance to make sure it's reasonable'
- {
- _hitset.remove(i);
- _fit = doFit(_hitset); //refit
- removeFromCandidates.add(i);
- keepgoing = true;
- break;
- }
+ String detector = ((BaseTrackerHitMC)i).getSimHits().get(0).getSubdetector().getName();
+ if (!det_hits.containsKey(i))
+ {
+ det_hits.put(detector,1);
}
- candidates.removeAll(removeFromCandidates);
- }
-
- //now, I think we have a track
- if (_hitset.size()>=_min_hits)
- {
- MCParticle majority_particle = findMajorityParticle(_hitset);
- double purity = findPurity(_hitset,majority_particle);
- double b_field = event.getDetector().getFieldMap().getField(_origin)[2];
-
- StandaloneOuterTrack track = new StandaloneOuterTrack(b_field,_fit,_hitset,majority_particle,purity);
-
- _tracks.add(track);
- _used.addAll(_hitset);
- purities.add(purity);
-
- //We love plots!
- generateHistogram("Number of hits per track",_hitset.size());
- generateHistogram("purity",purity);
- generateScatterPlot("purity vs. hits per track",purity,_hitset.size());
- generateScatterPlot("purity vs. chisql per nhits",purity,_fit.chisq()[1]/_hitset.size());
- generateScatterPlot("purity vs. chisqc per nhits",purity,_fit.chisq()[0]/_hitset.size());
- generateScatterPlot("purity vs. majority_particle momentum",purity,majority_particle.getMomentum().magnitude());
- generateScatterPlot("purity vs. z momentum",purity,majority_particle.getMomentum().z());
- generateScatterPlot("purity vs. theta",purity,Math.atan(_fit.parameters()[4]));
- generateScatterPlot("purity vs. phi",purity,_fit.parameters()[1]);
- generateScatterPlot("purity vs. z",purity,_fit.parameters()[3]);
- generateScatterPlot("purity vs. curvature",purity,_fit.parameters()[2]);
- generateScatterPlot("purity vs. dca",purity,_fit.parameters()[0]);
-
- HashMap<String,Integer> det_hits = new HashMap<String,Integer>();
-
-
- double max_distance = -1; //this will store the distance of the point farthest from the track.
- double total_distance = 0;
- double min_sep = 100000000;
- double total_sep = 0;
-
- for (TrackerHit i : _hitset) //calculate some info for plotting purposes
- {
- double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(i.getPosition()))));
- total_distance = total_distance + distance;
- if (distance>max_distance) max_distance=distance;
-
- double separation = _hit_separation.get(i).doubleValue();
- if (separation<min_sep) min_sep = separation;
- total_sep = total_sep + separation;
-
- String detector = ((BaseTrackerHitMC)i).getSimHits().get(0).getSubdetector().getName();
- if (!det_hits.containsKey(i))
- {
- det_hits.put(detector,1);
- }
-
- else
- {
- Integer already = det_hits.get(detector);
- det_hits.put(detector,already+1);
- }
+
+ else
+ {
+ Integer already = det_hits.get(detector);
+ det_hits.put(detector,already+1);
}
-
- generateScatterPlot("purity vs. number of subdetectors",purity,det_hits.keySet().size());
- generateScatterPlot("purity vs. max distance",purity,max_distance);
- generateScatterPlot("purity vs. average distance",purity,total_distance/((double)_hitset.size()));
- generateScatterPlot("purity vs. min seperation",purity,min_sep);
- generateScatterPlot("purity vs. average separation",purity,total_sep/((double)_hitset.size()));
}
- }
- }
- }
+
+ generateScatterPlot("purity vs. number of subdetectors",purity,det_hits.keySet().size());
+ generateScatterPlot("purity vs. max distance",purity,max_distance);
+ generateScatterPlot("purity vs. average distance",purity,total_distance/((double)_hitset.size()));
+ generateScatterPlot("purity vs. min seperation",purity,min_sep);
+ generateScatterPlot("purity vs. average separation",purity,total_sep/((double)_hitset.size()));
+ } // if
+ } //for 2
+ } //for 1
event.put(EventHeader.TRACKS, new ArrayList<StandaloneOuterTrack>(_tracks), Track.class, 0); //put tracks into event
@@ -422,6 +442,11 @@
public void setRequiredDistance(double d)
{_required_distance=d;}
+ /**
+ * Sets the maximum distance a hit may be from the final helix
+ *
+ * @param d double value for new required distance
+ */
public void setRequiredDistanceSecondPass(double d)
{_required_distance_second_pass=d;}
@@ -614,9 +639,17 @@
for (TrackerHit i : hits )
{
+
+ ArrayList<MCParticle> thisHitParticles = new ArrayList<MCParticle>();
for (MCParticle j : (((BaseTrackerHitMC) i).mcParticles()))
{
- if (particles.containsKey(j))
+ if (!thisHitParticles.contains(j)) thisHitParticles.add(j);
+ }
+
+ for (MCParticle j: thisHitParticles)
+ {
+ if (particles.containsKey(j))
+
{
particles.put(j,particles.get(j)+1);
}
@@ -650,8 +683,7 @@
for (TrackerHit i : hits)
{
- MCParticle hit_mp = findMajorityParticle(new ArrayList(Arrays.asList(new TrackerHit[]{i}))); //finds the particle contributing most to a given hit
- if (hit_mp.equals(mp))
+ if (((BaseTrackerHitMC)i).mcParticles().contains(mp)) // check to see if the hit contains the majority particle
{
num++;
}