lcsim/src/org/lcsim/contrib/tracking
diff -N AxialBarrelTrackFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AxialBarrelTrackFinder.java 22 Sep 2005 19:37:52 -0000 1.1
@@ -0,0 +1,631 @@
+/*
+ * AxialBarrelTrackFinder.java
+ *
+ * Created on August 18, 2005, 4:17 PM
+ *
+ */
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.util.Driver;
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.fit.circle.CircleFitter;
+
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.MultiLayerTracker;
+import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.geometry.subdetector.TrackerIDDecoder;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.*;
+
+import java.util.*;
+import java.lang.Math;
+
+/**
+ * AxialBarrelTrackFinder is a Driver that performs track pattern recognition
+ * using only the outer tracking layers: the VXD is ignored. This
+ * algorithm works from the outside inward, using all sets of isolated
+ * hits in three layers to find circles that pass close to the
+ * interaction point. When such a combination is found, the remaining
+ * layers are checked for nearby hits: first only hits that are
+ * isolated, but after refitting with all isolated hits, any hit close
+ * to the CircleFit is added. The resulting CircleFit objects are
+ * used to create StandaloneAxialBarrelTrack objects which are then stored in the
+ * event record. AxialBarrelTrackFinder currenltly only works in the barrel
+ * region but can be exteded to forward tracking in the future. Also,
+ * the algorithm currently uses unsmeared SimTrackerHits, but can be
+ * modified to use real TrackerHit objects when available.
+ *
+ * The code is currently very basic on not highly optimized. However,
+ * the various requirements and assumptions can be controlled by
+ * by the user: see the method descriptions. A number of
+ * histograms are also generated, which may be turned on at the user's
+ * descretion.
+ *
+ * @author Tim Nelson - SLAC
+ * @version %I%, %G%
+ * @see <a href="http://nicadd.niu.edu/cdsagenda//askArchive.php?base=agenda&categ=a0567&id=a0567s1t0/moreinfo">Standalone Tracking Talk - 9/2/05</a>
+ * @since 1.5
+ */
+public class AxialBarrelTrackFinder extends Driver
+{
+
+ // Data members
+ private String _input_hit_collection = "TkrBarrHits";
+ private String _output_hit_collection = "";
+ private double _module_length = 100.0;
+ private double _seedhit_isolation = 0.5;
+ private double _seed_ip_dca = 2.0;
+ private double _pass1_hit_dca = 0.5;
+ private double _pass1_hit_isolation = 1.0;
+ private double _pass2_hit_dca = 0.25;
+ private double _chisq_dof = 10.0;
+ private boolean _make_histograms = false;
+
+ private List<SimTrackerHit> _hits = null;
+ private Vector<List<SimTrackerHit>> _hitsbylayer = new Vector<List<SimTrackerHit>>();
+ private List<SimTrackerHit> _used_hits = new ArrayList<SimTrackerHit>();
+ private HashMap<SimTrackerHit,Double> _hit_separation = new HashMap<SimTrackerHit,Double>();
+
+ private LCMetaData _metadata = null; //event.getMetaData(hits);
+ private TrackerIDDecoder _decoder = null; //(TrackerIDDecoder) metadata.getIDDecoder();
+
+ private AIDA _aida = AIDA.defaultInstance();
+
+ private CircleFitter _fitter = new CircleFitter();
+ private CircleFit _fit = null;
+
+ /**
+ * Creates a new instance of AxialBarrelTrackFinder
+ */
+ public AxialBarrelTrackFinder()
+ {
+ }
+
+ /**
+ * Performs the event processing for standalone track finding. The steps are:
+ * <ul>
+ * <li> Seed generation from sets of three isolated hits
+ * <li> Addition of isolated hits in remaining layers, refitting after each addition
+ * <li> Addition of any nearby hits in remaining layers
+ * <li> Final fitting
+ * <li> Creation of StandaloneAxialBarrelTrack objects
+ * </ul>
+ *
+ * @param event The event header
+ */
+ protected void process(EventHeader event)
+ {
+
+ // Get detector and print out tracking subdetector names
+ Detector detector = event.getDetector();
+ Set<String> detnames = detector.getSubdetectorNames();
+
+ // Get z component of magnetic field
+ double[] ip = {0.,0.,0.};
+ double b_field = detector.getFieldMap().getField(ip)[2];
+
+ // Get the barrel tracker attributes
+ MultiLayerTracker barrel = (MultiLayerTracker)detector.getSubdetector("TrackerBarrel");
+ double[] radius = barrel.getInnerR();
+ double[] halflength = barrel.getOuterZ();
+ int nlayers = radius.length;
+
+ // Get the SimTrackerHits and metadata
+ _hits = event.getSimTrackerHits(_input_hit_collection);
+ _metadata = event.getMetaData(_hits);
+ _decoder = (TrackerIDDecoder) _metadata.getIDDecoder();
+
+ // Print out number of hits
+ if (_make_histograms) _aida.cloud1D("nHitsTotal").fill(_hits.size());
+
+ // Create/clear vector of SimTrackerHits for each layer
+ if (_hitsbylayer.size()==0)
+ {
+ for (int layer = 0; layer < nlayers; layer++)
+ {
+ _hitsbylayer.add(layer, new Vector<SimTrackerHit>());
+ }
+ }
+ else
+ {
+ for (int layer = 0; layer < nlayers; layer++)
+ {
+ _hitsbylayer.get(layer).clear();
+ }
+ }
+
+ for (SimTrackerHit hit : _hits)
+ {
+ _decoder.setID(hit.getCellID());
+ int layer = _decoder.getLayer();
+ _hitsbylayer.elementAt(layer).add(hit);
+ }
+
+ // Create map from hits to separation from nearest hit in layer
+ _hit_separation.clear();
+ for (int layer = 0; layer < nlayers; layer++)
+ {
+ List<SimTrackerHit> hits = _hitsbylayer.get(layer);
+ for (SimTrackerHit hit : hits)
+ {
+ _hit_separation.put(hit,rphiDistToNearest(hit,hits));
+ }
+ }
+
+ // Get the MCParticles
+ List<MCParticle> mc_particles = event.getMCParticles();
+ HashSet<MCParticle> all_layer_particles = new HashSet<MCParticle>();
+
+ int ntracks_all_layers = 0;
+ for (MCParticle mc_particle : mc_particles) {
+ boolean layers_hit[] = new boolean[nlayers];
+ for (SimTrackerHit hit : _hits) {
+ _decoder.setID(hit.getCellID());
+ int layer = _decoder.getLayer();
+ if (mc_particle == hit.getMCParticle())
+ {
+ layers_hit[layer] = true;
+ }
+ }
+
+ boolean all_layers = true;
+ for (int layer = 0; layer < nlayers ; layer++) {
+ all_layers = all_layers && layers_hit[layer];
+ }
+
+ if (all_layers) {
+ ntracks_all_layers++;
+ all_layer_particles.add(mc_particle);
+ if (_make_histograms) _aida.cloud1D("5-hit MCParticle momentum").fill(mc_particle.getMomentum().magnitude());
+ }
+ }
+
+ // Create set of found MCParticles
+ HashSet<MCParticle> found_particles = new HashSet<MCParticle>();
+
+ // Clear list of already used hits
+ _used_hits.clear();
+
+ // Count number of found tracks in this event
+ int nfound = 0;
+
+ // Create list of tracks
+ ArrayList<StandaloneAxialBarrelTrack> tracklist = new ArrayList<StandaloneAxialBarrelTrack>();
+
+ // Try all three-layer combinations from outside-in
+ for (int layer1 = nlayers-1; layer1 >= 0; layer1--)
+ {
+ for (int layer2 = layer1-1; layer2 >= 0; layer2--)
+ {
+ for (int layer3 = layer2-1; layer3 >= 0; layer3 --)
+ {
+
+ Vector<SimTrackerHit> hitset = new Vector<SimTrackerHit>();
+ Set<Integer> used_layers = new HashSet<Integer>();
+
+ List<SimTrackerHit> hits1 = _hitsbylayer.get(layer1);
+ List<SimTrackerHit> hits2 = _hitsbylayer.get(layer2);
+ List<SimTrackerHit> hits3 = _hitsbylayer.get(layer3);
+
+ double[] x = new double[5];
+ double[] y = new double[5];
+ double[] weight = new double[5];
+
+ // Try all combinations of hits
+ for (SimTrackerHit hit1 : hits1)
+ {
+ if (_hit_separation.get(hit1) < _seedhit_isolation) continue; // require well-separated hits
+ if (_used_hits.contains(hit1)) continue; // skip used hits
+ boolean zpos = (hit1.getPoint()[2] > 0.0); // find +/- z
+ double zrel = (hit1.getPoint()[2]/radius[layer1]);
+ x[0] = hit1.getPoint()[0];
+ y[0] = hit1.getPoint()[1];
+ weight[0] = 1000.0; // set to achieve reasonable chisq
+ for (SimTrackerHit hit2 : hits2)
+ {
+ if (_hit_separation.get(hit2) < _seedhit_isolation) continue; // require well-separated hits
+ if (_used_hits.contains(hit2)) continue; // skip used hits
+ if ((hit2.getPoint()[2] > 0.0) != zpos) continue; // require same side in z
+ if (Math.abs(hit2.getPoint()[2] - zrel*radius[layer2]) > _module_length/2.0) continue;
+ x[1] = hit2.getPoint()[0];
+ y[1] = hit2.getPoint()[1];
+ weight[1] = 1000.0; // set to achieve reasonable chisq
+ for (SimTrackerHit hit3 : hits3)
+ {
+ if (_hit_separation.get(hit3) < _seedhit_isolation) continue; // require well-separated hits
+ if (_used_hits.contains(hit3)) continue; // skip used hits
+ if ((hit3.getPoint()[2] > 0.0) != zpos) continue; // require same side in z
+ if (Math.abs(hit3.getPoint()[2] - zrel*radius[layer3]) > _module_length/2.0) continue;
+ x[2] = hit3.getPoint()[0];
+ y[2] = hit3.getPoint()[1];
+ weight[2] = 1000.0; // set to achieve reasonable chisq
+
+ // Perform circle fit with these three points, marking layers as used
+ int nhits = 3;
+
+ hitset.clear(); hitset.setSize(nlayers);
+ hitset.setElementAt(hit1,layer1);
+ hitset.setElementAt(hit2,layer2);
+ hitset.setElementAt(hit3,layer3);
+
+ used_layers.clear();
+ used_layers.add(layer1);
+ used_layers.add(layer2);
+ used_layers.add(layer3);
+
+ // perform fit to three points
+ boolean fitok = _fitter.fit(x, y, weight, nhits);
+ if (!fitok) System.out.println("Warning: seed fit failed!");
+ _fit = _fitter.getfit();
+
+ // If close to IP, loop over other layers from outside-in and attach unambiguous hits
+ if (Math.abs(_fit.dca())<_seed_ip_dca)
+ {
+
+ // Fill seed histograms
+ if (_make_histograms)
+ {
+ _aida.cloud1D("dca_seed").fill(_fit.dca());
+ _aida.cloud1D("curvature_seed").fill(_fit.curvature());
+ }
+
+ for (int layer = nlayers-1; layer >= 0; layer--)
+ {
+
+ if (used_layers.contains(layer)) continue;
+
+ // Add nearest hit if unambiguous -- check order of logic here
+ List<SimTrackerHit> hits = _hitsbylayer.get(layer);
+ double nearest_dist = 1000000.0;
+ double next_nearest_dist = 1000000.0;
+ SimTrackerHit nearest_hit = null;
+
+ for (SimTrackerHit hit : hits)
+ {
+ if (_used_hits.contains(hit)) continue; // skip used hits
+ if ((hit.getPoint()[2] > 0.0) != zpos) continue; // require same side in z
+ if (Math.abs(hit.getPoint()[2] - zrel*radius[layer]) > 50.0) continue;
+ _fit = _fitter.propagatefit(hit.getPoint()[0],hit.getPoint()[1]);
+ if (Math.abs(_fit.dca()) < nearest_dist)
+ {
+ next_nearest_dist = nearest_dist;
+ nearest_dist = Math.abs(_fit.dca());
+ nearest_hit = hit;
+ }
+ }
+ // Make histograms of nearest hits
+ if (_make_histograms)
+ {
+ if (nearest_dist < 10.0) {
+ _aida.cloud1D("dca_add1_nearest").fill(nearest_dist);
+ }
+ if (next_nearest_dist < 10.0)
+ {
+ _aida.cloud1D("dca_add1_nextnearest").fill(next_nearest_dist);
+ }
+ }
+
+ // require nearby hit that is well separated
+ if (nearest_dist < _pass1_hit_dca && next_nearest_dist > _pass1_hit_isolation)
+ {
+ used_layers.add(layer); // Mark layer used
+ hitset.setElementAt(nearest_hit,layer);
+ x[nhits] = nearest_hit.getPoint()[0];
+ y[nhits] = nearest_hit.getPoint()[1];
+ weight[nhits] = 1000.0; // set to achieve reasonable chisq
+ nhits++;
+
+ // Refit after adding each hit
+ fitok = _fitter.fit(x,y, weight, nhits);
+ if (!fitok) System.out.println("Warning: intermediate fit failed!");
+
+ }
+ }
+
+ // If still have unused layers, try for remaining hits
+ //----------------------------------------------------
+ if ( used_layers.size()<5 )
+ {
+
+ // Check layers from outside in
+ for (int layer = nlayers-1; layer >= 0; layer--)
+ {
+ if (used_layers.contains(layer)) continue;
+
+ // Find nearest hit
+ List<SimTrackerHit> hits = _hitsbylayer.get(layer);
+ double nearest_dist = 1000000.0; // 1km
+ double next_nearest_dist = 1000000.0; SimTrackerHit nearest_hit = null;
+ for (SimTrackerHit hit : hits)
+ {
+ if (_used_hits.contains(hit)) continue; // skip used hits
+ if ((hit.getPoint()[2] > 0.0) != zpos) continue; // require same side in z
+ if (Math.abs(hit.getPoint()[2] - zrel*radius[layer]) > 50.0) continue;
+ _fit = _fitter.propagatefit(hit.getPoint()[0],hit.getPoint()[1]);
+ if (Math.abs(_fit.dca()) < nearest_dist)
+ {
+ next_nearest_dist = nearest_dist;
+ nearest_dist = Math.abs(_fit.dca());
+ nearest_hit = hit;
+ }
+ }
+ // Make histograms of nearest hits
+ if (_make_histograms)
+ {
+ if (nearest_dist < 10.0) {
+ _aida.cloud1D("dca_add2_nearest").fill(nearest_dist);
+ }
+ if (next_nearest_dist < 10.0)
+ {
+ _aida.cloud1D("dca_add2_nextnearest").fill(next_nearest_dist);
+ }
+ }
+
+ // require nearby hit
+ if (nearest_dist < _pass2_hit_dca)
+ {
+ used_layers.add(layer); // Mark layer used
+ hitset.setElementAt(nearest_hit,layer);
+ x[nhits] = nearest_hit.getPoint()[0];
+ y[nhits] = nearest_hit.getPoint()[1];
+ weight[nhits] = 1000.0; // set to achieve reasonable chisq
+ nhits++;
+ }
+
+ }
+
+ } // if less than 5 layers
+
+ // Reset reference position
+ _fitter.setreferenceposition(0.0,0.0);
+
+ // If have successfully added hits, perform final fit
+ if (used_layers.size() >= 4)
+ {
+ fitok = _fitter.fit(x,y, weight, nhits);
+ if (!fitok) System.out.println("Warning: final fit failed!");
+
+ _fit = _fitter.getfit();
+
+ // If fit is good enough, fill final histograms
+ if (_fit.chisq()/nhits < _chisq_dof)
+ {
+ // Mark hits as used - find all MCParticles that contribute
+ HashSet<MCParticle> mc_particleset = new HashSet<MCParticle>();
+ for (SimTrackerHit hit : hitset)
+ {
+ // Find parent MCParticle if not null
+ if (hit == null) continue;
+
+ // Add to used hits
+ _used_hits.add(hit);
+// System.out.println("adding used hit!");
+
+ mc_particleset.add(hit.getMCParticle());
+
+ }
+
+ // Loop over MCParticles and find largest contribution
+ int nhits_max = 0;
+ MCParticle majority_particle = null;
+ for (MCParticle particle : mc_particleset)
+ {
+ int nhits_part = 0;
+ for (SimTrackerHit hit : hitset) {
+ if (hit == null) continue;
+ if (particle == hit.getMCParticle())
+ {
+ nhits_part++;
+ }
+ }
+ if (nhits_part > nhits_max)
+ {
+ majority_particle = particle;
+ nhits_max = nhits_part;
+ }
+ }
+
+ // increment found counter
+ if (all_layer_particles.contains(majority_particle) &&
+ !found_particles.contains(majority_particle))
+ {
+ found_particles.add(majority_particle);
+ nfound++;
+ // Make histogram of found particle momenta
+ if (_make_histograms)
+ {
+ if (nhits_max >= 3) {
+ _aida.cloud1D("Found MCParticle momentum")
+ .fill(majority_particle.getMomentum().magnitude());
+ }
+ }
+ }
+
+ // Make histograms of passing track qualities
+ if (_make_histograms)
+ {
+ _aida.cloud1D("chisq_pass").fill(_fit.chisq());
+ _aida.cloud1D("dca_pass").fill(_fit.dca());
+ _aida.cloud1D("nhits_pass").fill(used_layers.size());
+ _aida.cloud1D("purity").fill((double)nhits_max/(double)nhits);
+ }
+
+ // Create StandaloneAxialBarrelTrack
+ StandaloneAxialBarrelTrack track = new StandaloneAxialBarrelTrack(b_field, _module_length, _fit, hitset);
+ tracklist.add(track);
+
+ } // if chisq not crap
+ } // if more than 4 layers used
+
+
+ } // if 3-point dca < 0.5
+
+
+
+ } // hit3 loop
+ } // hit2 loop
+ } // hit1 loop
+
+ } // layer3 loop
+ } // layer2 loop
+ } // layer1 loop
+
+ // Make histograms of # tracks found and per-event efficiency
+ if (_make_histograms)
+ {
+ _aida.cloud1D("ntracks_found").fill(nfound);
+ if (ntracks_all_layers > 0) {
+ _aida.cloud1D("efficiency").fill((double)nfound/(double)ntracks_all_layers);
+ }
+ }
+
+ // Put found tracks into event
+ event.put(EventHeader.TRACKS, tracklist, Track.class, 0);
+
+ // If requested, put unused hits into event
+ if (_output_hit_collection != "")
+ {
+ List<SimTrackerHit> unused_hits = _hits;
+ for (SimTrackerHit hit : _used_hits)
+ {
+ unused_hits.remove(hit);
+ }
+ event.put(_output_hit_collection, unused_hits, SimTrackerHit.class, 0, "TkrBarrHits");
+ }
+
+ }
+
+
+ /**
+ *
+ * Generates plot of track-finding efficiency as a function of momentum whenever event processing
+ * is suspended.
+ *
+ */
+
+ protected void suspend()
+ {
+
+ if (_make_histograms)
+ {
+ IHistogramFactory hfactory = _aida.histogramFactory();
+
+ if (!_aida.cloud1D("Found MCParticle momentum").isConverted())
+ {
+ _aida.cloud1D("Found MCParticle momentum").convert(50, 0.0, 50.0);
+ }
+ if (!_aida.cloud1D("5-hit MCParticle momentum").isConverted())
+ {
+ _aida.cloud1D("5-hit MCParticle momentum").convert(50, 0.0, 50.0);
+ }
+
+ hfactory.divide("efficiency vs. momentum",
+ _aida.cloud1D("Found MCParticle momentum").histogram(),
+ _aida.cloud1D("5-hit MCParticle momentum").histogram());
+
+ }
+ }
+
+ // Public methods to control AxialBarrelTrackFinder
+ //======================================
+ /**
+ * Turn on histograms
+ */
+ public void makeHistograms() {_make_histograms = true;}
+ /**
+ * Set name of SimTrackerHit collection to make tracks from (default = "TkrBarrHits")
+ *
+ * @param input_hit_collection Name of input SimTrackerHit collection
+ *
+ */
+ public void setInputHits(String input_hit_collection) {_input_hit_collection = input_hit_collection;}
+ /**
+ * Set name of SimTrackerHit collection for unused hits (default = "" == not written out)
+ *
+ * @param output_hit_collection Name of output collection for unused SimTrackerHits
+ *
+ */
+ public void setOutputHits(String output_hit_collection) {_output_hit_collection = output_hit_collection;}
+ /**
+ * Set length of modules in z (default = 100.0)
+ *
+ * @param module_length Module length (mm)
+ *
+ */
+
+
+ public void setModuleLength(double module_length) {_module_length = module_length;}
+ /**
+ * Set Isolation requirement for hits used in seeds (default = 0.5)
+ *
+ * @param seedhit_isolation Isolation of hits used in seeds (mm)
+ *
+ */
+ public void setSeedhitIsolation(double seedhit_isolation) {_seedhit_isolation = seedhit_isolation;}
+ /**
+ * Set maximum DCA from IP for 3-hit seeds (default = 2.0)
+ *
+ * @param seed_ip_dca Maximum DCA for 3-hit seeds (mm)
+ *
+ */
+ public void setSeedIpDca(double seed_ip_dca) {_seed_ip_dca = seed_ip_dca;}
+ /**
+ * Set maximum distance to hits for first pass of hit addition (default = 0.5)
+ *
+ * @param pass1_hit_dca Maximum distance to additional hits in first pass (mm)
+ *
+ */
+ public void setPass1HitDca(double pass1_hit_dca) {_pass1_hit_dca = pass1_hit_dca;}
+ /**
+ * Set minimum distance to next-nearest hit for first pass of hit addition (default = 1.0)
+ *
+ * @param pass1_hit_isolation Minimum distance to next-nearest hit in first pass (mm)
+ *
+ */
+ public void setPass1HitIsolation(double pass1_hit_isolation) {_pass1_hit_isolation = pass1_hit_isolation;}
+ /**
+ * Set maximum distance to hits for second pass of hit addition (default = 0.25)
+ *
+ * @param pass2_hit_dca Maximum distance to additional hits in second pass (mm)
+ *
+ */
+ public void setPass2HitDca(double pass2_hit_dca) {_pass2_hit_dca = pass2_hit_dca;}
+ /**
+ * Set maximum chisq/nhits to create tracks (default = 10.0)
+ *
+ * @param chisq_dof Maximum chisq/nhits for tracks
+ *
+ */
+ public void setChisqDof(double chisq_dof) {_chisq_dof = chisq_dof;}
+
+ // Protected methods
+ //==================
+
+ /**
+ *
+ * Given a SimTrackerHit and a list of same, reports the r-phi distance to the nearest
+ * distinct hit in the list to the given hit. Used extensively in process().
+ *
+ * @param hit The SimTrackerHit of interest
+ * @param hits A list of SimTrackerHits to check for the nearest distinct hit
+ *
+ */
+ protected double rphiDistToNearest(SimTrackerHit hit, List<SimTrackerHit> hits)
+ {
+ double nearest_dist = 100000.0; // 1Km
+ for (SimTrackerHit hit2 : hits)
+ {
+ if (hit2 != hit)
+ {
+ double dx = hit2.getPoint()[0] - hit.getPoint()[0];
+ double dy = hit2.getPoint()[1] - hit.getPoint()[1];
+ nearest_dist = Math.min(Math.sqrt(dx*dx+dy*dy),nearest_dist);
+ }
+ }
+ return nearest_dist;
+ }
+
+}
\ No newline at end of file
lcsim/src/org/lcsim/contrib/tracking
diff -N StandaloneAxialBarrelTrack.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ StandaloneAxialBarrelTrack.java 22 Sep 2005 19:37:53 -0000 1.1
@@ -0,0 +1,403 @@
+/*
+ * StandaloneAxialBarrelTrack.java
+ *
+ * Created on August 29, 2005, 4:41 PM
+ *
+ * To change this template, choose Tools | Options and locate the template under
+ * the Source Creation and Management node. Right-click the template and choose
+ * Open. You can then make changes to the template in the Source Editor.
+ */
+
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.fit.circle.CircleFit;
+
+import java.util.*;
+
+/**
+ * Implements the Track interface for tracks found using only outer tracker information by TrackFinder
+ *
+ * @author tknelson
+ * @see AxialBarrelTrackFinder
+ *
+ */
+public class StandaloneAxialBarrelTrack implements Track
+{
+
+ // Data members
+ boolean _fit_success = false;
+ double _chi2 = 0.0;
+ int _ndof = 0;
+ double _dedx = 0.0;
+ double _dedx_error = 0.0;
+ double[][] _error_matrix = new double[5][5];
+ double[] _reference_point = new double[3];
+ int[] _hit_numbers = new int[0];
+ List<TrackerHit>_tracker_hits = new ArrayList<TrackerHit>();
+ List<SimTrackerHit>_sim_tracker_hits = new ArrayList<SimTrackerHit>();
+ double[] _track_parameters = new double[5];
+ double _b_field = 0.0;
+ double _module_length = 100.0;
+ List<Track> _tracks = new ArrayList<Track>();
+ boolean _isreferencepoint_pca = false;
+
+ // Static
+ static int _type = 2; // ???
+
+
+ // Constructors
+ //=============
+
+ /**
+ * Create a default instance of StandaloneAxialBarrelTrack
+ */
+ public StandaloneAxialBarrelTrack()
+ {
+ }
+
+ /**
+ * Create a StandaloneAxialBarrelTrack using information from TrackFinder
+ *
+ * @param b_field Magnetic field: required to generate momentum from curvature
+ * @param module_length Length of readout modules: required to estimate tan(lambda)
+ * @param fit The CircleFit for this track
+ * @param hits The list of SimTrackerHits used for this track
+ *
+ * @see org.lcsim.fit.circle.CircleFit
+ */
+ public StandaloneAxialBarrelTrack(double b_field, double module_length, CircleFit fit, List<SimTrackerHit> hits)
+ {
+ _b_field = b_field;
+ _module_length = module_length;
+
+ _fit_success = true;
+ _chi2 = fit.chisq();
+ _ndof = hits.size();
+
+ _reference_point[0] = fit.xref();
+ _reference_point[1] = fit.yref();
+ _reference_point[2] = 0.0;
+
+ _sim_tracker_hits = hits;
+
+ _track_parameters[0] = fit.dca();
+ _track_parameters[1] = fit.phi()-Math.PI; // who ordered this?
+ _track_parameters[2] = 10.0*fit.curvature(); // convert 1/mm to 1/cm
+ _track_parameters[3] = 0.0;
+ _track_parameters[4] = getTanLambda();
+
+ _isreferencepoint_pca = true;
+
+ _tracks.add(this);
+
+ }
+
+ // Methods needed to further initialize a StandaloneAxialBarrelTrack
+ //=======================================================
+
+ /**
+ * Set the dEdx for this track (default is 0.0)
+ */
+ public void setdEdx(double dedx)
+ {
+ _dedx = dedx;
+ }
+
+ /**
+ * Set the dEdx error for this track (default is 0.0)
+ */
+ public void setdEdxError(double dedx_error)
+ {
+ _dedx_error = dedx_error;
+ }
+
+ /**
+ * Set the tracks associated with this track
+ */
+ public void setTracks(List<Track> tracks)
+ {
+ _tracks = tracks;
+ }
+
+ // Implementation of Track interface
+ /**
+ * The charge of the particle creating this track
+ */
+ public int getCharge()
+ {
+ return (int)Math.signum(_track_parameters[2]);
+ }
+
+ /**
+ * The reference point for the track parameters
+ */
+ public double[] getReferencePoint()
+ {
+ return _reference_point;
+ }
+
+ /**
+ * The x coordinate of the reference point for the track parameters
+ */
+ public double getReferencePointX()
+ {
+ return _reference_point[0];
+ }
+
+ /**
+ * The y coordinate of the reference point for the track parameters
+ */
+ public double getReferencePointY()
+ {
+ return _reference_point[1];
+ }
+
+ /**
+ * The z coordinate of the reference point for the track parameters
+ */
+ public double getReferencePointZ()
+ {
+ return _reference_point[2];
+ }
+
+ /**
+ * Is this reference point the point of closest approach used to define the track parameters?
+ */
+ public boolean isReferencePointPCA()
+ {
+ return _isreferencepoint_pca;
+ }
+
+ /**
+ * The momentum of the particle corresponding to this track: depends upon b_field being set
+ * appropriately when constructing the StandaloneAxialBarrelTrack
+ */
+ public double[] getMomentum()
+ {
+ double[] momentum = new double[3];
+ double pt = (0.00029979*_b_field)/_track_parameters[2];
+ momentum[0] = pt*Math.cos(_track_parameters[1]);
+ momentum[1] = pt*Math.sin(_track_parameters[1]);
+ momentum[2] = pt*_track_parameters[4];
+ return momentum;
+ }
+
+ /**
+ * The x-component of the momentum of the particle corresponding to this track: depends upon
+ * b_field being appropriately when constructing the StandaloneAxialBarrelTrack
+ */
+ public double getPX()
+ {
+ return getMomentum()[0];
+ }
+
+ /**
+ * The y-component of the momentum of the particle corresponding to this track: depends upon
+ * b_field being appropriately when constructing the StandaloneAxialBarrelTrack
+ */
+ public double getPY()
+ {
+ return getMomentum()[1];
+ }
+
+ /**
+ * The z-component of the momentum of the particle corresponding to this track: depends upon
+ * b_field being appropriately when constructing the StandaloneAxialBarrelTrack
+ */
+ public double getPZ()
+ {
+ return getMomentum()[2];
+ }
+
+ /**
+ * Was track fit successful?
+ */
+ public boolean fitSuccess()
+ {
+ return _fit_success;
+ }
+
+ /**
+ * Retrieve the i-th track parameter
+ *
+ * @param i Index of track parameter to return
+ *
+ * @see org.lcsim.event.EventHeader.Track
+ *
+ */
+ public double getTrackParameter(int i)
+ {
+ return _track_parameters[i];
+ }
+
+ /**
+ * Retreive the vector of track parameters
+ *
+ * @see org.lcsim.event.EventHeader.Track
+ *
+ */
+ public double[] getTrackParameters()
+ {
+ return _track_parameters;
+ }
+
+ /**
+ * Retreive the (i,j) element of the error matrix of the track fit
+ * <p>
+ * <b>N.B. Currently returns 0.0</b>
+ */
+ public double getErrorMatrixElement(int i, int j)
+ {
+ return _error_matrix[i][j];
+ }
+
+ /**
+ * Retreive the error matrix of the track fit
+ * <p>
+ * <b>N.B. Currently returns an empty matrix</b>
+ */
+ public double[][] getErrorMatrix()
+ {
+ return _error_matrix;
+ }
+
+ /**
+ * Retreive the chisquared of the track fit
+ */
+ public double getChi2()
+ {
+ return _chi2;
+ }
+
+ /**
+ * Retrieve the number of degrees of freedom of the fit = nhits
+ */
+ public int getNDF()
+ {
+ return _ndof;
+ }
+
+ /**
+ * Retrieve the dEdx of the track
+ * <p>
+ * <b>N.B. Currently zero unless set explicitly</b>
+ */
+ public double getdEdx()
+ {
+ return _dedx;
+ }
+
+ /**
+ * Retrieve the dEdx error of the track: zero unless set explicitly
+ *<p>
+ *<b>N.B. Currently zero unless set explicitly</b>
+ */
+ public double getdEdxError()
+ {
+ return _dedx_error;
+ }
+
+ /**
+ * Get radius of innermost hit on track
+ */
+ public double getRadiusOfInnermostHit()
+ {
+ double min_radius = 1000000.0; // 1km
+ for (SimTrackerHit hit : _sim_tracker_hits){
+ min_radius = Math.min(Math.sqrt(hit.getPoint()[0]*hit.getPoint()[0] +
+ hit.getPoint()[1]*hit.getPoint()[1]), min_radius);
+ }
+ return min_radius;
+ }
+
+ /**
+ * Get subdetector hit numbers
+ *
+ * <b>N.B. Currently returns an empty vector</b>
+ */
+ public int[] getSubdetectorHitNumbers()
+ {
+ return _hit_numbers;
+ }
+
+ /**
+ * Get tracks associated with this track.
+ * <p>
+ * <b>Currently returns only this unless set explicitly</b>
+ */
+ public List<Track> getTracks()
+ {
+ return _tracks;
+ }
+
+ /**
+ * Get TrackerHits associated with this track.
+ * <p>
+ * <b>N.B. Since SimTrackerHits are currently used, this returns a blank list</b>
+ */
+ public List<TrackerHit> getTrackerHits()
+ {
+ return _tracker_hits;
+ }
+
+ /**
+ * Get type of track (type = 2 for StandaloneAxialBarrelTrack)
+ */
+ public int getType()
+ {
+ return _type;
+ }
+
+ // Other public methods
+ /**
+ * Get List of SimTrackerHits for this track
+ */
+ public List<SimTrackerHit> getSimTrackerHits()
+ {
+ return _sim_tracker_hits;
+ }
+
+ // Private methods
+ /**
+ * Calculate tan(lambda) using outermost hit radius and module length
+ */
+ private double getTanLambda() {
+ SimTrackerHit hit = getOutermostHit();
+ double x_outer = hit.getPoint()[0];
+ double y_outer = hit.getPoint()[1];
+ double z_outer_mod = hit.getPoint()[2]/_module_length;
+ double z_outer = Math.signum(z_outer_mod)*(Math.floor(Math.abs(z_outer_mod))+0.5)*_module_length;
+
+ double x_0 = _track_parameters[0]*Math.cos(_track_parameters[1]+Math.PI/2.0);
+ double y_0 = _track_parameters[0]*Math.sin(_track_parameters[1]+Math.PI/2.0);
+ double z_0 = _track_parameters[3];
+
+ double dist = Math.sqrt(Math.pow(x_outer-x_0, 2)+Math.pow(y_outer-y_0, 2));
+ double phi_path = 2.0*Math.asin(dist*(_track_parameters[2]/10.0)/2.0);
+ double pathlength = phi_path/(_track_parameters[2]/10.0);
+
+ return (z_outer-z_0)/pathlength;
+ }
+
+ /**
+ * Return outermost hit on track
+ */
+ private SimTrackerHit getOutermostHit()
+ {
+ double max_radius = 0.0;
+ SimTrackerHit outermost_hit = null;
+ for (SimTrackerHit hit : _sim_tracker_hits){
+ if (hit == null) continue;
+ double radius = Math.sqrt(hit.getPoint()[0]*hit.getPoint()[0] +
+ hit.getPoint()[1]*hit.getPoint()[1]);
+ if ( radius > max_radius)
+ {
+ max_radius = radius;
+ outermost_hit = hit;
+ }
+ }
+ return outermost_hit;
+ }
+
+}
\ No newline at end of file
