lcsim/src/org/lcsim/contrib/tracking
diff -N ZSeg.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ZSeg.java 18 Jul 2007 22:12:26 -0000 1.1
@@ -0,0 +1,128 @@
+/*
+ * ZSeg.java
+ *
+ * Created on July 9, 2007, 12:20 PM
+ *
+ * To change this template, choose Tools | Template Manager
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.tracking;
+
+import org.lcsim.event.SimTrackerHit;
+
+/**
+ * Segmentation of z-axis, segment length determined by user.
+ * Loops through list of hits to check all possible 3 hit combinations
+ * for consistency, using either extrapolation or interpolation.
+ * Takes hits from 2 layers and projects onto 3rd layer.
+ * Looks to see if 3rd hit falls within projected range.
+ * If inconsistent, returns false. Assumes only information known
+ * is which module the hit was on.
+ *
+ * @author Lori Stevens
+ */
+public class ZSeg {
+ //segment length, in mm
+ static int segs=300;
+
+ /** Creates a new instance of ZSeg */
+
+ public ZSeg() {
+ }
+
+ public static boolean zCheck(SimTrackerHit [] hits){
+ boolean zCheck = true;
+ //looping through hits to check for all possible 3 hit combinations
+ //as soon as one false is returned, false is returned for entire array
+ for (int i=0; i<hits.length; i++){
+ if (hits[i]==null) continue;
+
+ for (int j=0; j<hits.length; j++){
+ if (hits[j]==null || hits[j]==hits[i]) continue;
+
+ for (int k=0; k<hits.length; k++){
+ if(hits[k]==null || hits[k]==hits[i] || hits[k]==hits[j]) continue;
+
+ if (!zChecker(hits[i], hits[j], hits[k]))
+ return false;
+ }
+ }
+ }
+ return zCheck;
+ }
+
+ //checks whether or not 3rd hit falls within projected range
+ //includes check for whether to use interpolation or extrapolation
+ public static boolean zChecker(SimTrackerHit ht1, SimTrackerHit ht2, SimTrackerHit ht3){
+ boolean hitPass = false;
+
+ if (ht1!=null && ht2!=null && ht3 !=null){
+ double r1 = Math.sqrt(ht1.getPoint()[0]*ht1.getPoint()[0]+ht1.getPoint()[1]*ht1.getPoint()[1]);
+ double r2 = Math.sqrt(ht2.getPoint()[0]*ht2.getPoint()[0]+ht2.getPoint()[1]*ht2.getPoint()[1]);
+ double r3 = Math.sqrt(ht3.getPoint()[0]*ht3.getPoint()[0]+ht3.getPoint()[1]*ht3.getPoint()[1]);
+
+ double[] mod1 = moduleInfo(ht1.getPoint()[2], r1);
+ double[] mod2 = moduleInfo(ht2.getPoint()[2], r2);
+ double[] mod3 = moduleInfo(ht3.getPoint()[2], r3);
+
+ double z1 = zFinder(mod1[1], r1, mod2[0], r2, r3);
+ double z2 = zFinder(mod1[0], r1, mod2[1], r2, r3);
+
+ //extrapolation: used when hit is in layer on 1 side or other of 2 hits
+ //if any part of module lies between projected lines
+ //then hit passes extrapolation test
+ if ((r3 < r2 && r3 < r1) || (r3 > r2 && r3 > r1)){
+ //checking if minz lies within projected range
+ //only minz OR maxz needs to fall within range since
+ //we only have module info and not exact hit info
+ if ((mod3[0] >= z1 && mod3[0] <= z2) || (mod3[0] <= z1 && mod3[0] >= z2))
+ hitPass = true;
+ //checking if maxz lies within projected range
+ if ((mod3[1] >= z1 && mod3[1] <= z2) || (mod3[1] <= z1 && mod3[1] >= z2))
+ hitPass = true;
+ }
+ //interpolation: used when third hit is in layer between 2 other hits;
+ //if any part of module lies between projected lines
+ //then hit passes interpolation test
+ else if ((r3 < r1 && r3 > r2) || (r3 > r1 && r3 < r2)){
+ //checking if minz or maxz lies within projected range
+ if((mod3[0] >= z1 && mod3[0] <= z2) || (mod3[1] >= z1 && mod3[1] <= z2)){
+ hitPass = true;
+ }
+ }
+ }
+ return hitPass;
+ }
+
+ //only have information about which module the hit is on, not exact coordinate info
+ //method calculates min and max z of modules
+ public static double[] moduleInfo(double z, double r) {
+ int minz=0;
+ int maxz=0;
+ //first positive z segment has zmin=0
+ if (z>=0) {
+ minz = ((int)(z/segs))*segs;
+ maxz = minz+segs;
+ }
+ //neg numbers round in positive direction. ex. -2.2 becomes -2
+ if (z<0) {
+ maxz = ((int)(z/segs))*segs;
+ minz = maxz-segs;
+ if (maxz==minz){
+ maxz += segs;
+ }
+ }
+ double [ ] modInfo = {minz, maxz, r};
+
+ return modInfo;
+ }
+ //method used to draw projected lines onto layer at radius3
+ public static double zFinder(double z1, double r1, double z2, double r2, double r3) {
+ double m = (r1-r2)/(z1-z2);
+ double b = r1-m*z1;
+ double z3 = (r3-b)/m;
+ if (z1==z2) {z3=z1;}
+
+ return z3;
+ }
+}
lcsim/src/org/lcsim/contrib/tracking
diff -N AxialBarrelTrackFinderZ.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AxialBarrelTrackFinderZ.java 18 Jul 2007 22:12:26 -0000 1.1
@@ -0,0 +1,804 @@
+/*
+ * AxialBarrelTrackFinderZ.java
+ *
+ * Created on August 18, 2005, 4:17 PM
+ *
+ */
+
+package org.lcsim.contrib.tracking;
+
+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.detector.converter.compact.DeDetector;
+import org.lcsim.geometry.subdetector.MultiLayerTracker;
+import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.geometry.util.TrackerIDDecoder;
+import org.lcsim.util.aida.AIDA;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.aida.*;
+
+import java.util.*;
+import java.lang.Math;
+import org.lcsim.contrib.tracking.StandaloneAxialBarrelTrack1;
+
+/**
+ * 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 AxialBarrelTrackFinderZ 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; //0; //0.5;
+ private double _seed_ip_dca = 100.0; //900000; //100.0;
+ private double _pass1_hit_dca = 0.5; //900000; //0.5;
+ private double _pass1_hit_isolation = 1.0; //0; //1.0;
+ private double _pass2_hit_dca = 0.25; //900000; //0.25;
+ private double _chisq_dof = 10; //1.0E20; //10.0;
+ private boolean _make_histograms = false;
+ private int nAssocMin = 4;
+ private double _min_seed_pt = 0.5; //*Tyler
+ private double _max_phi_sep = (Math.PI)/2; //*Tyler
+
+
+ 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;
+
+ //zCheckOn turns z segmentation on and off //**************************
+ boolean zCheckOn = false;
+
+ /**
+ * Creates a new instance of AxialBarrelTrackFinderZ
+ */
+ public AxialBarrelTrackFinderZ()
+ {
+ }
+
+ /**
+ * 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 = Collections.synchronizedList(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<StandaloneAxialBarrelTrack1> tracklist = new ArrayList<StandaloneAxialBarrelTrack1>();
+ ArrayList<StandaloneAxialBarrelTrack1> faketracks = new ArrayList<StandaloneAxialBarrelTrack1>();
+ // 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[nlayers];
+ double[] y = new double[nlayers];
+ double[] weight = new double[nlayers];
+
+ // 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
+
+ //hit list for z segmentation *************************
+ SimTrackerHit [ ] hitList = new SimTrackerHit [5];
+ //adding first hit to list for z segmentation
+ hitList[0] = hit1;
+
+ for (SimTrackerHit hit2 : hits2)
+ {
+ if (_hit_separation.get(hit2) < _seedhit_isolation) continue; // require well-separated hits
+ if (_used_hits.contains(hit1)) break;
+ if (_used_hits.contains(hit2)) continue; // skip used hits
+
+ //if z segmentation is used, there is no **************************
+ //requirement for hit to be on same side in z
+ if (!zCheckOn) {
+ if ((hit2.getPoint()[2] > 0.0) != zpos) continue; //require same side in z
+ }
+ //adding 2nd hit to list for z segmentation
+ hitList[1] = hit2;
+
+ //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(hit2)) break;
+ if (_used_hits.contains(hit3)) continue; // skip used hits
+
+ //if z segmentation is used, there is no **************************
+ //requirement for hit to be on same side in z
+ if (!zCheckOn) {
+ if ((hit3.getPoint()[2] > 0.0) != zpos) continue; // require same side in z
+ }
+ if (zCheckOn) {
+ //adding 3rd hit to list for z segmentation
+ hitList[2] = hit3;
+ //z segmentation check n 3 hit seed
+ if (!ZSeg.zCheck(hitList))
+ continue;
+ }
+
+ //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;
+ Vector<SimTrackerHit> hitset = new Vector<SimTrackerHit>();
+ hitset.setSize(nlayers);
+ //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
+ //Require the fitted transverse momentum of the seed is of appropriate magnitude *Tyler
+ double pc = (0.00029979*b_field)/(_fit.curvature());
+ if (Math.abs(_fit.dca()) <_seed_ip_dca && Math.abs(pc) > _min_seed_pt )
+ {
+ // 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 z segmentation is used, there is no ********************************
+ //requirement for hit to be on same side in z
+ if (!zCheckOn) {
+ if ((hit.getPoint()[2] > 0.0) != zpos) continue; // require same side in z
+ }
+ if (zCheckOn){
+ //adding 4th hit to list for z segmentation check
+ hitList[3] = hit;
+ //z segmentation check on 4 hits in list; erases 4th hit if check fails
+ if(!ZSeg.zCheck(hitList)){
+ hitList[3]=null;
+ continue;
+ }
+ }
+
+ //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() < nlayers )
+ {
+ // 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 z segmentation is used, there is no ************************************
+ //requirement for hit to be on same side in z
+ if (!zCheckOn){
+ if ((hit.getPoint()[2] > 0.0) != zpos) continue; // require same side in z
+ }
+ if (zCheckOn){
+ //adding 5th hit to z segmentation list
+ hitList[4] = hit;
+ //check on all 5 hits; deletes hit from list if check fails
+ if(!ZSeg.zCheck(hitList)){
+ hitList[4]=null;
+ continue;
+ }
+ }
+
+ //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();
+ double[] phivalues2 = new double[5];
+ int k2 = 0;
+
+ //Create Array of Phi-hit values *Tyler
+ for (SimTrackerHit hit : hitset)
+ {
+ if(hit == null) continue;
+ k2++;
+ double[] hitpoint2 = hit.getPoint();
+ double hitphi2 = Math.atan2(hitpoint2[1],hitpoint2[0]);
+ phivalues2[k2-1] = hitphi2;
+ //rvalues[k-1] = Math.hypot(hitpoint1[0],hitpoint1[1]);
+ }
+
+ //If hit is null on layer, use the hitphi of the first hit *Tyler
+ if (phivalues2[4] == 0)
+ {
+ phivalues2[4] = phivalues2[0];
+ }
+
+ int b2=0;
+ for(int ps2=0; ps2<5; ps2++)
+ {
+
+ //Test hits to be sure all hits have a reasonable difference in phi *Tyler
+ for(int ts2=0; ts2<5; ts2++)
+ {
+ double phidifference2 = Math.abs(phivalues2[ps2]-phivalues2[ts2]);
+ if( (phidifference2 > _max_phi_sep) && (phidifference2 < 2*(Math.PI )-_max_phi_sep))
+ {
+ b2 = 1;
+ }
+ }
+ }
+
+ //If All hit phi values lie within specified phi *Tyler
+ if (b2 !=1)
+ {
+
+ // 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>();
+ int n = 0;
+ for (SimTrackerHit hit : hitset)
+ {
+ // Find parent MCParticle if not null
+ if (hit == null) continue;
+ // Add to used hits
+ n++;
+ _used_hits.add(hit);
+ //System.out.println("adding used hit!");
+
+ mc_particleset.add(hit.getMCParticle());
+
+ }
+
+ //System.out.println("Hitset : \n" + hitset + " \n Used Hits: \n" + _used_hits + "\n");
+ //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) && nhits_max >= nAssocMin)
+ {
+ 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());
+ }
+ }
+ }
+
+ double purity = (double)nhits_max/(double)nhits;
+
+ // 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 StandaloneAxialBarrelTrack1
+ StandaloneAxialBarrelTrack1 track = new StandaloneAxialBarrelTrack1(b_field, _module_length, _fit, hitset, majority_particle, purity);
+ tracklist.add(track);
+
+ } // 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
+
+ // 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;
+ int i = 0;
+ for (SimTrackerHit hit : _used_hits)
+ {
+ i++;
+ 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 AxialBarrelTrackFinderZ
+ //======================================
+
+
+ /**
+ * 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;}
+ /**
+ * Set minimum momentum required for the fitted track seed to be considered of appropriate magnitude
+ *
+ * @param min_seed_pt Minimum transverse momentum for track seeds
+ *
+ */
+ public void setMinSeedPt(double min_seed_pt)
+ {_min_seed_pt = min_seed_pt;} //*Tyler
+ /**
+ * Set maximum hit-phi seperation allowed for each hit in fitted track
+ *
+ * @param min_seed_pt Maximum hit-phi seperation
+ *
+ */
+ public void setMaxHitPhi(double max_phi_sep)
+ {_max_phi_sep = max_phi_sep;} //*Tyler
+
+
+
+
+
+
+
+ // 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;
+ }
+
+
+}