6 added files
lcsim/sandbox/LoriStevens
diff -N AxialBarrelTrackFinder1.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AxialBarrelTrackFinder1.java 27 Jul 2007 18:07:47 -0000 1.1
@@ -0,0 +1,802 @@
+/*
+ * AxialBarrelTrackFinder1.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.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 AxialBarrelTrackFinder1 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 AxialBarrelTrackFinder1
+ */
+ public AxialBarrelTrackFinder1()
+ {
+ }
+
+ /**
+ * 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 AxialBarrelTrackFinder1
+ //======================================
+
+
+ /**
+ * 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;
+ }
+
+
+}
lcsim/sandbox/LoriStevens
diff -N HelicalTrackFitter.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelicalTrackFitter.java 27 Jul 2007 18:07:47 -0000 1.1
@@ -0,0 +1,316 @@
+/*
+ * HelicalTrackFitter.java
+ *
+ * Created on March 25, 2006, 6:11 PM
+ *
+ * $Id: HelicalTrackFitter.java,v 1.1 2007/07/27 18:07:47 lstevens Exp $
+ */
+
+import hep.physics.matrix.SymmetricMatrix;
+import static java.lang.Math.atan2;
+import static java.lang.Math.abs;
+import static java.lang.Math.PI;
+import static java.lang.Math.cos;
+import static java.lang.Math.sin;
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.fit.circle.CircleFitter;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.line.SlopeInterceptLineFit;
+import org.lcsim.fit.line.SlopeInterceptLineFitter;
+import java.util.ArrayList;
+import java.util.Arrays;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.base.BaseTrackerHit;
+import java.util.Collections;
+import java.util.*;
+/**
+ * Fit a helix to a set of space points. First, a circle is fit to the x-y coordinates.
+ * A straight-line fit is then performed on the s-z coordinates.
+ *
+ * The r-phi and z coordinate measurements are assumed to be uncorrelated. A block
+ * diagonal covariance matrix is formed from the results of the circle and s-z fits,
+ * ignoring any correlations between these fits.
+ *
+ * The resulting track parameters follow the "L3 Convention" adopted by org.lcsim.
+ *
+ * Modified August 2006 by Richard Partridge
+ * @author Norman Graf
+ */
+
+public class HelicalTrackFitter //implements Comparator<TrackerHit>
+{
+ private CircleFitter _cfitter = new CircleFitter();
+ private SlopeInterceptLineFitter _lfitter = new SlopeInterceptLineFitter();
+
+ private HelicalTrackFit _fit;
+
+ private double[] _circleParameters = new double[3];
+
+ private ArrayList<TrackerHit> barrelHits = new ArrayList<TrackerHit>();
+ /**
+ * Creates a new instance of HelicalTrackFitter.
+ */
+ public HelicalTrackFitter()
+ {
+ }
+ /**
+ * Fit a helix to a set of space points. The space points should be ordered in
+ * increasing track length. The code will properly handle loopers provided
+ * that the change in direction between successive points is less than pi.
+ * @param x Array of x coordinates
+ * @param y Array of y coordinates
+ * @param z Array of z coordinates
+ * @param drphi Error in r-phi hit position
+ * @param dz Error in z coordinate. A negative value indicates that this point is to
+ * be excluded in the s-z fit.
+ * @param np Number of points
+ * @return True for successful fit
+ */
+
+ public boolean fit(List<TrackerHit> hits)
+ {
+ int np = hits.size();
+ int[] type = new int[np];
+ double[] wrphi = new double[np];
+ double[] wz = new double[np];
+ int nz=0;
+
+ order(hits);
+
+ for(int i=0; i<np; ++i){
+ type[i] = hits.get(i).getType();
+ //(1/(dxdx + dydy))
+ wrphi[i] = 1/(hits.get(i).getCovMatrix()[1]+hits.get(i).getCovMatrix()[3]);
+ //(1/dzdz)
+ wz[nz] = 1/hits.get(i).getCovMatrix()[6];
+ nz++;
+ }
+
+// return fitting(wrphi, wz, x, y, z, np, type);
+ return fitting(wrphi, wz, hits, type);
+
+ }
+
+ //overloaded method, fit for Cartesian coordinate arrays
+ /* public boolean fit(double[] x, double[] y, double[] z, double[] drphi, double[] dz, int np)
+
+ {
+ needSwap(x, y, z, drphi, dz);
+
+ // fit a circle in the x-y plane
+ double[] wrphi = new double[np];
+ double[] wz = new double[np];
+ int[] type = new int[np];
+ int nz=0;
+ for(int i=0; i<np; ++i)
+ {
+ if(dz[i]>0) type[i]=0;
+ if(dz[i]<0) type[i]=1;
+ wrphi[i] = 1/(drphi[i]*drphi[i]);
+ wz[nz] = 1/(dz[i]*dz[i]);
+ nz++;
+ }
+
+ return fitting(wrphi, wz, x, y, z, np, type);
+ }
+*/
+// public boolean fitting(double[] wrphi, double[] wz, double[] x, double[] y, double[] z, int np, int[] type)
+ public boolean fitting(double[] wrphi, double[] wz, List<TrackerHit> hits, int[] type)
+ {
+ double[] x = new double[hits.size()];
+ double[] y = new double[hits.size()];
+ double[] z = new double[hits.size()];
+ for(int i=0; i<hits.size(); i++)
+ {
+ x[i] = hits.get(i).getPosition()[0];
+ y[i] = hits.get(i).getPosition()[1];
+ z[i] = hits.get(i).getPosition()[2];
+ }
+ boolean success = _cfitter.fit(x, y, wrphi, hits.size());
+ if(!success) return false;
+
+ CircleFit cfit = _cfitter.getfit();
+ double radius = 1./cfit.curvature();
+ double phi0 = cfit.phi();
+ double dca = -cfit.dca();
+ double xc = (radius-dca)*sin(phi0);
+ double yc = -(radius-dca)*cos(phi0);
+ _circleParameters[0] = xc;
+ _circleParameters[1] = yc;
+ _circleParameters[2] = radius;
+ // fit a line in the s-z plane
+ // assumes points are in increasing order
+ double[] s = new double[hits.size()];
+ double[] sfit = new double[hits.size()];
+ double[] zfit = new double[hits.size()];
+ int nz = 0;
+ for(int i=0; i<hits.size(); i++)
+ {
+ if (i==0) {
+ s[0] = s(radius, xc, yc, x[0], y[0], 0., 0.);
+ }
+ else {
+ s[i] = s(radius, xc, yc, x[i], y[i], x[i-1], y[i-1]) + s[i-1];
+ }
+ //type 0 = 3D hit
+ if (type[i]==0)
+ {
+ sfit[nz] = s[i];
+ zfit[nz] = z[i];
+ nz++;
+ }
+ //z seg; type 1 = 2D hit; this still needs work
+ if (type[i]==1)
+ {
+ double[] pos = {x[i], y[i], z[i]};
+ barrelHits.add(new BaseTrackerHit(pos, null, 0, 0, 0));
+ if(barrelHits.size() >= 3){
+ // if(!ZSegTrackerHit.zCheck(barrelHits))
+ // return false;
+ }
+ }
+ }
+ success = _lfitter.fit(sfit, zfit, wz, nz);
+ if(!success) return false;
+ SlopeInterceptLineFit lfit = _lfitter.getFit();
+
+ double[] pars = new double[5];
+
+ SymmetricMatrix cov = new SymmetricMatrix(5);
+ double[] chisq = new double[2];
+ int[] ndf = new int[2];
+
+ chisq[0] = cfit.chisq();
+ chisq[1] = lfit.chisquared();
+
+ ndf[1] = lfit.ndf();
+ ndf[0] = hits.size() - 3;
+
+ // d0, xy impact parameter. Note: - sign due to opposite sign convention in CircleFitter
+ pars[0] = -cfit.dca();
+ // phi0
+ pars[1] = cfit.phi();
+ // omega signed curvature
+ pars[2] = cfit.curvature();
+ // z0
+ pars[3] = lfit.intercept();
+ // tan lambda
+ pars[4] = lfit.slope();
+
+ // fill in covariance matrix
+ cov.setElement(0,0, cfit.cov()[0]);
+ cov.setElement(0,1,-cfit.cov()[1]); // Need - sign due to sign convention in CircleFitter
+ cov.setElement(1,1, cfit.cov()[2]);
+ cov.setElement(0,2,-cfit.cov()[3]); // Need - sign due to sign convention in CircleFitter
+ cov.setElement(1,2, cfit.cov()[4]);
+ cov.setElement(2,2, cfit.cov()[5]);
+ //cov[0][3] = 0.;
+ //cov[1][3] = 0.;
+ //cov[2][3] = 0.;
+ cov.setElement(3,3, lfit.interceptUncertainty());
+ //cov[0][4] = 0.;
+ //cov[1][4] = 0.;
+ //cov[2][4] = 0.;
+ cov.setElement(3,4, lfit.covariance());
+ cov.setElement(4,4, lfit.slopeUncertainty());
+ _fit = new HelicalTrackFit(pars, cov, chisq, ndf);
+ _fit.setCircleParameters(_circleParameters);
+ return true;
+ }
+ /**
+ * Get the results of the most recent fit.
+ * @return HelicalTrackFit corresponding to most recent fit
+ */
+ public HelicalTrackFit getFit()
+ {
+ return _fit;
+ }
+ /**
+ * Get the circle paramaters (xc, yc, R) from the most recent fit.
+ * @return Circle parameters (xc, yc, R)
+ */
+ public double[] getCircleParameters()
+ {
+ return _circleParameters;
+ }
+
+ double s(double radius, double xcenter, double ycenter, double x1, double y1, double x2, double y2)
+ {
+ double phi1 = atan2( (y1-ycenter), (x1-xcenter) );
+ double phi2 = atan2( (y2-ycenter), (x2-xcenter) );
+ double dphi = abs(phi1-phi2);
+ if(dphi>PI) dphi = 2.*PI-dphi;
+ return abs(radius)*dphi;
+ }
+ public double[] sortDist(double[] x, double[] y, double[] z)
+ {
+ double[] dist = new double[x.length];
+ for(int i=0; i<dist.length; i++)
+ {
+ dist[i]=Math.sqrt(x[i]*x[i]+y[i]*y[i]+z[i]*z[i]);
+ }
+ Arrays.sort(dist);
+ return dist;
+ }
+ static void swap(double[] data, int first, int second)
+ {
+ double temp = data[first];
+ data[first] = data[second];
+ data[second] = temp;
+ }
+ public void needSwap(double[] x, double[] y, double[] z, double[] drphi, double[] dz)
+ {
+ double[] dist = sortDist(x, y, z);
+ double[] ztest = z.clone();
+ Arrays.sort(ztest);
+
+ for(int i=0; i<dist.length; i++){
+ //looking for index of hit with min dist from origin
+ if(dist[0]!=Math.sqrt(x[i]*x[i]+y[i]*y[i]+z[i]*z[i])) continue;
+ double zmin = z[i];
+ //ordering points according to z coordinate
+ for(int j=1; j<dist.length; j++){
+ for(int k=(dist.length-1); k>=j; k--){
+ //if z closest to origin is a minimum, order from least to greatest
+ //if z closest to origin is a maximum, order from greatest to least
+ if((zmin==ztest[0] && z[k-1] > z[k]) || (zmin==ztest[dist.length-1] && z[k-1] < z[k])){
+ swap(x, k, k-1);
+ swap(y, k, k-1);
+ swap(z, k, k-1);
+ swap(drphi, k, k-1);
+ swap(dz, k, k-1);
+ }
+ }
+ }
+ break;
+ }
+ }
+
+ public void order(List<TrackerHit> hits)
+ {
+ Comparator<TrackerHit> z = new Comp();
+ Collections.sort(hits, z);
+ double[] hit1 = hits.get(0).getPosition();
+ double[] hit2 = hits.get(hits.size()-1).getPosition();
+ double dist1 = Math.sqrt(hit1[0]*hit1[0]+hit1[1]*hit1[1]+hit1[2]*hit1[2]);
+ double dist2 = Math.sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]+hit2[2]*hit2[2]);
+ if(dist1 > dist2){
+ Collections.reverse(hits);
+ }
+ }
+}
+class Comp implements Comparator<TrackerHit>
+{
+ public int compare(TrackerHit hit1, TrackerHit hit2)
+ {
+ double z1 = hit1.getPosition()[2];
+ double z2 = hit2.getPosition()[2];
+
+ if (z1 < z2)
+ return -1;
+ else if (z1 > z2)
+ return 1;
+ else
+ return 0;
+ }
+}
lcsim/sandbox/LoriStevens
diff -N HelixSwimmer.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelixSwimmer.java 27 Jul 2007 18:07:47 -0000 1.1
@@ -0,0 +1,149 @@
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.constants.Constants;
+import org.lcsim.event.Track;
+import org.lcsim.spacegeom.SpacePoint;
+//import org.lcsim.util.swim.Helix;
+import org.lcsim.util.swim.Trajectory;
+import org.lcsim.util.swim.Line;
+
+import static java.lang.Math.abs;
+import static java.lang.Math.signum;
+
+/**
+ * A simple helix smimmer for use in org.lcsim. Uses standard lcsim units
+ * Tesla, mm, GeV. This swimmer works for charged and neutral tracks.
+ * <p>
+ * For more info on swimming see <a href="doc-files/transport.pdf">this paper</a>
+ * by Paul Avery.
+ * @author tonyj
+ * @version $Id: HelixSwimmer.java,v 1.1 2007/07/27 18:07:47 lstevens Exp $
+ */
+public class HelixSwimmer
+{
+ public class SpatialParameters
+ {
+ public double px;
+ public double py;
+ public double pz;
+ public int charge;
+ public boolean isInvalid = true;
+ }
+ private double field;
+ private Trajectory trajectory;
+ private SpatialParameters spatialParms;
+ private Track track;
+ /** Creates a new instance of HelixSwimmer
+ * @param B field strength in Tesla; uniform, solenoidal, directed along z-axis
+ */
+ public HelixSwimmer(double B)
+ {
+ field = B*Constants.fieldConversion;
+ spatialParms = new SpatialParameters();
+ }
+ /**
+ * Sets parameters for helix swimmmer.
+ *
+ * @param p 3-momentum (px,py,pz)
+ * @param r0 initial position(x0,y0,z0)
+ * @param iq charge iq = q/|e| = +1/0/-1
+ */
+ public void setTrack(Hep3Vector p, Hep3Vector r0, int iq)
+ {
+ double temp = p.x()*p.x()+p.y()*p.y();
+ double pmom = Math.sqrt(temp + p.z()*p.z());
+ double sin_lambda = p.z()/pmom;
+ double phi = Math.atan2(p.y(),p.x());
+ double lambda = Math.asin(sin_lambda);
+// no bug here
+// System.out.println("temp "+temp+" pmom "+pmom+" sin_lambda "+sin_lambda+" phi "+phi+" lambda "+lambda);
+
+ if (iq != 0 && field != 0)
+ {
+ double pt = Math.sqrt(temp);
+ double radius = pt/(iq*field);
+ trajectory = new Helix(r0,radius,phi,lambda);
+// System.out.println("r0 "+r0+" radius "+radius+" phi "+phi+" lambda "+lambda);
+// System.out.println("pt "+pt+" iq "+iq+" field "+field);
+ }
+ else
+ {
+ trajectory = new Line(r0, phi, lambda);
+// System.out.println("r0 "+r0+" phi "+phi+" lambda "+lambda);
+ }
+ spatialParms.isInvalid = true;
+ }
+
+ public void setTrack(Track t)
+ {
+ double omega = t.getTrackParameter(2);
+ double z0 = t.getTrackParameter(3);
+ double phi0 = t.getTrackParameter(1);
+ double lambda = Math.atan(t.getTrackParameter(4));
+ double d0 = t.getTrackParameter(0);
+ double[] ref = t.getReferencePoint();
+ // origin of the circle that is the x-y projection of the helix
+ Hep3Vector origin = new BasicHep3Vector( ref[0] -d0 * Math.sin(phi0), ref[1] + d0 * Math.cos(phi0), ref[2] + z0);
+ trajectory = new Helix(origin,1/omega,phi0,lambda);
+ spatialParms.isInvalid = true;
+ track=t;
+ }
+ public SpacePoint getPointAtDistance(double alpha)
+ {
+ if (trajectory == null) throw new RuntimeException("Trajectory not set");
+ return trajectory.getPointAtDistance(alpha);
+ }
+ public double getDistanceToRadius(double r)
+ {
+ if (trajectory == null) throw new RuntimeException("Trajectory not set");
+ System.out.println("distance "+trajectory.getDistanceToInfiniteCylinder(r));
+ return trajectory.getDistanceToInfiniteCylinder(r);
+ }
+ public double getDistanceToZ(double z)
+ {
+ if (trajectory == null) throw new RuntimeException("Trajectory not set");
+ double result = trajectory.getDistanceToZPlane(z);
+ if (result<0) result = trajectory.getDistanceToZPlane(-z);
+ return result;
+ }
+ public double getDistanceToCylinder(double r,double z)
+ {
+ double x1 = getDistanceToRadius(r); //problem here
+ double x2 = getDistanceToZ(z);
+ System.out.println("radius "+r);
+ System.out.println("x1 "+x1+" x2 "+x2); //Lori
+// System.out.println("min "+Math.min(x1,x2));
+// if(Double.isNaN(x1)) System.out.println("x2 "+x2); else System.out.println("minim "+Math.min(x1,x2));
+ return Double.isNaN(x1) ? x2 : Math.min(x1,x2);
+ }
+
+
+ /**
+ * Returns the distance along the trajectory to get to the point of closest approach
+ * @param point The point to swim as close as possible to
+ * @return the length parameter by how much the trajectory has to be swum
+ */
+ public double getDistanceToPoint(Hep3Vector point) {
+ return trajectory.getDistanceToPoint(point);
+ }
+
+ /**
+ * Calculates px, py, pz, iq from the Track Parameters and the B field.
+ * @return a Parameter object with the new parameters
+ */
+ // This unfortunately has to go here, because the Track parameters can't be used to get the momentum.
+ // The B field is needed in addition.
+ public SpatialParameters getSpatialParameters() {
+ if (spatialParms.isInvalid) {
+ double omega = track.getTrackParameter(2);
+ double Pt = abs((1./omega) * field);
+ spatialParms.px = Pt * Math.cos(track.getTrackParameter(1));
+ spatialParms.py = Pt * Math.sin(track.getTrackParameter(1));
+ spatialParms.pz = Pt * track.getTrackParameter(4);
+ spatialParms.charge = (int) signum(omega);
+ }
+ return spatialParms;
+ }
+
+ public Trajectory getTrajectory() { return trajectory; }
+}
lcsim/sandbox/LoriStevens
diff -N ZSeg.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ZSeg.java 27 Jul 2007 18:07:47 -0000 1.1
@@ -0,0 +1,132 @@
+
+/*
+ * 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.
+ */
+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 double segs;
+
+ /** Creates a new instance of ZSeg */
+
+ public ZSeg(){
+ segs = 100.0;
+ }
+
+ public ZSeg(double moduleLength) {
+ segs = moduleLength;
+ }
+
+ 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) {
+ double minz=0.0;
+ double maxz=0.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/sandbox/LoriStevens
diff -N ZSegTrackerHit.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ZSegTrackerHit.java 27 Jul 2007 18:07:47 -0000 1.1
@@ -0,0 +1,133 @@
+
+/*
+ * 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.
+ */
+import org.lcsim.event.TrackerHit;
+import java.util.ArrayList;
+
+/**
+ * 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 ZSegTrackerHit {
+ //segment length, in mm
+ static double segs;
+
+ /** Creates a new instance of ZSeg */
+
+ public ZSegTrackerHit(){
+ segs = 100.0;
+ }
+
+ public ZSegTrackerHit(double moduleLength) {
+ segs = moduleLength;
+ }
+
+ public static boolean zCheck(ArrayList<TrackerHit> 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 (TrackerHit hit1: hits){
+ if (hit1==null) continue;
+
+ for (TrackerHit hit2: hits){
+ if (hit2==null || hit2==hit1) continue;
+
+ for (TrackerHit hit3: hits){
+ if(hit3==null || hit3==hit1 || hit3==hit2) continue;
+
+ if (!zChecker(hit1, hit2, hit3))
+ return false;
+ }
+ }
+ }
+ return zCheck;
+ }
+
+ //checks whether or not 3rd hit falls within projected range
+ //includes check for whether to use interpolation or extrapolation
+ private static boolean zChecker(TrackerHit ht1, TrackerHit ht2, TrackerHit ht3){
+ boolean hitPass = false;
+
+ if (ht1!=null && ht2!=null && ht3 !=null){
+ double r1 = Math.sqrt(ht1.getPosition()[0]*ht1.getPosition()[0]+ht1.getPosition()[1]*ht1.getPosition()[1]);
+ double r2 = Math.sqrt(ht2.getPosition()[0]*ht2.getPosition()[0]+ht2.getPosition()[1]*ht2.getPosition()[1]);
+ double r3 = Math.sqrt(ht3.getPosition()[0]*ht3.getPosition()[0]+ht3.getPosition()[1]*ht3.getPosition()[1]);
+
+ double[] mod1 = moduleInfo(ht1.getPosition()[2], r1);
+ double[] mod2 = moduleInfo(ht2.getPosition()[2], r2);
+ double[] mod3 = moduleInfo(ht3.getPosition()[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) {
+ double minz=0.0;
+ double maxz=0.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
+ private 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/sandbox/LoriStevens
diff -N ZSeg_temp.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ZSeg_temp.java 27 Jul 2007 18:07:47 -0000 1.1
@@ -0,0 +1,181 @@
+ public ZSeg() {
+ }
+
+ public boolean hitLooper(SimTrackerHit [] hits){
+ boolean zCheck = true;
+ //looping through hits to check for all possible 3 hit combinations
+ 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;
+ }
+ public boolean hitLooper(SpacePoint [] sp){
+ boolean zCheck = true;
+ //looping through hits to check for all possible 3 hit combinations
+ for (int i=0; i<sp.length; i++){
+ if (sp[i]==null) continue;
+
+ for (int j=0; j<sp.length; j++){
+ if (sp[j]==null || sp[j]==sp[i]) continue;
+
+ for (int k=0; k<sp.length; k++){
+ if(sp[k]==null || sp[k]==sp[i] || sp[k]==sp[j]) continue;
+
+ if (!zChecker(sp[i], sp[j], sp[k]))
+ return false;
+ }
+ }
+ }
+ return zCheck;
+ }
+ public boolean hitLooper(double[] x, double[] y, double[] z){
+ boolean zCheck = true;
+ SpacePoint[] sp = new SpacePoint[x.length];
+ for(int i=0; i<x.length; i++)
+ {
+ sp[i] = new SpacePoint(new BasicHep3Vector(x[i], y[i], z[i]));
+ }
+ //looping through hits to check for all possible 3 hit combinations
+ for (int i=0; i<sp.length; i++){
+ if (sp[i]==null) continue;
+
+ for (int j=0; j<sp.length; j++){
+ if (sp[j]==null || sp[j]==sp[i]) continue;
+
+ for (int k=0; k<sp.length; k++){
+ if(sp[k]==null || sp[k]==sp[i] || sp[k]==sp[j]) continue;
+
+ if (!zChecker(sp[i], sp[j], sp[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 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;
+ }
+ //checks whether or not 3rd hit falls within projected range
+ //includes check for whether to use interpolation or extrapolation
+ public boolean zChecker(SpacePoint sp1, SpacePoint sp2, SpacePoint sp3){
+ boolean hitPass = false;
+
+ if (sp1!=null && sp2!=null && sp3 !=null){
+ double r1 = Math.sqrt(sp1.x()*sp1.x()+sp1.y()*sp1.y());
+ double r2 = Math.sqrt(sp2.x()*sp2.x()+sp2.y()*sp2.y());
+ double r3 = Math.sqrt(sp3.x()*sp3.x()+sp3.y()*sp3.y());
+
+ double[] mod1 = moduleInfo(sp1.z(), r1);
+ double[] mod2 = moduleInfo(sp2.z(), r2);
+ double[] mod3 = moduleInfo(sp3.z(), 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 double[] moduleInfo(double z, double r) {
+ int minz=0, 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;
+ }
+}
CVSspam 0.2.8
