lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -N TrackerHitTrackFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackerHitTrackFinder.java 19 Jul 2007 23:16:32 -0000 1.1
@@ -0,0 +1,425 @@
+/*
+ * TrackerHitTrackFinder.java
+ *
+ * Created on July 18, 2007, 1:53 PM
+ *
+ * To change this template, choose Tools | Template Manager
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.CosminDeaconu;
+
+import hep.physics.vec.BasicHep3Vector;
+import java.util.Arrays;
+import java.util.Collection;
+import java.util.Collections;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.List;
+import org.lcsim.contrib.tracking.TrackerHitCheater;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseTrackerHitMC;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.util.Driver;
+import java.util.ArrayList;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author Cosmin Deaconu
+ */
+
+//===============================================================================//
+// Driver to find helical tracks from a collection of TrackerHits
+public class TrackerHitTrackFinder extends Driver{
+
+ private ArrayList<TrackerHit> _input_hit_collection = null;
+ private HashMap<TrackerHit,Double> _hit_separation = new HashMap<TrackerHit,Double>();
+
+ private ArrayList<TrackerHit> _hitset = new ArrayList<TrackerHit>();
+ private ArrayList<StandaloneOuterTrack> _tracks = new ArrayList<StandaloneOuterTrack>();
+ private ArrayList<TrackerHit> _used = new ArrayList<TrackerHit>();
+
+ private HelicalTrackFitter _fitter = new HelicalTrackFitter();
+ private HelicalTrackFit _fit = null;
+
+ double chisq_l = 10;
+ double chisq_c = 10;
+
+ int min_hits = 4;
+ int max_hits = 8;
+
+ private AIDA _aida = AIDA.defaultInstance();
+ private boolean _make_histograms = true;
+ private String _id = "";
+
+ double required_distance = 0.3;
+ double required_separation = 0.5;
+
+ boolean require_different_layers = false;
+
+ double drphi = 0.03;
+ double dz = 100;
+
+ double[] origin = {0.,0.,0.};
+
+
+ /** Creates a new instance of TrackerHitTrackFinder */
+ public TrackerHitTrackFinder() {
+
+ }
+
+ protected void process(EventHeader event) {
+
+ // addSimTrackerHitInput((ArrayList<SimTrackerHit>) event.getSimTrackerHits("TkrBarrHits"));
+ addSimTrackerHitInput((ArrayList<SimTrackerHit>) event.getSimTrackerHits("TkrForwardHits"));
+ addSimTrackerHitInput((ArrayList<SimTrackerHit>) event.getSimTrackerHits("TkrEndcapHits"));
+
+
+ if (_input_hit_collection==null)
+ {
+ System.out.println("Warning: TrackerHitTrackFinder has no input hits.");
+ return;
+ }
+
+
+ DistanceCompare<TrackerHit> c = new DistanceCompare<TrackerHit>();
+ Collections.sort(_input_hit_collection,c); // sort from outside in... yeah this may take a while if there are a lot of hits
+
+
+ // get hit separation
+ for (TrackerHit hit : _input_hit_collection)
+ {
+ _hit_separation.put(hit,findNearest(hit,_input_hit_collection));
+ }
+
+ for (TrackerHit hit1 : _input_hit_collection)
+ {
+ for (TrackerHit hit2 : _input_hit_collection)
+ {
+ if (hit1==hit2) continue;
+ if (require_different_layers && sameLayer(hit1,hit2)) continue;
+
+ for (TrackerHit hit3 : _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;
+
+ _hitset.clear();
+ _hitset.add(hit1);
+ _hitset.add(hit2);
+ _hitset.add(hit3);
+
+ _fit = doFit(_hitset);
+
+ makeHistogram("chisql_1",_fit.chisq()[1]);
+ //check line chisquare to see if it's reasonable'
+ if (_fit.chisq()[1]/3.0>chisq_l) continue;
+ makeHistogram("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)
+ {
+ if (_hitset.contains(new_hit)) continue; //make sure not already used
+
+
+ //check separation/distance
+ if (Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(new_hit.getPosition())))) > 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;
+
+ 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);
+ }
+
+ makeHistogram("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);
+
+ makeHistogram("chisql_3",_fit.chisq()[1]);
+ makeHistogram("chisqc_1",_fit.chisq()[0]);
+
+ int while_loop_iterations = 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)
+ {
+ if (_hitset.size()<min_hits) break;
+
+ int size = _hitset.size();
+ double farthest_distance = -1;
+
+ int removeMe=-1;
+
+ for (int i = 0; i<candidates.size(); i++)
+ {
+ double distance = Math.abs(_fit.getHelix().getSignedClosestDifferenceToPoint(new SpacePoint(new BasicHep3Vector(candidates.get(i).getPosition()))));
+
+ if (distance > farthest_distance)
+ {
+ farthest_distance = distance;
+ removeMe = i;
+ }
+ }
+
+ if(!_hitset.remove(candidates.get(removeMe))) System.out.println("Infinite Loop! Diagnostic information: "+removeMe+" "+candidates.size());
+ candidates.remove(removeMe);
+
+ while_loop_iterations++;
+
+ //fit again
+ _fit = doFit(_hitset);
+ }
+
+ makeHistogram("while loop iter",while_loop_iterations);
+
+ //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);
+
+ makeHistogram("Number of hits per track",_hitset.size());
+ makeHistogram("Purity",purity);
+ }
+ }
+ }
+ }
+ event.put(EventHeader.TRACKS, _tracks, Track.class, 0); //put tracks into event
+ makeHistogram("n_tracks",_tracks.size());
+ purgeLists(); //do some cleanup
+ }
+
+
+
+//-------------------------Public Controller Methods--------------------------------------------------------------
+ public void addInput(ArrayList<TrackerHit> hits)
+ {
+ if (_input_hit_collection==null)
+ {
+ _input_hit_collection=hits;
+ }
+
+ else
+ {
+ _input_hit_collection.addAll(hits);
+ }
+ }
+
+ public void addSimTrackerHitInput(ArrayList<SimTrackerHit> hits)
+ {
+ TrackerHitCheater _cheat = new TrackerHitCheater();
+ addInput ((ArrayList<TrackerHit>) _cheat.makeTrackerHits(hits));
+ }
+
+//--------------------------Private Helper Methods-----------------------------------------------------
+
+ private double findNearest(TrackerHit hit, List<TrackerHit> hits)
+ {
+ double returnme = 100000000.0;
+
+ for (TrackerHit i : hits)
+ {
+ double d = distance(((BaseTrackerHitMC)hit).getPosition(),((BaseTrackerHitMC)i).getPosition());
+ if (d < returnme)
+ {
+ returnme = d;
+ }
+ }
+ return returnme;
+ }
+
+ private double distance(double[] a, double[] b) //this is probably implemented elsewhere, but whatevs
+ {
+ return Math.sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])+(a[2]-b[2])*(a[2]-b[2]));
+ }
+
+ private double[] extractCoordinate(int coord, ArrayList<TrackerHit> h)
+ {
+ double[] returnme = new double[h.size()];
+
+ for (int i =0; i<h.size(); i++)
+ {
+ returnme[i] = (h.get(i).getPosition())[coord];
+ }
+
+ return returnme;
+ }
+
+
+ private double[] get_drphi(int nhits) //returns an array of the right size with the right drphi value
+ {
+ double[] returnme = new double[nhits];
+ for (int i = 0; i<nhits;i++)
+ {returnme[i]=drphi;}
+ return returnme;
+ }
+
+ private double[] get_dz(int nhits) //similar to above, but with dz
+ {
+ double[] returnme = new double[nhits];
+ for (int i = 0; i<nhits;i++)
+ {returnme[i]=dz;}
+
+ return returnme;
+ }
+
+ private boolean sameLayer(TrackerHit a, TrackerHit b)
+ {
+ int a_layer = ((BaseTrackerHitMC) a).getSimHits().get(0).getLayer();
+ int b_layer = ((BaseTrackerHitMC) b).getSimHits().get(0).getLayer();
+
+ String a_detector = ((BaseTrackerHitMC) a).getSimHits().get(0).getSubdetector().getName();
+ String b_detector = ((BaseTrackerHitMC) b).getSimHits().get(0).getSubdetector().getName();
+
+ return (a_layer==b_layer && a_detector.equals(b_detector));
+ }
+
+ private HelicalTrackFit doFit(ArrayList<TrackerHit> hitset)
+ {
+ boolean fitok = _fitter.fit(extractCoordinate(0,hitset),extractCoordinate(1,hitset),extractCoordinate(2,hitset),get_drphi(hitset.size()),get_dz(hitset.size()),hitset.size());
+ if (!fitok){System.out.println("WARNING: Couldn't get fit");}
+ return _fitter.getFit();
+ }
+
+ private MCParticle findMajorityParticle(ArrayList<TrackerHit> hits)
+ {
+ HashMap<MCParticle,Integer> particles = new HashMap<MCParticle,Integer>();
+
+ for (TrackerHit i : hits )
+ {
+ for (MCParticle j : (((BaseTrackerHitMC) i).mcParticles()))
+ {
+ if (particles.containsKey(j))
+ {
+ particles.put(j,particles.get(j)+1);
+ }
+
+ else
+ {
+ particles.put(j,1);
+ }
+ }
+ }
+
+ MCParticle currentMax = null;
+ int max_value = -1;
+
+
+ for (MCParticle i : particles.keySet())
+ {
+ if (particles.get(i).intValue()>max_value)
+ {
+ max_value = particles.get(i).intValue();
+ currentMax = i;
+ }
+ }
+
+ return currentMax;
+ }
+
+ private double findPurity(ArrayList<TrackerHit> hits, MCParticle mp)
+ {
+ int denom = hits.size();
+ int num = 0;
+
+ 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))
+ {
+ num++;
+ }
+ }
+ return ((double) num)/((double) denom);
+ }
+
+ private void makeHistogram(String name, double value)
+ {
+ if (_make_histograms)
+ {
+ _aida.cloud1D(name+_id).fill(value);
+ }
+
+ }
+
+ private void purgeLists()
+ {
+ _input_hit_collection.clear();
+ _hitset.clear();
+ _hit_separation.clear();
+ _used.clear();
+ _tracks.clear();
+ }
+}
+
+//comparator class used to sort points from outermost to innermost
+class DistanceCompare<TrackerHit> implements Comparator<TrackerHit>
+{
+ public int compare(Object a, Object b)
+ {
+ double[] origin = {0.0,0.0,0.0};
+
+ double[] a_pos = ((BaseTrackerHitMC)a).getPosition();
+ double[] b_pos = ((BaseTrackerHitMC)b).getPosition();
+
+ double a_distance = distance(a_pos,origin);
+ double b_distance = distance(b_pos,origin);
+
+ if (a_distance>b_distance) return -1;
+ else if (b_distance>a_distance) return 1;
+ else return 0;
+ }
+
+ private double distance(double[] a, double[] b)
+ {
+ return Math.sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])+(a[2]-b[2])*(a[2]-b[2]));
+ }
+}
+
+
lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -u -r1.2 -r1.3
--- OuterTrackFinder.java 19 Jul 2007 00:05:14 -0000 1.2
+++ OuterTrackFinder.java 19 Jul 2007 23:16:32 -0000 1.3
@@ -153,7 +153,7 @@
}
- /**HelicalTrackFit
+ /* process
* Performs the event processing for standalone track finding. The steps are:
* <ul>
* <li> Seed generation from sets of three isolated hits
@@ -282,7 +282,7 @@
}
}
- //System.out.println(layer_converter.toString());
+ //System.out.println(layer_converter.toString());
nlayers = used.size(); //set nlayers to the number of layers actually used
}
@@ -890,16 +890,10 @@
} // if chisq not crap
} //if hitphi is appropriate
} // if more than 4 layers used
-
-
} // if 3-point dca < 0.5
-
-
-
} // hit3 loop
} // hit2 loop
} // hit1 loop
-
} // layer3 loop
} // layer2 loop
} // layer1 loop
@@ -1152,10 +1146,11 @@
public void useDefaultBarrelSettings() //~Cosmin
{
_useAllLayers=true;
- setInputHits("TkrBarrelHits");
+ setInputHits("TkrBarrHits");
setSubdetectors(new String[]{"TrackerBarrel"});
+ dz = ((double)_module_length)/(Math.sqrt(12.0));
_phi_cut = true;
- _h_chisq_dof[1]=10000000000000000000000000000000000000000.0; // emulate circle fit by setting ridiculously high line chisq cutoff for helical fit.
+ //_h_chisq_dof[1]=1000000000000000000000000000000000000000000.0; // emulate circle fit by setting ridiculously high line chisq cutoff for helical fit.
_zseg=true;
}
@@ -1169,6 +1164,12 @@
{_phi_cut = usePhiCut;}
+ public void set_dz(double new_dz)
+ {dz = new_dz;};
+
+ public void set_drphi(double new_drphi)
+ {drphi = new_drphi;}
+
/**
* Enable/Disable using zSegmentation
*
lcsim/src/org/lcsim/contrib/CosminDeaconu
diff -u -r1.1 -r1.2
--- StandaloneOuterTrack.java 18 Jul 2007 19:01:06 -0000 1.1
+++ StandaloneOuterTrack.java 19 Jul 2007 23:16:32 -0000 1.2
@@ -94,7 +94,7 @@
_reference_point[2] = 0.0;
_tracker_hits = hits;
- double[] _ip = {0.,0.,0.};
+ double[] _ip = {0.,0.,0.};
_track_parameters[0] = fit.parameters()[0];
_track_parameters[1] = fit.parameters()[1]-Math.PI; // who ordered this?
_track_parameters[2] = fit.parameters()[2]; // convert 1/mm to 1/cm
@@ -110,6 +110,35 @@
}
+ public StandaloneOuterTrack(double b_field, HelicalTrackFit fit, List<TrackerHit> hits, MCParticle majority_particle, double purity)
+ {
+ _b_field = b_field;
+
+ _fit_success = true;
+ _h_chi2 = fit.chisq();
+ _helical_fit=true;
+ _ndof = hits.size();
+
+ _reference_point[0] = 0.; //yeah, this is bad
+ _reference_point[1] = 0.; //yeah, this is bad ~/Cosmin
+ _reference_point[2] = 0.0;
+
+ _tracker_hits = hits;
+ double[] _ip = {0.,0.,0.};
+ _track_parameters[0] = fit.parameters()[0];
+ _track_parameters[1] = fit.parameters()[1]-Math.PI; // who ordered this?
+ _track_parameters[2] = fit.parameters()[2]; // convert 1/mm to 1/cm
+ _track_parameters[3] = fit.parameters()[3];
+ _track_parameters[4] = fit.parameters()[4];
+
+ _isreferencepoint_pca = true;
+
+ _majority_particle = majority_particle;
+ _purity = purity;
+
+ _tracks.add(this);
+
+ }
// Methods needed to further initialize a StandaloneAxialBarrelTrack